Stupid xyz2om2 converter in fortran (cuz it was slow to do 1M files with a python converter)

xyz2om2.f90

! PUBLIC DOMAIN LICENSE 2017 BY ANDERS S. CHRISTENSEN
!
! I WROTE THIS BECAUSE I WAS BORED - I DON'T RECOMMEND
! WRITING FILE PARSERS IN FORTRAN BECAUSE IT IS NOT
! PRODUCTIVE.

program convert
    
    implicit none

    integer :: n

    character :: title

    double precision :: x, y, z
    double precision :: x1, y1, z1
    double precision :: x2, y2, z2
    double precision :: x3, y3, z3
    
    character(len=2) :: atom_type
    character(len=2) :: atom_type1
    character(len=2) :: atom_type2
    character(len=2) :: atom_type3

    double precision, dimension(3,3) :: xyz
    double precision, dimension(3) :: ba, bc
    double precision :: cosine_angle, angle

    double precision, parameter :: pi = 4 * atan(1.0d0)

    integer :: i
    integer :: nargs

    CHARACTER(len=32) :: arg
    CHARACTER(len=32) :: arg2

    integer :: opt

    integer :: optarg

    opt = 0

    ! 0 = Single-point energy
    ! 1 = Single-point gradient
    ! 2 = optimization

    optarg = 0

    nargs = iargc()
    if (nargs == 2) then
        call getarg(2, arg2)
        read(arg2, *) optarg
    endif


    if (optarg == 0) then
        ! Single-point
        write(*,*) "OM2 1SCF MULLIK PRECISE"
    else if (optarg == 1) then
        ! gradient
        write(*,*) "OM2 PRECISE jop=-2"
    else if (optarg == 2) then
        ! optimization
        write(*,*) "OM2 PRECISE"
        opt = 1
    endif

    write(*,*) 
    write(*,*) 

    CALL getarg(1, arg)
    open(unit=1, file=arg)
    read(1,*) n
    read(1,'(a)') title


    if (n > 3) then
        
        do i=1, n
            read(1,*) atom_type, x, y, z
            write(*,"(A3,F20.12,I3F20.12,I3,F20.12,I3)") atom_type, x, opt, y, opt,  z, opt
        enddo

    else if (n == 3) then

        read(1,*) atom_type1, xyz(1,1), xyz(2,1), xyz(3,1)
        read(1,*) atom_type2, xyz(1,2), xyz(2,2), xyz(3,2)
        read(1,*) atom_type3, xyz(1,3), xyz(2,3), xyz(3,3)

        ba = xyz(:,2) - xyz(:,1)
        bc = xyz(:,2) - xyz(:,3)
        
        cosine_angle = dot_product(ba, bc) / (norm2(ba) * norm2(bc))
        angle = acos(cosine_angle) / pi * 180.0d0


        write(*,"(A3)") atom_type1
        write(*,"(A3,F20.12,I3)") atom_type2, norm2(ba), opt
        write(*,"(A3,F20.12,I3F20.12,I3)") atom_type3, norm2(bc), opt, angle, opt

    else if (n == 2) then
        
        read(1,*) atom_type1, xyz(1,1), xyz(2,1), xyz(3,1)
        read(1,*) atom_type2, xyz(1,2), xyz(2,2), xyz(3,2)

        ba = xyz(:,2) - xyz(:,1)
        
        write(*,"(A3)") atom_type1
        write(*,"(A3,F20.12,I3)") atom_type2, norm2(ba), opt

    else if (n == 1) then
        
        read(1,*) atom_type1, xyz(1,1), xyz(2,1), xyz(3,1)
        
        write(*,"(A3)") atom_type1

    endif


    close(1)


end program convert