! This program reads the standard geometry and normal modes
! from a frequency calculation by Gaussian 09 program package.
! The formula to do the distortion is
! x=x0+sqrt(hcut/(m*omega))*normal coordinates*factor
! where x0 is the standard optimized geometry in angstrom units
! hcut in atomic units is 1
! m is output in Gaussian program is amu. This has to be 
! converted into atomic units. The conversion factor is 
! given as c1. omega is output in Gaussian program is cm-1. 
! The conversion factor from cm-1 to au is c2. 
! The normal coordinates are in atomic units (in bohr). This has 
! to be converted into angstrom. The conversion factor is c3.
! The factor will be a small number between -1 to 1.


program distort
implicit none
integer, parameter :: rdp=selected_real_kind(15,100)      
integer, parameter :: i4=selected_int_kind(4)
integer (kind=i4)  :: i, ierr,its,ins,idummy,num,j,k,kend,k1,k2
real (kind=rdp)    :: x,y,z,factor
character (len=100) :: linebuf
character (len=6),allocatable, dimension(:) :: modelabel
character (len=21)  :: stanor='Standard orientation:'
character (len=19)    :: norcor='normal coordinates:'
character (len=2), allocatable, dimension(:)   :: atsymb
integer (kind=i4), allocatable, dimension(:)   :: atnum,modenum
real (kind=rdp), allocatable, dimension(:,:)  :: xyz,q,xyzplus,xyzmin,runit
real (kind=rdp), allocatable, dimension(:)  :: freq,mass,force
character (len=25) :: chardummy
character (len=80) :: fileout
real (kind=rdp), parameter :: cm_1_to_au= 4.5563352529120d-6  ! cm-1 unit to atomic unit
real (kind=rdp), parameter :: amu_to_au = 1.82289d3    ! amu to au unit
real (kind=rdp), parameter :: bohr_to_ang = 0.52917706 ! bohr to angrstong conversion

open(unit=10,file='acet-ortho-freq-iop.log',status='old',iostat=ierr,action='read')

fileinp: if (ierr == 0) then
    write(*,*)'The Gaussian frequency file exits'
    write(*,*)'I will proceed for further computation'
    exit fileinp
else
    write(*,*)'The Gaussian frequency file does not exits'
endif fileinp


outerdo: do
read(10,'(a100)',iostat=ierr)linebuf
if (ierr /= 0) then
exit outerdo
else
its=index(linebuf,stanor)
if (its /= 0 ) then
write(*,*)linebuf
exit outerdo
endif
endif
enddo outerdo

do i=1,4
read(10,'(a100)')linebuf
enddo

num=0
size: do
read(10,*,iostat=ierr)idummy,idummy,idummy,x,y,z
if (ierr /= 0 ) then
exit 
else
num=num+1
endif
enddo size

allocate (atnum(num),xyz(num,3),modelabel(3*num-6),freq(3*num-6))
allocate (mass(3*num-6),force(3*num-6),q(3*num,3*num-6),modenum(3*num-6))
allocate (atsymb(num),runit(3*num,3*num-6),xyzplus(num,3),xyzmin(num,3))

do i=1,num+1
backspace(10)
enddo

do i=1,num
read(10,*)idummy,atnum(i),idummy,(xyz(i,j),j=1,3)
enddo

innerdo: do
read(10,'(a120)',iostat=ierr)linebuf
if(ierr/=0) then
exit innerdo
endif
ins=index(linebuf,norcor)        
if (ins /=0) then
exit innerdo
endif
enddo innerdo

do k=1,3*num-6,5
if (k+4 > 3*num-6) then
kend=3*num-6
else
kend=k+4        
endif
read(10,*)(modenum(i),i=k,kend)
read(10,*)(modelabel(i),i=k,kend)
read(10,90)chardummy,(freq(i),i=k,kend)
read(10,90)chardummy,(mass(i),i=k,kend)
read(10,90)chardummy,(force(i),i=k,kend)
do i=1,2
read(10,'(a100)')linebuf
enddo
do i=1,3*num
read(10,*)idummy,idummy,idummy,(q(i,j),j=k,kend)
enddo
enddo

94 format(5(f8.5,2X))
write(*,*)'Enter the factor'
read(*,*)factor

write(*,*)'Enter the output file name'
read(*,'(a80)')fileout

open(unit=11,file=trim(fileout),status='new',action='write',iostat=ierr)
open(unit=12,file='f.out',status='new',action='write',iostat=ierr)

if(ierr /= 0) then
  Write(*,*)'The output file already exists. Please sort it out and re-execute'
else
  Write(*,*)'THe output will be written to ',fileout
