Skip to content

Commit 982f2d5

Browse files
jayponderjayponder
committed
add VRESPA and BRESPA options for MTS MD integrators
1 parent 6d90a1e commit 982f2d5

File tree

18 files changed

+378
-106
lines changed

18 files changed

+378
-106
lines changed

‎doc/tinker-guide.pdf‎

61 Bytes
Binary file not shown.

‎example/aplussmall.key‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ barostat MONTECARLO
88

99
neighbor-list
1010
list-buffer 0.2
11-
vdw-cutoff 9.0
11+
vdw-cutoff 8.5
1212
ewald
1313
ewald-cutoff 7.0
1414

‎example/watersmall.key‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ barostat MONTECARLO
88

99
neighbor-list
1010
list-buffer 0.2
11-
vdw-cutoff 9.0
11+
vdw-cutoff 8.5
1212
ewald
1313
ewald-cutoff 7.0
1414

‎interface/c/tinker/detail/moldyn.hh‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,8 @@ extern "C" {
88
extern double* TINKER_MOD(moldyn, v);
99
extern double* TINKER_MOD(moldyn, a);
1010
extern double* TINKER_MOD(moldyn, aalt);
11+
extern double* TINKER_MOD(moldyn, aslow);
12+
extern double* TINKER_MOD(moldyn, afast);
1113
#ifdef __cplusplus
1214
}
1315
#endif

‎interface/c/tinker/routines.h‎

Lines changed: 8 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -54,6 +54,8 @@ void attach_();
5454
// baoab.f
5555
void baoab_(int* istep, double* dt);
5656
#define tinker_f_baoab baoab_
57+
void obabo_(int* istep, double* dt);
58+
#define tinker_f_obabo obabo_
5759
void oprep_(double* dt, double* vfric, double* vrand);
5860
#define tinker_f_oprep oprep_
5961

@@ -1926,10 +1928,10 @@ void mdinit_(double* dt);
19261928
// mdrest.f
19271929
void mdrest_(int* istep);
19281930
#define tinker_f_mdrest mdrest_
1929-
void rgdrest_();
1930-
#define tinker_f_rgdrest rgdrest_
19311931
void xyzrest_();
19321932
#define tinker_f_xyzrest xyzrest_
1933+
void rgdrest_();
1934+
#define tinker_f_rgdrest rgdrest_
19331935

19341936
// mdsave.f
19351937
void mdsave_(int* istep, double* dt, double* epot, double* eksum);
@@ -2390,8 +2392,10 @@ void replica_(double* cutoff);
23902392
#define tinker_f_replica replica_
23912393

23922394
// respa.f
2393-
void respa_(int* istep, double* dt);
2394-
#define tinker_f_respa respa_
2395+
void vrespa_(int* istep, double* dt);
2396+
#define tinker_f_vrespa vrespa_
2397+
void brespa_(int* istep, double* dt);
2398+
#define tinker_f_brespa brespa_
23952399
void gradfast_(double* energy, double* derivs);
23962400
#define tinker_f_gradfast gradfast_
23972401
void gradslow_(double* energy, double* derivs);

‎interface/cpp/tinker/detail/moldyn.hh‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,14 +6,20 @@ namespace tinker { namespace moldyn {
66
extern double*& v;
77
extern double*& a;
88
extern double*& aalt;
9+
extern double*& aslow;
10+
extern double*& afast;
911

1012
#ifdef TINKER_FORTRAN_MODULE_CPP
1113
extern "C" double* TINKER_MOD(moldyn, v);
1214
extern "C" double* TINKER_MOD(moldyn, a);
1315
extern "C" double* TINKER_MOD(moldyn, aalt);
16+
extern "C" double* TINKER_MOD(moldyn, aslow);
17+
extern "C" double* TINKER_MOD(moldyn, afast);
1418

1519
double*& v = TINKER_MOD(moldyn, v);
1620
double*& a = TINKER_MOD(moldyn, a);
1721
double*& aalt = TINKER_MOD(moldyn, aalt);
22+
double*& aslow = TINKER_MOD(moldyn, aslow);
23+
double*& afast = TINKER_MOD(moldyn, afast);
1824
#endif
1925
} }

‎interface/cpp/tinker/routines.h‎

Lines changed: 8 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -54,6 +54,8 @@ void attach_();
5454
// baoab.f
5555
void baoab_(int* istep, double* dt);
5656
#define tinker_f_baoab baoab_
57+
void obabo_(int* istep, double* dt);
58+
#define tinker_f_obabo obabo_
5759
void oprep_(double* dt, double* vfric, double* vrand);
5860
#define tinker_f_oprep oprep_
5961

@@ -1926,10 +1928,10 @@ void mdinit_(double* dt);
19261928
// mdrest.f
19271929
void mdrest_(int* istep);
19281930
#define tinker_f_mdrest mdrest_
1929-
void rgdrest_();
1930-
#define tinker_f_rgdrest rgdrest_
19311931
void xyzrest_();
19321932
#define tinker_f_xyzrest xyzrest_
1933+
void rgdrest_();
1934+
#define tinker_f_rgdrest rgdrest_
19331935

