1c 2c 3c ################################################### 4c ## COPYRIGHT (C) 1990 by Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################### 7c 8c ############################################################# 9c ## ## 10c ## subroutine beeman -- Beeman molecular dynamics step ## 11c ## ## 12c ############################################################# 13c 14c 15c "beeman" performs a single molecular dynamics time step 16c via the Beeman multistep recursion formula; uses original 17c coefficients or Bernie Brooks' "Better Beeman" values 18c 19c literature references: 20c 21c D. Beeman, "Some Multistep Methods for Use in Molecular 22c Dynamics Calculations", Journal of Computational Physics, 23c 20, 130-139 (1976) 24c 25c B. R. Brooks, "Algorithms for Molecular Dynamics at Constant 26c Temperature and Pressure", DCRT Report, NIH, April 1988 27c 28c 29 subroutine beeman (istep,dt) 30 use atomid 31 use atoms 32 use freeze 33 use ielscf 34 use mdstuf 35 use moldyn 36 use polar 37 use units 38 use usage 39 implicit none 40 integer i,j,k 41 integer istep 42 real*8 dt,dt_x,factor 43 real*8 etot,eksum,epot 44 real*8 temp,pres 45 real*8 part1,part2 46 real*8 dt_2,term 47 real*8 ekin(3,3) 48 real*8 stress(3,3) 49 real*8, allocatable :: xold(:) 50 real*8, allocatable :: yold(:) 51 real*8, allocatable :: zold(:) 52 real*8, allocatable :: derivs(:,:) 53c 54c 55c set time values and coefficients for Beeman integration 56c 57 factor = dble(bmnmix) 58 dt_x = dt / factor 59 part1 = 0.5d0*factor + 1.0d0 60 part2 = part1 - 2.0d0 61c 62c perform dynamic allocation of some local arrays 63c 64 allocate (xold(n)) 65 allocate (yold(n)) 66 allocate (zold(n)) 67 allocate (derivs(3,n)) 68c 69c store the current atom positions, then find half-step 70c velocities and full-step positions via Beeman recursion 71c 72 do i = 1, nuse 73 k = iuse(i) 74 do j = 1, 3 75 v(j,k) = v(j,k) + (part1*a(j,k)-aalt(j,k))*dt_x 76 end do 77 xold(k) = x(k) 78 yold(k) = y(k) 79 zold(k) = z(k) 80 x(k) = x(k) + v(1,k)*dt 81 y(k) = y(k) + v(2,k)*dt 82 z(k) = z(k) + v(3,k)*dt 83 end do 84c 85c apply Verlet half-step updates for any auxiliary dipoles 86c 87 if (use_ielscf) then 88 dt_2 = 0.5d0 * dt 89 do i = 1, nuse 90 k = iuse(i) 91 do j = 1, 3 92 vaux(j,k) = vaux(j,k) + aaux(j,k)*dt_2 93 vpaux(j,k) = vpaux(j,k) + apaux(j,k)*dt_2 94 uaux(j,k) = uaux(j,k) + vaux(j,k)*dt 95 upaux(j,k) = upaux(j,k) + vpaux(j,k)*dt 96 end do 97 end do 98 end if 99c 100c get constraint-corrected positions and half-step velocities 101c 102 if (use_rattle) call rattle (dt,xold,yold,zold) 103c 104c get the potential energy and atomic forces 105c 106 call gradient (epot,derivs) 107c 108c make half-step temperature and pressure corrections 109c 110 call temper2 (dt,temp) 111 call pressure2 (epot,temp) 112c 113c use Newton's second law to get the next accelerations; 114c find the full-step velocities using the Beeman recursion 115c 116 do i = 1, nuse 117 k = iuse(i) 118 do j = 1, 3 119 aalt(j,k) = a(j,k) 120 a(j,k) = -ekcal * derivs(j,k) / mass(k) 121 v(j,k) = v(j,k) + (part2*a(j,k)+aalt(j,k))*dt_x 122 end do 123 end do 124c 125c apply Verlet full-step updates for any auxiliary dipoles 126c 127 if (use_ielscf) then 128 term = 2.0d0 / (dt*dt) 129 do i = 1, nuse 130 k = iuse(i) 131 do j = 1, 3 132 aaux(j,k) = term * (uind(j,k)-uaux(j,k)) 133 apaux(j,k) = term * (uinp(j,k)-upaux(j,k)) 134 vaux(j,k) = vaux(j,k) + aaux(j,k)*dt_2 135 vpaux(j,k) = vpaux(j,k) + apaux(j,k)*dt_2 136 end do 137 end do 138 end if 139c 140c perform deallocation of some local arrays 141c 142 deallocate (xold) 143 deallocate (yold) 144 deallocate (zold) 145 deallocate (derivs) 146c 147c find the constraint-corrected full-step velocities 148c 149 if (use_rattle) call rattle2 (dt) 150c 151c make full-step temperature and pressure corrections 152c 153 call temper (dt,eksum,ekin,temp) 154 call pressure (dt,epot,ekin,temp,pres,stress) 155c 156c total energy is sum of kinetic and potential energies 157c 158 etot = eksum + epot 159c 160c compute statistics and save trajectory for this step 161c 162 call mdstat (istep,dt,etot,epot,eksum,temp,pres) 163 call mdsave (istep,dt,epot,eksum) 164 call mdrest (istep) 165 return 166 end 167