Skip to content

Commit cbe6b45

Browse files
jayponderjayponder
committed
alter MD integrators to better enforce constraints for NPT
1 parent 415e8b6 commit cbe6b45

File tree

7 files changed

+350
-76
lines changed

7 files changed

+350
-76
lines changed

‎other/obabo.f‎

Lines changed: 173 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,173 @@
1+
c
2+
c
3+
c ###############################################################
4+
c ## COPYRIGHT (C) 2016 by Charles Matthews & Ben Leimkuhler ##
5+
c ## COPYRIGHT (C) 2017 by Jay William Ponder ##
6+
c ## All Rights Reserved ##
7+
c ###############################################################c
8+
c
9+
c ############################################################
10+
c ## ##
11+
c ## subroutine obabo -- OBABO stochastic dynamics step ##
12+
c ## ##
13+
c ############################################################
14+
c
15+
c
16+
c "obabo" implements a constrained stochastic dynamics time
17+
c step using the geodesic OBABO scheme
18+
c
19+
c literature reference:
20+
c
21+
c B. Leimkuhler and C. Matthews, "Efficient Molecular Dynamics
22+
c Using Geodesic Integration and Solvent-Solute Splitting",
23+
c Proceedings of the Royal Society A, 472, 20160138 (2016)
24+
c
25+
c note this code likely needs to have the constraint portions
26+
c reworked, similar to the BAOAB routine in the main source
27+
c
28+
c
29+
subroutine obabo (istep,dt)
30+
use atomid
31+
use atoms
32+
use freeze
33+
use moldyn
34+
use units
35+
use usage
36+
implicit none
37+
integer i,j,k
38+
integer istep
39+
integer nrattle
40+
real*8 dt,dt_2,dtr
41+
real*8 etot,epot
42+
real*8 eksum
43+
real*8 temp,pres
44+
real*8 ekin(3,3)
45+
real*8 stress(3,3)
46+
real*8, allocatable :: xold(:)
47+
real*8, allocatable :: yold(:)
48+
real*8, allocatable :: zold(:)
49+
real*8, allocatable :: vfric(:)
50+
real*8, allocatable :: vrand(:,:)
51+
real*8, allocatable :: derivs(:,:)
52+
c
53+
c
54+
c set some time values for the dynamics integration
55+
c
56+
dt_2 = 0.5d0 * dt
57+
nrattle = 1
58+
dtr = dt / dble(nrattle)
59+
c
60+
c perform dynamic allocation of some local arrays
61+
c
62+
allocate (xold(n))
63+
allocate (yold(n))
64+
allocate (zold(n))
65+
allocate (derivs(3,n))
66+
allocate (vfric(n))
67+
allocate (vrand(3,n))
68+
c
69+
c update velocities with frictional and random components
70+
c
71+
call oprep (dt_2,vfric,vrand)
72+
do i = 1, nuse
73+
k = iuse(i)
74+
do j = 1, 3
75+
v(j,k) = v(j,k)*vfric(k) + vrand(j,k)
76+
end do
77+
end do
78+
if (use_rattle) call rattle2 (dt_2)
79+
c
80+
c find half-step velocities via the Verlet recursion
81+
c
82+
do i = 1, nuse
83+
k = iuse(i)
84+
do j = 1, 3
85+
v(j,k) = v(j,k) + a(j,k)*dt_2
86+
end do
87+
end do
88+
c
89+
c find the constraint-corrected full-step velocities
90+
c
91+
if (use_rattle) call rattle2 (dt)
92+
c
93+
c take first A step according to the BAOAB sequence
94+
c
95+
do j = 1, nrattle
96+
do i = 1, nuse
97+
k = iuse(i)
98+
xold(k) = x(k)
99+
yold(k) = y(k)
100+
zold(k) = z(k)
101+
x(k) = x(k) + v(1,k)*dtr
102+
y(k) = y(k) + v(2,k)*dtr
103+
z(k) = z(k) + v(3,k)*dtr
104+
end do
105+
if (use_rattle) call rattle (dtr,xold,yold,zold)
106+
if (use_rattle) call rattle2 (dtr)
107+
end do
108+
c
109+
c find the constraint-corrected full-step velocities
110+
c
111+
if (use_rattle) call rattle2 (dt)
112+
c
113+
c get the potential energy and atomic forces
114+
c
115+
call gradient (epot,derivs)
116+
c
117+
c use Newton's second law to get the next accelerations;
118+
c find the full-step velocities using the Verlet recursion
119+
c
120+
do i = 1, nuse
121+
k = iuse(i)
122+
do j = 1, 3
123+
a(j,k) = -ekcal * derivs(j,k) / mass(k)
124+
v(j,k) = v(j,k) + a(j,k)*dt_2
125+
end do
126+
end do
127+
c
128+
c find the constraint-corrected full-step velocities
129+
c
130+
if (use_rattle) call rattle2 (dt)
131+
c
132+
c update velocities with frictional and random components
133+
c
134+
call oprep (dt_2,vfric,vrand)
135+
do i = 1, nuse
136+
k = iuse(i)
137+
do j = 1, 3
138+
v(j,k) = v(j,k)*vfric(k) + vrand(j,k)
139+
end do
140+
end do
141+
if (use_rattle) call rattle2 (dt_2)
142+
c
143+
c perform deallocation of some local arrays
144+
c
145+
deallocate (xold)
146+
deallocate (yold)
147+
deallocate (zold)
148+
deallocate (derivs)
149+
deallocate (vfric)
150+
deallocate (vrand)
151+
c
152+
c find the constraint-corrected full-step velocities
153+
c
154+
if (use_rattle) call rattle2 (dt)
155+
c
156+
c compute the kinetic energy and control the pressure;
157+
c half-step kinetic energy gives better temperature control
158+
c
159+
c call kinetic (eksum,ekin,temp)
160+
call pressure (dt,epot,ekin,temp,pres,stress)
161+
call pressure2 (epot,temp)
162+
c
163+
c total energy is sum of kinetic and potential energies
164+
c
165+
etot = eksum + epot
166+
c
167+
c compute statistics and save trajectory for this step
168+
c
169+
call mdstat (istep,dt,etot,epot,eksum,temp,pres)
170+
call mdsave (istep,dt,epot,eksum)
171+
call mdrest (istep)
172+
return
173+
end