19341936
// mdsave.f
19351937
void mdsave_(int* istep, double* dt, double* epot, double* eksum);
@@ -2390,8 +2392,10 @@ void replica_(double* cutoff);
23902392
#define tinker_f_replica replica_
23912393

23922394
// respa.f
2393-
void respa_(int* istep, double* dt);
2394-
#define tinker_f_respa respa_
2395+
void vrespa_(int* istep, double* dt);
2396+
#define tinker_f_vrespa vrespa_
2397+
void brespa_(int* istep, double* dt);
2398+
#define tinker_f_brespa brespa_
23952399
void gradfast_(double* energy, double* derivs);
23962400
#define tinker_f_gradfast gradfast_
23972401
void gradslow_(double* energy, double* derivs);

‎source/anneal.f‎

Lines changed: 8 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -238,8 +238,10 @@ program anneal
238238
call ghmcstep (istep,dt)
239239
else if (integrate .eq. 'RIGIDBODY') then
240240
call rgdstep (istep,dt)
241-
else if (integrate .eq. 'RESPA') then
242-
call respa (istep,dt)
241+
else if (integrate .eq. 'VRESPA') then
242+
call vrespa (istep,dt)
243+
else if (integrate .eq. 'BRESPA') then
244+
call brespa (istep,dt)
243245
else
244246
call beeman (istep,dt)
245247
end if
@@ -295,8 +297,10 @@ program anneal
295297
call ghmcstep (istep,dt)
296298
else if (integrate .eq. 'RIGIDBODY') then
297299
call rgdstep (istep,dt)
298-
else if (integrate .eq. 'RESPA') then
299-
call respa (istep,dt)
300+
else if (integrate .eq. 'VRESPA') then
301+
call vrespa (istep,dt)
302+
else if (integrate .eq. 'BRESPA') then
303+
call brespa (istep,dt)
300304
else
301305
call beeman (istep,dt)
302306
end if

‎source/baoab.f‎

