c c c ################################################### c ## COPYRIGHT (C) 1990 by Jay William Ponder ## c ## All Rights Reserved ## c ################################################### c c ############################################################### c ## ## c ## subroutine mdinit -- initialize a dynamics trajectory ## c ## ## c ############################################################### c c c "mdinit" initializes the velocities and accelerations c for a molecular dynamics trajectory, including restarts c c subroutine mdinit (dt) use atomid use atoms use bath use bound use couple use extfld use files use freeze use group use inform use iounit use keys use mdstuf use molcul use moldyn use mpole use output use potent use rgddyn use rigid use stodyn use units use usage implicit none integer i,j,istep integer idyn,lext integer size,next integer ndummy integer freeunit integer trimtext real*8 dt,e,ekt,qterm real*8 maxwell,speed real*8 amass,gmass real*8 vec(3) real*8, allocatable :: derivs(:,:) logical exist character*7 ext character*20 keyword character*20 text character*240 dynfile character*240 record character*240 string c c c set default parameters for the dynamics trajectory c integrate = 'BEEMAN' bmnmix = 8 nrespa = max(1,nint(dt/0.0005d0)) nfree = 0 ndummy = 0 irest = -1 iprint = 100 use_wrap = .true. velsave = .false. frcsave = .false. uindsave = .false. friction = 0.5d0 if (use_solv) friction = 91.0d0 use_sdarea = .false. c c set default values for temperature and pressure control c thermostat = 'BUSSI' tautemp = -1.0d0 collide = 0.1d0 do i = 1, maxnose vnh(i) = 0.0d0 qnh(i) = 0.0d0 gnh(i) = 0.0d0 end do barostat = 'BUSSI' prestyp = 'ISOTROPIC' taupres = -1.0d0 compress = 0.000046d0 vbar = 0.0d0 qbar = 0.0d0 gbar = 0.0d0 voltrial = 25 volmove = 100.0d0 volscale = 'MOLECULAR' c c check for keywords containing any altered parameters c do i = 1, nkey next = 1 record = keyline(i) call gettext (record,keyword,next) call upcase (keyword) string = record(next:240) if (keyword(1:11) .eq. 'INTEGRATOR ') then call getword (record,integrate,next) call upcase (integrate) if (integrate .eq. 'RESPA') integrate = 'BRESPA' else if (keyword(1:14) .eq. 'BEEMAN-MIXING ') then read (string,*,err=10,end=10) bmnmix else if (keyword(1:12) .eq. 'RESPA-INNER ') then read (string,*,err=10,end=10) nrespa else if (keyword(1:16) .eq. 'DEGREES-FREEDOM ') then read (string,*,err=10,end=10) nfree else if (keyword(1:15) .eq. 'REMOVE-INERTIA ') then read (string,*,err=10,end=10) irest else if (keyword(1:12) .eq. 'UNWRAP-COORDS ') then use_wrap = .false. else if (keyword(1:14) .eq. 'SAVE-VELOCITY ') then velsave = .true. else if (keyword(1:11) .eq. 'SAVE-FORCE ') then frcsave = .true. else if (keyword(1:13) .eq. 'SAVE-INDUCED ') then uindsave = .true. else if (keyword(1:9) .eq. 'FRICTION ') then read (string,*,err=10,end=10) friction else if (keyword(1:17) .eq. 'FRICTION-SCALING ') then use_sdarea = .true. else if (keyword(1:11) .eq. 'THERMOSTAT ') then call getword (record,text,next) call upcase (text) if (text(1:5) .eq. 'BUSSI') thermostat = 'BUSSI' if (text(1:9) .eq. 'BERENDSEN') thermostat = 'BERENDSEN' if (text(1:4) .eq. 'NOSE') thermostat = 'NOSE-HOOVER' if (text(1:8) .eq. 'ANDERSEN') thermostat = 'ANDERSEN' else if (keyword(1:16) .eq. 'TAU-TEMPERATURE ') then read (string,*,err=10,end=10) tautemp else if (keyword(1:10) .eq. 'COLLISION ') then read (string,*,err=10,end=10) collide else if (keyword(1:9) .eq. 'BAROSTAT ') then call getword (record,text,next) call upcase (text) if (text(1:5) .eq. 'BUSSI') barostat = 'BUSSI' if (text(1:9) .eq. 'BERENDSEN') barostat = 'BERENDSEN' if (text(1:4) .eq. 'NOSE') barostat = 'NOSE-HOOVER' if (text(1:10) .eq. 'MONTECARLO') barostat = 'MONTECARLO' else if (keyword(1:13) .eq. 'TAU-PRESSURE ') then read (string,*,err=10,end=10) taupres else if (keyword(1:9) .eq. 'PRESSURE ') then call getword (record,text,next) call upcase (text) if (text(1:3) .eq. 'ISO') prestyp = 'ISOTROPIC' if (text(1:5) .eq. 'ANISO') prestyp = 'ANISO' if (text(1:4) .eq. 'SEMI') prestyp = 'SEMIISO' else if (keyword(1:9) .eq. 'COMPRESS ') then read (string,*,err=10,end=10) compress else if (keyword(1:13) .eq. 'VOLUME-TRIAL ') then read (string,*,err=10,end=10) voltrial else if (keyword(1:12) .eq. 'VOLUME-MOVE ') then read (string,*,err=10,end=10) volmove else if (keyword(1:13) .eq. 'VOLUME-SCALE ') then call getword (record,text,next) call upcase (text) if (text(1:9) .eq. 'MOLECULAR') volscale = 'MOLECULAR' if (text(1:6) .eq. 'ATOMIC') volscale = 'ATOMIC' else if (keyword(1:9) .eq. 'PRINTOUT ') then read (string,*,err=10,end=10) iprint end if 10 continue end do c c check for use of induced dipole prediction methods c if (use_polar) call predict c c make sure all atoms or groups have a nonzero mass c if (integrate .eq. 'RIGIDBODY') then do i = 1, ngrp if (grpmass(i) .le. 0.0d0) then grpmass(i) = 1.0d0 if (igrp(1,i) .le. igrp(2,i)) then totmass = totmass + 1.0d0 if (verbose) then write (iout,20) i 20 format (/,' MDINIT -- Warning, Mass of Group', & i6,' Set to 1.0 for Dynamics') end if end if end if end do else do i = 1, n if (atomic(i) .le. 0) ndummy = ndummy + 1 if (use(i) .and. mass(i).le.0.0d0) then mass(i) = 1.0d0 if (atomic(i) .gt. 0) totmass = totmass + 1.0d0 if (verbose) then write (iout,30) i 30 format (/,' MDINIT -- Warning, Mass of Atom', & i6,' Set to 1.0 for Dynamics') end if end if end do end if c c decide whether to remove center of mass motion; note c should be applied by default to stochastic dynamics c dorest = .true. if (irest .eq. 0) then dorest = .false. else if (irest .lt. 0) then if (nuse .ne. n) dorest = .false. c if (integrate .eq. 'BAOAB') dorest = .false. c if (integrate .eq. 'OBABO') dorest = .false. c if (integrate .eq. 'SRESPA') dorest = .false. c if (integrate .eq. 'STOCHASTIC') dorest = .false. if (integrate .eq. 'GHMC') dorest = .false. if (isothermal .and. thermostat.eq.'ANDERSEN') & dorest = .false. end if if (dorest) then if (irest .lt. 0) irest = 100 else irest = 0 end if c c enforce use of velocity Verlet with Andersen thermostat c if (thermostat .eq. 'ANDERSEN') then if (integrate .eq. 'BEEMAN') integrate = 'VERLET' if (integrate .eq. 'BRESPA') integrate = 'VRESPA' end if c c couple Nose-Hoover thermostat and barostat with integrator c if (integrate .eq. 'NOSE-HOOVER') then thermostat = 'NOSE-HOOVER' barostat = 'NOSE-HOOVER' else if (thermostat.eq.'NOSE-HOOVER' .and. & barostat.eq.'NOSE-HOOVER') then integrate = 'NOSE-HOOVER' end if c c apply default values for thermostat and barostat coupling c if (tautemp .lt. 0.0d0) then tautemp = 0.2d0 if (thermostat .eq. 'NOSE-HOOVER') tautemp = 1.0d0 end if if (taupres .lt. 0.0d0) then taupres = 2.0d0 if (barostat .eq. 'NOSE-HOOVER') taupres = 10.0d0 if (prestyp .eq. 'ANISO') taupres = 10.0d0 end if c c check for options not allowed with Monte Carlo barostat c if (barostat .eq. 'MONTECARLO') then if (use_freeze .and. volscale.eq.'ATOMIC') then write (iout,40) 40 format (/,' MDINIT -- No Atom-Based Monte Carlo', & ' Barostat with Constraints') call fatal else if (prestyp .eq. 'SEMIISO') then write (iout,50) 50 format (/,' MDINIT -- No Monte Carlo Barostat', & ' with Semi-Isotropic Pressure') call fatal end if end if c c perform dynamic allocation of some global arrays c if (integrate .eq. 'RIGIDBODY') then if (.not. allocated(xcmo)) allocate (xcmo(n)) if (.not. allocated(ycmo)) allocate (ycmo(n)) if (.not. allocated(zcmo)) allocate (zcmo(n)) if (.not. allocated(vcm)) allocate (vcm(3,ngrp)) if (.not. allocated(wcm)) allocate (wcm(3,ngrp)) if (.not. allocated(lm)) allocate (lm(3,ngrp)) if (.not. allocated(vc)) allocate (vc(3,ngrp)) if (.not. allocated(wc)) allocate (wc(3,ngrp)) if (.not. allocated(linear)) allocate (linear(ngrp)) else if (.not. allocated(v)) allocate (v(3,n)) if (.not. allocated(a)) allocate (a(3,n)) if (.not. allocated(aalt)) allocate (aalt(3,n)) if (.not. allocated(aslow)) allocate (aslow(3,n)) if (.not. allocated(afast)) allocate (afast(3,n)) end if c c set the number of degrees of freedom for the system c if (nfree .eq. 0) then if (integrate .eq. 'RIGIDBODY') then call grpline nfree = 6 * ngrp do i = 1, ngrp size = igrp(2,i) - igrp(1,i) + 1 if (size .eq. 1) nfree = nfree - 3 if (linear(i)) nfree = nfree - 1 end do else nfree = 3 * (nuse-ndummy) end if if (use_freeze) then nfree = nfree - nrat do i = 1, nratx nfree = nfree - kratx(i) end do nfree = nfree - 3*nwat end if if (isothermal .and. thermostat.ne.'ANDERSEN' & .and. integrate.ne.'BAOAB' & .and. integrate.ne.'OBABO' & .and. integrate.ne.'SRESPA' & .and. integrate.ne.'STOCHASTIC' & .and. integrate.ne.'GHMC') then if (.not. use_exfld) then if (use_bounds) then if (integrate.ne.'RIGIDBODY' .and. ngrp.ne.0) then nfree = nfree - 6*(ngrp+1) else nfree = nfree - 3 end if else nfree = nfree - 6 end if end if end if end if c c check for a nonzero number of degrees of freedom c if (nfree .lt. 0) nfree = 0 if (debug) then write (iout,60) nfree 60 format (/,' Number of Degrees of Freedom for Dynamics :',i10) end if if (nfree .eq. 0) then write (iout,70) 70 format (/,' MDINIT -- No Degrees of Freedom for Dynamics') call fatal end if c c set masses for Nose-Hoover thermostat and barostat c if (thermostat .eq. 'NOSE-HOOVER') then ekt = gasconst * kelvin qterm = ekt * tautemp * tautemp do j = 1, maxnose if (qnh(j) .eq. 0.0d0) qnh(j) = qterm end do qnh(1) = dble(nfree) * qnh(1) end if if (barostat .eq. 'NOSE-HOOVER') then ekt = gasconst * kelvin qterm = ekt * taupres * taupres qbar = dble(nfree+1) * qterm end if c c perform dynamic allocation of some local arrays c allocate (derivs(3,n)) c c try to restart using prior velocities and accelerations c dynfile = filename(1:leng)//'.dyn' call version (dynfile,'old') inquire (file=dynfile,exist=exist) if (exist) then call gradient (e,derivs) idyn = freeunit () open (unit=idyn,file=dynfile,status='old') rewind (unit=idyn) call readdyn (idyn) close (unit=idyn) write (iout,80) dynfile(1:trimtext(dynfile)) 80 format (/,' Restarting Molecular Dynamics Using : ',a) c c set translational velocities for rigid body dynamics c else if (integrate .eq. 'RIGIDBODY') then call gradient (e,derivs) do i = 1, ngrp gmass = grpmass(i) speed = maxwell (gmass,kelvin) call ranvec (vec) do j = 1, 3 vcm(j,i) = speed * vec(j) wcm(j,i) = 0.0d0 lm(j,i) = 0.0d0 end do end do if (nuse .eq. n) then istep = 0 call mdrest (istep) end if c c set velocities and fast/slow accelerations for RESPA methods c else if (integrate.eq.'VRESPA' .or. integrate.eq.'BRESPA' & .or. integrate.eq.'SRESPA') then call gradslow (e,derivs) do i = 1, n amass = mass(i) if (use(i) .and. amass.ne.0.0d0) then speed = maxwell (amass,kelvin) call ranvec (vec) do j = 1, 3 v(j,i) = speed * vec(j) a(j,i) = -ekcal * derivs(j,i) / mass(i) aslow(j,i) = a(j,i) end do else do j = 1, 3 v(j,i) = 0.0d0 a(j,i) = 0.0d0 aslow(j,i) = 0.0d0 end do end if end do call gradfast (e,derivs) do i = 1, n amass = mass(i) if (use(i) .and. amass.ne.0.0d0) then do j = 1, 3 aalt(j,i) = -ekcal * derivs(j,i) / amass afast(j,i) = aalt(j,i) end do else do j = 1, 3 aalt(j,i) = 0.0d0 afast(j,i) = aalt(j,i) end do end if end do if (nuse .eq. n) then istep = 0 call mdrest (istep) end if c c set velocities and accelerations for Cartesian dynamics c else call gradient (e,derivs) do i = 1, n amass = mass(i) if (use(i) .and. amass.ne.0.0d0) then speed = maxwell (amass,kelvin) call ranvec (vec) do j = 1, 3 v(j,i) = speed * vec(j) a(j,i) = -ekcal * derivs(j,i) / amass aalt(j,i) = a(j,i) end do else do j = 1, 3 v(j,i) = 0.0d0 a(j,i) = 0.0d0 aalt(j,i) = 0.0d0 end do end if end do if (nuse .eq. n) then istep = 0 call mdrest (istep) end if end if c c perform deallocation of some local arrays c deallocate (derivs) c c check for any prior dynamics coordinate sets c i = 0 exist = .true. do while (exist) i = i + 1 lext = 3 call numeral (i,ext,lext) dynfile = filename(1:leng)//'.'//ext(1:lext) inquire (file=dynfile,exist=exist) if (.not.exist .and. i.lt.100) then lext = 2 call numeral (i,ext,lext) dynfile = filename(1:leng)//'.'//ext(1:lext) inquire (file=dynfile,exist=exist) end if if (.not.exist .and. i.lt.10) then lext = 1 call numeral (i,ext,lext) dynfile = filename(1:leng)//'.'//ext(1:lext) inquire (file=dynfile,exist=exist) end if end do nprior = i - 1 return end