1c 2c 3c ################################################### 4c ## COPYRIGHT (C) 1990 by Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################### 7c 8c ################################################################# 9c ## ## 10c ## program dynamic -- run molecular or stochastic dynamics ## 11c ## ## 12c ################################################################# 13c 14c 15c "dynamic" computes a molecular or stochastic dynamics trajectory 16c in one of the standard statistical mechanical ensembles and using 17c any of several possible integration methods 18c 19c 20 program dynamic 21 use atoms 22 use bath 23 use bndstr 24 use bound 25 use inform 26 use iounit 27 use keys 28 use mdstuf 29 use potent 30 use stodyn 31 use usage 32 implicit none 33 integer i,istep,nstep 34 integer mode,next 35 real*8 dt,dtsave 36 logical exist 37 character*20 keyword 38 character*240 record 39 character*240 string 40c 41c 42c set up the structure and molecular mechanics calculation 43c 44 call initial 45 call getxyz 46 call mechanic 47c 48c initialize the temperature, pressure and coupling baths 49c 50 kelvin = 0.0d0 51 atmsph = 0.0d0 52 isothermal = .false. 53 isobaric = .false. 54c 55c check for keywords containing any altered parameters 56c 57 integrate = 'BEEMAN' 58 do i = 1, nkey 59 next = 1 60 record = keyline(i) 61 call gettext (record,keyword,next) 62 call upcase (keyword) 63 string = record(next:240) 64 if (keyword(1:11) .eq. 'INTEGRATOR ') then 65 call getword (record,integrate,next) 66 call upcase (integrate) 67 end if 68 end do 69c 70c initialize the simulation length as number of time steps 71c 72 nstep = -1 73 call nextarg (string,exist) 74 if (exist) read (string,*,err=10,end=10) nstep 75 10 continue 76 do while (nstep .lt. 0) 77 write (iout,20) 78 20 format (/,' Enter the Number of Dynamics Steps to be', 79 & ' Taken : ',$) 80 read (input,30,err=40) nstep 81 30 format (i10) 82 if (nstep .lt. 0) nstep = 0 83 40 continue 84 end do 85c 86c get the length of the dynamics time step in picoseconds 87c 88 dt = -1.0d0 89 call nextarg (string,exist) 90 if (exist) read (string,*,err=50,end=50) dt 91 50 continue 92 do while (dt .lt. 0.0d0) 93 write (iout,60) 94 60 format (/,' Enter the Time Step Length in Femtoseconds', 95 & ' [1.0] : ',$) 96 read (input,70,err=80) dt 97 70 format (f20.0) 98 if (dt .le. 0.0d0) dt = 1.0d0 99 80 continue 100 end do 101 dt = 0.001d0 * dt 102c 103c enforce bounds on thermostat and barostat coupling times 104c 105 tautemp = max(tautemp,dt) 106 taupres = max(taupres,dt) 107c 108c set the time between trajectory snapshot coordinate saves 109c 110 dtsave = -1.0d0 111 call nextarg (string,exist) 112 if (exist) read (string,*,err=90,end=90) dtsave 113 90 continue 114 do while (dtsave .lt. 0.0d0) 115 write (iout,100) 116 100 format (/,' Enter Time between Saves in Picoseconds', 117 & ' [0.1] : ',$) 118 read (input,110,err=120) dtsave 119 110 format (f20.0) 120 if (dtsave .le. 0.0d0) dtsave = 0.1d0 121 120 continue 122 end do 123 iwrite = nint(dtsave/dt) 124c 125c get choice of statistical ensemble for periodic system 126c 127 if (use_bounds) then 128 mode = -1 129 call nextarg (string,exist) 130 if (exist) read (string,*,err=130,end=130) mode 131 130 continue 132 do while (mode.lt.1 .or. mode.gt.4) 133 write (iout,140) 134 140 format (/,' Available Statistical Mechanical Ensembles :', 135 & //,4x,'(1) Microcanonical (NVE)', 136 & /,4x,'(2) Canonical (NVT)', 137 & /,4x,'(3) Isoenthalpic-Isobaric (NPH)', 138 & /,4x,'(4) Isothermal-Isobaric (NPT)', 139 & //,' Enter the Number of the Desired Choice', 140 & ' [1] : ',$) 141 read (input,150,err=160) mode 142 150 format (i10) 143 if (mode .le. 0) mode = 1 144 160 continue 145 end do 146 if (integrate.eq.'BUSSI' .or. integrate.eq.'NOSE-HOOVER' 147 & .or. integrate.eq.'GHMC') then 148 if (mode .ne. 4) then 149 mode = 4 150 write (iout,170) 151 170 format (/,' Switching to NPT Ensemble as Required', 152 & ' by Chosen Integrator') 153 end if 154 end if 155 if (mode.eq.2 .or. mode.eq.4) then 156 isothermal = .true. 157 kelvin = -1.0d0 158 call nextarg (string,exist) 159 if (exist) read (string,*,err=180,end=180) kelvin 160 180 continue 161 do while (kelvin .le. 0.0d0) 162 write (iout,190) 163 190 format (/,' Enter the Desired Temperature in Degrees', 164 & ' K [298] : ',$) 165 read (input,200,err=210) kelvin 166 200 format (f20.0) 167 if (kelvin .le. 0.0d0) kelvin = 298.0d0 168 210 continue 169 end do 170 end if 171 if (mode.eq.3 .or. mode.eq.4) then 172 isobaric = .true. 173 atmsph = -1.0d0 174 call nextarg (string,exist) 175 if (exist) read (string,*,err=220,end=220) atmsph 176 220 continue 177 do while (atmsph .eq. -1.0d0) 178 write (iout,230) 179 230 format (/,' Enter the Desired Pressure in Atm', 180 & ' [1.0] : ',$) 181 read (input,240,err=250) atmsph 182 240 format (f20.0) 183 if (atmsph .eq. -1.0d0) atmsph = 1.0d0 184 250 continue 185 end do 186 end if 187 end if 188c 189c use constant energy or temperature for nonperiodic system 190c 191 if (.not. use_bounds) then 192 mode = -1 193 call nextarg (string,exist) 194 if (exist) read (string,*,err=260,end=260) mode 195 260 continue 196 do while (mode.lt.1 .or. mode.gt.2) 197 write (iout,270) 198 270 format (/,' Available Simulation Control Modes :', 199 & //,4x,'(1) Constant Total Energy Value (E)', 200 & /,4x,'(2) Constant Temperature via Thermostat (T)', 201 & //,' Enter the Number of the Desired Choice', 202 & ' [1] : ',$) 203 read (input,280,err=290) mode 204 280 format (i10) 205 if (mode .le. 0) mode = 1 206 290 continue 207 end do 208 if (mode .eq. 2) then 209 isothermal = .true. 210 kelvin = -1.0d0 211 call nextarg (string,exist) 212 if (exist) read (string,*,err=300,end=300) kelvin 213 300 continue 214 do while (kelvin .le. 0.0d0) 215 write (iout,310) 216 310 format (/,' Enter the Desired Temperature in Degrees', 217 & ' K [298] : ',$) 218 read (input,320,err=330) kelvin 219 320 format (f20.0) 220 if (kelvin .le. 0.0d0) kelvin = 298.0d0 221 330 continue 222 end do 223 end if 224 end if 225c 226c perform the setup functions needed to run dynamics 227c 228 call mdinit 229c 230c print out a header line for the dynamics computation 231c 232 if (integrate .eq. 'VERLET') then 233 write (iout,340) 234 340 format (/,' Molecular Dynamics Trajectory via', 235 & ' Velocity Verlet Algorithm') 236 else if (integrate .eq. 'STOCHASTIC') then 237 write (iout,350) 238 350 format (/,' Stochastic Dynamics Trajectory via', 239 & ' Velocity Verlet Algorithm') 240 else if (integrate .eq. 'BAOAB') then 241 write (iout,360) 242 360 format (/,' Constrained Stochastic Dynamics Trajectory', 243 & ' via BAOAB Algorithm') 244 else if (integrate .eq. 'BUSSI') then 245 write (iout,370) 246 370 format (/,' Molecular Dynamics Trajectory via', 247 & ' Bussi-Parrinello NPT Algorithm') 248 else if (integrate .eq. 'NOSE-HOOVER') then 249 write (iout,380) 250 380 format (/,' Molecular Dynamics Trajectory via', 251 & ' Nose-Hoover NPT Algorithm') 252 else if (integrate .eq. 'GHMC') then 253 write (iout,390) 254 390 format (/,' Stochastic Dynamics Trajectory via', 255 & ' Generalized Hybrid Monte Carlo') 256 else if (integrate .eq. 'RIGIDBODY') then 257 write (iout,400) 258 400 format (/,' Molecular Dynamics Trajectory via', 259 & ' Rigid Body Algorithm') 260 else if (integrate .eq. 'RESPA') then 261 write (iout,410) 262 410 format (/,' Molecular Dynamics Trajectory via', 263 & ' r-RESPA MTS Algorithm') 264 else 265 write (iout,420) 266 420 format (/,' Molecular Dynamics Trajectory via', 267 & ' Modified Beeman Algorithm') 268 end if 269 flush (iout) 270c 271c integrate equations of motion to take a time step 272c 273 do istep = 1, nstep 274 if (integrate .eq. 'VERLET') then 275 call verlet (istep,dt) 276 else if (integrate .eq. 'STOCHASTIC') then 277 call sdstep (istep,dt) 278 else if (integrate .eq. 'BAOAB') then 279 call baoab (istep,dt) 280 else if (integrate .eq. 'BUSSI') then 281 call bussi (istep,dt) 282 else if (integrate .eq. 'NOSE-HOOVER') then 283 call nose (istep,dt) 284 else if (integrate .eq. 'GHMC') then 285 call ghmcstep (istep,dt) 286 else if (integrate .eq. 'RIGIDBODY') then 287 call rgdstep (istep,dt) 288 else if (integrate .eq. 'RESPA') then 289 call respa (istep,dt) 290 else 291 call beeman (istep,dt) 292 end if 293 end do 294c 295c perform any final tasks before program exit 296c 297 call final 298 end 299