Lines changed: 32 additions & 53 deletions
Original file line numberDiff line numberDiff line change
@@ -68,9 +68,9 @@ subroutine baoab (istep,dt)
6868
allocate (xold(n))
6969
allocate (yold(n))
7070
allocate (zold(n))
71-
allocate (derivs(3,n))
7271
allocate (vfric(n))
7372
allocate (vrand(3,n))
73+
allocate (derivs(3,n))
7474
c
7575
c use a first B step to find the half-step velocities
7676
c
@@ -96,20 +96,15 @@ subroutine baoab (istep,dt)
9696
if (use_rattle) then
9797
call rattle (dtr,xold,yold,zold)
9898
call rattle2 (dtr)
99+
do i = 1, 3
100+
do k = 1, 3
101+
vir(k,i) = 0.0d0
102+
end do
103+
end do
99104
end if
100105
end do
101106
c
102-
c initialize virial for stochastic and constraint forces
103-
c
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
111-
c
112-
c next an O step to get frictional and random components
107+
c use an O step to get frictional and random components
113108
c
114109
call oprep (dt,vfric,vrand)
115110
do i = 1, nuse
@@ -169,22 +164,19 @@ subroutine baoab (istep,dt)
169164
v(j,k) = v(j,k) + a(j,k)*dt_2
170165
end do
171166
end do
172-
c
173-
c find the constraint-corrected full-step velocities
174-
c
175167
if (use_rattle) then
176-
do i = 1, nuse
177-
k = iuse(i)
178-
xold(k) = x(k)
179-
yold(k) = y(k)
180-
zold(k) = z(k)
181-
end do
182168
call rattle2 (dt)
183169
do i = 1, 3
184170
do j = 1, 3
185171
vir(j,i) = vir(j,i) + virrat(j,i)
186172
end do
187173
end do
174+
do i = 1, nuse
175+
k = iuse(i)
176+
xold(k) = x(k)
177+
yold(k) = y(k)
178+
zold(k) = z(k)
179+
end do
188180
end if
189181
c
190182
c compute full-step kinetic energy and pressure correction;
@@ -203,9 +195,9 @@ subroutine baoab (istep,dt)
203195
deallocate (xold)
204196
deallocate (yold)
205197
deallocate (zold)
206-
deallocate (derivs)
207198
deallocate (vfric)
208199
deallocate (vrand)
200+
deallocate (derivs)
209201
c
210202
c total energy is sum of kinetic and potential energies
211203
c
@@ -281,21 +273,11 @@ subroutine obabo (istep,dt)
281273
allocate (xold(n))
282274
allocate (yold(n))
283275
allocate (zold(n))
284-
allocate (derivs(3,n))
285276
allocate (vfric(n))
286277
allocate (vrand(3,n))
278+
allocate (derivs(3,n))
287279
c
288-
c initialize virial for stochastic and constraint forces
289-
c
290-
if (use_rattle) then
291-
do i = 1, 3
292-
do j = 1, 3
293-
vir(j,i) = 0.0d0
294-
end do
295-
end do
296-
end if
297-
c
298-
c update velocities with frictional and random components
280+
c take first O step for frictional and random components
299281
c
300282
call oprep (dt_2,vfric,vrand)
301283
do i = 1, nuse
@@ -305,6 +287,11 @@ subroutine obabo (istep,dt)
305287
end do
306288
end do
307289
if (use_rattle) then
290+
do i = 1, 3
291+
do j = 1, 3
292+
vir(j,i) = 0.0d0
293+
end do
294+
end do
308295
call rattle2 (dt_2)
309296
do i = 1, 3
310297
do j = 1, 3
@@ -314,7 +301,7 @@ subroutine obabo (istep,dt)
314301
end do
315302
end if
316303
c
317-
c find half-step velocities via the Verlet recursion
304+
c use a first B step to find the half-step velocities
318305
c
319306
do i = 1, nuse
320307
k = iuse(i)
@@ -323,7 +310,7 @@ subroutine obabo (istep,dt)
323310
end do
324311
end do
325312
c
326-
c take first A step according to the BAOAB sequence
313+
c take full-step A step according to the BAOAB sequence
327314
c
328315
do j = 1, nrattle
329316
do i = 1, nuse
@@ -355,8 +342,7 @@ subroutine obabo (istep,dt)
355342
c
356343
c call kinetic (eksum,ekin,temp)
357344
c
358-
c use Newton's second law to get the next accelerations;
359-
c find the full-step velocities using the Verlet recursion
345+
c second B step for accelerations and full-step velocities
360346
c
361347
do i = 1, nuse
362348
k = iuse(i)
@@ -365,12 +351,9 @@ subroutine obabo (istep,dt)
365351
v(j,k) = v(j,k) + a(j,k)*dt_2
366352
end do
367353
end do
368-
c
369-
c find the constraint-corrected full-step velocities
370-
c
371354
if (use_rattle) call rattle2 (dt)
372355
c
373-
c update velocities with frictional and random components
356+
c use second O step for frictional and random components
374357
c
375358
call oprep (dt_2,vfric,vrand)
376359
do i = 1, nuse
@@ -379,23 +362,19 @@ subroutine obabo (istep,dt)
379362
v(j,k) = v(j,k)*vfric(k) + vrand(j,k)
380363
end do
381364
end do
382-
if (use_rattle) call rattle2 (dt_2)
383-
c
384-
c find the constraint-corrected full-step velocities
385-
c
386365
if (use_rattle) then
366+
call rattle2 (dt_2)
367+
do i = 1, 3
368+
do j = 1, 3
369+
vir(j,i) = vir(j,i) + virrat(j,i)
370+
end do
371+
end do
387372
do i = 1, nuse
388373
k = iuse(i)
389374
xold(k) = x(k)
390375
yold(k) = y(k)
391376
zold(k) = z(k)
392377
end do
393-
call rattle2 (dt)
394-
do i = 1, 3
395-
do j = 1, 3
396-
vir(j,i) = vir(j,i) + virrat(j,i)
397-
end do
398-
end do
399378
end if
400379
c
401380
c compute the kinetic energy and control the pressure
@@ -413,9 +392,9 @@ subroutine obabo (istep,dt)
413392
deallocate (xold)
414393
deallocate (yold)
415394
deallocate (zold)
416-
deallocate (derivs)
417395
deallocate (vfric)
418396
deallocate (vrand)
397+
deallocate (derivs)
419398
c
420399
c total energy is sum of kinetic and potential energies
421400
c

‎source/beeman.f‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -39,7 +39,7 @@ subroutine beeman (istep,dt)
3939
implicit none
4040
integer i,j,k
4141
integer istep
42-
real*8 dt,dt_x,factor
42+
real*8 dt,dtx,dmix
4343
real*8 etot,eksum,epot
4444
real*8 temp,pres
4545
real*8 part1,part2
@@ -54,10 +54,10 @@ subroutine beeman (istep,dt)
5454
c
5555
c set time values and coefficients for Beeman integration
5656
c
57-
factor = dble(bmnmix)
58-
dt_x = dt / factor
59-
part1 = 0.5d0*factor + 1.0d0
57+
dmix = dble(bmnmix)
58+
part1 = 0.5d0*dmix + 1.0d0
6059
part2 = part1 - 2.0d0
60+
dtx = dt / dmix
6161
c
6262
c perform dynamic allocation of some local arrays
6363
c
@@ -72,7 +72,7 @@ subroutine beeman (istep,dt)
7272
do i = 1, nuse
7373
k = iuse(i)
7474
do j = 1, 3
75-
v(j,k) = v(j,k) + (part1*a(j,k)-aalt(j,k))*dt_x
75+
v(j,k) = v(j,k) + (part1*a(j,k)-aalt(j,k))*dtx
7676
end do
7777
xold(k) = x(k)
7878
yold(k) = y(k)
@@ -118,7 +118,7 @@ subroutine beeman (istep,dt)
118118
do j = 1, 3
119119
aalt(j,k) = a(j,k)
120120
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
121+
v(j,k) = v(j,k) + (part2*a(j,k)+aalt(j,k))*dtx
122122
end do
123123
end do
124124
c

0 commit comments

Comments
 (0)