Skip to content

Instantly share code, notes, and snippets.

@machsix

machsix/readin.f90

Created April 17, 2018 05:13
Show Gist options
  • Save machsix/d0853834fa7422c76db60074d497f92f to your computer and use it in GitHub Desktop.
Save machsix/d0853834fa7422c76db60074d497f92f to your computer and use it in GitHub Desktop.
Download ZIP
Raw
readin.f90
subroutine readin
!
use param, only: mnmol, mnat, mnsurf, mngauss, mntmod, mntor, &
mnfrag, mnmol, mntraj, mnfrag, mnoutcome, &
mnbond, pi, amutoau, autoatm, autofs, autoev, kb, &
autoang
use c_output, only: idebug, minprinticon, minprinttrapz, maxprint, &
lwrite, nprint
use c_dd, only: lreadin
use c_struct, only: cell, nvibmods, bdfmthr, rijco, nattype, xx0, &
mm0, symbolm, icharge, imulti, iatom, natom, &
nat, nunit, natunit, symbol0, nrotmods, &
rotaxis, momiprin, mmofag, symbol0, mmm, xxm, &
charge, nmol, ldiat, icell, ldiat, nfrag, &
indatom, rij, cvec, xunit, linearmol, bdbrthr, &
invcvec, indmm, natoms
use c_sys, only: tflag, vfreq, taup, taut, sysfreq, iadjtemp, &
rampfact, nequil, n_ramp, nrampstep, e0param, &
nsurft, cparam, methflag, repflag, qcpack, &
itempzpe, itherm, dimen, p0, intflag, reactcol, &
hstep0, temp0, e0param, nramp, vol0, ntraj, &
trajlist, periodic, maxtraj, eps, potflag, &
ieqtherm, iadjpress
use c_initial, only: deltate, lfixte, ifreinitrot, ifnoang, &
demax, demin, ipicktrj, pickthr, repflgi, &
jjad, equilibrium, vibene, temp0im, ewell, &
rotsel, initx, initp, vibsel, ivibdist, &
ivibtype, vibdistn, ifeqwithrot, nsurfi, nsurf0, &
ifr0eff, vvad, itosave, nsave, ademod, stateselect, &
b_fix, bandwd, icompp, itunn, drr, ilinben, ibend, &
istr, ntunn, r0effect, r0collision, te0fixed, &
admode, nlinben, nimptor, nbend, nstr, ntor, etrans, &
qnmvib0, compp, comxx, icomxx, transelect, rran, &
eroti, qnmrot0, verte, ecolad, bparam, arrad, rrad0, &
itor, b0impact, escatad, imptor, &
bmaxad, bminad, phmode,irottemp,iscaleenergy,irotfreeze
use atomdata, only: atmradco, atomsymbols
use c_traj, only: frusmeth, stodecoflag, lmodpop, ralpha, ntrapz, &
ltunn, rbin, icfrag, heatcap, totenergy, potential,&
kinetic, cohesive, nstat, nreinit, lzerotrans, lzerocom,&
cre, cim, vol, erot, tuntype, ngauss, nvers, prodzpe, &
freqlim, tinyrho, eta
use c_term, only: nbondform, nbondbreak, agmerge, agdefrag, &
attrans, determ, n_frags, n_turn, termflag, &
trmax, tibreak, trmin, dihlim, dihd, ntorsion, &
tiform, t_gradmax, t_stime, t_nstep, n_atoffr, &
termsteps, t_gradmag
implicit none
include 'c_mopac.f90'
! inattmp: number of duplicated atom types in an AG
! symtmp: temp storage of atom symbols
! mmtmp: temp storage of atom weights
! relmu: Reduce mass between AGs
! rran(1) rran(1) is the diameter of the sphere
! ran random number
! ilwrite temp storage of lwrite information
! ismb temp storage of atom index
! atomdiatom: if a run is atom-diatom calculation. To minimize the
! modification to the old program, this is just a local variable
! trtype(o) = Bond threshold provided by the user, input in A,
! immediately converted to bohr.
! titype = atom pairs whose bond distance threshold will be changed
! by the user (atom symmbol as input)
! nbondtype = Number of bond types whose bond threshold will be changed
! initr0col = User provides r0collision for reative collision run
!
integer :: im,ii,i,j,k,l,ismb(mnoutcome,2),itotline,endreccon, &
outflag,ilwrite(1000),inattmp,i1,i2,nbondtype, &
ifkeyexist,ierror,ifblank,nmodel
integer :: ranseed
double precision :: radmin(mnmol),radmax(mnmol),mmtmp,maxtmp,mintmp, &
comxxm(3),volm(3,3),detmtbt,tmp, &
relmu,trtype(mnoutcome),ran
double precision :: small
logical :: initr0col,atomdiatom,rotst(mnmol), &
roteng(mnmol),ireadchg,find_real
character :: records*200,symbtmp*2,chartmp*2,lowercase*200
character*2 :: titype(mnoutcome,2)
! Default values set here
include 'defaults.f90'
! inititial values
b0impact=0.d0
initr0col = .false.
do i=1,mnat
mm0(i)=0.d0
do j=1,3
xx0(j,i)=0.d0
enddo
enddo
small = 1.d-10
!
! Read inputfile into records line by line
! lines with the FIRST!!! character as # will be deemed as comment line
! and will be ignored
! open(10,file='input.txt',status='OLD')
! file 12, store the input without comments
open(12,status='scratch')
! file 13, temporary file used to read in variables
open(13,status='scratch')
!
write(6,*) 'Your inputfile:'
do while(.true.)
read(5,'(a)',end=100) records
if (records(1:1) .ne. '#') then
records = lowercase(records)
write(12,'(a)') records
write(6,'(a)') records
endif
enddo
!
100 continue
! close(10)
write(12,'(a)') 'thend0'
rewind(12)
!
! Read in values
! After (length of keywords) + 1 '=' is the value we want.
! Write the rest of the line records to a temporary tmptxt
! and then read the value we want from the beginning of tmptxt
! This is the best way I can think of in Fortran 77
! The input can be written in two formats:
! 1. $system var1=1 var2=2 $end
! 2. $system var1=1
! var2=2
! $end
!
read(12,'(a)',end=9999) records
! ---------------------------------------------------------
! System control keywords and values
do 51 while(.true.)
if (index(records,'$control') .GT. 0) then
endreccon=0
do while(endreccon .ne. 9999)
! Reaction collision
if (index(records,'reactcol') .GT. 0) reactcol = .true.
! Potential used
ifkeyexist=index(records,'potflag')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('potflag')+1:)
rewind(13)
read(13,*) potflag
rewind(13)
endif
! QC package used
ifkeyexist=index(records,'qcpack')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('qcpack')+1:)
rewind(13)
read(13,*) qcpack
rewind(13)
endif
! Temperature desired to run the simulation
ifkeyexist=index(records,'temp0')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('temp0')+1:)
rewind(13)
read(13,*) temp0
rewind(13)
! Default temperature for each AG is temp0
do i=1,mnmol
temp0im(i) = temp0
enddo
endif
! Pressure
ifkeyexist=index(records,'pressure')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('pressure')+1:)
rewind(13)
read(13,*) p0
rewind(13)
p0=p0/autoatm
endif
! Diabatic or adiabatic propagation
! if (index(records,'diabatic') .GT. 0) repflag = 1
if (index(records,'diabatic').GT.0.AND.index(records, &
'adiabatic').EQ.0) repflag = 1
if (index(records,'periodic') .GT. 0) periodic = .true.
if (index(records,'nve') .GT. 0) itherm=0
if (index(records,'nvt') .GT. 0) itherm=1
if (index(records,'npt') .GT. 0) itherm=2
! Add zero vibration: at 0 K, the AG has a KE of 0.5 ZPE
if (index(records,'izpetemp') .GT. 0) itempzpe=.true.
! Integration method
! Bulirsch-Stoer integrator method
if (index(records,'bulstoint') .GT. 0) intflag=0
! 4th order Runge-Kutta integrator method
if (index(records,'runkut4') .GT. 0) intflag=1
! Simple Verlet algorithm
if (index(records,'verlet0') .GT. 0) intflag=2
if (index(records,'vverlet') .GT. 0) intflag=3
if (index(records,'beeman') .GT. 0) intflag=4
if (index(records,'liouville') .GT. 0) intflag=5
if (index(records,'bulstointhack') .GT. 0) intflag=6
if (index(records,'fixedte0') .GT. 0) lfixte=.true.
! Initial Total energy to be fixed
ifkeyexist=index(records,'te0fixed')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('te0fixed')+1:)
rewind(13)
read(13,*) te0fixed
rewind(13)
te0fixed=te0fixed/autoev
lfixte=.true.
endif
! The width of the TE can be varied
ifkeyexist=index(records,'deltate')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('deltate')+1:)
rewind(13)
read(13,*) deltate
rewind(13)
deltate=deltate/autoev
endif
! hstep0 initial time step
ifkeyexist=index(records,'hstep0')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('hstep0')+1:)
rewind(13)
read(13,*) hstep0
rewind(13)
hstep0=hstep0/autofs
endif
! eps if intflag=1
ifkeyexist=index(records,'eps')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('eps')+1:)
rewind(13)
read(13,*) eps
rewind(13)
endif
! Random number seed
ifkeyexist=index(records,'ranseed')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('ranseed')+1:)
rewind(13)
read(13,*) ranseed
rewind(13)
endif
! Atom-diatom calculation
if (index(records,'atomdiatom') .GT. 0) atomdiatom = .true.
! Meet the end of this section
if (index(records,'$end') .GT. 0) then
endreccon=9999
else
read(12,'(a)',end=9999) records
endif
enddo
read(12,'(a)',end=9999) records
endif
!
! ---------------------------------------------------------
! Surface information
if (index(records,'$surface') .GT. 0) then
endreccon=0
do while(endreccon .ne. 9999)
! Initial electronic surface
ifkeyexist=index(records,'nsurf0')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('nsurf0')+1:)
rewind(13)
read(13,*) nsurf0
rewind(13)
endif
! Initial conditions prepared from distribution of NURFI
ifkeyexist=index(records,'nsurfi')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('nsurfi')+1:)
rewind(13)
read(13,*) nsurfi
rewind(13)
endif
! Initial conditions prepared from distribution of
! NURFI,RepFlagi
ifkeyexist=index(records,'repflgi')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('repflgi')+1:)
rewind(13)
read(13,*) repflgi
rewind(13)
endif
! Total number of electronic surfaces
ifkeyexist=index(records,'nsurft')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('nsurft')+1:)
rewind(13)
read(13,*) nsurft
rewind(13)
endif
! Surface propagation method
ifkeyexist=index(records,'methflag')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('methflag')+1:)
rewind(13)
read(13,*) methflag
rewind(13)
endif
! Analysis of mode populations ON or OFF
if (index(records,'lmodpop') .GT. 0) lmodpop = .true.
! Coefficient for the determination of the amount of vibration we want
! to keep to calculate mode populations
ifkeyexist=index(records,'ralpha')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('ralpha')+1:)
rewind(13)
read(13,*) ralpha
rewind(13)
endif
! Methods to treat frustrated hops
ifkeyexist=index(records,'frusmeth')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('frusmeth')+1:)
rewind(13)
read(13,*) frusmeth
rewind(13)
endif
! Use of a Projected Hessian or TRAPZ-like methods ?
ifkeyexist=index(records,'ntrapz')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('ntrapz')+1:)
rewind(13)
read(13,*) ntrapz
rewind(13)
endif
if ((ntrapz.eq.2.or.ntrapz.eq.3).and.methflag.ne.1.and. &
methflag.ne.5) then
write(6,*) 'TRAPZ can only be used with FS, FSTU or FSTU/SD'
stop
endif
! Which version of TRAPZ ?
ifkeyexist=index(records,'nvers')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('nvers')+1:)
rewind(13)
read(13,*) nvers
rewind(13)
endif
! Criterion to avoid entering TRAPZ-like methods
ifkeyexist=index(records,'freqlim')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('freqlim')+1:)
rewind(13)
read(13,*) freqlim
rewind(13)
endif
! Value of the product ZPE for the mTRAPZ* method (nvers = 2)
! prodzpe must be expressed in eV in the input file
ifkeyexist=index(records,'prodzpe')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('prodzpe')+1:)
rewind(13)
read(13,*) prodzpe
rewind(13)
prodzpe = prodzpe/autoev
if (nvers.ne.2) then
write(6,*) 'mTRAPZ* is not selected: prodzpe is not used'
write(6,*)
endif
endif
! Use of stochastic decoherence?
ifkeyexist=index(records,'stodecoflag')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('stodecoflag')+1:)
rewind(13)
read(13,*) stodecoflag
rewind(13)
endif
if (stodecoflag.eq.1.and.methflag.ne.1.and. &
methflag.ne.5) then
write(6,*) 'Stochastic decoherence can only be used', &
' with FS or FSTU'
stop
endif
! Cparam: For CSDM and SCDM methods
ifkeyexist=index(records,'cparam')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('cparam')+1:)
rewind(13)
read(13,*) cparam
rewind(13)
endif
! E0param: For CSDM and SCDM methods
ifkeyexist=index(records,'e0param')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('e0param')+1:)
rewind(13)
read(13,*) e0param
rewind(13)
endif
! tinyrho to avoid dividing by zero in DM methods
ifkeyexist=index(records,'tinyrho')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('tinyrho')+1:)
rewind(13)
read(13,*) tinyrho
rewind(13)
endif
! Meet the end of this section
if (index(records,'$end') .GT. 0) then
endreccon=9999
else
read(12,'(a)',end=9999) records
endif
enddo
read(12,'(a)',end=9999) records
endif
!
! ---------------------------------------------------------
! Reaction container information (cell)
if (index(records,'$cell') .GT. 0) then
endreccon=0
do while(endreccon .ne. 9999)
! cubic
if (index(records,'cubic') .GT. 0) then
read(12,*) cell(1)
icell=1
endif
! Cuboid
if (index(records,'cuboid') .GT. 0) then
read(12,*) (cell(j),j=1,3)
icell=2
endif
! Triclinic
if (index(records,'arbitrary') .GT. 0) then
read(12,*) (cell(j),j=1,6)
icell=3
endif
! Triclinic
if (index(records,'isolate') .GT. 0) then
icell=0
endif
! Meet the end of this section
if (index(records,'$end') .GT. 0) then
endreccon=9999
else
read(12,'(a)',end=9999) records
endif
enddo
! Move to next line
read(12,'(a)',end=9999) records
endif
!
! ---------------------------------------------------------
! Reaction collision parameters
if (index(records,'$rxcollision') .GT. 0) then
endreccon=0
reactcol=.true.
do while(endreccon .ne. 9999)
! State select
if (index(records,'equilibrium') .GT. 0) stateselect=.false.
if (index(records,'stateselect') .GT. 0) stateselect=.true.
ifkeyexist=index(records,'b0impact')
! Impact parameter. If provided, it would be fixed
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('b0impact')+1:)
rewind(13)
read(13,*) b0impact
rewind(13)
b0impact=b0impact/autoang
b_fix = .true.
endif
! If user provide r0collision value, at the same time it means program
! will not decide r0collision itself, thus initr0col = 1
! Initial distance information between the two AGs
ifkeyexist=index(records,'r0collision')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('r0collision')+1:)
rewind(13)
read(13,*) r0collision
rewind(13)
r0collision=r0collision/autoang
initr0col = .true.
endif
! Effective collision radius between two reactants
ifkeyexist=index(records,'r0effectcol')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('r0effectcol')+1:)
rewind(13)
read(13,*) r0effect
rewind(13)
r0effect=r0effect/autoang
ifr0eff=.true.
endif
! Meet the end of this section
if (index(records,'$end') .GT. 0) then
endreccon=9999
else
read(12,'(a)',end=9999) records
endif
enddo
! Move to next line
read(12,'(a)',end=9999) records
endif
!
! ---------------------------------------------------------
! Termination condition
if (index(records,'$termcon') .GT. 0) then
endreccon=0
do while(endreccon .ne. 9999)
! Termination Flag
ifkeyexist=index(records,'termflag')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('termflag')+1:)
rewind(13)
read(13,*) termflag
rewind(13)
endif
! For fragment monitoring (termflag=4) only
if (index(records,'noagdefrag') .GT. 0) agdefrag=.false.
if (index(records,'noagmerge') .GT. 0) agmerge=.false.
if (index(records,'noattrans') .GT. 0) attrans=.false.
! Number of products to be monitored in unimolecular
! fragmentation monintor condition (termflag=7)
ifkeyexist=index(records,'n_frags')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('n_frags')+1:)
rewind(13)
read(13,*) n_frags
rewind(13)
endif
! For termflag =3 4 5 only
ifkeyexist=index(records,'termsteps')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('termsteps')+1:)
rewind(13)
read(13,*) termsteps
rewind(13)
endif
! Maximum number of steps per trajectory
ifkeyexist=index(records,'t_nstep')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('t_nstep')+1:)
rewind(13)
read(13,*) t_nstep
rewind(13)
! if t_nstep <0, unlimited, set to a very large value
if (t_nstep.lt.0) t_nstep=1000000000
endif
! Simulation time
ifkeyexist=index(records,'t_stime')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('t_stime')+1:)
rewind(13)
read(13,*) t_stime
rewind(13)
t_stime = t_stime/autofs
endif
ifkeyexist=index(records,'t_gradmag')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('t_gradmag')+1:)
rewind(13)
read(13,*) t_gradmag
rewind(13)
t_gradmag=t_gradmag*autoang/autoev
endif
ifkeyexist=index(records,'t_gradmax')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('t_gradmax')+1:)
rewind(13)
read(13,*) t_gradmax
rewind(13)
t_gradmax=t_gradmax*autoang/autoev
endif
! Total energy threshold in termflag=4 if te < determ, term
! Currently not used
! ifkeyexist=index(records,'determ')
! if (ifkeyexist .GT. 0) then
! write(13,*) records(ifkeyexist+len('determ')+1:)
! rewind(13)
! read(13,*) determ
! rewind(13)
! endif
! used for reactcol, if reactant turn n_turn times, reactive
ifkeyexist=index(records,'n_turn')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('n_turn')+1:)
rewind(13)
read(13,*) n_turn
rewind(13)
endif
! Bond forming threshold
ifkeyexist=index(records,'bondfmthre')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('bondfmthre')+1:)
rewind(13)
read(13,*) bdfmthr
rewind(13)
endif
! Bond broken threshold
ifkeyexist=index(records,'bondbrthre')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('bondbrthre')+1:)
rewind(13)
read(13,*) bdbrthr
rewind(13)
endif
ifkeyexist=index(records,'nbondform')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('nbondform')+1:)
rewind(13)
read(13,*) nbondform
rewind(13)
! These pair of distances must be written in seperated
! lines line by line.
do j=1,nbondform
! Read from the next line below the line of nbondform
read(12,'(a)',end=9999) records
write(13,*) records
rewind(13)
read(13,*) tiform(j,1),tiform(j,2),trmin(j)
rewind(13)
enddo
read(12,'(a)',end=9999) records
endif
ifkeyexist=index(records,'nbondbreak')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('nbondbreak')+1:)
rewind(13)
read(13,*) nbondbreak
rewind(13)
! These pair of distances must be written in separated
! lines, line by line.
do j=1,nbondbreak
! Read from the next line below the line of nbondbreak
read(12,'(a)',end=9999) records
write(13,*) records
rewind(13)
read(13,*) tibreak(j,1),tibreak(j,2),trmax(j)
rewind(13)
enddo
read(12,'(a)',end=9999) records
endif
ifkeyexist=index(records,'nbondtype')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('nbondtype')+1:)
rewind(13)
read(13,*) nbondtype
rewind(13)
! These pair of distances must be written in separate
! lines, line by line.
do j=1,nbondtype
! Read from the next line below the line of nbondtype
read(12,'(a)',end=9999) records
write(13,*) records
rewind(13)
read(13,*) titype(j,1),titype(j,2),trtype(j)
trtype(j)=trtype(j)/autoang
do i=1,111
if ( titype(j,1).eq.atomsymbols(i) ) ismb(j,1)=i
if ( titype(j,2).eq.atomsymbols(i) ) ismb(j,2)=i
enddo
rijco( ismb(j,1),ismb(j,2) ) = trtype(j)
rijco( ismb(j,2),ismb(j,1) ) = trtype(j)
rewind(13)
enddo
read(12,'(a)',end=9999) records
endif
! The fragments to be monitored. termflag=7 only
ifkeyexist=index(records,'n_atomtypes')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('n_atomtypes')+1:)
rewind(13)
read(13,*) j
rewind(13)
n_atoffr=0
do i=1,j
! Read from the next line below the line of n_atomtypes
read(12,*) k,nattype(mnmol+1,k)
n_atoffr=n_atoffr+nattype(mnmol+1,k)
enddo
read(12,'(a)',end=9999) records
endif
!
! JZ torsion angles
ifkeyexist=index(records,'ntorsion')
if (ifkeyexist .gt. 0) then
write(13,*) records(ifkeyexist+len('ntorsion')+1:)
rewind(13)
read(13,*) ntorsion
rewind(13)
do j = 1, ntorsion
read(12,'(a)',end=9999) records
write(13,*) records
rewind(13)
read(13,*) dihd(1,j),dihd(2,j),dihd(3,j),dihd(4,j), &
dihlim(1,j),dihlim(2,j)
rewind(13)
if(dihlim(1,j).lt.0.d0.or.dihlim(2,j).lt.0.d0) then !RMP
write(6,*) "The values that defined the limits of the" & !RMP
,"torsional angles must be positive." !RMP
stop
endif
enddo
read(12,'(a)',end=9999) records
endif
! Meet the end of this section
if (index(records,'$end') .GT. 0) then
endreccon=9999
else
read(12,'(a)',end=9999) records
endif
enddo
! Move to next line
read(12,'(a)',end=9999) records
endif
!
! Trajectory control information
if (index(records,'$traject') .GT. 0) then
endreccon=0
do while(endreccon .ne. 9999)
ifkeyexist=index(records,'ntraj')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('ntraj')+1:)
rewind(13)
read(13,*) ntraj
rewind(13)
if (ntraj.gt.mntraj) then
write(6,111) 'NTRAJ = ',ntraj,' exceeds the maximum', &
' number of trajectories the program set: ',mntraj
stop
endif
! default: maxtraj = ntraj
maxtraj = ntraj
endif
ifkeyexist=index(records,'tflag1')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('tflag1')+1:)
rewind(13)
read(13,*) tflag(1)
rewind(13)
if (tflag(1).lt.0 .or. tflag(1).gt.5) then
write(6,*) 'TFLAG1 != 0, 1, 2, 3, 4, 5 is not supported'
stop
endif
endif
! For termflag=7 only. Equilibrium to provide initial Mom and
! Structure
if (index(records,'noequilibrium') .GT. 0) stateselect=.true.
! Save coord. and momentum from itosave'th step
ifkeyexist=index(records,'itosave')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('itosave')+1:)
rewind(13)
read(13,*) itosave
rewind(13)
endif
! Save coord. and momentum every nsave steps
! ifkeyexist=index(records,'nsave')
! if (ifkeyexist .GT. 0) then
! write(13,*) records(ifkeyexist+len('nsave')+1:)
! rewind(13)
! read(13,*) nsave
! rewind(13)
! endif
! ensemble used for equilibrium run (react col and termflag=7)
ifkeyexist=index(records,'ieqtherm')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('ieqtherm')+1:)
rewind(13)
read(13,*) ieqtherm
rewind(13)
if (ieqtherm.lt.0 .or. ieqtherm.gt.2) then
write(6,*) 'ERROR: IEQTHERM = 0,1,2 only'
stop
endif
endif
! Temperature adjusting method
ifkeyexist=index(records,'iadjtemp')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('iadjtemp')+1:)
rewind(13)
read(13,*) iadjtemp
rewind(13)
if (iadjtemp.lt.0 .or. iadjtemp.gt.3) then
write(6,*) 'ERROR: IADJTEMP = 0,1,2,3 only'
stop
endif
endif
! Pressure adjusting method
ifkeyexist=index(records,'iadjpress')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('iadjpress')+1:)
rewind(13)
read(13,*) iadjpress
rewind(13)
endif
! Coupling parameter between the bath and the system
ifkeyexist=index(records,'taut')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('taut')+1:)
rewind(13)
read(13,*) taut
rewind(13)
endif
ifkeyexist=index(records,'taup')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('taup')+1:)
rewind(13)
read(13,*) taup
rewind(13)
endif
!C 'Mass' used for Nose-Hoover thermostat
! ifkeyexist=index(records,'qnose1')
! if (ifkeyexist .GT. 0) then
! write(13,*) records(ifkeyexist+len('qnose1')+1:)
! rewind(13)
! read(13,*) nh(1,1)
! rewind(13)
! endif
! ifkeyexist=index(records,'qnose2')
! if (ifkeyexist .GT. 0) then
! write(13,*) records(ifkeyexist+len('qnose2')+1:)
! rewind(13)
! read(13,*) nh(1,2)
! rewind(13)
! endif
! Charactoristic vibration frequency of the system (1/fs) used
! for Nose-Hoover thermostat
ifkeyexist=index(records,'sysfreq')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('sysfreq')+1:)
rewind(13)
read(13,*) sysfreq
rewind(13)
endif
! Collision frequency used by Andersen thermostat
ifkeyexist=index(records,'vfreq')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('vfreq')+1:)
rewind(13)
read(13,*) vfreq
rewind(13)
endif
! If remove CoM angular before entering DRIVER
if (index(records,'keepangmom') .GT. 0) ifnoang = .false.
if (index(records,'withrotineq') .GT.0) ifeqwithrot= .true.
if (index(records,'noreinitrot') .GT.0) ifreinitrot= .false.
if (index(records,'scaleenergy') .GT.0) iscaleenergy = .true.
! Remove angmom every nreinit steps
ifkeyexist=index(records,'nreinit')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('nreinit')+1:)
rewind(13)
read(13,*) nreinit
rewind(13)
endif
! If reset the CoM to origin
if (index(records,'nozerocom') .GT. 0) lzerocom = .false.
! If reset the CoM to origin
if (index(records,'keeptrans') .GT. 0) lzerotrans = .false.
! trajectory picking method
ifkeyexist=index(records,'ipicktrj')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('ipicktrj')+1:)
rewind(13)
read(13,*) ipicktrj
rewind(13)
if (ipicktrj.lt.0 .or. ipicktrj.gt.2) then
write(6,*) 'ERROR: only 3 choices for ipicktrj: 0,1,2'
stop
endif
endif
! Probability of choose trajectories in Equil run
ifkeyexist=index(records,'pickthr')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('pickthr')+1:)
rewind(13)
read(13,*) pickthr
rewind(13)
endif
! TE threshold (low) to pick trajectories in Equil run
ifkeyexist=index(records,'demin')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('demin')+1:)
rewind(13)
read(13,*) demin
rewind(13)
endif
! TE threshold (upper) to pick trajectories in Equil run
ifkeyexist=index(records,'demax')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('demax')+1:)
rewind(13)
read(13,*) demax
rewind(13)
endif
! change temperature at ramptime
ifkeyexist=index(records,'rampstep0')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('rampstep0')+1:)
rewind(13)
read(13,*) nrampstep
rewind(13)
endif
ifkeyexist=index(records,'rampfact')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('rampfact')+1:)
rewind(13)
read(13,*) rampfact
rewind(13)
endif
! ramping time interval
ifkeyexist=index(records,'n_ramp')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('n_ramp')+1:)
rewind(13)
read(13,*) n_ramp
rewind(13)
endif
! equilibrate for tequil time after a ramping interval
ifkeyexist=index(records,'nequilramp')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('nequilramp')+1:)
rewind(13)
read(13,*) nequil
rewind(13)
endif
! ramping times
ifkeyexist=index(records,'nramp')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('nramp')+1:)
rewind(13)
read(13,*) nramp
rewind(13)
endif
! restarting trajectory information
if (index(records,'restart') .GT. 0) then
tflag(2) = 1
! trajectory list must be from the next line of restart
read(12,*) (trajlist(k),k=1,ntraj)
maxtraj = trajlist(1)
do k=1,ntraj
maxtraj = max(maxtraj,trajlist(k))
enddo
read(12,'(a)',end=9999) records
endif
! End of records
if (index(records,'$end') .GT. 0) then
endreccon=9999
else
read(12,'(a)',end=9999) records
endif
enddo
! Move to next line
read(12,'(a)',end=9999) records
endif
! JZ tunneling
if(index(records,'$tunneling') .GT.0) then
ntunn = 0
endreccon=0
ltunn = .true.
do while(endreccon.ne.9999)
ifkeyexist=index(records,'eta')
if(ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('eta')+1:)
rewind(13)
read(13,*) eta
rewind(13)
if (eta.lt.0.d0.or.eta.gt.1.d0) then !RMP
write(6,*) "The value of eta must be ",& !RMP
"between 0 and 1." !RMP
stop !RMP
endif !RMP
endif
! ifkeyexist=index(records,'probthr')
! if(ifkeyexist .GT. 0) then
! write(13,*) records(ifkeyexist+len('probthr')+1:)
! rewind(13)
! read(13,*) actint_thr
! rewind(13)
! actint_thr = -0.5d0*log(actint_thr)
! endif
ifkeyexist=index(records,'ngauss')
if(ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('ngauss')+1:)
rewind(13)
read(13,*) ngauss
rewind(13)
if (ngauss.gt.512) then !RMP
write(6,*) "Number of Gaussian-quadrature nodes must be ", &
"no more than 512."
stop
endif
endif
! ifkeyexist=index(records,'ntunn_mode') ! RMP remove this part correspond to the normal modes
! if (ifkeyexist .gt. 0) then
! write(13,*) records(ifkeyexist+len('ntunn_mode')+1:)
! rewind(13)
! read(13,*) ntunn_mode
! rewind(13)
! read(12,'(a)',end=9999) records
! write(13,*) records
! rewind(13)
! read(13,*) (itmod(j),j=1,ntunn_mode)
! rewind(13)
! endif
! streching internal coordinates
ifkeyexist=index(records,'nstr')
if (ifkeyexist .gt. 0) then
write(13,*) records(ifkeyexist+len('nstr')+1:)
rewind(13)
read(13,*) nstr
rewind(13)
k = 0
do i = 1, nstr
k = k + 1
read(12,'(a)',end=9999) records
write(13,*) records
rewind(13)
if(index(records,'*').gt.0) then
ntunn = ntunn + 1
itunn(k) = 1
tuntype = 1
drr(k) = 1d0
read(13,*) (istr(j,i),j=1,2)
elseif(index(records,'+').gt.0) then
itunn(k) = 1
tuntype = 2
read(13,*) (istr(j,i),j=1,2),chartmp,drr(k)
elseif(index(records,'-').gt.0) then
itunn(k) = -1
tuntype = 2
read(13,*) (istr(j,i),j=1,2),chartmp,drr(k)
drr(k) = -drr(k)
else
read(13,*) (istr(j,i),j=1,2)
endif
rewind(13)
enddo
endif
! bond angles
ifkeyexist=index(records,'nbend')
if (ifkeyexist .gt. 0) then
write(13,*) records(ifkeyexist+len('nbend')+1:)
rewind(13)
read(13,*) nbend
rewind(13)
do i = 1, nbend
k = k + 1
! RMP we should double check all the possible options in the input file
read(12,'(a)',end=9999) records
if(index(records,'*').gt.0) then
ntunn = ntunn + 1
itunn(k) = 1
tuntype = 1
drr(k) = 1d0
! elseif(index(records,'+').gt.0) then ! RMP it is not needed
! itunn(k) = 1
! tuntype = 2
! elseif(index(records,'-').gt.0) then
! itunn(k) = -1
! tuntype = 2
endif
write(13,*) records
rewind(13)
read(13,*) (ibend(j,i),j=1,3)
rewind(13)
enddo
endif
! RMP
! linear bendings
!integer :: nlinben,ilinben(3,3*mnat) in c_initial.f
ifkeyexist=index(records,'linben')
if (ifkeyexist .gt. 0) then
write(13,*) records(ifkeyexist+len('linben')+1:)
read(13,*) nlinben
rewind(13)
do i = 1, nlinben
k = k + 1
read(12,'(a)',end=9999) records
if(index(records,'*').gt.0) then
ntunn = ntunn + 1
itunn(k) = 1
tuntype = 1
drr(k) = 1d0
! elseif(index(records,'+').gt.0) then !RMP is not needed
! itunn(k) = 1
! tuntype = 2
! elseif(index(records,'-').gt.0) then
! itunn(k) = -1
! tuntype = 2
endif
write(13,*) records
rewind(13)
read(13,*) (ilinben(j,i),j=1,3)
rewind(13)
enddo
endif
! dihedral angles
ifkeyexist=index(records,'ntor')
if (ifkeyexist .gt. 0) then
write(13,*) records(ifkeyexist+len('ntor')+1:)
rewind(13)
read(13,*) ntor
rewind(13)
do i = 1, ntor
k = k + 1
read(12,'(a)',end=9999) records
if(index(records,'*').gt.0) then
ntunn = ntunn + 1
itunn(k) = 1
tuntype = 1
drr(k) = 1d0
! elseif(index(records,'+').gt.0) then !RMP is not needed
! itunn(k) = 1
! tuntype = 2
! elseif(index(records,'-').gt.0) then
! itunn(k) = -1
! tuntype = 2
endif
write(13,*) records
rewind(13)
read(13,*) (itor(j,i),j=1,4)
rewind(13)
enddo
endif
! out of plane / improper torsion
ifkeyexist=index(records,'nimptor')
if (ifkeyexist .gt. 0) then
write(13,*) records(ifkeyexist+len('nimptor')+1:)
rewind(13)
read(13,*) nimptor
rewind(13)
do i = 1, nimptor
k = k + 1
read(12,'(a)',end=9999) records
if(index(records,'*').gt.0) then
ntunn = ntunn + 1
itunn(k) = 1
tuntype = 1
drr(k) = 1d0
! elseif(index(records,'+').gt.0) then ! RMP is not needed
! itunn(k) = 1
! tuntype = 2
! elseif(index(records,'-').gt.0) then
! itunn(k) = -1
! tuntype = 2
endif
write(13,*) records
rewind(13)
read(13,*) (imptor(j,i),j=1,4)
rewind(13)
enddo
endif
if(index(records,'$end') .GT. 0) then
endreccon=9999
else
read(12,'(a)',end=9999) records
endif
enddo
! Move to next line
read(12,'(a)',end=9999) records
endif
!
! JZ MOPAC
if(index(records,'$mopac') .GT.0) then
potflag = -2
endreccon=0
read(12,'(a)',end=9999) records
do while(endreccon.ne.9999)
if(index(records,'$end') .GT. 0) then
endreccon=9999
else
keywrd = records
write(6,*) 'keywrd',keywrd
read(12,'(a)',end=9999) records
endif
enddo
! Move to next line
read(12,'(a)',end=9999) records
endif
!
!
! Output writing list
if (index(records,'$output') .GT. 0) then
endreccon=0
do while(endreccon .ne. 9999)
! Print information every nprint steps
ifkeyexist=index(records,'nprint')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('nprint')+1:)
rewind(13)
read(13,*) nprint
rewind(13)
endif
! maxprint controls if print information in detail
if (index(records,'maxprint') .GT. 0) maxprint = .true.
! minprinticon removes printing of selection of initial conditions
if (index(records,'minprinticon') .GT. 0) &
minprinticon = .true.
! minprinttrapz removes intermediate information from TRAPZ
if (index(records,'minprinttrapz') .GT. 0) &
minprinttrapz = .true.
! idebug controls if print information in detail
if (index(records,'idebug') .GT. 0) then
idebug = .true.
maxprint = .true.
endif
! outputflag
ifkeyexist=index(records,'outflag')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('outflag')+1:)
rewind(13)
read(13,*) outflag
rewind(13)
if (outflag.eq.0) then
do j=1,99
lwrite(j)=.true.
enddo
elseif (outflag.gt.0 .and. outflag.le.99) then
do j=1,1000
lwrite(j)=.false.
enddo
! readin output writting list
read(12,*,iostat=ierror) (ilwrite(k),k=1,outflag)
if (ierror.ne.0) then
write(6,222) 'ERROR: output list reading error.', &
' Last read: ',ilwrite(k-1)
stop
endif
do j=1,outflag
lwrite(ilwrite(j)) = .true.
if (lwrite(48)) lwrite(47)=.true.
enddo
else
write(6,*) 'OUTFLAG= ',outflag,' is not supported'
stop
endif
endif
! End of records
if (index(records,'$end') .GT. 0) then
endreccon=9999
else
read(12,'(a)',end=9999) records
endif
enddo
! Move to next line
read(12,'(a)',end=9999) records
endif
!
! ---------------------------------------------------------
! Atom-diatom section
if (index(records,'$atomdiatom') .GT. 0) then
! if input for $atomdiatom appears, it is a special triatom run
atomdiatom = .true.
endreccon=0
do while(endreccon .ne. 9999)
! Total energy in eV
ifkeyexist=index(records,'escatad')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('escatad')+1:)
rewind(13)
read(13,*) escatad
rewind(13)
endif
! Initial vibrational quantum number for the diatom
ifkeyexist=index(records,'vvad')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('vvad')+1:)
rewind(13)
read(13,*) vvad
rewind(13)
endif
! Initial rotational quantum number for the diatom
ifkeyexist=index(records,'jjad')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('jjad')+1:)
rewind(13)
read(13,*) jjad
rewind(13)
endif
! Initial atom-diatom separation in Angstrom
ifkeyexist=index(records,'rrad0')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('rrad0')+1:)
rewind(13)
read(13,*) rrad0
rewind(13)
endif
! Initial molecular arrangement, 1: AB+C, 2: BC+A, 3: AC+B
ifkeyexist=index(records,'arrad')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('arrad')+1:)
rewind(13)
read(13,*) arrad
rewind(13)
endif
! RMP14
! Impact parameter
ifkeyexist=index(records,'bparam')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('bparam')+1:)
rewind(13)
read(13,*) bparam
bparam=bparam/autoang
rewind(13)
endif
ifkeyexist=index(records,'bmaxad')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('bmaxad')+1:)
rewind(13)
read(13,*) bmaxad ! in angstrom
bmaxad=bmaxad/autoang
rewind(13)
endif
ifkeyexist=index(records,'bminad')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('bminad')+1:)
rewind(13)
read(13,*) bminad ! in angstrom
bminad=bminad/autoang
rewind(13)
endif
! RMP14
! Atom diatom Modes:
ifkeyexist=index(records,'mode')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('mode')+1:)
rewind(13)
read(13,*) admode
if(admode.lt.0.or.admode.gt.3) then
write(6,*) "Error in Atom-diatom mode"
write(6,*) "Currently only mode=1,2,3 are available."
stop
endif
rewind(13)
endif
! Atom diatom Phase Modes:
ifkeyexist=index(records,'phmode')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('phmode')+1:)
rewind(13)
read(13,*) phmode
if(phmode.lt.1.or.phmode.gt.2) then
write(6,*) "Error in Atom-diatom phase mode"
write(6,*) "Currently only mode=1,2 are available."
stop
endif
rewind(13)
endif
! RMP14
! Atom diatom Modes:
ifkeyexist=index(records,'emode')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('emode')+1:)
rewind(13)
read(13,*) ademod
if(ademod.lt.1.or.ademod.gt.2) then
write(6,*) "Error in Atom-diatom energy mode"
write(6,*) "The available options are:"
write(6,*) " emode = 1 Scattering energy by input"
write(6,*) " emode = 2 Collision energy by input"
stop
endif
rewind(13)
endif
! RMP14
! Atom diatom Modes:
ifkeyexist=index(records,'ecoli')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('ecoli')+1:)
rewind(13)
read(13,*) ecolad
rewind(13)
endif
! RMP14
! Controls how the diatomic potential is computed
ifkeyexist=index(records,'flagdiat')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('flagdiat')+1:)
rewind(13)
ldiat=.true.
endif
! Meet the end of this section
if (index(records,'$end') .GT. 0) then
endreccon=9999
else
read(12,'(a)',end=9999) records
endif
enddo
! Move to next line
read(12,'(a)',end=9999) records
endif
!
! Results analysis control deck
if (index(records,'$analysis') .GT. 0) then
endreccon=0
do while(endreccon .ne. 9999)
ifkeyexist=index(records,'nstat')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('nstat')+1:)
rewind(13)
read(13,*) nstat
rewind(13)
endif
if (index(records,'cohesive') .GT. 0) cohesive=.true.
if (index(records,'kinetic') .GT. 0) kinetic=.true.
if (index(records,'potential') .GT. 0) potential=.true.
if (index(records,'totenergy') .GT. 0) totenergy=.true.
if (index(records,'fragcon') .GT. 0) icfrag=.true.
ifkeyexist=index(records,'heatcapacity')
if (ifkeyexist .GT. 0) then
heatcap=.true.
totenergy=.true.
endif
! Bin distanc used for the bin
ifkeyexist=index(records,'rbin')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('rbin')+1:)
rewind(13)
read(13,*) rbin
rewind(13)
rbin = rbin/autoang
endif
! End of records
if (index(records,'$end') .GT. 0) then
endreccon=9999
else
read(12,'(a)',end=9999) records
endif
enddo
! Move to next line
read(12,'(a)',end=9999) records
endif
! ---------------------------------------------------------
! Coordinates and other data for every AG goes here. Here the format is
! relatively fixed. No need to use endreccon to judge if the section is
! finished
if (index(records,'$data') .GT. 0) then
! Number of molecules or atom groups (AGs)
read(12,'(a)',end=9999) records
ifkeyexist=index(records,'nmol')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('nmol')+1:)
rewind(13)
read(13,*) nmol
rewind(13)
if (nmol.lt.1 .or. nmol.gt.mnmol) then
write(6,222) 'ERROR: number of AGs cant be less than 1', &
' or larger than ',mnmol
stop
elseif (atomdiatom .and. nmol.ne.1) then
write(6,108) 'ERROR: atom-diatom calculation cant', &
' be used with multi AGs (nmol>1)'
stop
elseif (reactcol .and. nmol.ne.2) then
write(6,108) 'ERROR: currently reactive collision can ', &
'only be used with two AGs (nmol=2)'
stop
else
continue
endif
! Move to the next line after reading nmol
read(12,'(a)',end=9999) records
endif
!
! Data for each AG, loop over nmol AGs, had to be formated at present
! version of ants
!
nat = 0
do im=1,nmol
! Natom, initx, initp
! Input for natom, initx, and initp should be in one line
! User can choose not to provide initx and initp, but must for natom
ifkeyexist=index(records,'natom')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('natom')+1:)
rewind(13)
read(13,*) natom(im)
rewind(13)
if (natom(im).lt.1 .or. natom(im).gt.mnat) then
write(6,222)'ERROR: number of atoms cant be less than', &
' 1 or larger than ',mnat
stop
endif
endif
! read in charge information for each atoms or not
if (index(records,'readinchg') .GT. 0) ireadchg=.true.
! initx information
ifkeyexist=index(records,'initx')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('initx')+1:)
rewind(13)
read(13,*) initx(im)
rewind(13)
if (initx(im).lt.0 .or. initx(im).gt.5) then
write(6,333) 'ERROR: initial conditions flag = ', &
initx(im),' is not supported'
stop
endif
if (initx(im).eq.4) then
if (natom(im).ne.3.and.nmol.ne.1) then
write(6,108) 'INITx=3 must be used with a triatomic', &
' molecule, and only a single AG'
write(6,*) '.ie. natom=3 and nmol=1'
stop
endif
endif
if ((initx(im).eq.3 .or. initx(im).eq.5) .and. nmol.ne.1) then
write(6,444) 'ERROR: randomly generate box only for ', &
'single AG run'
stop
endif
endif
ifkeyexist=index(records,'initp')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('initp')+1:)
rewind(13)
read(13,*) initp(im)
rewind(13)
if (initp(im).lt.-1 .or. initp(im).gt.1) then
write(6,333) 'ERROR: initial conditions flag = ', &
initp(im),' is not supported'
stop
endif
endif
ifkeyexist=index(records,'temp0im')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('temp0im')+1:)
rewind(13)
read(13,*) temp0im(im)
rewind(13)
endif
! ifkeyexist=index(records,'vibene')
! if (ifkeyexist .GT. 0) then
! write(13,*) records(ifkeyexist+len('vibene')+1:)
! rewind(13)
! read(13,*) vibene(im)
! rewind(13)
! vibene(im) = vibene(im)/autoev
! endif
ifkeyexist=index(records,'verte')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('verte')+1:)
rewind(13)
read(13,*) verte(im)
rewind(13)
verte(im) = verte(im)/autoev
endif
ifkeyexist=index(records,'bandwd')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('bandwd')+1:)
rewind(13)
read(13,*) bandwd(im)
rewind(13)
bandwd(im) = bandwd(im)/autoev
endif
! Initial relative translation energy between two AGs
! (unit: kT)
ifkeyexist=index(records,'ereltrans')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('ereltrans')+1:)
rewind(13)
read(13,*) etrans(im)
rewind(13)
etrans(im) = etrans(im)*kb*temp0im(im)
transelect = .true.
endif
! Initial vibrational state information
ifkeyexist=index(records,'vibselect')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('vibselect')+1:)
rewind(13)
read(13,*) vibsel(im)
rewind(13)
endif
! coordinate & momentum distribution of vibration
ifkeyexist=index(records,'vibdist')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('vibdist')+1:)
rewind(13)
read(13,*) ivibdist(im)
rewind(13)
endif
! Initial geometry is a minimum or a saddle point
ifkeyexist=index(records,'vibtype')
if (ifkeyexist .GT. 0) then
write(13,*) records(ifkeyexist+len('vibtype')+1:)
rewind(13)
read(13,*) ivibtype(im)
rewind(13)
endif
! Mode populations were implemented assuming that populations are obtained
! from a quasiclassical distribution. So vibdist must equal 0 when lmodpop
! is true.
if (lmodpop) ivibdist(im) = 0
! Set values for initx and initp according to the specification in the
! keyword sections
if (atomdiatom) then
initx(im) = 4
initp(im) = -1
endif
! If the user choose to do a normal mode analysis,
! initp should be set as -1
if ( vibsel(im).ne.0 .and. vibsel(im).ne.3) then
initp(im) = -1
endif
if ( vibsel(im).eq.3) then
itempzpe=.true.
endif
!
! iatom(im): the index of the first atom of AG im in all the atoms,
! starts from 1.
! nat: total number of atoms in a run, sum over natom(im)
iatom(im) = nat + 1
nat = nat+natom(im)
if (nat.gt.mnat) then
write(6,222) 'ERROR: the number of all atoms cant be', &
' larger than',mnat
stop
endif
! Since charges, multi, coord. are written in a formatted format, it is
! no need to use file 13 to transfer data
! Charge and multiplicity, in one line, move to next line after reading
read(12,*) icharge(im),imulti(im)
! Coordinates (if any) and other information related to coordinates
! INITx controls the initial coordinates
! read different info depending on the value of INITx(IM)
if (initx(im).eq.0) then
do i =1,natom(im)
ii=iatom(im)+i-1
! read in atomic symbol, mass, and coords
if (ireadchg) then
read(12,*) symbolm(i),charge(i),mmm(i), &
(xxm(k,i),k=1,3)
else
read(12,*) symbolm(i),mmm(i),(xxm(k,i),k=1,3)
endif
! convert from amu to au before calculating principal axes of rotation
mmm(i) = mmm(i) * amutoau
do k =1,3
! convert from angstrom to bohr
xxm(k,i) = xxm(k,i)/autoang
! put them to global storage
enddo
mm0(ii) = mmm(i)
symbol0(ii) = symbolm(i)
enddo
! Judge if this molecule is linear or not
natoms=natom(im)
call comgen(xxm,mmm,natoms,comxxm,mmofag(im))
! Get principal axes of rotation, get rotaxis by calling
! rotprin
call rotprin(xxm,mmm,natoms,momiprin,rotaxis, &
linearmol(im))
! Rot to the principal axes of rotation
! not for periodic condition
if (.not.periodic) then
call rottran(rotaxis,xxm,natoms)
endif
do i=1,natom(im)
ii=iatom(im)+i-1
do k =1,3
xx0(k,ii) = xxm(k,i)
enddo
enddo
if (linearmol(im).eq.0) then
nvibmods(im) = 0
nrotmods(im) = 0
elseif (linearmol(im).eq.1) then
nvibmods(im) = 3*natom(im) - 5
nrotmods(im) = 2
else
nvibmods(im) = 3*natom(im) - 6
nrotmods(im) = 3
endif
elseif (initx(im).eq.1 .or. initx(im).eq.2 &
.or. initx(im).eq.3) then
if (natom(im).le.1) then
write(6,108) 'ERROR: Only one atom, no need to', &
' generate coordinates by random method.'
stop
endif
if (initx(im).eq.1) read(12,*) rran(1,im)
if (initx(im).eq.2.) read(12,*) (rran(j,im),j=1,3)
! read in atomic symbol, and mass. User uses, e.g. 1 H 2 C to
! input, this is easy if AG contains a large number of atoms
! inattmp to store the number of one atom type
! j counts the atoms in this AG, sum of inattmp
j=0
do while(j.lt.natom(im))
read(12,*) inattmp,symbtmp,mmtmp
do k=iatom(im)+j,iatom(im)+j+inattmp-1
symbol0(k) = symbtmp
mm0(k) = mmtmp
enddo
j=j+inattmp
enddo
if (j.ne.natom(im)) then
write(6,108) 'ERROR: sum of atom numbers does not', &
' match the total number of atoms for this AG'
stop
endif
! random mol.
elseif (initx(im).eq.5) then
read(12,*) (rran(j,im),j=1,3)
read(12,*) nunit,natunit
natom(im)=nunit*natunit
do j=1,natunit
read(12,*) symbtmp,mmtmp,(xunit(k,j),k=1,3)
do k=1,3
xunit(k,j)=xunit(k,j)/autoang
enddo
do i=1,nunit
ii=(i-1)*natunit+j
symbol0(ii) = symbtmp
mm0(ii) = mmtmp
enddo
enddo
! Atom-diatom calculation
elseif (initx(im).eq.4) then
do i = iatom(im), iatom(im)+natom(im)-1
! read in atomic symbol, and mass
read(12,*) symbol0(i), mm0(i)
enddo
else
write(6,333) 'ERROR: initial conditions initx = ', &
initx(im),' is not supported'
! this endif ends information on initx
endif
! Move to next line after readin structures
read(12,'(a)',end=9999) records
!
! Read in vibration quantum numbers the user specify.
if (index(records,'vibstate').gt.0) then
vibsel(im) = 1
if (natom(im).eq.1) then
write(6,108) 'ERROR: there is no vibrational mode ', &
'for just one atom'
stop
endif
read(12,*,iostat=ierror) (qnmvib0(k,im), k=1,nvibmods(im))
if (ierror.ne.0) then
write(6,222) 'ERROR: reading vibrational quantum ', &
'numbers for AG ',im
stop
endif
! RVM
if (ivibtype(im).eq.1) then
write(6,*)"Initial geometry is a saddle point"
write(6,*)"Unbound mode will be assigned an energy of ", &
qnmvib0(nvibmods(im),im)," eV"
qnmvib0(nvibmods(im),im)=qnmvib0(nvibmods(im),im)/autoev
else
write(6,*)"Initial geometry is a local minimum"
endif
! Move to next line after readin vib modes
read(12,'(a)',end=9999) records
endif
! Read in the vibrational energy for each mode
if (index(records,'vibene').gt.0) then
if (natom(im).eq.1) then
write(6,108) 'ERROR: there is no vibrational mode ', &
'for just one atom '
stop
endif
if( vibsel(im).eq.4.or.vibsel(im).eq.6) then
read(12,*,iostat=ierror) tmp
vibene(:,im) = tmp
elseif(vibsel(im).eq.5.or.vibsel(im).eq.7) then
read(12,*,iostat=ierror) (vibene(k,im), k=1,nvibmods(im))
else
write(6,*)'ERROR: vibene should be used together ', &
'with vibselect = 4, 5, 6, or 7'
endif
vibene(:,im)=vibene(:,im)/autoev
if (ierror.ne.0) then
write(6,222)'ERROR: reading vibrational energy ', &
'for AG ',im
stop
endif
! Move to next line after readin vibene
read(12,'(a)',end=9999) records
endif
! Read in vibdist for each individual mode
if (index(records,'vibdistn').gt.0) then
if (natom(im).eq.1) then
write(6,108) 'ERROR: there is no vibrational mode ', &
'for just one atom '
stop
endif
if(ivibdist(im).ne.9) then
write(6,*) 'ERROR: to use vibdistn, the vibdist must',' be 9'
stop
endif
read(12,*,iostat=ierror) (vibdistn(k,im), k=1,nvibmods(im))
! Move to next line after readin vibene
read(12,'(a)',end=9999) records
endif
! Readin rotational quantum numbers if rotselect=1
! and only for linear molecules
! Force sample rotational energy based on temp(im)
if (index(records,'rottemp') .gt. 0 )then
irottemp = .true.
read(12,'(a)',end=9999) records
else
irottemp = .false.
endif
if (index(records,'rotfreeze') .gt. 0 )then
irotfreeze = .true.
read(12,'(a)',end=9999) records
else
irotfreeze = .false.
endif
if (index(records,'rotstate').gt.0 ) then
rotsel(im) = .true.
rotst(im) = .true.
if (linearmol(im).eq.1) then
read(12,*,iostat=ierror) qnmrot0(1,im)
if (ierror.ne.0) then
write(6,222) 'ERROR: reading rotational quantum ', &
'numbers for AG ',im
stop
endif
eroti(1,im)=dble(qnmrot0(1,im)*(qnmrot0(1,im)+1))/ &
( 2.d0*momiprin(3) )
erot(im) = eroti(1,im)
elseif (linearmol(im).eq.0) then
write(6,108) 'ERROR: there is no rotation for', &
' just one atom'
stop
else
write(6,108) 'ERROR: Cant deal with nonlinear ', &
'molecule for stateselect rotation'
stop
endif
! Move to next line after readin rot modes
read(12,'(a)',end=9999) records
endif
if (index(records,'rotenerg').gt.0 ) then
rotsel(im) = .true.
roteng(im) = .true.
if (linearmol(im).ne.0) then
read(12,*,iostat=ierror) (eroti(k,im), k=1,nrotmods(im))
if (ierror.ne.0) then
write(6,222) 'ERROR: reading rotational energies ', &
'for AG ',im
stop
endif
else
write(6,108) 'ERROR: there is no rotation for', &
' just one atom'
stop
endif
erot(im)=0.d0
do k=1,nrotmods(im)
eroti(k,im) = eroti(k,im)*kb*temp0im(im)
erot(im) = erot(im) + eroti(k,im)
enddo
! Move to next line after readin rot energies
read(12,'(a)',end=9999) records
endif
ireadchg = .false.
! this enddo ends the loop over AG data
enddo
! Move to next line
read(12,'(a)',end=9999) records
endif
!
! ---------------------------------------------------------
! Initial orientation
if (index(records,'$comxx') .GT. 0) then
! coordinates start from next line of $comxx
icomxx = .true.
do im=1,nmol
read(12,*) (comxx(k,im),k=1,3)
enddo
! Meet the end of this section, formated input, no need to find
! $end. first line to read $end second line to readin the line
! after $end
read(12,'(a)',end=9999) records
! Move to next line after reading $end
read(12,'(a)',end=9999) records
!
! ---------------------------------------------------------
! Initial translation momentum
endif
if (index(records,'$compp') .GT. 0) then
! coordinates start from next line of $compp
icompp = .true.
transelect = .true.
do im=1,nmol
read(12,*) (compp(k,im),k=1,3)
enddo
! Meet the end of this section, formated input, no need to find
! $end. first line to read $end second line to readin the line
! after $end
read(12,'(a)',end=9999) records
read(12,'(a)',end=9999) records
endif
! If a line contains only spaces (blank line), then move to next line
! How about a line composed of characters all useless to the program?
ifblank=0
do j=1,74
if (records(j:j).eq.' ') then
ifblank=ifblank+1
endif
enddo
if (ifblank.eq.74) then
read(12,'(a)',end=9999) records
endif
if (index(records,'thend0') .GT. 0) goto 9999
51 continue
9999 continue
close(12)
close(13)
!
! Check for bad combinations
if (reactcol) then
if (atomdiatom) then
write(6,108) 'ERROR: cant run reaction collision and ', &
'special atom-diatom at the same time'
stop
endif
if (termflag.eq.2) then
write(6,108) 'ERROR: TERMFLAG = 2 cant be used with', &
' react collision run'
stop
endif
if (tflag(1).ne.0) then
write(6,108) 'ERROR: tflag1 != 0 cant be used with', &
' react collision run'
stop
endif
if (termflag.eq.7) then
write(6,108) 'ERROR: Monitoring unimolecular fragmentation', &
' should not be used with react collision run'
stop
endif
endif
if (termflag.eq.6 .and. (.not.reactcol)) then
write(6,108) 'ERROR: TERMFLAG = 6 works only for', &
' reactive collision run'
stop
endif
if (termflag.eq.7 .and. nmol.ne.1) then
write(6,444) 'ERROR: Monitoring fragmentation ', &
'works only for 1 AG run (nmol=1)'
stop
endif
if (bdfmthr.gt.bdbrthr) then
write(6,*) 'ERROR: Bondbrthre must be larger than bondfmthre'
stop
endif
do im=1,nmol
if (vibsel(im).ne.0) then
if ((initx(im).ge.1 .and. initx(im).le.3) .or. initx(im).eq.5) then
write(6,108) 'ERROR: random structure cant do ', &
'normal mode analysis'
stop
endif
endif
if (vibsel(im).lt.4.and.vibene(1,im).gt.1d-10) then
write(6,*) 'ERROR: Fixed energy for each mode must use '
write(6,*) 'vibsel=4, 5, 6, or 7'
endif
enddo
if (lfixte) then
write(6,*) 'Choosing fixed initial total energy'
write(6,555) 'Warning: currently only works for initial ', &
'conditions prepared with normod analysis', &
' without other motions (trans. & rot.)'
lzerotrans = .true.
ifnoang = .true.
if (nmol.ne.1) then
write(6,*) 'Currently works only for single AG simulation'
stop
endif
if (deltate.gt.small) then
write(6,666) 'TE0 fixed to between ',(te0fixed-deltate)*autoev, &
' eV and ',(te0fixed+deltate),' eV'
else
write(6,*) 'TE0 fixed to ',(te0fixed-deltate)*autoev,' eV'
endif
endif
!
! ---------------------------------------------------------
! All printout information for the input goes here
! Potential
write(6,'(1x,74a1)') ('-',i=1,74)
write(6,*) 'Potential information'
write(6,'(1x,74a1)') ('-',i=1,74)
write(6,*) 'POTFLAG = ',potflag
call prepot
if (potflag.eq.0) then
write(6,*) 'Using the POTLIB MM-HO-1 interface'
if (nsurft.ne.1) then
write(6,108) 'POTLIB MM-HO-1 is for SINGLE surface only, ', &
'NSURFT should be 1'
stop
endif
elseif (potflag.eq.1) then
write(6,*) 'Using the POTLIB 3-2V interface'
if (nat.ne.3) then
write(6,*) 'POTLIB 3-2V works only for 3 atoms'
stop
endif
if (nsurft.ne.2) then
write(6,*)'NSURFT = ',nsurft,' is not allowed for POTFLAG = 1'
stop
endif
elseif (potflag.eq.2) then
write(6,*) 'Using the POTLIB MM-HE-1 interface'
if (nsurft.ne.1) then
write(6,108) 'POTLIB MM-HE-1 is for SINGLE surface only, ', &
'NSURFT should be 1'
stop
endif
elseif (potflag.eq.3) then
write(6,*)'Using the NH3 interface'
if (nat.ne.4) then
write(6,*) 'NH3 potential works only for 4 atoms'
stop
endif
if (nsurft.ne.2) then
write(6,*)'NSURFT = ',nsurft,' is not allowed for NH3 potential'
stop
endif
elseif (potflag.eq.4) then
write(6,*) 'Using the 4-XS interface'
if (nsurft.ne.1) then
write(6,108) 'POTLIB 4-XS is for SINGLE surface only, ', &
'NSURFT should be 1'
stop
endif
if (nat.gt.4) then
write(6,*) 'POTLIB 4-XS is 4 body potential, but nat= ',nat
stop
endif
elseif (potflag.eq.5) then
write(6,*) 'Using the POTLIB MM-HE-1 interface'
if (nsurft.ne.1) then
write(6,108) 'POTLIB MM-HE-1 is for SINGLE surface only, ', &
'NSURFT should be 1'
stop
endif
elseif (potflag.eq.6) then
write(6,*) 'Using the HBr interface'
elseif (potflag.eq.7) then
write(6,*) 'Using the CH2BrCl interface'
elseif (potflag.eq.8) then
write(6,*) 'Using the HN2 interface'
elseif (potflag.eq.9) then
write(6,*) 'Using the HX2 interface'
elseif (potflag.eq.10) then
write(6,*) 'Using the Phenol interface'
elseif (potflag.eq.11) then
write(6,*) 'Using the N4 interface'
elseif (potflag.eq.12) then
write(6,*) 'Using the HOCO interface'
elseif (potflag.eq.-1) then
if(qcpack .eq. 1) then
write(6,*) ' Direct dynamics with Gaussian package'
elseif(qcpack .eq. 2) then
write(6,*) ' Direct dynamics with Molpro package '
else
write(6,*) ' Direct dynamics with unsupported QC package '
stop
endif
elseif (potflag.eq.-2) then
write(6,*) ' Direct dynamics with MOPAC-mn'
else
write(6,*) 'POTFLAG = ',potflag,' is not supported'
stop
endif
! Adiabatic or non-adiabatic
write(6,*) 'REPFLAG= ',repflag
if (repflag.eq.0) then
write(6,*) 'The adiabatic representation will be used.'
else
write(6,*) 'The diabatic representation will be used.'
endif
! Integrator information
write(6,'(1x,74a1)') ('-',i=1,74)
write(6,*) 'Integrator information'
write(6,'(1x,74a1)') ('-',i=1,74)
write(6,*) 'INTFLAG = ',intflag
if (intflag.eq.0) then
write(6,*) 'Using the Bulirsch-Stoer integrator'
write(6,*) 'with initial time step ',hstep0*autofs,' fs'
write(6,*) 'and with integrator tolerance ',eps
elseif (intflag.eq.1) then
write(6,*)'Using the 4th-order Runge-Kutta integrator'
write(6,*) 'with time step ',hstep0*autofs,' fs'
elseif (intflag.eq.2) then
write(6,*)'Using the simplest Verlet integrator'
write(6,*) 'with time step ',hstep0*autofs,' fs'
elseif (intflag.eq.3) then
write(6,*)'Using the velocity Verlet integrator'
write(6,*) 'with time step ',hstep0*autofs,' fs'
elseif (intflag.eq.4) then
write(6,*)'Using the Beeman integrator'
write(6,*) 'with time step ',hstep0*autofs,' fs'
elseif (intflag.eq.5) then
write(6,*)'Using the Liouville approach of velocity Verlet', &
' integrator'
write(6,*) 'with time step ',hstep0*autofs,' fs'
elseif (intflag.eq.6) then
write(6,*)'Using the Bulirsch-Stoer integrator with Hack''s ', &
'method'
write(6,*) 'with initial time step ',hstep0*autofs,' fs'
write(6,*) 'and with integrator tolerance ',eps
if (intflag.eq.6.and.nsurft.eq.1) then
write(6,*) 'BS integrator with Hack''s method is only', &
' meaningful for multisurface calculations'
write(6,*) 'Choose bulstoint in the input instead'
stop
endif
if (intflag.eq.6.and.methflag.eq.2) then
write(6,*) 'BS integrator with Hack''s method is not', &
' meaningful for SE calculations'
write(6,*) 'Choose bulstoint in the input instead'
stop
endif
else
write(6,*)'INTFLAG = ',intflag,' is not supported'
stop
endif
if ( (itherm.ne.0 .or. ieqtherm.ne.0) .and. iadjtemp.eq.3 .and. &
intflag.ge.2 .and. intflag.le.4) then
write(6,444) 'Nose-Hoover thermostat cant be used with simple-', &
', velocity- and Beeman- Verlet integrator'
write(6,*)'Use Liouville approach to Velocity Verlet'
stop
endif
!
! Some initial conditions related with thermostat
if (reactcol .or.atomdiatom) ifnoang=.false.
! CoM motion and angular momentum for anderson thermostat (iadjtemp=2)
! are not conserved. So do not zero them.
if (itherm.ne.0 .and. iadjtemp.eq.2) then
! lzerotrans = .false.
! ifnoang = .false.
! ifeqwithrot=.true.
if (lzerotrans) then
write(6,444) 'WARNING: You chose to remove CoM in Anderson ', &
'thermostat, which is not conserved.'
endif
if (ifnoang) then
write(6,444) 'WARNING: You chose to keep angular momentrum ', &
'in Anderson thermostat, which is not conserved.'
endif
endif
! Total number of freedom for the whole run
! Periodic running information
write(6,*)
write(6,'(1x,74a1)') ('-',i=1,74)
write(6,*) 'Periodic and cell information'
write(6,'(1x,74a1)') ('-',i=1,74)
write(6,*) 'Doing a ',nmol,' AG calculation'
if (icell.eq.0) then
write(6,*) 'Simulation for an isolated AG in vacuum'
! periodic condition
periodic = .false.
! NVT can not be used
if (itherm.eq.2) then
write(6,444) 'ERROR: Cannot run isolated cluster in vacuum', &
' with NVT ensemble'
stop
endif
else
if (icell.eq.1) then
write(6,777) 'The simulation is in a cubic cell with a = ', &
cell(1),' A'
cell(2)=cell(1)
cell(3)=cell(1)
elseif (icell.eq.2) then
write(6,888)'The simulation is in a cuboid cell with a b c = ', &
(cell(j),j=1,3),' A'
elseif (icell.eq.3) then
write(6,999) 'The simulation is in an arbitrary cell with', &
' a b c = ',(cell(j),j=1,3),' A'
write(6,888) 'and with alpha, beta, gamma = ',&
(cell(j),j=4,6),' Degree'
endif
! Convert the cell parameters to cell vectors, with
! a (cell(1)) parallel to x axis. cvec(3,3) is a vector matrix,
! i.e. the first column is the a vector, second b vector
do j=1,3
cell(j) = cell(j)/autoang
enddo
if (icell.eq.1 .or. icell.eq.2) then
do j=4,6
cell(j) = 90.0d0
enddo
endif
do j=4,6
if (cell(j).gt.180.d0) then
write(6,*) 'ERROR: Angles cannot be greater than 180 degree'
write(6,*) 'Angles: ',(cell(i),i=4,6)
stop
endif
cell(j) = pi*cell(j)/180.d0
enddo
cvec(1,1) = cell(1)
cvec(2,1) = 0.0d0
cvec(3,1) = 0.0d0
cvec(1,2) = cell(2)*dcos(cell(6))
cvec(2,2) = cell(2)*dsin(cell(6))
cvec(3,2) = 0.0d0
cvec(1,3) = cell(3)*dcos(cell(5))
! Following equations are incorrect
cvec(2,3) = (cell(2)*cell(3)*dcos(cell(4))-cvec(1,2)*cvec(1,3))/cvec(2,2)
cvec(3,3) = dsqrt(cell(3)**2-cvec(1,3)**2-cvec(2,3)**2)
! Volume of the cell det(cvec) = a(bxc)
vol0 = detmtbt(cvec)
call invmtbt(cvec,invcvec)
vol=vol0
write(6,777) 'The initial volume of the cell = ',&
vol0*autoang**3,' A^3'
endif
write(6,*)
! Random seed
write(6,'(1x,74a1)') ('-',i=1,74)
write(6,*) 'Random number generator information'
write(6,'(1x,74a1)') ('-',i=1,74)
! RMP14
call seed(ranseed)
! if (ranseed.eq.0) then
!C If user not to provide ranseed, the program will determine
!C one automatically by calling make_sprng_seed()
! ranseed = make_sprng_seed()
! endif
!
! write(6,*) 'Random number seed is ',ranseed
! write(6,*)
! Surface information
write(6,'(1x,74a1)') ('-',i=1,74)
write(6,*) 'Electronic surfaces'
write(6,'(1x,74a1)') ('-',i=1,74)
write(6,*) 'Initial electronic suface = ',nsurf0
if (nsurfi.eq.0) then
nsurfi=nsurf0
repflgi=repflag
else
write(6,222) 'Initial condition prepared from electronic ', &
' surface ',nsurfi
if (repflgi.eq.-1) repflgi=repflag
endif
if ((nsurfi.ne.nsurf0 .or. repflgi.ne.repflgi) &
.and. deltate.lt.small) then
write(6,444) 'WARNING: Frank-Condon effect cannot fix the total', &
' initial energy exactly'
endif
write(6,*) 'Number of electronic surfaces = ',nsurft
if (methflag.eq.0) then
write(6,*) 'METHFLAG = 0: Single surface propagation'
elseif (methflag.eq.1) then
write(6,108) 'METHFLAG = 1: Electronically nonadiabatic',&
' calculation with'
write(6,*) 'Tully''s fewest switches method'
elseif (methflag.eq.2) then
write(6,108) 'METHFLAG = 2: Electronically nonadiabatic', &
' calculation with'
write(6,*) 'Semiclassical Ehrenfest method'
elseif (methflag.eq.3) then
write(6,108) 'METHFLAG = 3: Electronically nonadiabatic', &
' calculation with'
write(6,*) 'Self consistent decay of mixing method'
elseif (methflag.eq.4) then
write(6,108) 'METHFLAG = 4: Electronically nonadiabatic', &
' calculation with'
write(6,108) 'Coherent switching with decay of mixing method'
elseif (methflag.eq.5) then
write(6,108) 'METHFLAG = 5: Electronically nonadiabatic', &
' calculation with'
write(6,108) 'Fewest switches with time uncertainty method'
else
write(6,*) 'METHFLAG = ',methflag,' is not allowed'
stop
endif
! RVM Information for Stochastic Decoherence and TRAPZ
if (stodecoflag.eq.1) then
write(6,*) 'STODECOFLAG = 1: Using Stochastic Decoherence'
elseif (stodecoflag.lt.0.or.stodecoflag.gt.1) then
write(6,*) 'stodecoflag= ',stodecoflag,' is not allowed'
stop
endif
if (ntrapz.eq.0) then
write(6,108) 'Using traditional Hessian for initial normal', &
' mode analysis'
write(6,*) 'TRAPZ is not used for the dynamics'
elseif (ntrapz.eq.1) then
write(6,108) 'Using Projected Hessian for initial normal', &
' mode analysis'
write(6,*) 'TRAPZ is not used for the dynamics'
elseif(ntrapz.eq.2) then
write(6,108) 'Using traditional Hessian for initial normal', &
' mode analysis'
write(6,*) 'Using TRAPZ for the dynamics'
elseif(ntrapz.eq.3) then
write(6,108) 'Using Projected Hessian for initial normal', &
' mode analysis'
write(6,*) 'Using TRAPZ for the dynamics'
elseif(ntrapz.le.0.or.ntrapz.gt.3) then
write(6,*) 'ntrapz= ',ntrapz,' is not allowed'
stop
endif
! RVM Information about version of TRAPZ used
if (nvers.eq.0) then
write(6,*) 'Using general TRAPZ method'
elseif (nvers.eq.1) then
write(6,108) 'Using the mTRAPZ method'
elseif(nvers.eq.2) then
write(6,108) 'Using the mTRAPZ* method'
elseif(nvers.lt.0.or.nvers.gt.2) then
write(6,*) 'nvers= ',nvers,' is not allowed for TRAPZ'
stop
endif
do i=1,nsurft*2
cre(i) = 0.d0
cim(i) = 0.d0
enddo
cre(nsurf0) = 1.d0
write(6,*)
!
! Find out the corresponding index of the atom in the periodic table
! Model atoms types are added by J. Zheng May 2013
nmodel = 0
do i=1,nat
find_real = .false.
do j=1,111
if ( symbol0(i).eq.atomsymbols(j)) then
indatom(i)=j
find_real=.true.
endif
enddo
if(.not.find_real) then
nmodel = nmodel + 1
indatom(i)= nmodel
endif
enddo
!
! Count the number of atoms which are the same for different AG
do im=1,nmol
do j=1,121
nattype(im,j)=0
enddo
enddo
do im=1,nmol
do i=1,natom(im)
ii=indatom(iatom(im)+i-1)
nattype(im,ii)=nattype(im,ii)+1
enddo
enddo
!
! Atomic group data (initx, initp, coordinates)
write(6,'(1x,74a1)') ('-',i=1,74)
write(6,*) 'Atom group (AG) specification'
write(6,'(1x,74a1)') ('-',i=1,74)
! Atom-Diatom information
if (atomdiatom .or. initx(1).eq.4) then
write(6,*) 'Doing a special atom-diatom simulation'
write(6,910) ' # symb','mass (amu)'
! atom-diatom has only one AG, so always begin with 1
mmofag(im)=0.d0
do i = 1, natom(1)
write(6,920) i,symbol0(i),mm0(i)
! Convert from amu to au
mm0(i)=mm0(i)*amutoau
mmofag(im)=mmofag(im)+mm0(i)
enddo
if (ademod .eq. 1) then
write(6,*) 'Total energy = ',escatad,' eV'
else
write(6,*) 'Collision energy = ',ecolad,' eV'
endif
write(6,*) 'Initial vib. quantum number v = ',vvad
write(6,*) 'Initial rot. quantum number j = ',jjad
write(6,*) 'Initial distance = ',rrad0
if (arrad.eq.1) then
write(6,*) 'Initial arrangement = AB + C'
elseif (arrad.eq.2) then
write(6,*) 'Initial arrangement = BC + A'
elseif (arrad.eq.3) then
write(6,*) 'Initial arrangement = AC + B'
else
write(6,*) 'ERROR: arrad = ',arrad,' is not supported'
stop
endif
ecolad=ecolad/autoev
escatad=escatad/autoev
rrad0=rrad0/autoang
! None atom-diatom AG information
else
! start loop over AG
do im=1,nmol
write(6,'(1x,74a1)') ('-',i=1,74)
write(6,*) 'AG ',im
write(6,'(1x,74a1)') ('-',i=1,74)
write(6,200) 'This AG has ',natom(im),' atoms (', &
'atoms ',iatom(im),' to ',iatom(im)+natom(im)-1,')'
write(6,201) 'Initial conditions prepared using INITx = ',&
initx(im),' and INITp = ',initp(im)
! information on coordinates and initx
! Print out initial coordinates and save them to global storage
if (initx(im).eq.0) then
write(6,910)' # symb','mass (amu)','x (A)','y (A)','z (A)'
! do i = iatom(im), iatom(im)+natom(im)-1
do i = 1,natom(im)
ii=iatom(im)+i-1
write(6,920) i,symbol0(ii),mm0(ii)/amutoau, &
(xx0(k,ii)*autoang,k=1,3)
enddo
elseif ((initx(im).ge.1 .and. initx(im).le.3) .or. &
initx(im).eq.5) then
maxtmp=0.d0
mintmp=200.d0
do i=iatom(im),iatom(im)+natom(im)-1
do j=i,iatom(im)+natom(im)-1
maxtmp = max(maxtmp,rijco(indatom(i),indatom(j)))
mintmp = min(mintmp,rijco(indatom(i),indatom(j)))
enddo
enddo
if (initx(im).eq.1) then
write(6,930) 'Generate random spherical clusters', &
' of radius ',rran(1,im),' A'
! Judge if this radius is reasonable
radmin(im)=mintmp*dble(natom(im))**(1.d0/3.d0)/2.d0
radmax(im)=maxtmp*dble(natom(im))**(1.d0/3.d0)/2.d0
if (radmin(im)*autoang.gt.rran(1,im) ) then
write(6,108) 'WARNING: The radius given maybe too ', &
'small for the program to get a random structure'
endif
if (2.0d0*radmax(im)*autoang.lt.rran(1,im) ) then
write(6,108) 'WARNING: The radius given maybe too ', &
'large to get a spheric cluster'
endif
elseif (initx(im).eq.2) then
write(6,940) 'Generate random atoms in cubic with a,b,c ', &
(rran(i,im),i=1,3),' A'
k=int((rran(1,im)*rran(2,im)*rran(3,im))/ &
(pi*((mintmp*autoang/2.d0)**3.d0)*4.d0/3.d0))
if (k.lt.natom(im)) then
write(6,108) 'WARNING: volumn of the cuboid specified', &
' maybe too small'
endif
elseif (initx(im).eq.3) then
write(6,108) 'Generate random clusters with arbitrary ', &
'shape using default bond threshold'
elseif (initx(im).eq.5) then
write(6,940) 'Generate random molecules in a box of: ', &
(rran(i,im),i=1,3),' A'
endif
do i=1,3
rran(i,im) = rran(i,im)/autoang
enddo
write(6,910) ' # symb','mass (amu)'
! Calculate total weight of each AG
mmofag(im)=0.d0
do i = iatom(im), iatom(im)+natom(im)-1
write(6,920) i,symbol0(i),mm0(i)
! Convert from amu to au
mm0(i)=mm0(i)*amutoau
mmofag(im)=mmofag(im)+mm0(i)
enddo
else
write(6,*) 'ERROR: INITx = ',initx(im),' is not supported'
stop
endif
! information on initp
if (initp(im).eq.0) then
write(6,*) 'INITp = 0: Initial momenta set to zero'
write(6,*)
elseif (initp(im).eq.1) then
write(6,930) 'Target temperature for random thermal ', &
'distribution is ',temp0im(im),' K'
write(6,*)
elseif (initp(im).eq.-1) then
write(6,*) 'INITp =-1: Initial momenta determined by INITx'
write(6,*)
else
write(6,*) 'ERROR: INITp = ',initp(im),' is not supported'
stop
endif
if (vibsel(im).eq.0) then
write(6,*) 'Vibrational state determined otherwise'
elseif (vibsel(im).eq.1) then
write(6,*) 'User provide vibrational states for each AG'
write(6,*) 'Initial vibrational quantum numbers:'
! RVM
if(ivibtype(im).eq.1) then
write(6,'(1x,15f4.1/)') (qnmvib0(j,im),j=1,nvibmods(im)-1)
else
write(6,'(1x,15f4.1/)') (qnmvib0(j,im),j=1,nvibmods(im))
endif
elseif (vibsel(im).eq.2) then
write(6,*) 'Vibrational state randomly selected'
elseif (vibsel(im).eq.3) then
write(6,*) 'Do frequency analysis'
elseif (vibsel(im).eq.4) then
write(6,*) 'Using the same amount of energy for each ' &
,'vibrational mode'
elseif (vibsel(im).eq.5) then
write(6,*) 'Using the individual amount of energy for each ' &
,'vibrational mode'
elseif (vibsel(im).eq.6) then
write(6,*) 'Each vibrational mode has energy as ' &
,' min{(0.5*h*nu_m, vibene_1} '
elseif (vibsel(im).eq.7) then
write(6,*) 'Each vibrational mode has energy as ' &
, ' min{(0.5*h*nu_m, vibene_m} '
else
write(6,*) 'ERROR: VIBSELECT = ',vibsel(im),' not supported'
stop
endif
if (vibsel(im).ne.0) then
if (ivibdist(im).eq.0) then
write(6,444) 'Using classical HO distribution to prepare ', &
'initial coordinates and momentum'
elseif (ivibdist(im).eq.1) then
write(6,444) 'Using HO distribution to prepare ', &
'initial coordinates and momentum'
elseif (ivibdist(im).eq.2) then
write(6,444) 'Using Wigner distribution to prepare ',&
'initial coordinates and momentum'
elseif (ivibdist(im).eq.9) then
write(6,444) 'Each mode using different distributions to ',&
'prepare initial coordinates and momentum'
j = nvibmods(im)
if(ivibtype(im).eq.1) j = nvibmods(im)-1
do i = 1, j
if(vibdistn(i,im).lt.0) then
write(6,*) 'The distribution method for mode', i, &
' is not given. Check vibdistn keyword...'
endif
enddo
else
write(6,*) 'ERROR: VIBDIST = ',ivibdist(im),' not supported'
stop
endif
endif
if (rotsel(im)) then
write(6,*) 'User provide rotational state (energy)'
else
write(6,*) 'Rotational state (energy) randomly decided'
endif
if (rotst(im)) then
write(6,*) 'Initial rotational quantum numbers:'
write(6,'(1x,3i4)') (qnmrot0(j,im),j=1,nrotmods(im))
endif
if (roteng(im)) then
write(6,*) 'Initial rotational energies (eV) of AG ',im
write(6,'(1x,3f12.5)') (eroti(j,im)*autoev,j=1,nrotmods(im))
endif
enddo ! End of loop over AG
ENDIF ! End of AG section
!
! Orientation of AGs
write(6,'(1x,74a1)') ('-',i=1,74)
write(6,*) 'AG initial orientation and translation motion'
write(6,'(1x,74a1)') ('-',i=1,74)
if (nmol.gt.1) then
write(6,*) 'Multi AG simulation'
do im=1,nmol
write(6,*) 'For AG ',im
if (icomxx) then
write(6,940) 'Center of mass coords (x,y,z) (A) = ', &
(comxx(k,im),k=1,3)
else
write(6,940) 'Center of mass coords determined otherwise'
endif
if (icompp) then
write(6,940) 'Center of mass motion (x,y,z) (au) = ', &
(compp(k,im),k=1,3)
! Calculate translational energy
etrans(im)=0.d0
do k=1,3
etrans(im)=etrans(im) + compp(k,im)**2/(2.d0*mmofag(im))
enddo
else
write(6,940) 'Center of mass motion determined otherwise'
endif
do k=1,3
comxx(k,im) = comxx(k,im)/autoang
enddo
enddo
do im=1,nmol
write(6,960) 'Translational energy for AG ',im, &
' is ',etrans(im)*autoev,' eV'
enddo
else
write(6,*) 'Single AG simulation. No orientation'
! Single AG calculation
do k=1,3
comxx(k,1) = 0.d0
compp(k,1) = 0.d0
enddo
endif
!
! Reaction collision information
if (reactcol) then
write(6,'(1x,74a1)') ('-',i=1,74)
write(6,*) 'Reaction collision information'
write(6,'(1x,74a1)') ('-',i=1,74)
if (.not.stateselect) then
write(6,*) 'Do an equilibrium run before collision'
write(6,*) 'The probability to pick a configuration is',pickthr
else
write(6,*) 'Do a state select reaction collision'
write(6,*) 'States specified by the user will be fixed'
endif
if (ipicktrj.ne.0) then
if (demin.eq.-1.d50 .and. demax.eq.1.d50) then
write(6,*) 'No constraint on the energy of the points picked'
elseif (demin.ne.-1.d50) then
write(6,*) 'Pick those with energy >= than ',demin,' eV'
demin=demin/autoev
elseif (demax.ne.1.d50) then
write(6,*) 'Pick those with energy < than ',demax,' eV'
demax=demax/autoev
endif
endif
if (transelect) then
write(6,*) 'User specify relative translation energy'
else
write(6,*) 'Randomly determined relative translation energy'
endif
write(6,*) 'Information on initial distance between 2 AGs'
if (.not.initr0col) then
write(6,108) 'The initial distance between 2 AGs will', &
' be calcualted based on covalent ', &
'radii of the atoms in the AG.'
write(6,108) 'Attention: atomic radii is just a very', &
' rough estimation of the size of an atom'
! Call a subroutine to provide r0collision here. If random cluster, only
! calculate for sphere cluster. It is unlikely that one would like to
! run rxcollision using none sphere random clusters
r0collision=0.d0
do im=1,nmol
natoms=natom(im)
if (natoms.eq.1) then
radmax(im) = 2.2d0*atmradco(indatom(iatom(im)))/autoang
else
if (initx(im).eq.0) then
do i=1,natoms
ii=iatom(im)+i-1
do j=1,3
xxm(j,i) = xx0(j,ii)
enddo
indmm(i)=indatom(ii)
enddo
call rijmatr
! Find out the largest distance(bohr) in rijmatrix
maxtmp=0.d0
do i=1,natoms
do j=i+1,natoms
maxtmp=dmax1(maxtmp,rij(i,j))
if (rij(i,j).eq.maxtmp) then
i1=indmm(i)
i2=indmm(j)
endif
enddo
enddo
radmax(im)=2.2d0*(maxtmp+rijco(i1,i2))/2.d0
elseif (initx(im).eq.1) then
! provided in printing $data out
radmax(im)=rran(1,im)
else
write(6,108) 'ERROR: Program will only decide ', &
'r0collision for INITx = 0 or 1'
stop
endif
endif
if (idebug) then
write(6,*) 'radmax(im)',radmax(im)*autoang
endif
r0collision=r0collision+radmax(im)
enddo
else
write(6,*) 'User provide initial collision distance'
endif
write(6,950) 'Initial distance R0Collision = ', &
r0collision*autoang,' A'
if (b_fix) then
write(6,950) 'Initial impact parameter b fixed at ', &
b0impact*autoang,' A'
else
write(6,*) 'Integrate over all impact parameter'
endif
if (r0collision.le.b0impact) then
write(6,*) 'ERROR: R0Collision shorter than b0impact '
stop
endif
! endif of reaction collision
ENDIF
!
! Trajectory information
write(6,'(1x,74a1)') ('-',i=1,74)
write(6,*) 'Trajectory control flags'
write(6,'(1x,74a1)') ('-',i=1,74)
write(6,*) 'Running ',ntraj,' trajectories'
if (tflag(1).eq.0) then
write(6,*) 'No special options for tflag1'
elseif (tflag(1).eq.1) then
write(6,108)'TFLAG(1)= 1: Momenta set to zero at every ', &
'step, i.e., steepest descent minimization'
itherm=0
elseif (tflag(1).eq.2) then
write(6,*) 'TFLAG(1) = 2: Temperature ramping'
write(6,333) 'The momentum will be ramped at step ', &
nrampstep,' by '
write(6,202) 'a factor of RAMPFACT = ',rampfact,' and ', &
nramp,' times'
write(6,333) 'The momentum will be ramped for ', &
n_ramp,' steps and '
write(6,333) 'followed by an equilibrium run of ', &
nequil,' steps each time'
tmp=dble(nrampstep+nramp*(n_ramp+nequil))*hstep0
if (t_stime.lt.tmp) then
write(6,203)'Terminating time too short, it must be longer', &
' than ',tmp*autofs,' fs'
stop
endif
elseif (tflag(1).eq.3.or.tflag(1).eq.4.or.tflag(1).eq.5) then
if (tflag(1).eq.3) then
write(6,*) 'TFLAG(1) = 3: simulated heating/cooling'
elseif (tflag(1).eq.4) then
write(6,444) 'TFLAG(1) = 4: simulated heating/cooling ', &
'followed by additional geometry optimization using BFGS'
elseif (tflag(1).eq.5) then
write(6,204) 'TFLAG(1) = 4: simulated heating/cooling ', &
'do geometry optimization using BFGS at steps with ', &
'a frequency of ',pickthr
endif
if (tflag(1).eq.4 .or. tflag(1).eq.5) then
write(6,108) 'Geometry optimization will terminate ', &
'when the RMS of the gradient'
write(6,*) ' < ',t_gradmag*autoev/autoang,' eV/A'
write(6,108) 'and when the maximum of the gradient'
write(6,*) ' < ',t_gradmax*autoev/autoang,' eV/A'
endif
tmp = rampfact/(dble(n_ramp+nequil)*hstep0*autofs/1000.d0)
if (rampfact.gt.0.d0) then
write(6,*) 'The heating rate is ',tmp,' K/ps '
elseif (rampfact.lt.0.d0) then
write(6,*) 'The cooling rate is ',tmp,' K/ps '
endif
write(6,777) 'The time interval for the heating/cooling is ', &
dble(n_ramp+nequil)*hstep0*autofs,' fs '
write(6,333) 'Properties are averaged over ', &
nequil,' steps for each heating/cooling step'
tmp = temp0 + dble(nramp)*rampfact
write(6,*) 'The ending temperature would be ',tmp,' K'
if (tmp.lt.0.d0) then
write(6,205)'NRAMP is too large, it must be smaller', &
' than ',-int(temp0/rampfact),' times'
stop
endif
tmp = dble(nrampstep+nramp*(n_ramp+nequil))*hstep0
if (t_stime.lt.tmp) then
write(6,203)'T_STIME is too short, it must be longer', &
' than ',tmp*autofs,' fs'
stop
endif
else
write(6,*) 'ERROR: tflag1 = ',tflag(1),' not supported'
stop
endif
! Set nreinit to a very large number if do not zero CoM, CoM motion
! and angular momentum
if ((.not.lzerocom).and.(.not.lzerotrans).and.(.not.ifnoang))&
nreinit = 1000000000
if (lzerotrans) then
write(6,*) 'Remove CoM motion every step'
lzerocom = .true.
endif
if (lzerocom) then
write(6,*) 'Move the CoM back to origin every step'
endif
if (nmol.eq.1 .and. ifnoang) ifreinitrot=.false.
if ((.not.reactcol).and.(initx(1).ne.4)) &
write(333,*) 'Remove overall rotational momentum every ',nreinit, ' steps'
if (tflag(2).eq.1) then
write(6,*) 'TFLAG(2) = 1: Restart trajectories'
write(6,*) 'Restarting trajectories: ',(trajlist(k),k=1,ntraj)
endif
if (lzerotrans) write(6,444) 'Translation momentum will be ', &
'removed before entering DRIVE.'
if (ifnoang) write(6,444) 'Overall CoM angular momentum will be ', &
'removed before entering DRIVE.'
if (ifeqwithrot) write(6,444) 'Overall CoM angular momentum is ', &
'kept in equilibrium run.'
if (ifreinitrot)write(6,444) 'Overall CoM angular momentum will ', &
'reset after equilibrium run.'
write(6,*)
!
! Termination condition
write(6,'(1x,74a1)') ('-',i=1,74)
write(6,*) 'Termination condition information'
write(6,'(1x,74a1)') ('-',i=1,74)
if (termflag.eq.0) then
write(6,108)'TERMFLAG = 0: Fixed number of steps per ', &
'trajectory'
elseif (termflag.eq.1) then
write(6,777) 'TERMFLAG = 1: Trajectories will run for ', &
t_stime*autofs,' fs'
t_nstep=1000000000
elseif (termflag.eq.2) then
write(6,108) 'TERMFLAG = 2: Trajectories will terminate ', &
'when the RMS of the gradient'
write(6,*) ' is less than ',t_gradmag*autoev/autoang,' eV/A'
write(6,108) 'and when the maximum of the gradient'
write(6,*) ' < ',t_gradmax*autoev/autoang,' eV/A'
elseif (termflag.eq.3) then
write(6,*) 'Outcome label Bond Distance (A)'
if (nbondform.gt.0) then
write(6,108) 'TERMFLAG = 3: Bond forming termination ', &
'conditions'
do j=1,nbondform
write(6,900) j,tiform(j,1),tiform(j,2),trmin(j)
trmin(j)=trmin(j)/autoang
enddo
endif
if (nbondbreak.gt.0) then
write(6,108) 'TERMFLAG = 3: Bond breaking termination ', &
'conditions'
do j=1,nbondbreak
write(6,900) j,tibreak(j,1),tibreak(j,2),trmax(j)
trmax(j)=trmax(j)/autoang
enddo
endif
elseif (termflag.eq.4) then
write(6,108) 'TERMFLAG = 4: Monitor fragments'
if (agdefrag) write(6,*) 'Monitor AG fragmentation'
if (agmerge) then
write(6,*) 'Monitor AG merge'
write(6,777) 'Total energy less than ',determ, &
' eV is deemed as merge'
determ=determ/autoev
endif
if (attrans) write(6,*) 'Monitor atom/group transfer/exchange'
elseif (termflag.eq.5) then
write(6,108) 'TERMFLAG = 5: Monitor any bond breaking or', &
' forming'
elseif (termflag.eq.6) then
write(6,108) 'TERMFLAG = 6: Average on sticking time '
elseif (termflag.eq.7) then
write(6,*) 'Monitor unimolecular fragmentation probability'
t_nstep=1000000000
write(6,*)'Maximum simulation time is ',t_stime*autofs,' fs'
if (n_frags.eq.0) then
write(6,*) 'Number of products is unlimited'
elseif (n_frags.ge.2) then
write(6,*) 'Monitor ',n_frags,' product fragmentation'
else
write(6,*) 'ERROR: n_frags cant be 1 or negative'
stop
endif
if (initx(im).eq.0) then
do i=1,natom(1)
do j=1,3
xxm(j,i) = xx0(j,i)
enddo
indmm(i)=indatom(i)
enddo
call rijmatr
call frag
if (nfrag.ne.1) then
write(6,444) 'ERROR: Monitoring fragmentation ', &
'all atoms in the AG must be bonded one another'
stop
endif
endif
if (.not.stateselect) then
write(6,*) 'Equilibrium run to provide initial conditions'
write(6,*) 'Probability to pick a configuration is ',pickthr
endif
if (ipicktrj.ne.0) then
if (demin.eq.-1.d50 .and. demax.eq.1.d50) then
write(6,*) 'No constraint on the energy of the points picked'
elseif (demin.ne.-1.d50) then
write(6,*) 'Pick those with energy >= than ',demin,' eV'
demin=demin/autoev
elseif (demax.ne.1.d50) then
write(6,*) 'Pick those with energy < than ',demax,' eV'
demax=demax/autoev
endif
endif
elseif (termflag .eq. 8) then
write(6,108) 'TERMFLAG = 8: Trajectories will terminate ', &
'when specified dihedral angles rotate'
else
write(6,*) 'TERMFLAG = ',termflag,' is not supported'
stop
endif
if (t_nstep.ne.1000000000) then
write(6,206) 'Maximum number of steps per trajectory is ', &
t_nstep
else
write(6,108) 'Maximum number of steps per trajectory is', &
' unlimited'
endif
if (nstat.gt.t_nstep) then
write(6,108) 'Steps used to average properties', &
' cant be smaller than the total number of steps'
stop
endif
if (nbondtype.gt.0) then
write(6,*) 'Default bond forming/breaking threshold changed'
write(6,*) 'Outcome label Bond Distance (A)'
do j=1,nbondtype
write(6,905) j,titype(j,1),titype(j,2),trtype(j)*autoang
enddo
else
write(6,444) 'Default bond breaking/forming threshold used:', &
' 2.5/1.5 times of default atomic covalent bond distance'
endif
!
! Output information
write(6,*)
write(6,'(1x,74a1)') ('-',i=1,74)
write(6,*) 'Output write to units:'
write(6,'(1x,74a1)') ('-',i=1,74)
if (outflag.eq.0) then
write(6,333) 'OUTFLAG= ',outflag, &
' : Write output to all units'
else
write(6,333) 'OUTFLAG= ',outflag,': Write output to units', &
' specified by user:'
write(6,*) (ilwrite(j),j=1,outflag)
endif
if (maxprint) write(6,*) 'Print out maximum information'
if (minprinticon) write(6,*) 'Minimum printout of initial', &
' conditions'
if (minprinttrapz) write(6,*) 'Skipping intermediate info', &
' in TRAPZ'
if (idebug) write(6,108) 'Print out debug information. ',&
'WARNING: Very lengthy'
write(6,*) 'Print general information every ',nprint,' steps'
!
! Some values control the program
ii=mnmol+1
mmofag(ii)=0.d0
nvibmods(ii)=0
nrotmods(ii)=0
do im=1,nmol
if (natom(im).eq.1) then
nvibmods(im)=0
nrotmods(im)=0
elseif (natom(im).eq.2) then
nvibmods(im)=1
nrotmods(im)=1
else
if (nvibmods(im).eq.0 ) nvibmods(im)=3*natom(im)-6
if (nrotmods(im).eq.0 ) nrotmods(im)=3
endif
mmofag(ii)=mmofag(ii)+mmofag(im)
nvibmods(ii)=nvibmods(ii)+nvibmods(im)
nrotmods(ii)=nrotmods(ii)+nrotmods(im)
enddo
!
!
! Ensemble information
! For reactive collision run, use nve ensemble instead
if (itherm.ne.0) then
if (reactcol) &
write(6,444)'WARNING!!: using NVT or NPT ensemble for reactive', &
' collision run'
if (initx(1).eq.4) &
write(6,444) 'WARNING!!: using NVT or NPT ensemble for atom-', &
'diatom run'
endif
if (itherm.eq.0) then
write(6,*) 'NVE ensemble'
elseif (itherm.eq.1) then
write(6,*) 'NVT ensemble'
else
write(6,*) 'NPT ensemble'
write(6,*) 'Pressure is: ',p0*autoatm,' atm'
if (icell.eq.1) then
tmp=(dble(nat)*kb*temp0/p0)**(1.d0/3.d0)
write(6,*) 'Desired cell should have a box of ',tmp*autoang
endif
endif
if (itherm.ne.0 .or. ieqtherm.ne.0) then
tmp=(10.d0/autoang)**2*dsqrt(8.d0*pi*(mmofag(ii)+39.948d0*amutoau) &
/(mmofag(ii)*39.948d0*amutoau*kb*temp0))*p0
write(6,*) 'Collision frequency is ',tmp/autofs,' 1/fs'
if (taut.ne.0.d0) then
taut = taut/autofs
else
! taut = 1.d0/tmp
taut = hstep0/0.0025d0
endif
if (taup.ne.0.d0) then
taup = taup/(autofs*autoatm)
else
taup = hstep0*p0/0.0025d0
endif
if (vfreq.ne.0.d0) then
vfreq = vfreq*autofs
else
! vfreq = 0.1d0/hstep0
vfreq = tmp
endif
if (itempzpe) then
write(6,*) 'You ve chosen to count ZPE into temperature'
write(6,444)'The real average temperature should belower than', &
' the calculated temperature by a constant value'
endif
if (iadjtemp.eq.0) then
write(6,203) 'Berendsen thermostat. Coupling parameter tau ', &
'is : ',taut*autofs,' fs'
elseif (iadjtemp.eq.1) then
write(6,444) 'Simplet temperature control scheme: Scaling ',&
'momentum by square root of T0/T'
elseif (iadjtemp.eq.2) then
write(6,777) 'Andersen thermostat. Collision frequency is ',&
vfreq/autofs,' 1/fs'
elseif (iadjtemp.eq.3) then
write(6,203) 'Nose-Hoover thermostat. system vibration ', &
'frequency is: ',sysfreq,' 1/fs'
! Convert to vibration period ( au )
sysfreq = 1.d0/(sysfreq*autofs)
! write(6,*) 'Nose-Hoover thermostat. Q2 used is ',nh(1,2),&
! ' eV fs**2'
! nh(1,1)=nh(1,1)/(autoev*autofs**2)
! nh(1,2)=nh(1,2)/(autoev*autofs**2)
else
write(6,444) 'ERROR: Other thermostat have not been ', &
'implemented yet'
stop
endif
if (itherm.eq.2) then
if (iadjpress.eq.0) then
write(6,203) 'Berendsen barostat. Coupling parameter tau ', &
'is : ',taup*autofs*autoatm,' fs atm'
else
write(6,444) 'ERROR: Other barostat have not been ', &
'implemented yet'
stop
endif
endif
endif
write(6,*)
! All write formats come here
900 format(1x,i5,15x,i2,'-',i2,f15.5)
905 format(1x,i5,15x,a2,'-',a2,f15.5)
910 format(1x,a10,4a12)
920 format(1x,i5,a5,4f12.5)
930 format(1x,a,a,f12.5,1x,a)
940 format(1x,a,3f12.5,a)
950 format(1x,a,f12.5,a)
960 format(1x,a,i4,a,f12.5,a)
108 format(1x,10a)
111 format(1x,a,i10,2a,i10)
222 format(1x,2a,i10)
333 format(1x,a,i10,a)
444 format(1x,2a)
555 format(1x,3a)
666 format(1x,a,e16.8,a,e16.8,a)
777 format(1x,a,e16.8,a)
888 format(1x,a,3e16.8,a)
999 format(1x,2a,3e16.8,a)
200 format(1x,a,i10,2a,i10,a,i10)
201 format(1x,a,i10,a,i10)
202 format(1x,a,e16.8,a,i15,a)
203 format(1x,2a,e16.8,a)
204 format(1x,3a,e16.8)
205 format(1x,2a,i10,a)
206 format(1x,a,i15)
END subroutine readin
!
function lowercase (string)
!
! Function which takes a string of 200 characters and converts the
! upper case letters in the string to lower case letters
! This function is taken from POLYRATE with a slight modification for
! character length.
!
character*200 :: string, lowercase
character*1 :: xlett
!
do 10 i = 1, 200
xlett = string(i:i)
itry = ichar (xlett)
if (xlett .ge. 'A' .and. xlett .le. 'Z') then
itry = itry + 32
string (i:i) = char (itry)
endif
10 continue
!
lowercase = string
end function lowercase
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment