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