1c 2c 3c ################################################### 4c ## COPYRIGHT (C) 2011 by Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################### 7c 8c ############################################################# 9c ## ## 10c ## subroutine respa -- r-RESPA molecular dynamics step ## 11c ## ## 12c ############################################################# 13c 14c 15c "respa" performs a single multiple time step molecular dynamics 16c step using the reversible reference system propagation algorithm 17c (r-RESPA) via a Verlet core with the potential split into fast- 18c and slow-evolving portions 19c 20c literature references: 21c 22c D. D. Humphreys, R. A. Friesner and B. J. Berne, "A Multiple- 23c Time-Step Molecular Dynamics Algorithm for Macromolecules", 24c Journal of Physical Chemistry, 98, 6885-6892 (1994) 25c 26c X. Qian and T. Schlick, "Efficient Multiple-Time-Step Integrators 27c with Distance-Based Force Splitting for Particle-Mesh-Ewald 28c Molecular Dynamics Simulations", Journal of Chemical Physics, 29c 115, 4019-4029 (2001) 30c 31c 32 subroutine respa (istep,dt) 33 use atomid 34 use atoms 35 use freeze 36 use ielscf 37 use mdstuf 38 use moldyn 39 use polar 40 use units 41 use usage 42 use virial 43 implicit none 44 integer i,j,k,m 45 integer istep 46 integer nalt 47 real*8 dt,dt_2 48 real*8 dta,dta_2 49 real*8 epot,etot 50 real*8 eksum,eps 51 real*8 temp,pres 52 real*8 ealt,dalt 53 real*8 term 54 real*8 ekin(3,3) 55 real*8 stress(3,3) 56 real*8 viralt(3,3) 57 real*8, allocatable :: xold(:) 58 real*8, allocatable :: yold(:) 59 real*8, allocatable :: zold(:) 60 real*8, allocatable :: derivs(:,:) 61c 62c 63c set some time values for the dynamics integration 64c 65 eps = 0.00000001d0 66 dalt = arespa 67 nalt = int(dt/(dalt+eps)) + 1 68 dalt = dble(nalt) 69 dt_2 = 0.5d0 * dt 70 dta = dt / dalt 71 dta_2 = 0.5d0 * dta 72c 73c store the current atom positions, then find half-step 74c velocities via velocity Verlet recursion 75c 76 do i = 1, nuse 77 m = iuse(i) 78 do j = 1, 3 79 v(j,m) = v(j,m) + a(j,m)*dt_2 80 end do 81 end do 82c 83c initialize virial from fast-evolving potential energy terms 84c 85 do i = 1, 3 86 do j = 1, 3 87 viralt(j,i) = 0.0d0 88 end do 89 end do 90c 91c perform dynamic allocation of some local arrays 92c 93 allocate (xold(n)) 94 allocate (yold(n)) 95 allocate (zold(n)) 96 allocate (derivs(3,n)) 97c 98c find fast-evolving velocities and positions via Verlet recursion 99c 100 do k = 1, nalt 101 do i = 1, nuse 102 m = iuse(i) 103 do j = 1, 3 104 v(j,m) = v(j,m) + aalt(j,m)*dta_2 105 end do 106 xold(m) = x(m) 107 yold(m) = y(m) 108 zold(m) = z(m) 109 x(m) = x(m) + v(1,m)*dta 110 y(m) = y(m) + v(2,m)*dta 111 z(m) = z(m) + v(3,m)*dta 112 end do 113 if (use_rattle) call rattle (dta,xold,yold,zold) 114c 115c get the fast-evolving potential energy and atomic forces 116c 117 call gradfast (ealt,derivs) 118c 119c use Newton's second law to get fast-evolving accelerations; 120c update fast-evolving velocities using the Verlet recursion 121c 122 do i = 1, nuse 123 m = iuse(i) 124 do j = 1, 3 125 aalt(j,m) = -ekcal * derivs(j,m) / mass(m) 126 v(j,m) = v(j,m) + aalt(j,m)*dta_2 127 end do 128 end do 129 if (use_rattle) call rattle2 (dta) 130c 131c find average virial from fast-evolving potential terms 132c 133 do i = 1, 3 134 do j = 1, 3 135 viralt(j,i) = viralt(j,i) + vir(j,i)/dalt 136 end do 137 end do 138 end do 139c 140c apply Verlet half-step updates for any auxiliary dipoles 141c 142 if (use_ielscf) then 143 do i = 1, nuse 144 m = iuse(i) 145 do j = 1, 3 146 vaux(j,m) = vaux(j,m) + aaux(j,m)*dt_2 147 vpaux(j,m) = vpaux(j,m) + apaux(j,m)*dt_2 148 uaux(j,m) = uaux(j,m) + vaux(j,m)*dt 149 upaux(j,m) = upaux(j,m) + vpaux(j,m)*dt 150 end do 151 end do 152 end if 153c 154c get the slow-evolving potential energy and atomic forces 155c 156 call gradslow (epot,derivs) 157 epot = epot + ealt 158c 159c make half-step temperature and pressure corrections 160c 161 call temper2 (dt,temp) 162 call pressure2 (epot,temp) 163c 164c use Newton's second law to get the slow accelerations; 165c find full-step velocities using velocity Verlet recursion 166c 167 do i = 1, nuse 168 m = iuse(i) 169 do j = 1, 3 170 a(j,m) = -ekcal * derivs(j,m) / mass(m) 171 v(j,m) = v(j,m) + a(j,m)*dt_2 172 end do 173 end do 174c 175c apply Verlet full-step updates for any auxiliary dipoles 176c 177 if (use_ielscf) then 178 term = 2.0d0 / (dt*dt) 179 do i = 1, nuse 180 m = iuse(i) 181 do j = 1, 3 182 aaux(j,m) = term * (uind(j,m)-uaux(j,m)) 183 apaux(j,m) = term * (uinp(j,m)-upaux(j,m)) 184 vaux(j,m) = vaux(j,m) + aaux(j,m)*dt_2 185 vpaux(j,m) = vpaux(j,m) + apaux(j,m)*dt_2 186 end do 187 end do 188 end if 189c 190c perform deallocation of some local arrays 191c 192 deallocate (xold) 193 deallocate (yold) 194 deallocate (zold) 195 deallocate (derivs) 196c 197c find the constraint-corrected full-step velocities 198c 199 if (use_rattle) call rattle2 (dt) 200c 201c increment total virial from sum of fast and slow parts 202c 203 do i = 1, 3 204 do j = 1, 3 205 vir(j,i) = vir(j,i) + viralt(j,i) 206 end do 207 end do 208c 209c make full-step temperature and pressure corrections 210c 211 call temper (dt,eksum,ekin,temp) 212 call pressure (dt,epot,ekin,temp,pres,stress) 213c 214c total energy is sum of kinetic and potential energies 215c 216 etot = eksum + epot 217c 218c compute statistics and save trajectory for this step 219c 220 call mdstat (istep,dt,etot,epot,eksum,temp,pres) 221 call mdsave (istep,dt,epot,eksum) 222 call mdrest (istep) 223 return 224 end 225c 226c 227c ################################################################## 228c ## ## 229c ## subroutine gradfast -- fast energy & gradient components ## 230c ## ## 231c ################################################################## 232c 233c 234c "gradfast" calculates the potential energy and first derivatives 235c for the fast-evolving local valence potential energy terms 236c 237c 238 subroutine gradfast (energy,derivs) 239 use limits 240 use potent 241 implicit none 242 real*8 energy 243 real*8 derivs(3,*) 244 logical save_vdw,save_repuls 245 logical save_disp,save_charge 246 logical save_chgdpl,save_dipole 247 logical save_mpole,save_polar 248 logical save_chgtrn,save_rxnfld 249 logical save_solv,save_list 250c 251c 252c save the original state of slow-evolving potentials 253c 254 save_vdw = use_vdw 255 save_repuls = use_repuls 256 save_disp = use_disp 257 save_charge = use_charge 258 save_chgdpl = use_chgdpl 259 save_dipole = use_dipole 260 save_mpole = use_mpole 261 save_polar = use_polar 262 save_chgtrn = use_chgtrn 263 save_rxnfld = use_rxnfld 264 save_solv = use_solv 265 save_list = use_list 266c 267c turn off slow-evolving nonbonded potential energy terms 268c 269 use_vdw = .false. 270 use_repuls = .false. 271 use_disp = .false. 272 use_charge = .false. 273 use_chgdpl = .false. 274 use_dipole = .false. 275 use_mpole = .false. 276 use_polar = .false. 277 use_chgtrn = .false. 278 use_rxnfld = .false. 279 use_solv = .false. 280 use_list = .false. 281c 282c get energy and gradient for fast-evolving potential terms 283c 284 call gradient (energy,derivs) 285c 286c restore the original state of slow-evolving potentials 287c 288 use_vdw = save_vdw 289 use_repuls = save_repuls 290 use_disp = save_disp 291 use_charge = save_charge 292 use_chgdpl = save_chgdpl 293 use_dipole = save_dipole 294 use_mpole = save_mpole 295 use_polar = save_polar 296 use_chgtrn = save_chgtrn 297 use_rxnfld = save_rxnfld 298 use_solv = save_solv 299 use_list = save_list 300 return 301 end 302c 303c 304c ################################################################## 305c ## ## 306c ## subroutine gradslow -- slow energy & gradient components ## 307c ## ## 308c ################################################################## 309c 310c 311c "gradslow" calculates the potential energy and first derivatives 312c for the slow-evolving nonbonded potential energy terms 313c 314c 315 subroutine gradslow (energy,derivs) 316 use potent 317 implicit none 318 real*8 energy 319 real*8 derivs(3,*) 320 logical save_bond,save_angle 321 logical save_strbnd,save_urey 322 logical save_angang,save_opbend 323 logical save_opdist,save_improp 324 logical save_imptor,save_tors 325 logical save_pitors,save_strtor 326 logical save_angtor,save_tortor 327 logical save_geom,save_metal 328 logical save_extra 329c 330c 331c save the original state of fast-evolving potentials 332c 333 save_bond = use_bond 334 save_angle = use_angle 335 save_strbnd = use_strbnd 336 save_urey = use_urey 337 save_angang = use_angang 338 save_opbend = use_opbend 339 save_opdist = use_opdist 340 save_improp = use_improp 341 save_imptor = use_imptor 342 save_tors = use_tors 343 save_pitors = use_pitors 344 save_strtor = use_strtor 345 save_angtor = use_angtor 346 save_tortor = use_tortor 347 save_geom = use_geom 348 save_metal = use_metal 349 save_extra = use_extra 350c 351c turn off fast-evolving valence potential energy terms 352c 353 use_bond = .false. 354 use_angle = .false. 355 use_strbnd = .false. 356 use_urey = .false. 357 use_angang = .false. 358 use_opbend = .false. 359 use_opdist = .false. 360 use_improp = .false. 361 use_imptor = .false. 362 use_tors = .false. 363 use_pitors = .false. 364 use_strtor = .false. 365 use_angtor = .false. 366 use_tortor = .false. 367 use_geom = .false. 368 use_metal = .false. 369 use_extra = .false. 370c 371c get energy and gradient for slow-evolving potential terms 372c 373 call gradient (energy,derivs) 374c 375c restore the original state of fast-evolving potentials 376c 377 use_bond = save_bond 378 use_angle = save_angle 379 use_strbnd = save_strbnd 380 use_urey = save_urey 381 use_angang = save_angang 382 use_opbend = save_opbend 383 use_opdist = save_opdist 384 use_improp = save_improp 385 use_imptor = save_imptor 386 use_tors = save_tors 387 use_pitors = save_pitors 388 use_strtor = save_strtor 389 use_angtor = save_angtor 390 use_tortor = save_tortor 391 use_geom = save_geom 392 use_metal = save_metal 393 use_extra = save_extra 394 return 395 end 396