‎source/baoab.f‎

Lines changed: 74 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,7 @@ subroutine baoab (istep,dt)
3434
use moldyn
3535
use units
3636
use usage
37+
use virial
3738
implicit none
3839
integer i,j,k
3940
integer istep
@@ -42,8 +43,10 @@ subroutine baoab (istep,dt)
4243
real*8 etot,epot
4344
real*8 eksum
4445
real*8 temp,pres
46+
real*8 drattle
4547
real*8 ekin(3,3)
4648
real*8 stress(3,3)
49+
real*8 virrat(3,3)
4750
real*8, allocatable :: xold(:)
4851
real*8, allocatable :: yold(:)
4952
real*8, allocatable :: zold(:)
@@ -54,8 +57,10 @@ subroutine baoab (istep,dt)
5457
c
5558
c set some time values for the dynamics integration
5659
c
57-
dt_2 = 0.5d0 * dt
5860
nrattle = 1
61+
if (use_rattle) nrattle = 5
62+
drattle = dble(nrattle)
63+
dt_2 = 0.5d0 * dt
5964
dtr = dt_2 / dble(nrattle)
6065
c
6166
c perform dynamic allocation of some local arrays
@@ -64,23 +69,19 @@ subroutine baoab (istep,dt)
6469
allocate (yold(n))
6570
allocate (zold(n))
6671
allocate (derivs(3,n))
67-
allocate (vfric(n))
72+
allocate (vfric(n))
6873
allocate (vrand(3,n))
6974
c
70-
c find half-step velocities via the Verlet recursion
75+
c use a first B step to find the half-step velocities
7176
c
7277
do i = 1, nuse
7378
k = iuse(i)
7479
do j = 1, 3
7580
v(j,k) = v(j,k) + a(j,k)*dt_2
76-
end do
81+
end do
7782
end do
7883
c
79-
c find the constraint-corrected full-step velocities
80-
c
81-
if (use_rattle) call rattle2 (dt)
82-
c
83-
c take first A step according to the BAOAB sequence
84+
c take the first A step to get the half-step positions
8485
c
8586
do j = 1, nrattle
8687
do i = 1, nuse
@@ -92,15 +93,23 @@ subroutine baoab (istep,dt)
9293
y(k) = y(k) + v(2,k)*dtr
9394
z(k) = z(k) + v(3,k)*dtr
9495
end do
95-
if (use_rattle) call rattle (dtr,xold,yold,zold)
96-
if (use_rattle) call rattle2 (dtr)
97-
end do
96+
if (use_rattle) then
97+
call rattle (dtr,xold,yold,zold)
98+
call rattle2 (dtr)
99+
end if
100+
end do
98101
c
99-
c find the constraint-corrected full-step velocities
102+
c initialize virial for stochastic and constraint forces
100103
c
101-
if (use_rattle) call rattle2 (dt)
104+
if (use_rattle) then
105+
do i = 1, 3
106+
do j = 1, 3
107+
vir(j,i) = 0.0d0
108+
end do
109+
end do
110+
end if
102111
c
103-
c update velocities with frictional and random components
112+
c next an O step to get frictional and random components
104113
c
105114
call oprep (dt,vfric,vrand)
106115
do i = 1, nuse
@@ -109,9 +118,17 @@ subroutine baoab (istep,dt)
109118
v(j,k) = v(j,k)*vfric(k) + vrand(j,k)
110119
end do
111120
end do
112-
if (use_rattle) call rattle2 (dt)
121+
if (use_rattle) then
122+
call rattle2 (dt)
123+
do i = 1, 3
124+
do j = 1, 3
125+
virrat(j,i) = vir(j,i)
126+
vir(j,i) = 0.0d0
127+
end do
128+
end do
129+
end if
113130
c
114-
c take second A step according to the BAOAB sequence
131+
c take a second A step to get the full-step positions
115132
c
116133
do j = 1, nrattle
117134
do i = 1, nuse
@@ -123,9 +140,17 @@ subroutine baoab (istep,dt)
123140
y(k) = y(k) + v(2,k)*dtr
124141
z(k) = z(k) + v(3,k)*dtr
125142
end do
126-
if (use_rattle) call rattle (dtr,xold,yold,zold)
127-
if (use_rattle) call rattle2 (dtr)
128-
end do
143+
if (use_rattle) then
144+
call rattle (dtr,xold,yold,zold)
145+
call rattle2 (dtr)
146+
do i = 1, 3
147+
do k = 1, 3
148+
virrat(k,i) = virrat(k,i) + vir(k,i)/drattle
149+
vir(k,i) = 0.0d0
150+
end do
151+
end do
152+
end if
153+
end do
129154
c
130155
c get the potential energy and atomic forces
131156
c
@@ -136,8 +161,7 @@ subroutine baoab (istep,dt)
136161
call kinetic (eksum,ekin,temp)
137162
call pressure2 (epot,temp)
138163
c
139-
c use Newton's second law to get the next accelerations;
140-
c find the full-step velocities using the Verlet recursion
164+
c second B step for accelerations and full-step velocities
141165
c
142166
do i = 1, nuse
143167
k = iuse(i)
@@ -147,28 +171,45 @@ subroutine baoab (istep,dt)
147171
end do
148172
end do
149173
c
150-
c perform deallocation of some local arrays
151-
c
152-
deallocate (xold)
153-
deallocate (yold)
154-
deallocate (zold)
155-
deallocate (derivs)
156-
deallocate (vfric)
157-
deallocate (vrand)
158-
c
159174
c find the constraint-corrected full-step velocities
160175
c
161-
if (use_rattle) call rattle2 (dt)
176+
if (use_rattle) then
177+
do i = 1, nuse
178+
k = iuse(i)
179+
xold(k) = x(k)
180+
yold(k) = y(k)
181+
zold(k) = z(k)
182+
end do
183+
call rattle2 (dt)
184+
do i = 1, 3
185+
do j = 1, 3
186+
vir(j,i) = vir(j,i) + virrat(j,i)
187+
end do
188+
end do
189+
end if
162190
c
163191
c compute the kinetic energy and control the pressure;
164192
c half-step kinetic energy gives better temperature control
165193
c
166194
c call kinetic (eksum,ekin,temp)
167195
call pressure (dt,epot,ekin,temp,pres,stress)
168196
c
197+
c final constraint step to enforce position convergence
198+
c
199+
if (use_rattle) call shake (xold,yold,zold)
200+
c
201+
c perform deallocation of some local arrays
202+
c
203+
deallocate (xold)
204+
deallocate (yold)
205+
deallocate (zold)
206+
deallocate (derivs)
207+
deallocate (vfric)
208+
deallocate (vrand)
209+
c
169210
c total energy is sum of kinetic and potential energies
170211
c
171-
etot = eksum + epot
212+
etot = eksum + epot
172213
c
173214
c compute statistics and save trajectory for this step
174215
c

0 commit comments

Comments
 (0)