endif
do i=1,num 
if(atnum(i)==6) then
atsymb(i)='C'
endif
if(atnum(i)==1) then
atsymb(i)='H'
endif
if(atnum(i)==8) then
atsymb(i)='O'
endif
if(atnum(i)==7) then
atsymb(i)='N'
endif
if(atnum(i)==5) then
atsymb(i)='B'
endif
if(atnum(i)==9) then
atsymb(i)='F'
endif

enddo

do i=1,3*num-6
do j=1,3*num
runit(j,i)=bohr_to_ang*sqrt(1_rdp/(amu_to_au*cm_1_to_au))*sqrt(1.0_rdp/(mass(i)*freq(i)))*q(j,i)
enddo

95 format(2x,f8.6)

xyzplus=0.0
do k=1,num
do k1=1,3
k2=3*(k-1)+k1
xyzplus(k,k1)=xyz(k,k1)+runit(k2,i)*factor
enddo
enddo

write(11,*)'$CONTRL SCFTYP=MCSCF RUNTYP=ENERGY UNITS=ANGS ISPHER=1'
write(11,*)'MPLEVL=2 $END'
write(11,*)'$SYSTEM MWORDS=1500 MEMDDI=15000 $END'
write(11,*)'$BASIS GBASIS=CCD $END'
write(11,*)'$GUESS GUESS=MOREAD NORB=602 $END'
write(11,*)'$MCSCF CISTEP=GUGA MAXIT=120 FULLNR=.T. FORS=.T.'
write(11,*)'DIABAT=.T. $END'
write(11,*)'$DRT GROUP=C1 STSYM=A FORS=.T. NMCC=127 NDOC=2 NVAL=2 $END'
write(11,*)'$GUGDIA NSTATE=8 ITERMX=120 $END'
write(11,*)'$GUGDM2 WSTATE(1)=1,1,1,1,1,1,1,1 $END'
write(11,*)'$MRMP MRPT=MCQDPT $END'
write(11,*)'$MCQDPT KSTATE(1)=1,1,1,1,1,1,1,1 XZERO=.T. EDSHFT=0.04 $END'
write(11,*)'$DIABAT SLCTTH=0.2 NMLAP=4 REFMOS=.T. REFGRP=.T. $END'
write(11,*)''
write(11,*)'$DATA'
write(11,*)'DISTORTION ALONG NORMAL COORDINATE ',i,'BY',factor
write(11,*)'C1'
do j=1, num
write(11,99)atsymb(j),real(atnum(j)),(xyzplus(j,k),k=1,3)
enddo
write(11,*)'$END'
write(11,*)''

xyzmin=0.0
do k=1,num
do k1=1,3
k2=3*(k-1)+k1
xyzmin(k,k1)=xyz(k,k1)-runit(k2,i)*factor
enddo
enddo

write(11,*)'$CONTRL SCFTYP=MCSCF RUNTYP=ENERGY UNITS=ANGS ISPHER=1'
write(11,*)'MPLEVL=2 $END'
write(11,*)'$SYSTEM MWORDS=1500 MEMDDI=15000 $END'
write(11,*)'$BASIS GBASIS=CCD $END'
write(11,*)'$GUESS GUESS=MOREAD NORB=602 $END'
write(11,*)'$MCSCF CISTEP=GUGA MAXIT=120 FULLNR=.T. FORS=.T.'
write(11,*)'DIABAT=.T. $END'
write(11,*)'$DRT GROUP=C1 STSYM=A FORS=.T. NMCC=127 NDOC=2 NVAL=2 $END'
write(11,*)'$GUGDIA NSTATE=8 ITERMX=120 $END'
write(11,*)'$GUGDM2 WSTATE(1)=1,1,1,1,1,1,1,1 $END'
write(11,*)'$MRMP MRPT=MCQDPT $END'
write(11,*)'$MCQDPT KSTATE(1)=1,1,1,1,1,1,1,1 XZERO=.T. EDSHFT=0.04 $END'
write(11,*)'$DIABAT SLCTTH=0.2 NMLAP=4 REFMOS=.T. REFGRP=.T. $END'
write(11,*)''
write(11,*)'$DATA'
write(11,*)'DISTORTION ALONG NORMAL COORDINATE ',i,'BY',factor
write(11,*)'C1'
do j=1, num
write(11,99)atsymb(j),real(atnum(j)),(xyzmin(j,k),k=1,3)
enddo
write(11,*)'$END'
write(11,*)''
enddo
do i=1,3*num-6
write(12,*)freq(i)
enddo
99 format(a4,2x,f5.2,2x,3(f12.6,2x))
90 format(a24,5(f8.4,2x))
91 format(5a5)
end program distort