1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * GROwing Monsters And Cloning Shrimps
54 #include "gmx_random.h"
68 #include "gmx_wallcycle.h"
70 /*For debugging, start at v(-dt/2) for velolcity verlet -- uncomment next line */
71 /*#define STARTFROMDT2*/
91 /* The random state */
97 gmx_sd_sigma_t
*sdsig
;
102 typedef struct gmx_update
107 /* Variables for the deform algorithm */
108 gmx_large_int_t deformref_step
;
109 matrix deformref_box
;
113 void store_rvec(rvec
*from
, rvec
*to
, int n
) {
116 copy_rvec(from
[i
],to
[i
]);
120 static void do_update_md(int start
,int nrend
,double dt
,
121 t_grp_tcstat
*tcstat
,t_grp_acc
*gstat
,double nh_vxi
[],
122 rvec accel
[],ivec nFreeze
[],real invmass
[],
123 unsigned short ptype
[],unsigned short cFREEZE
[],
124 unsigned short cACC
[],unsigned short cTC
[],
125 rvec x
[],rvec xprime
[],rvec v
[],
132 real vn
,vv
,va
,vb
,vnrel
;
138 /* Update with coupling to extended ensembles, used for
139 * Nose-Hoover and Parrinello-Rahman coupling
140 * Nose-Hoover uses the reversible leap-frog integrator from
141 * Holian et al. Phys Rev E 52(3) : 2338, 1995
143 for(n
=start
; n
<nrend
; n
++)
158 lg
= tcstat
[gt
].lambda
;
162 rvec_sub(v
[n
],gstat
[ga
].u
,vrel
);
166 if((ptype
[n
] != eptVSite
) && (ptype
[n
] != eptShell
) && !nFreeze
[gf
][d
])
168 vnrel
= (lg
*vrel
[d
] + dt
*(imass
*f
[n
][d
] - 0.5*vxi
*vrel
[d
]
169 - iprod(M
[d
],vrel
)))/(1 + 0.5*vxi
*dt
);
170 /* do not scale the mean velocities u */
171 vn
= gstat
[ga
].u
[d
] + accel
[ga
][d
]*dt
+ vnrel
;
173 xprime
[n
][d
] = x
[n
][d
]+vn
*dt
;
178 xprime
[n
][d
] = x
[n
][d
];
185 /* Classic version of update, used with berendsen coupling */
186 for(n
=start
; n
<nrend
; n
++)
188 w_dt
= invmass
[n
]*dt
;
201 lg
= tcstat
[gt
].lambda
;
206 if((ptype
[n
] != eptVSite
) && (ptype
[n
] != eptShell
) && !nFreeze
[gf
][d
])
208 vv
= lg
*vn
+ f
[n
][d
]*w_dt
;
210 /* do not scale the mean velocities u */
212 va
= vv
+ accel
[ga
][d
]*dt
;
213 vb
= va
+ (1.0-lg
)*u
;
215 xprime
[n
][d
] = x
[n
][d
]+vb
*dt
;
220 xprime
[n
][d
] = x
[n
][d
];
227 static void do_update_vv_vel(int start
,int nrend
,double dt
,
228 t_grp_tcstat
*tcstat
,t_grp_acc
*gstat
,
229 rvec accel
[],ivec nFreeze
[],real invmass
[],
230 unsigned short ptype
[],
231 unsigned short cFREEZE
[],
232 unsigned short cACC
[],unsigned short cTC
[],
233 rvec v
[],rvec f
[],bool bExtended
, real veta
, real alpha
)
238 real vn
,vv
,va
,vb
,vnrel
;
243 g
= 0.25*dt
*veta
*alpha
;
245 mv2
= series_sinhx(g
);
247 for(n
=start
; n
<nrend
; n
++)
249 w_dt
= invmass
[n
]*dt
;
262 lg
= tcstat
[gt
].lambda
;
263 rvec_sub(v
[n
],gstat
[ga
].u
,vrel
);
267 if((ptype
[n
] != eptVSite
) && (ptype
[n
] != eptShell
) && !nFreeze
[gf
][d
])
271 vnrel
= mv1
*(mv1
*vrel
[d
] + 0.5*w_dt
*mv2
*f
[n
][d
]);
272 /* do not scale the mean velocities u: MRS: may not be done correctly currently? */
273 vn
= gstat
[ga
].u
[d
] + accel
[ga
][d
]*dt
+ vnrel
;
278 vv
= lg
*v
[n
][d
] + 0.5*f
[n
][d
]*w_dt
;
279 /* do not scale the mean velocities u MRS: may not be done correctly currently? */
281 va
= vv
+ 0.5*accel
[ga
][d
]*dt
;
282 vb
= va
+ (1.0-lg
)*u
;
292 } /* do_update_vv_vel */
294 static void do_update_vv_pos(int start
,int nrend
,double dt
,
295 t_grp_tcstat
*tcstat
,t_grp_acc
*gstat
,
296 rvec accel
[],ivec nFreeze
[],real invmass
[],
297 unsigned short ptype
[],
298 unsigned short cFREEZE
[],
299 unsigned short cACC
[],unsigned short cTC
[],
300 rvec x
[],rvec xprime
[],rvec v
[],
301 rvec f
[],bool bExtended
, real veta
, real alpha
)
306 real vn
,vv
,va
,vb
,vnrel
;
309 double g
,mr1
,mr2
,mv1
,mv2
;
314 mr2
= series_sinhx(g
);
322 for(n
=start
; n
<nrend
; n
++) {
330 if ((ptype
[n
] != eptVSite
) && (ptype
[n
] != eptShell
) && !nFreeze
[gf
][d
])
332 xprime
[n
][d
] = mr1
*(mr1
*x
[n
][d
]+mr2
*dt
*v
[n
][d
]);
336 xprime
[n
][d
] = x
[n
][d
];
340 }/* do_update_vv_pos */
342 static void do_update_visc(int start
,int nrend
,double dt
,
343 t_grp_tcstat
*tcstat
,real invmass
[],double nh_vxi
[],
344 unsigned short ptype
[],unsigned short cTC
[],
345 rvec x
[],rvec xprime
[],rvec v
[],
346 rvec f
[],matrix M
,matrix box
,real
358 fac
= 2*M_PI
/(box
[ZZ
][ZZ
]);
361 /* Update with coupling to extended ensembles, used for
362 * Nose-Hoover and Parrinello-Rahman coupling
364 for(n
=start
; n
<nrend
; n
++) {
370 lg
= tcstat
[gt
].lambda
;
371 cosz
= cos(fac
*x
[n
][ZZ
]);
373 copy_rvec(v
[n
],vrel
);
385 if((ptype
[n
] != eptVSite
) && (ptype
[n
] != eptShell
))
387 vn
= (lg
*vrel
[d
] + dt
*(imass
*f
[n
][d
] - 0.5*vxi
*vrel
[d
]
388 - iprod(M
[d
],vrel
)))/(1 + 0.5*vxi
*dt
);
391 vn
+= vc
+ dt
*cosz
*cos_accel
;
394 xprime
[n
][d
] = x
[n
][d
]+vn
*dt
;
398 xprime
[n
][d
] = x
[n
][d
];
405 /* Classic version of update, used with berendsen coupling */
406 for(n
=start
; n
<nrend
; n
++)
408 w_dt
= invmass
[n
]*dt
;
413 lg
= tcstat
[gt
].lambda
;
414 cosz
= cos(fac
*x
[n
][ZZ
]);
420 if((ptype
[n
] != eptVSite
) && (ptype
[n
] != eptShell
))
425 /* Do not scale the cosine velocity profile */
426 vv
= vc
+ lg
*(vn
- vc
+ f
[n
][d
]*w_dt
);
427 /* Add the cosine accelaration profile */
428 vv
+= dt
*cosz
*cos_accel
;
432 vv
= lg
*(vn
+ f
[n
][d
]*w_dt
);
435 xprime
[n
][d
] = x
[n
][d
]+vv
*dt
;
440 xprime
[n
][d
] = x
[n
][d
];
447 static gmx_stochd_t
*init_stochd(FILE *fplog
,t_inputrec
*ir
)
456 /* Initiate random number generator for langevin type dynamics,
457 * for BD, SD or velocity rescaling temperature coupling.
459 sd
->gaussrand
= gmx_rng_init(ir
->ld_seed
);
461 ngtc
= ir
->opts
.ngtc
;
463 if (ir
->eI
== eiBD
) {
464 snew(sd
->bd_rf
,ngtc
);
465 } else if (EI_SD(ir
->eI
)) {
467 snew(sd
->sdsig
,ngtc
);
470 for(n
=0; n
<ngtc
; n
++) {
471 sdc
[n
].gdt
= ir
->delta_t
/ir
->opts
.tau_t
[n
];
472 sdc
[n
].eph
= exp(sdc
[n
].gdt
/2);
473 sdc
[n
].emh
= exp(-sdc
[n
].gdt
/2);
474 sdc
[n
].em
= exp(-sdc
[n
].gdt
);
475 if (sdc
[n
].gdt
>= 0.05) {
476 sdc
[n
].b
= sdc
[n
].gdt
*(sdc
[n
].eph
*sdc
[n
].eph
- 1)
477 - 4*(sdc
[n
].eph
- 1)*(sdc
[n
].eph
- 1);
478 sdc
[n
].c
= sdc
[n
].gdt
- 3 + 4*sdc
[n
].emh
- sdc
[n
].em
;
479 sdc
[n
].d
= 2 - sdc
[n
].eph
- sdc
[n
].emh
;
482 /* Seventh order expansions for small y */
483 sdc
[n
].b
= y
*y
*y
*y
*(1/3.0+y
*(1/3.0+y
*(17/90.0+y
*7/9.0)));
484 sdc
[n
].c
= y
*y
*y
*(2/3.0+y
*(-1/2.0+y
*(7/30.0+y
*(-1/12.0+y
*31/1260.0))));
485 sdc
[n
].d
= y
*y
*(-1+y
*y
*(-1/12.0-y
*y
/360.0));
488 fprintf(debug
,"SD const tc-grp %d: b %g c %g d %g\n",
489 n
,sdc
[n
].b
,sdc
[n
].c
,sdc
[n
].d
);
496 void get_stochd_state(gmx_update_t upd
,t_state
*state
)
498 gmx_rng_get_state(upd
->sd
->gaussrand
,state
->ld_rng
,state
->ld_rngi
);
501 void set_stochd_state(gmx_update_t upd
,t_state
*state
)
503 gmx_rng_set_state(upd
->sd
->gaussrand
,state
->ld_rng
,state
->ld_rngi
[0]);
506 gmx_update_t
init_update(FILE *fplog
,t_inputrec
*ir
)
512 if (ir
->eI
== eiBD
|| EI_SD(ir
->eI
) || ir
->etc
== etcVRESCALE
)
514 upd
->sd
= init_stochd(fplog
,ir
);
523 static void do_update_sd1(gmx_stochd_t
*sd
,
524 int start
,int homenr
,double dt
,
525 rvec accel
[],ivec nFreeze
[],
526 real invmass
[],unsigned short ptype
[],
527 unsigned short cFREEZE
[],unsigned short cACC
[],
528 unsigned short cTC
[],
529 rvec x
[],rvec xprime
[],rvec v
[],rvec f
[],
531 int ngtc
,real tau_t
[],real ref_t
[])
543 if (homenr
> sd
->sd_V_nalloc
)
545 sd
->sd_V_nalloc
= over_alloc_dd(homenr
);
546 srenew(sd
->sd_V
,sd
->sd_V_nalloc
);
548 gaussrand
= sd
->gaussrand
;
550 for(n
=0; n
<ngtc
; n
++)
553 /* The mass is encounted for later, since this differs per atom */
554 sig
[n
].V
= sqrt(2*kT
*(1 - sdc
[n
].em
));
557 for(n
=start
; n
<start
+homenr
; n
++)
559 ism
= sqrt(invmass
[n
]);
575 if((ptype
[n
] != eptVSite
) && (ptype
[n
] != eptShell
) && !nFreeze
[gf
][d
])
577 sd_V
= ism
*sig
[gt
].V
*gmx_rng_gaussian_table(gaussrand
);
579 v
[n
][d
] = v
[n
][d
]*sdc
[gt
].em
580 + (invmass
[n
]*f
[n
][d
] + accel
[ga
][d
])*tau_t
[gt
]*(1 - sdc
[gt
].em
)
583 xprime
[n
][d
] = x
[n
][d
] + v
[n
][d
]*dt
;
588 xprime
[n
][d
] = x
[n
][d
];
594 static void do_update_sd2(gmx_stochd_t
*sd
,bool bInitStep
,
595 int start
,int homenr
,
596 rvec accel
[],ivec nFreeze
[],
597 real invmass
[],unsigned short ptype
[],
598 unsigned short cFREEZE
[],unsigned short cACC
[],
599 unsigned short cTC
[],
600 rvec x
[],rvec xprime
[],rvec v
[],rvec f
[],
602 int ngtc
,real tau_t
[],real ref_t
[],
607 /* The random part of the velocity update, generated in the first
608 * half of the update, needs to be remembered for the second half.
620 if (homenr
> sd
->sd_V_nalloc
)
622 sd
->sd_V_nalloc
= over_alloc_dd(homenr
);
623 srenew(sd
->sd_V
,sd
->sd_V_nalloc
);
626 gaussrand
= sd
->gaussrand
;
630 for (n
=0; n
<ngtc
; n
++)
633 /* The mass is encounted for later, since this differs per atom */
634 sig
[n
].V
= sqrt(kT
*(1-sdc
[n
].em
));
635 sig
[n
].X
= sqrt(kT
*sqr(tau_t
[n
])*sdc
[n
].c
);
636 sig
[n
].Yv
= sqrt(kT
*sdc
[n
].b
/sdc
[n
].c
);
637 sig
[n
].Yx
= sqrt(kT
*sqr(tau_t
[n
])*sdc
[n
].b
/(1-sdc
[n
].em
));
641 for (n
=start
; n
<start
+homenr
; n
++)
643 ism
= sqrt(invmass
[n
]);
663 if((ptype
[n
] != eptVSite
) && (ptype
[n
] != eptShell
) && !nFreeze
[gf
][d
])
669 sd_X
[n
][d
] = ism
*sig
[gt
].X
*gmx_rng_gaussian_table(gaussrand
);
671 Vmh
= sd_X
[n
][d
]*sdc
[gt
].d
/(tau_t
[gt
]*sdc
[gt
].c
)
672 + ism
*sig
[gt
].Yv
*gmx_rng_gaussian_table(gaussrand
);
673 sd_V
[n
-start
][d
] = ism
*sig
[gt
].V
*gmx_rng_gaussian_table(gaussrand
);
675 v
[n
][d
] = vn
*sdc
[gt
].em
676 + (invmass
[n
]*f
[n
][d
] + accel
[ga
][d
])*tau_t
[gt
]*(1 - sdc
[gt
].em
)
677 + sd_V
[n
-start
][d
] - sdc
[gt
].em
*Vmh
;
679 xprime
[n
][d
] = x
[n
][d
] + v
[n
][d
]*tau_t
[gt
]*(sdc
[gt
].eph
- sdc
[gt
].emh
);
684 /* Correct the velocities for the constraints.
685 * This operation introduces some inaccuracy,
686 * since the velocity is determined from differences in coordinates.
689 (xprime
[n
][d
] - x
[n
][d
])/(tau_t
[gt
]*(sdc
[gt
].eph
- sdc
[gt
].emh
));
691 Xmh
= sd_V
[n
-start
][d
]*tau_t
[gt
]*sdc
[gt
].d
/(sdc
[gt
].em
-1)
692 + ism
*sig
[gt
].Yx
*gmx_rng_gaussian_table(gaussrand
);
693 sd_X
[n
][d
] = ism
*sig
[gt
].X
*gmx_rng_gaussian_table(gaussrand
);
695 xprime
[n
][d
] += sd_X
[n
][d
] - Xmh
;
704 xprime
[n
][d
] = x
[n
][d
];
711 static void do_update_bd(int start
,int nrend
,double dt
,
713 real invmass
[],unsigned short ptype
[],
714 unsigned short cFREEZE
[],unsigned short cTC
[],
715 rvec x
[],rvec xprime
[],rvec v
[],
716 rvec f
[],real friction_coefficient
,
717 int ngtc
,real tau_t
[],real ref_t
[],
718 real
*rf
,gmx_rng_t gaussrand
)
720 /* note -- these appear to be full step velocities . . . */
726 if (friction_coefficient
!= 0)
728 invfr
= 1.0/friction_coefficient
;
729 for(n
=0; n
<ngtc
; n
++)
731 rf
[n
] = sqrt(2.0*BOLTZ
*ref_t
[n
]/(friction_coefficient
*dt
));
736 for(n
=0; n
<ngtc
; n
++)
738 rf
[n
] = sqrt(2.0*BOLTZ
*ref_t
[n
]);
741 for(n
=start
; (n
<nrend
); n
++)
751 for(d
=0; (d
<DIM
); d
++)
753 if((ptype
[n
]!=eptVSite
) && (ptype
[n
]!=eptShell
) && !nFreeze
[gf
][d
])
755 if (friction_coefficient
!= 0) {
756 vn
= invfr
*f
[n
][d
] + rf
[gt
]*gmx_rng_gaussian_table(gaussrand
);
760 /* NOTE: invmass = 1/(mass*friction_constant*dt) */
761 vn
= invmass
[n
]*f
[n
][d
]*dt
762 + sqrt(invmass
[n
])*rf
[gt
]*gmx_rng_gaussian_table(gaussrand
);
765 xprime
[n
][d
] = x
[n
][d
]+vn
*dt
;
771 xprime
[n
][d
] = x
[n
][d
];
777 static void dump_it_all(FILE *fp
,const char *title
,
778 int natoms
,rvec x
[],rvec xp
[],rvec v
[],rvec f
[])
783 fprintf(fp
,"%s\n",title
);
784 pr_rvecs(fp
,0,"x",x
,natoms
);
785 pr_rvecs(fp
,0,"xp",xp
,natoms
);
786 pr_rvecs(fp
,0,"v",v
,natoms
);
787 pr_rvecs(fp
,0,"f",f
,natoms
);
792 static void calc_ke_part_normal(rvec v
[], t_grpopts
*opts
,t_mdatoms
*md
,
793 gmx_ekindata_t
*ekind
,t_nrnb
*nrnb
,bool bEkinAveVel
,
796 int start
=md
->start
,homenr
=md
->homenr
;
797 int g
,d
,n
,m
,ga
=0,gt
=0;
800 t_grp_tcstat
*tcstat
=ekind
->tcstat
;
801 t_grp_acc
*grpstat
=ekind
->grpstat
;
804 /* three main: VV with AveVel, vv with AveEkin, leap with AveEkin. Leap with AveVel is also
805 an option, but not supported now. Additionally, if we are doing iterations.
806 bEkinAveVel: If TRUE, we sum into ekin, if FALSE, into ekinh.
807 bSavEkinOld: If TRUE (in the case of iteration = bIterate is TRUE), we don't copy over the ekinh_old.
808 If FALSE, we overrwrite it.
811 /* group velocities are calculated in update_ekindata and
812 * accumulated in acumulate_groups.
813 * Now the partial global and groups ekin.
815 for(g
=0; (g
<opts
->ngtc
); g
++)
819 copy_mat(tcstat
[g
].ekinh
,tcstat
[g
].ekinh_old
);
822 clear_mat(tcstat
[g
].ekinf
);
824 clear_mat(tcstat
[g
].ekinh
);
827 tcstat
[g
].ekinscalef_nhc
= 1.0; /* need to clear this -- logic is complicated! */
830 ekind
->dekindl_old
= ekind
->dekindl
;
833 for(n
=start
; (n
<start
+homenr
); n
++)
843 hm
= 0.5*md
->massT
[n
];
845 for(d
=0; (d
<DIM
); d
++)
847 v_corrt
[d
] = v
[n
][d
] - grpstat
[ga
].u
[d
];
849 for(d
=0; (d
<DIM
); d
++)
851 for (m
=0;(m
<DIM
); m
++)
853 /* if we're computing a full step velocity, v_corrt[d] has v(t). Otherwise, v(t+dt/2) */
856 tcstat
[gt
].ekinf
[m
][d
]+=hm
*v_corrt
[m
]*v_corrt
[d
];
860 tcstat
[gt
].ekinh
[m
][d
]+=hm
*v_corrt
[m
]*v_corrt
[d
];
864 if (md
->nMassPerturbed
&& md
->bPerturbed
[n
])
866 dekindl
-= 0.5*(md
->massB
[n
] - md
->massA
[n
])*iprod(v_corrt
,v_corrt
);
869 ekind
->dekindl
= dekindl
;
870 inc_nrnb(nrnb
,eNR_EKIN
,homenr
);
873 static void calc_ke_part_visc(matrix box
,rvec x
[],rvec v
[],
874 t_grpopts
*opts
,t_mdatoms
*md
,
875 gmx_ekindata_t
*ekind
,
876 t_nrnb
*nrnb
, bool bEkinAveVel
, bool bSaveEkinOld
)
878 int start
=md
->start
,homenr
=md
->homenr
;
882 t_grp_tcstat
*tcstat
=ekind
->tcstat
;
883 t_cos_acc
*cosacc
=&(ekind
->cosacc
);
888 for(g
=0; g
<opts
->ngtc
; g
++)
890 copy_mat(ekind
->tcstat
[g
].ekinh
,ekind
->tcstat
[g
].ekinh_old
);
891 clear_mat(ekind
->tcstat
[g
].ekinh
);
893 ekind
->dekindl_old
= ekind
->dekindl
;
895 fac
= 2*M_PI
/box
[ZZ
][ZZ
];
898 for(n
=start
; n
<start
+homenr
; n
++)
904 hm
= 0.5*md
->massT
[n
];
906 /* Note that the times of x and v differ by half a step */
907 /* MRS -- would have to be changed for VV */
908 cosz
= cos(fac
*x
[n
][ZZ
]);
909 /* Calculate the amplitude of the new velocity profile */
910 mvcos
+= 2*cosz
*md
->massT
[n
]*v
[n
][XX
];
912 copy_rvec(v
[n
],v_corrt
);
913 /* Subtract the profile for the kinetic energy */
914 v_corrt
[XX
] -= cosz
*cosacc
->vcos
;
915 for (d
=0; (d
<DIM
); d
++)
917 for (m
=0; (m
<DIM
); m
++)
919 /* if we're computing a full step velocity, v_corrt[d] has v(t). Otherwise, v(t+dt/2) */
922 tcstat
[gt
].ekinf
[m
][d
]+=hm
*v_corrt
[m
]*v_corrt
[d
];
926 tcstat
[gt
].ekinh
[m
][d
]+=hm
*v_corrt
[m
]*v_corrt
[d
];
930 if(md
->nPerturbed
&& md
->bPerturbed
[n
])
932 dekindl
-= 0.5*(md
->massB
[n
] - md
->massA
[n
])*iprod(v_corrt
,v_corrt
);
935 ekind
->dekindl
= dekindl
;
936 cosacc
->mvcos
= mvcos
;
938 inc_nrnb(nrnb
,eNR_EKIN
,homenr
);
941 void calc_ke_part(t_state
*state
,t_grpopts
*opts
,t_mdatoms
*md
,
942 gmx_ekindata_t
*ekind
,t_nrnb
*nrnb
, bool bEkinAveVel
, bool bSaveEkinOld
)
944 if (ekind
->cosacc
.cos_accel
== 0)
946 calc_ke_part_normal(state
->v
,opts
,md
,ekind
,nrnb
,bEkinAveVel
,bSaveEkinOld
);
950 calc_ke_part_visc(state
->box
,state
->x
,state
->v
,opts
,md
,ekind
,nrnb
,bEkinAveVel
,bSaveEkinOld
);
954 void init_ekinstate(ekinstate_t
*ekinstate
,t_inputrec
*ir
)
956 ekinstate
->ekin_n
= ir
->opts
.ngtc
;
957 snew(ekinstate
->ekinh
,ekinstate
->ekin_n
);
958 snew(ekinstate
->ekinf
,ekinstate
->ekin_n
);
959 snew(ekinstate
->ekinh_old
,ekinstate
->ekin_n
);
960 snew(ekinstate
->ekinscalef_nhc
,ekinstate
->ekin_n
);
961 snew(ekinstate
->ekinscaleh_nhc
,ekinstate
->ekin_n
);
962 snew(ekinstate
->vscale_nhc
,ekinstate
->ekin_n
);
963 ekinstate
->dekindl
= 0;
964 ekinstate
->mvcos
= 0;
967 void update_ekinstate(ekinstate_t
*ekinstate
,gmx_ekindata_t
*ekind
)
971 for(i
=0;i
<ekinstate
->ekin_n
;i
++)
973 copy_mat(ekind
->tcstat
[i
].ekinh
,ekinstate
->ekinh
[i
]);
974 copy_mat(ekind
->tcstat
[i
].ekinf
,ekinstate
->ekinf
[i
]);
975 copy_mat(ekind
->tcstat
[i
].ekinh_old
,ekinstate
->ekinh_old
[i
]);
976 ekinstate
->ekinscalef_nhc
[i
] = ekind
->tcstat
[i
].ekinscalef_nhc
;
977 ekinstate
->ekinscaleh_nhc
[i
] = ekind
->tcstat
[i
].ekinscaleh_nhc
;
978 ekinstate
->vscale_nhc
[i
] = ekind
->tcstat
[i
].vscale_nhc
;
981 copy_mat(ekind
->ekin
,ekinstate
->ekin_total
);
982 ekinstate
->dekindl
= ekind
->dekindl
;
983 ekinstate
->mvcos
= ekind
->cosacc
.mvcos
;
987 void restore_ekinstate_from_state(t_commrec
*cr
,
988 gmx_ekindata_t
*ekind
,ekinstate_t
*ekinstate
)
994 for(i
=0;i
<ekinstate
->ekin_n
;i
++)
996 copy_mat(ekinstate
->ekinh
[i
],ekind
->tcstat
[i
].ekinh
);
997 copy_mat(ekinstate
->ekinf
[i
],ekind
->tcstat
[i
].ekinf
);
998 copy_mat(ekinstate
->ekinh_old
[i
],ekind
->tcstat
[i
].ekinh_old
);
999 ekind
->tcstat
[i
].ekinscalef_nhc
= ekinstate
->ekinscalef_nhc
[i
];
1000 ekind
->tcstat
[i
].ekinscaleh_nhc
= ekinstate
->ekinscaleh_nhc
[i
];
1001 ekind
->tcstat
[i
].vscale_nhc
= ekinstate
->vscale_nhc
[i
];
1004 copy_mat(ekinstate
->ekin_total
,ekind
->ekin
);
1006 ekind
->dekindl
= ekinstate
->dekindl
;
1007 ekind
->cosacc
.mvcos
= ekinstate
->mvcos
;
1008 n
= ekinstate
->ekin_n
;
1013 gmx_bcast(sizeof(n
),&n
,cr
);
1016 gmx_bcast(DIM
*DIM
*sizeof(ekind
->tcstat
[i
].ekinh
[0][0]),
1017 ekind
->tcstat
[i
].ekinh
[0],cr
);
1018 gmx_bcast(DIM
*DIM
*sizeof(ekind
->tcstat
[i
].ekinf
[0][0]),
1019 ekind
->tcstat
[i
].ekinf
[0],cr
);
1020 gmx_bcast(DIM
*DIM
*sizeof(ekind
->tcstat
[i
].ekinh_old
[0][0]),
1021 ekind
->tcstat
[i
].ekinh_old
[0],cr
);
1023 gmx_bcast(sizeof(ekind
->tcstat
[i
].ekinscalef_nhc
),
1024 &(ekind
->tcstat
[i
].ekinscalef_nhc
),cr
);
1025 gmx_bcast(sizeof(ekind
->tcstat
[i
].ekinscaleh_nhc
),
1026 &(ekind
->tcstat
[i
].ekinscaleh_nhc
),cr
);
1027 gmx_bcast(sizeof(ekind
->tcstat
[i
].vscale_nhc
),
1028 &(ekind
->tcstat
[i
].vscale_nhc
),cr
);
1030 gmx_bcast(DIM
*DIM
*sizeof(ekind
->ekin
[0][0]),
1033 gmx_bcast(sizeof(ekind
->dekindl
),&ekind
->dekindl
,cr
);
1034 gmx_bcast(sizeof(ekind
->cosacc
.mvcos
),&ekind
->cosacc
.mvcos
,cr
);
1038 void set_deform_reference_box(gmx_update_t upd
,gmx_large_int_t step
,matrix box
)
1040 upd
->deformref_step
= step
;
1041 copy_mat(box
,upd
->deformref_box
);
1044 static void deform(gmx_update_t upd
,
1045 int start
,int homenr
,rvec x
[],matrix box
,matrix
*scale_tot
,
1046 const t_inputrec
*ir
,gmx_large_int_t step
)
1048 matrix bnew
,invbox
,mu
;
1052 elapsed_time
= (step
+ 1 - upd
->deformref_step
)*ir
->delta_t
;
1054 for(i
=0; i
<DIM
; i
++)
1056 for(j
=0; j
<DIM
; j
++)
1058 if (ir
->deform
[i
][j
] != 0)
1061 upd
->deformref_box
[i
][j
] + elapsed_time
*ir
->deform
[i
][j
];
1065 /* We correct the off-diagonal elements,
1066 * which can grow indefinitely during shearing,
1067 * so the shifts do not get messed up.
1069 for(i
=1; i
<DIM
; i
++)
1071 for(j
=i
-1; j
>=0; j
--)
1073 while (bnew
[i
][j
] - box
[i
][j
] > 0.5*bnew
[j
][j
])
1075 rvec_dec(bnew
[i
],bnew
[j
]);
1077 while (bnew
[i
][j
] - box
[i
][j
] < -0.5*bnew
[j
][j
])
1079 rvec_inc(bnew
[i
],bnew
[j
]);
1083 m_inv_ur0(box
,invbox
);
1085 mmul_ur0(box
,invbox
,mu
);
1087 for(i
=start
; i
<start
+homenr
; i
++)
1089 x
[i
][XX
] = mu
[XX
][XX
]*x
[i
][XX
]+mu
[YY
][XX
]*x
[i
][YY
]+mu
[ZZ
][XX
]*x
[i
][ZZ
];
1090 x
[i
][YY
] = mu
[YY
][YY
]*x
[i
][YY
]+mu
[ZZ
][YY
]*x
[i
][ZZ
];
1091 x
[i
][ZZ
] = mu
[ZZ
][ZZ
]*x
[i
][ZZ
];
1095 /* The transposes of the scaling matrices are stored,
1096 * so we need to do matrix multiplication in the inverse order.
1098 mmul_ur0(*scale_tot
,mu
,*scale_tot
);
1102 static void combine_forces(int nstlist
,
1103 gmx_constr_t constr
,
1104 t_inputrec
*ir
,t_mdatoms
*md
,t_idef
*idef
,
1105 t_commrec
*cr
,gmx_large_int_t step
,t_state
*state
,
1106 int start
,int nrend
,
1107 rvec f
[],rvec f_lr
[],
1112 /* f contains the short-range forces + the long range forces
1113 * which are stored separately in f_lr.
1116 if (constr
!= NULL
&& !(ir
->eConstrAlg
== econtSHAKE
&& ir
->epc
== epcNO
))
1118 /* We need to constrain the LR forces separately,
1119 * because due to the different pre-factor for the SR and LR
1120 * forces in the update algorithm, we can not determine
1121 * the constraint force for the coordinate constraining.
1122 * Constrain only the additional LR part of the force.
1124 /* MRS -- need to make sure this works with trotter integration -- the constraint calls may not be right.*/
1125 constrain(NULL
,FALSE
,FALSE
,constr
,idef
,ir
,NULL
,cr
,step
,0,md
,
1126 state
->x
,f_lr
,f_lr
,state
->box
,state
->lambda
,NULL
,
1127 NULL
,NULL
,nrnb
,econqForce
,ir
->epc
==epcMTTK
,state
->veta
,state
->veta
);
1130 /* Add nstlist-1 times the LR force to the sum of both forces
1131 * and store the result in forces_lr.
1134 for(i
=start
; i
<nrend
; i
++)
1136 for(d
=0; d
<DIM
; d
++)
1138 f_lr
[i
][d
] = f
[i
][d
] + nm1
*f_lr
[i
][d
];
1143 void update_extended(FILE *fplog
,
1144 gmx_large_int_t step
,
1145 t_inputrec
*inputrec
,
1148 gmx_ekindata_t
*ekind
,
1149 matrix pcoupl_mu
, /* pressure coupling */
1151 gmx_wallcycle_t wcycle
,
1157 bool bExtended
,bNH
,bPR
,bTrotter
,bLastStep
,bLog
=FALSE
,bEner
=FALSE
,bCouple
;
1160 int start
,homenr
,nrend
,i
,n
,m
,g
;
1163 homenr
= md
->homenr
;
1164 nrend
= start
+homenr
;
1166 /* if using vv, we do this elsewhere in the code */
1167 if ((IR_NVT_TROTTER(inputrec
)) || (IR_NPT_TROTTER(inputrec
)))
1172 /* We should only couple after a step where energies were determined */
1174 bCouple
= (inputrec
->nstcalcenergy
== 1 ||
1175 do_per_step(step
+inputrec
->nstcalcenergy
-1,
1176 inputrec
->nstcalcenergy
));
1177 dtc
= inputrec
->nstcalcenergy
*inputrec
->delta_t
;
1179 bNH
= (inputrec
->etc
== etcNOSEHOOVER
);
1180 bPR
= (inputrec
->epc
== epcPARRINELLORAHMAN
);
1181 bExtended
= (bNH
|| bPR
);
1183 dt
= inputrec
->delta_t
;
1185 clear_mat(pcoupl_mu
);
1186 for(i
=0; i
<DIM
; i
++)
1188 pcoupl_mu
[i
][i
] = 1.0;
1195 switch (inputrec
->etc
)
1200 berendsen_tcoupl(&(inputrec
->opts
),ekind
,dtc
);
1203 nosehoover_tcoupl(&(inputrec
->opts
),ekind
,dtc
,
1204 state
->nosehoover_xi
,state
->nosehoover_vxi
,MassQ
);
1207 vrescale_tcoupl(&(inputrec
->opts
),ekind
,dtc
,
1208 state
->therm_integral
,upd
->sd
->gaussrand
);
1212 switch (inputrec
->epc
)
1214 /* We can always pcoupl, even if we did not sum the energies
1215 * the previous step, since state->pres_prev is only updated
1216 * when the energies have been summed.
1220 case (epcBERENDSEN
):
1223 berendsen_pcoupl(fplog
,step
,inputrec
,dtc
,state
->pres_prev
,state
->box
,
1227 case (epcPARRINELLORAHMAN
):
1228 parrinellorahman_pcoupl(fplog
,step
,inputrec
,dtc
,state
->pres_prev
,
1229 state
->box
,state
->box_rel
,state
->boxv
,
1230 M
,pcoupl_mu
,bInitStep
);
1238 /* Set the T scaling lambda to 1 to have no scaling */
1239 for(i
=0; (i
<inputrec
->opts
.ngtc
); i
++)
1241 ekind
->tcstat
[i
].lambda
= 1.0;
1246 static rvec
*get_xprime(const t_state
*state
,gmx_update_t upd
)
1248 if (state
->nalloc
> upd
->xp_nalloc
)
1250 upd
->xp_nalloc
= state
->nalloc
;
1251 srenew(upd
->xp
,upd
->xp_nalloc
);
1257 void update_constraints(FILE *fplog
,
1258 gmx_large_int_t step
,
1259 real
*dvdlambda
, /* FEP stuff */
1260 t_inputrec
*inputrec
, /* input record and box stuff */
1261 gmx_ekindata_t
*ekind
,
1265 rvec force
[], /* forces on home particles */
1268 tensor vir
, /* tensors for virial and ekin, needed for computing */
1271 gmx_wallcycle_t wcycle
,
1273 gmx_constr_t constr
,
1279 bool bExtended
,bTrotter
,bLastStep
,bLog
=FALSE
,bEner
=FALSE
,bDoConstr
=FALSE
;
1282 int start
,homenr
,nrend
,i
,n
,m
,g
,d
;
1284 rvec
*vbuf
,*xprime
=NULL
;
1286 if (constr
) {bDoConstr
=TRUE
;}
1287 if (bFirstHalf
&& !EI_VV(inputrec
->eI
)) {bDoConstr
=FALSE
;}
1289 /* for now, SD update is here -- though it really seems like it
1290 should be reformulated as a velocity verlet method, since it has two parts */
1293 homenr
= md
->homenr
;
1294 nrend
= start
+homenr
;
1296 dt
= inputrec
->delta_t
;
1301 * APPLY CONSTRAINTS:
1304 * When doing PR pressure coupling we have to constrain the
1305 * bonds in each iteration. If we are only using Nose-Hoover tcoupling
1306 * it is enough to do this once though, since the relative velocities
1307 * after this will be normal to the bond vector
1312 /* clear out constraints before applying */
1313 clear_mat(vir_part
);
1315 xprime
= get_xprime(state
,upd
);
1317 bLastStep
= (step
== inputrec
->init_step
+inputrec
->nsteps
);
1318 bLog
= (do_per_step(step
,inputrec
->nstlog
) || bLastStep
|| (step
< 0));
1319 bEner
= (do_per_step(step
,inputrec
->nstenergy
) || bLastStep
);
1320 /* Constrain the coordinates xprime */
1321 wallcycle_start(wcycle
,ewcCONSTR
);
1322 if (EI_VV(inputrec
->eI
) && bFirstHalf
)
1324 constrain(NULL
,bLog
,bEner
,constr
,idef
,
1325 inputrec
,ekind
,cr
,step
,1,md
,
1326 state
->x
,state
->v
,state
->v
,
1327 state
->box
,state
->lambda
,dvdlambda
,
1328 NULL
,bCalcVir
? &vir_con
: NULL
,nrnb
,econqVeloc
,
1329 inputrec
->epc
==epcMTTK
,state
->veta
,vetanew
);
1333 constrain(NULL
,bLog
,bEner
,constr
,idef
,
1334 inputrec
,ekind
,cr
,step
,1,md
,
1335 state
->x
,xprime
,NULL
,
1336 state
->box
,state
->lambda
,dvdlambda
,
1337 state
->v
,bCalcVir
? &vir_con
: NULL
,nrnb
,econqCoord
,
1338 inputrec
->epc
==epcMTTK
,state
->veta
,state
->veta
);
1340 wallcycle_stop(wcycle
,ewcCONSTR
);
1344 dump_it_all(fplog
,"After Shake",
1345 state
->natoms
,state
->x
,xprime
,state
->v
,force
);
1349 if (inputrec
->eI
== eiSD2
)
1351 /* A correction factor eph is needed for the SD constraint force */
1352 /* Here we can, unfortunately, not have proper corrections
1353 * for different friction constants, so we use the first one.
1355 for(i
=0; i
<DIM
; i
++)
1357 for(m
=0; m
<DIM
; m
++)
1359 vir_part
[i
][m
] += upd
->sd
->sdc
[0].eph
*vir_con
[i
][m
];
1365 m_add(vir_part
,vir_con
,vir_part
);
1369 pr_rvecs(debug
,0,"constraint virial",vir_part
,DIM
);
1375 if ((inputrec
->eI
== eiSD2
) && !(bFirstHalf
))
1377 xprime
= get_xprime(state
,upd
);
1379 /* The second part of the SD integration */
1380 do_update_sd2(upd
->sd
,FALSE
,start
,homenr
,
1381 inputrec
->opts
.acc
,inputrec
->opts
.nFreeze
,
1382 md
->invmass
,md
->ptype
,
1383 md
->cFREEZE
,md
->cACC
,md
->cTC
,
1384 state
->x
,xprime
,state
->v
,force
,state
->sd_X
,
1385 inputrec
->opts
.ngtc
,inputrec
->opts
.tau_t
,
1386 inputrec
->opts
.ref_t
,FALSE
);
1387 inc_nrnb(nrnb
, eNR_UPDATE
, homenr
);
1391 /* Constrain the coordinates xprime */
1392 wallcycle_start(wcycle
,ewcCONSTR
);
1393 constrain(NULL
,bLog
,bEner
,constr
,idef
,
1394 inputrec
,NULL
,cr
,step
,1,md
,
1395 state
->x
,xprime
,NULL
,
1396 state
->box
,state
->lambda
,dvdlambda
,
1397 NULL
,NULL
,nrnb
,econqCoord
,FALSE
,0,0);
1398 wallcycle_stop(wcycle
,ewcCONSTR
);
1402 /* We must always unshift after updating coordinates; if we did not shake
1403 x was shifted in do_force */
1405 if (!(bFirstHalf
)) /* in the first half of vv, no shift. */
1407 if (graph
&& (graph
->nnodes
> 0))
1409 unshift_x(graph
,state
->box
,state
->x
,upd
->xp
);
1410 if (TRICLINIC(state
->box
))
1412 inc_nrnb(nrnb
,eNR_SHIFTX
,2*graph
->nnodes
);
1416 inc_nrnb(nrnb
,eNR_SHIFTX
,graph
->nnodes
);
1418 copy_rvecn(upd
->xp
,state
->x
,start
,graph
->start
);
1419 copy_rvecn(upd
->xp
,state
->x
,graph
->start
+graph
->nnodes
,nrend
);
1423 copy_rvecn(upd
->xp
,state
->x
,start
,nrend
);
1426 dump_it_all(fplog
,"After unshift",
1427 state
->natoms
,state
->x
,upd
->xp
,state
->v
,force
);
1429 /* ############# END the update of velocities and positions ######### */
1432 void update_box(FILE *fplog
,
1433 gmx_large_int_t step
,
1434 t_inputrec
*inputrec
, /* input record and box stuff */
1438 rvec force
[], /* forces on home particles */
1442 gmx_wallcycle_t wcycle
,
1447 bool bExtended
,bTrotter
,bLastStep
,bLog
=FALSE
,bEner
=FALSE
;
1450 int start
,homenr
,nrend
,i
,n
,m
,g
;
1454 homenr
= md
->homenr
;
1455 nrend
= start
+homenr
;
1458 (inputrec
->etc
== etcNOSEHOOVER
) ||
1459 (inputrec
->epc
== epcPARRINELLORAHMAN
) ||
1460 (inputrec
->epc
== epcMTTK
);
1462 dt
= inputrec
->delta_t
;
1466 /* now update boxes */
1467 switch (inputrec
->epc
) {
1470 case (epcBERENDSEN
):
1471 berendsen_pscale(inputrec
,pcoupl_mu
,state
->box
,state
->box_rel
,
1472 start
,homenr
,state
->x
,md
->cFREEZE
,nrnb
);
1474 case (epcPARRINELLORAHMAN
):
1475 /* The box velocities were updated in do_pr_pcoupl in the update
1476 * iteration, but we dont change the box vectors until we get here
1477 * since we need to be able to shift/unshift above.
1483 state
->box
[i
][m
] += dt
*state
->boxv
[i
][m
];
1486 preserve_box_shape(inputrec
,state
->box_rel
,state
->box
);
1488 /* Scale the coordinates */
1489 for(n
=start
; (n
<start
+homenr
); n
++)
1491 tmvmul_ur0(pcoupl_mu
,state
->x
[n
],state
->x
[n
]);
1495 switch (inputrec
->epct
)
1497 case (epctISOTROPIC
):
1498 /* DIM * eta = ln V. so DIM*eta_new = DIM*eta_old + DIM*dt*veta =>
1499 ln V_new = ln V_old + 3*dt*veta => V_new = V_old*exp(3*dt*veta) =>
1500 Side length scales as exp(veta*dt) */
1502 msmul(state
->box
,exp(state
->veta
*dt
),state
->box
);
1504 /* Relate veta to boxv. veta = d(eta)/dT = (1/DIM)*1/V dV/dT.
1505 o If we assume isotropic scaling, and box length scaling
1506 factor L, then V = L^DIM (det(M)). So dV/dt = DIM
1507 L^(DIM-1) dL/dt det(M), and veta = (1/L) dL/dt. The
1508 determinant of B is L^DIM det(M), and the determinant
1509 of dB/dt is (dL/dT)^DIM det (M). veta will be
1510 (det(dB/dT)/det(B))^(1/3). Then since M =
1511 B_new*(vol_new)^(1/3), dB/dT_new = (veta_new)*B(new). */
1513 msmul(state
->box
,state
->veta
,state
->boxv
);
1523 if ((!(IR_NPT_TROTTER(inputrec
))) && scale_tot
)
1525 /* The transposes of the scaling matrices are stored,
1526 * therefore we need to reverse the order in the multiplication.
1528 mmul_ur0(*scale_tot
,pcoupl_mu
,*scale_tot
);
1531 if (DEFORM(*inputrec
))
1533 deform(upd
,start
,homenr
,state
->x
,state
->box
,scale_tot
,inputrec
,step
);
1536 dump_it_all(fplog
,"After update",
1537 state
->natoms
,state
->x
,upd
->xp
,state
->v
,force
);
1540 void update_coords(FILE *fplog
,
1541 gmx_large_int_t step
,
1542 t_inputrec
*inputrec
, /* input record and box stuff */
1545 rvec
*f
, /* forces on home particles */
1549 gmx_ekindata_t
*ekind
,
1551 gmx_wallcycle_t wcycle
,
1555 t_commrec
*cr
, /* these shouldn't be here -- need to think about it */
1557 gmx_constr_t constr
,
1560 bool bExtended
,bNH
,bPR
,bTrotter
,bLastStep
,bLog
=FALSE
,bEner
=FALSE
;
1562 real
*imass
,*imassin
;
1565 int start
,homenr
,nrend
,i
,j
,d
,n
,m
,g
;
1566 int blen0
,blen1
,iatom
,jatom
,nshake
,nsettle
,nconstr
,nexpand
;
1569 rvec
*vcom
,*xcom
,*vall
,*xall
,*xin
,*vin
,*forcein
,*fall
,*xpall
,*xprimein
,*xprime
;
1572 /* Running the velocity half does nothing except for velocity verlet */
1573 if (UpdatePart
== etrtVELOCITY
)
1575 if (!EI_VV(inputrec
->eI
)) {return;} /*THIS IS FOR STARTING AT v(t=0)*/
1579 homenr
= md
->homenr
;
1580 nrend
= start
+homenr
;
1582 xprime
= get_xprime(state
,upd
);
1584 dt
= inputrec
->delta_t
;
1587 /* We need to update the NMR restraint history when time averaging is used */
1588 if (state
->flags
& (1<<estDISRE_RM3TAV
))
1590 update_disres_history(fcd
,&state
->hist
);
1592 if (state
->flags
& (1<<estORIRE_DTAV
))
1594 update_orires_history(fcd
,&state
->hist
);
1597 bNH
= inputrec
->etc
== etcNOSEHOOVER
;
1598 bPR
= ((inputrec
->epc
== epcPARRINELLORAHMAN
) || (inputrec
->epc
== epcMTTK
));
1600 bExtended
= bNH
|| bPR
;
1602 if (bDoLR
&& inputrec
->nstlist
> 1 && !EI_VV(inputrec
->eI
)) /* get this working with VV? */
1604 /* Store the total force + nstlist-1 times the LR force
1605 * in forces_lr, so it can be used in a normal update algorithm
1606 * to produce twin time stepping.
1608 /* is this correct in the new construction? MRS */
1609 combine_forces(inputrec
->nstlist
,constr
,inputrec
,md
,idef
,cr
,step
,state
,
1610 start
,nrend
,f
,f_lr
,nrnb
);
1618 /* ############# START The update of velocities and positions ######### */
1620 dump_it_all(fplog
,"Before update",
1621 state
->natoms
,state
->x
,xprime
,state
->v
,force
);
1623 switch (inputrec
->eI
) {
1625 if (ekind
->cosacc
.cos_accel
== 0) {
1626 /* use normal version of update */
1627 do_update_md(start
,nrend
,dt
,
1628 ekind
->tcstat
,ekind
->grpstat
,state
->nosehoover_vxi
,
1629 inputrec
->opts
.acc
,inputrec
->opts
.nFreeze
,md
->invmass
,md
->ptype
,
1630 md
->cFREEZE
,md
->cACC
,md
->cTC
,
1631 state
->x
,xprime
,state
->v
,force
,M
,
1636 do_update_visc(start
,nrend
,dt
,
1637 ekind
->tcstat
,md
->invmass
,state
->nosehoover_vxi
,
1638 md
->ptype
,md
->cTC
,state
->x
,xprime
,state
->v
,force
,M
,
1639 state
->box
,ekind
->cosacc
.cos_accel
,ekind
->cosacc
.vcos
,bNH
,bPR
);
1643 do_update_sd1(upd
->sd
,start
,homenr
,dt
,
1644 inputrec
->opts
.acc
,inputrec
->opts
.nFreeze
,
1645 md
->invmass
,md
->ptype
,
1646 md
->cFREEZE
,md
->cACC
,md
->cTC
,
1647 state
->x
,xprime
,state
->v
,force
,state
->sd_X
,
1648 inputrec
->opts
.ngtc
,inputrec
->opts
.tau_t
,inputrec
->opts
.ref_t
);
1651 /* The SD update is done in 2 parts, because an extra constraint step
1654 do_update_sd2(upd
->sd
,bInitStep
,start
,homenr
,
1655 inputrec
->opts
.acc
,inputrec
->opts
.nFreeze
,
1656 md
->invmass
,md
->ptype
,
1657 md
->cFREEZE
,md
->cACC
,md
->cTC
,
1658 state
->x
,xprime
,state
->v
,force
,state
->sd_X
,
1659 inputrec
->opts
.ngtc
,inputrec
->opts
.tau_t
,inputrec
->opts
.ref_t
,
1663 do_update_bd(start
,nrend
,dt
,
1664 inputrec
->opts
.nFreeze
,md
->invmass
,md
->ptype
,
1665 md
->cFREEZE
,md
->cTC
,
1666 state
->x
,xprime
,state
->v
,force
,
1668 inputrec
->opts
.ngtc
,inputrec
->opts
.tau_t
,inputrec
->opts
.ref_t
,
1669 upd
->sd
->bd_rf
,upd
->sd
->gaussrand
);
1673 alpha
= 1.0 + DIM
/((double)inputrec
->opts
.nrdf
[0]); /* assuming barostat coupled to group 0. */
1674 switch (UpdatePart
) {
1675 case (etrtVELOCITY
):
1676 do_update_vv_vel(start
,nrend
,dt
,
1677 ekind
->tcstat
,ekind
->grpstat
,
1678 inputrec
->opts
.acc
,inputrec
->opts
.nFreeze
,
1679 md
->invmass
,md
->ptype
,
1680 md
->cFREEZE
,md
->cACC
,md
->cTC
,
1681 state
->v
,force
,bExtended
,state
->veta
,alpha
);
1683 case (etrtPOSITION
):
1684 do_update_vv_pos(start
,nrend
,dt
,
1685 ekind
->tcstat
,ekind
->grpstat
,
1686 inputrec
->opts
.acc
,inputrec
->opts
.nFreeze
,
1687 md
->invmass
,md
->ptype
,
1688 md
->cFREEZE
,md
->cACC
,md
->cTC
,
1689 state
->x
,xprime
,state
->v
,force
,
1690 bExtended
,state
->veta
,alpha
);
1695 gmx_fatal(FARGS
,"Don't know how to update coordinates");
1701 void correct_ekin(FILE *log
,int start
,int end
,rvec v
[],rvec vcm
,real mass
[],
1702 real tmass
,tensor ekin
)
1705 * This is a debugging routine. It should not be called for production code
1707 * The kinetic energy should calculated according to:
1708 * Ekin = 1/2 m (v-vcm)^2
1709 * However the correction is not always applied, since vcm may not be
1710 * known in time and we compute
1711 * Ekin' = 1/2 m v^2 instead
1712 * This can be corrected afterwards by computing
1713 * Ekin = Ekin' + 1/2 m ( -2 v vcm + vcm^2)
1715 * Ekin = Ekin' - m v vcm + 1/2 m vcm^2
1722 /* Local particles */
1725 /* Processor dependent part. */
1727 for(i
=start
; (i
<end
); i
++)
1731 for(j
=0; (j
<DIM
); j
++)
1737 svmul(1/tmass
,vcm
,vcm
);
1738 svmul(0.5,vcm
,hvcm
);
1740 for(j
=0; (j
<DIM
); j
++)
1742 for(k
=0; (k
<DIM
); k
++)
1744 dekin
[j
][k
] += vcm
[k
]*(tm
*hvcm
[j
]-mv
[j
]);
1747 pr_rvecs(log
,0,"dekin",dekin
,DIM
);
1748 pr_rvecs(log
,0," ekin", ekin
,DIM
);
1749 fprintf(log
,"dekin = %g, ekin = %g vcm = (%8.4f %8.4f %8.4f)\n",
1750 trace(dekin
),trace(ekin
),vcm
[XX
],vcm
[YY
],vcm
[ZZ
]);
1751 fprintf(log
,"mv = (%8.4f %8.4f %8.4f)\n",
1752 mv
[XX
],mv
[YY
],mv
[ZZ
]);