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
[],
127 gmx_bool bNH
,gmx_bool bPR
)
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
[],unsigned short cACC
[],
233 gmx_bool bExtended
, real veta
, real alpha
)
238 real u
,vn
,vv
,va
,vb
,vnrel
;
244 g
= 0.25*dt
*veta
*alpha
;
246 mv2
= series_sinhx(g
);
253 for(n
=start
; n
<nrend
; n
++)
255 w_dt
= invmass
[n
]*dt
;
267 if((ptype
[n
] != eptVSite
) && (ptype
[n
] != eptShell
) && !nFreeze
[gf
][d
])
269 v
[n
][d
] = mv1
*(mv1
*v
[n
][d
] + 0.5*(w_dt
*mv2
*f
[n
][d
]))+0.5*accel
[ga
][d
]*dt
;
277 } /* do_update_vv_vel */
279 static void do_update_vv_pos(int start
,int nrend
,double dt
,
280 t_grp_tcstat
*tcstat
,t_grp_acc
*gstat
,
281 rvec accel
[],ivec nFreeze
[],real invmass
[],
282 unsigned short ptype
[],
283 unsigned short cFREEZE
[],
284 rvec x
[],rvec xprime
[],rvec v
[],
285 rvec f
[],gmx_bool bExtended
, real veta
, real alpha
)
295 mr2
= series_sinhx(g
);
303 for(n
=start
; n
<nrend
; n
++) {
311 if ((ptype
[n
] != eptVSite
) && (ptype
[n
] != eptShell
) && !nFreeze
[gf
][d
])
313 xprime
[n
][d
] = mr1
*(mr1
*x
[n
][d
]+mr2
*dt
*v
[n
][d
]);
317 xprime
[n
][d
] = x
[n
][d
];
321 }/* do_update_vv_pos */
323 static void do_update_visc(int start
,int nrend
,double dt
,
324 t_grp_tcstat
*tcstat
,real invmass
[],double nh_vxi
[],
325 unsigned short ptype
[],unsigned short cTC
[],
326 rvec x
[],rvec xprime
[],rvec v
[],
327 rvec f
[],matrix M
,matrix box
,real
329 gmx_bool bNH
,gmx_bool bPR
)
339 fac
= 2*M_PI
/(box
[ZZ
][ZZ
]);
342 /* Update with coupling to extended ensembles, used for
343 * Nose-Hoover and Parrinello-Rahman coupling
345 for(n
=start
; n
<nrend
; n
++) {
351 lg
= tcstat
[gt
].lambda
;
352 cosz
= cos(fac
*x
[n
][ZZ
]);
354 copy_rvec(v
[n
],vrel
);
366 if((ptype
[n
] != eptVSite
) && (ptype
[n
] != eptShell
))
368 vn
= (lg
*vrel
[d
] + dt
*(imass
*f
[n
][d
] - 0.5*vxi
*vrel
[d
]
369 - iprod(M
[d
],vrel
)))/(1 + 0.5*vxi
*dt
);
372 vn
+= vc
+ dt
*cosz
*cos_accel
;
375 xprime
[n
][d
] = x
[n
][d
]+vn
*dt
;
379 xprime
[n
][d
] = x
[n
][d
];
386 /* Classic version of update, used with berendsen coupling */
387 for(n
=start
; n
<nrend
; n
++)
389 w_dt
= invmass
[n
]*dt
;
394 lg
= tcstat
[gt
].lambda
;
395 cosz
= cos(fac
*x
[n
][ZZ
]);
401 if((ptype
[n
] != eptVSite
) && (ptype
[n
] != eptShell
))
406 /* Do not scale the cosine velocity profile */
407 vv
= vc
+ lg
*(vn
- vc
+ f
[n
][d
]*w_dt
);
408 /* Add the cosine accelaration profile */
409 vv
+= dt
*cosz
*cos_accel
;
413 vv
= lg
*(vn
+ f
[n
][d
]*w_dt
);
416 xprime
[n
][d
] = x
[n
][d
]+vv
*dt
;
421 xprime
[n
][d
] = x
[n
][d
];
428 static gmx_stochd_t
*init_stochd(FILE *fplog
,t_inputrec
*ir
)
437 /* Initiate random number generator for langevin type dynamics,
438 * for BD, SD or velocity rescaling temperature coupling.
440 sd
->gaussrand
= gmx_rng_init(ir
->ld_seed
);
442 ngtc
= ir
->opts
.ngtc
;
446 snew(sd
->bd_rf
,ngtc
);
448 else if (EI_SD(ir
->eI
))
451 snew(sd
->sdsig
,ngtc
);
454 for(n
=0; n
<ngtc
; n
++)
456 if (ir
->opts
.tau_t
[n
] > 0)
458 sdc
[n
].gdt
= ir
->delta_t
/ir
->opts
.tau_t
[n
];
459 sdc
[n
].eph
= exp(sdc
[n
].gdt
/2);
460 sdc
[n
].emh
= exp(-sdc
[n
].gdt
/2);
461 sdc
[n
].em
= exp(-sdc
[n
].gdt
);
465 /* No friction and noise on this group */
471 if (sdc
[n
].gdt
>= 0.05)
473 sdc
[n
].b
= sdc
[n
].gdt
*(sdc
[n
].eph
*sdc
[n
].eph
- 1)
474 - 4*(sdc
[n
].eph
- 1)*(sdc
[n
].eph
- 1);
475 sdc
[n
].c
= sdc
[n
].gdt
- 3 + 4*sdc
[n
].emh
- sdc
[n
].em
;
476 sdc
[n
].d
= 2 - sdc
[n
].eph
- sdc
[n
].emh
;
481 /* Seventh order expansions for small y */
482 sdc
[n
].b
= y
*y
*y
*y
*(1/3.0+y
*(1/3.0+y
*(17/90.0+y
*7/9.0)));
483 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))));
484 sdc
[n
].d
= y
*y
*(-1+y
*y
*(-1/12.0-y
*y
/360.0));
487 fprintf(debug
,"SD const tc-grp %d: b %g c %g d %g\n",
488 n
,sdc
[n
].b
,sdc
[n
].c
,sdc
[n
].d
);
495 void get_stochd_state(gmx_update_t upd
,t_state
*state
)
497 gmx_rng_get_state(upd
->sd
->gaussrand
,state
->ld_rng
,state
->ld_rngi
);
500 void set_stochd_state(gmx_update_t upd
,t_state
*state
)
502 gmx_rng_set_state(upd
->sd
->gaussrand
,state
->ld_rng
,state
->ld_rngi
[0]);
505 gmx_update_t
init_update(FILE *fplog
,t_inputrec
*ir
)
511 if (ir
->eI
== eiBD
|| EI_SD(ir
->eI
) || ir
->etc
== etcVRESCALE
)
513 upd
->sd
= init_stochd(fplog
,ir
);
522 static void do_update_sd1(gmx_stochd_t
*sd
,
523 int start
,int homenr
,double dt
,
524 rvec accel
[],ivec nFreeze
[],
525 real invmass
[],unsigned short ptype
[],
526 unsigned short cFREEZE
[],unsigned short cACC
[],
527 unsigned short cTC
[],
528 rvec x
[],rvec xprime
[],rvec v
[],rvec f
[],
530 int ngtc
,real tau_t
[],real ref_t
[])
542 if (homenr
> sd
->sd_V_nalloc
)
544 sd
->sd_V_nalloc
= over_alloc_dd(homenr
);
545 srenew(sd
->sd_V
,sd
->sd_V_nalloc
);
547 gaussrand
= sd
->gaussrand
;
549 for(n
=0; n
<ngtc
; n
++)
552 /* The mass is encounted for later, since this differs per atom */
553 sig
[n
].V
= sqrt(2*kT
*(1 - sdc
[n
].em
));
556 for(n
=start
; n
<start
+homenr
; n
++)
558 ism
= sqrt(invmass
[n
]);
574 if((ptype
[n
] != eptVSite
) && (ptype
[n
] != eptShell
) && !nFreeze
[gf
][d
])
576 sd_V
= ism
*sig
[gt
].V
*gmx_rng_gaussian_table(gaussrand
);
578 v
[n
][d
] = v
[n
][d
]*sdc
[gt
].em
579 + (invmass
[n
]*f
[n
][d
] + accel
[ga
][d
])*tau_t
[gt
]*(1 - sdc
[gt
].em
)
582 xprime
[n
][d
] = x
[n
][d
] + v
[n
][d
]*dt
;
587 xprime
[n
][d
] = x
[n
][d
];
593 static void do_update_sd2(gmx_stochd_t
*sd
,gmx_bool bInitStep
,
594 int start
,int homenr
,
595 rvec accel
[],ivec nFreeze
[],
596 real invmass
[],unsigned short ptype
[],
597 unsigned short cFREEZE
[],unsigned short cACC
[],
598 unsigned short cTC
[],
599 rvec x
[],rvec xprime
[],rvec v
[],rvec f
[],
601 int ngtc
,real tau_t
[],real ref_t
[],
606 /* The random part of the velocity update, generated in the first
607 * half of the update, needs to be remembered for the second half.
619 if (homenr
> sd
->sd_V_nalloc
)
621 sd
->sd_V_nalloc
= over_alloc_dd(homenr
);
622 srenew(sd
->sd_V
,sd
->sd_V_nalloc
);
625 gaussrand
= sd
->gaussrand
;
629 for (n
=0; n
<ngtc
; n
++)
632 /* The mass is encounted for later, since this differs per atom */
633 sig
[n
].V
= sqrt(kT
*(1-sdc
[n
].em
));
634 sig
[n
].X
= sqrt(kT
*sqr(tau_t
[n
])*sdc
[n
].c
);
635 sig
[n
].Yv
= sqrt(kT
*sdc
[n
].b
/sdc
[n
].c
);
636 sig
[n
].Yx
= sqrt(kT
*sqr(tau_t
[n
])*sdc
[n
].b
/(1-sdc
[n
].em
));
640 for (n
=start
; n
<start
+homenr
; n
++)
642 ism
= sqrt(invmass
[n
]);
662 if((ptype
[n
] != eptVSite
) && (ptype
[n
] != eptShell
) && !nFreeze
[gf
][d
])
668 sd_X
[n
][d
] = ism
*sig
[gt
].X
*gmx_rng_gaussian_table(gaussrand
);
670 Vmh
= sd_X
[n
][d
]*sdc
[gt
].d
/(tau_t
[gt
]*sdc
[gt
].c
)
671 + ism
*sig
[gt
].Yv
*gmx_rng_gaussian_table(gaussrand
);
672 sd_V
[n
-start
][d
] = ism
*sig
[gt
].V
*gmx_rng_gaussian_table(gaussrand
);
674 v
[n
][d
] = vn
*sdc
[gt
].em
675 + (invmass
[n
]*f
[n
][d
] + accel
[ga
][d
])*tau_t
[gt
]*(1 - sdc
[gt
].em
)
676 + sd_V
[n
-start
][d
] - sdc
[gt
].em
*Vmh
;
678 xprime
[n
][d
] = x
[n
][d
] + v
[n
][d
]*tau_t
[gt
]*(sdc
[gt
].eph
- sdc
[gt
].emh
);
683 /* Correct the velocities for the constraints.
684 * This operation introduces some inaccuracy,
685 * since the velocity is determined from differences in coordinates.
688 (xprime
[n
][d
] - x
[n
][d
])/(tau_t
[gt
]*(sdc
[gt
].eph
- sdc
[gt
].emh
));
690 Xmh
= sd_V
[n
-start
][d
]*tau_t
[gt
]*sdc
[gt
].d
/(sdc
[gt
].em
-1)
691 + ism
*sig
[gt
].Yx
*gmx_rng_gaussian_table(gaussrand
);
692 sd_X
[n
][d
] = ism
*sig
[gt
].X
*gmx_rng_gaussian_table(gaussrand
);
694 xprime
[n
][d
] += sd_X
[n
][d
] - Xmh
;
703 xprime
[n
][d
] = x
[n
][d
];
710 static void do_update_bd(int start
,int nrend
,double dt
,
712 real invmass
[],unsigned short ptype
[],
713 unsigned short cFREEZE
[],unsigned short cTC
[],
714 rvec x
[],rvec xprime
[],rvec v
[],
715 rvec f
[],real friction_coefficient
,
716 int ngtc
,real tau_t
[],real ref_t
[],
717 real
*rf
,gmx_rng_t gaussrand
)
719 /* note -- these appear to be full step velocities . . . */
725 if (friction_coefficient
!= 0)
727 invfr
= 1.0/friction_coefficient
;
728 for(n
=0; n
<ngtc
; n
++)
730 rf
[n
] = sqrt(2.0*BOLTZ
*ref_t
[n
]/(friction_coefficient
*dt
));
735 for(n
=0; n
<ngtc
; n
++)
737 rf
[n
] = sqrt(2.0*BOLTZ
*ref_t
[n
]);
740 for(n
=start
; (n
<nrend
); n
++)
750 for(d
=0; (d
<DIM
); d
++)
752 if((ptype
[n
]!=eptVSite
) && (ptype
[n
]!=eptShell
) && !nFreeze
[gf
][d
])
754 if (friction_coefficient
!= 0) {
755 vn
= invfr
*f
[n
][d
] + rf
[gt
]*gmx_rng_gaussian_table(gaussrand
);
759 /* NOTE: invmass = 1/(mass*friction_constant*dt) */
760 vn
= invmass
[n
]*f
[n
][d
]*dt
761 + sqrt(invmass
[n
])*rf
[gt
]*gmx_rng_gaussian_table(gaussrand
);
764 xprime
[n
][d
] = x
[n
][d
]+vn
*dt
;
770 xprime
[n
][d
] = x
[n
][d
];
776 static void dump_it_all(FILE *fp
,const char *title
,
777 int natoms
,rvec x
[],rvec xp
[],rvec v
[],rvec f
[])
782 fprintf(fp
,"%s\n",title
);
783 pr_rvecs(fp
,0,"x",x
,natoms
);
784 pr_rvecs(fp
,0,"xp",xp
,natoms
);
785 pr_rvecs(fp
,0,"v",v
,natoms
);
786 pr_rvecs(fp
,0,"f",f
,natoms
);
791 static void calc_ke_part_normal(rvec v
[], t_grpopts
*opts
,t_mdatoms
*md
,
792 gmx_ekindata_t
*ekind
,t_nrnb
*nrnb
,gmx_bool bEkinAveVel
,
793 gmx_bool bSaveEkinOld
)
795 int start
=md
->start
,homenr
=md
->homenr
;
796 int g
,d
,n
,m
,ga
=0,gt
=0;
799 t_grp_tcstat
*tcstat
=ekind
->tcstat
;
800 t_grp_acc
*grpstat
=ekind
->grpstat
;
803 /* three main: VV with AveVel, vv with AveEkin, leap with AveEkin. Leap with AveVel is also
804 an option, but not supported now. Additionally, if we are doing iterations.
805 bEkinAveVel: If TRUE, we sum into ekin, if FALSE, into ekinh.
806 bSavEkinOld: If TRUE (in the case of iteration = bIterate is TRUE), we don't copy over the ekinh_old.
807 If FALSE, we overrwrite it.
810 /* group velocities are calculated in update_ekindata and
811 * accumulated in acumulate_groups.
812 * Now the partial global and groups ekin.
814 for(g
=0; (g
<opts
->ngtc
); g
++)
818 copy_mat(tcstat
[g
].ekinh
,tcstat
[g
].ekinh_old
);
821 clear_mat(tcstat
[g
].ekinf
);
823 clear_mat(tcstat
[g
].ekinh
);
826 tcstat
[g
].ekinscalef_nhc
= 1.0; /* need to clear this -- logic is complicated! */
829 ekind
->dekindl_old
= ekind
->dekindl
;
832 for(n
=start
; (n
<start
+homenr
); n
++)
842 hm
= 0.5*md
->massT
[n
];
844 for(d
=0; (d
<DIM
); d
++)
846 v_corrt
[d
] = v
[n
][d
] - grpstat
[ga
].u
[d
];
848 for(d
=0; (d
<DIM
); d
++)
850 for (m
=0;(m
<DIM
); m
++)
852 /* if we're computing a full step velocity, v_corrt[d] has v(t). Otherwise, v(t+dt/2) */
855 tcstat
[gt
].ekinf
[m
][d
]+=hm
*v_corrt
[m
]*v_corrt
[d
];
859 tcstat
[gt
].ekinh
[m
][d
]+=hm
*v_corrt
[m
]*v_corrt
[d
];
863 if (md
->nMassPerturbed
&& md
->bPerturbed
[n
])
865 dekindl
-= 0.5*(md
->massB
[n
] - md
->massA
[n
])*iprod(v_corrt
,v_corrt
);
868 ekind
->dekindl
= dekindl
;
869 inc_nrnb(nrnb
,eNR_EKIN
,homenr
);
872 static void calc_ke_part_visc(matrix box
,rvec x
[],rvec v
[],
873 t_grpopts
*opts
,t_mdatoms
*md
,
874 gmx_ekindata_t
*ekind
,
875 t_nrnb
*nrnb
, gmx_bool bEkinAveVel
, gmx_bool bSaveEkinOld
)
877 int start
=md
->start
,homenr
=md
->homenr
;
881 t_grp_tcstat
*tcstat
=ekind
->tcstat
;
882 t_cos_acc
*cosacc
=&(ekind
->cosacc
);
887 for(g
=0; g
<opts
->ngtc
; g
++)
889 copy_mat(ekind
->tcstat
[g
].ekinh
,ekind
->tcstat
[g
].ekinh_old
);
890 clear_mat(ekind
->tcstat
[g
].ekinh
);
892 ekind
->dekindl_old
= ekind
->dekindl
;
894 fac
= 2*M_PI
/box
[ZZ
][ZZ
];
897 for(n
=start
; n
<start
+homenr
; n
++)
903 hm
= 0.5*md
->massT
[n
];
905 /* Note that the times of x and v differ by half a step */
906 /* MRS -- would have to be changed for VV */
907 cosz
= cos(fac
*x
[n
][ZZ
]);
908 /* Calculate the amplitude of the new velocity profile */
909 mvcos
+= 2*cosz
*md
->massT
[n
]*v
[n
][XX
];
911 copy_rvec(v
[n
],v_corrt
);
912 /* Subtract the profile for the kinetic energy */
913 v_corrt
[XX
] -= cosz
*cosacc
->vcos
;
914 for (d
=0; (d
<DIM
); d
++)
916 for (m
=0; (m
<DIM
); m
++)
918 /* if we're computing a full step velocity, v_corrt[d] has v(t). Otherwise, v(t+dt/2) */
921 tcstat
[gt
].ekinf
[m
][d
]+=hm
*v_corrt
[m
]*v_corrt
[d
];
925 tcstat
[gt
].ekinh
[m
][d
]+=hm
*v_corrt
[m
]*v_corrt
[d
];
929 if(md
->nPerturbed
&& md
->bPerturbed
[n
])
931 dekindl
-= 0.5*(md
->massB
[n
] - md
->massA
[n
])*iprod(v_corrt
,v_corrt
);
934 ekind
->dekindl
= dekindl
;
935 cosacc
->mvcos
= mvcos
;
937 inc_nrnb(nrnb
,eNR_EKIN
,homenr
);
940 void calc_ke_part(t_state
*state
,t_grpopts
*opts
,t_mdatoms
*md
,
941 gmx_ekindata_t
*ekind
,t_nrnb
*nrnb
, gmx_bool bEkinAveVel
, gmx_bool bSaveEkinOld
)
943 if (ekind
->cosacc
.cos_accel
== 0)
945 calc_ke_part_normal(state
->v
,opts
,md
,ekind
,nrnb
,bEkinAveVel
,bSaveEkinOld
);
949 calc_ke_part_visc(state
->box
,state
->x
,state
->v
,opts
,md
,ekind
,nrnb
,bEkinAveVel
,bSaveEkinOld
);
953 void init_ekinstate(ekinstate_t
*ekinstate
,const t_inputrec
*ir
)
955 ekinstate
->ekin_n
= ir
->opts
.ngtc
;
956 snew(ekinstate
->ekinh
,ekinstate
->ekin_n
);
957 snew(ekinstate
->ekinf
,ekinstate
->ekin_n
);
958 snew(ekinstate
->ekinh_old
,ekinstate
->ekin_n
);
959 snew(ekinstate
->ekinscalef_nhc
,ekinstate
->ekin_n
);
960 snew(ekinstate
->ekinscaleh_nhc
,ekinstate
->ekin_n
);
961 snew(ekinstate
->vscale_nhc
,ekinstate
->ekin_n
);
962 ekinstate
->dekindl
= 0;
963 ekinstate
->mvcos
= 0;
966 void update_ekinstate(ekinstate_t
*ekinstate
,gmx_ekindata_t
*ekind
)
970 for(i
=0;i
<ekinstate
->ekin_n
;i
++)
972 copy_mat(ekind
->tcstat
[i
].ekinh
,ekinstate
->ekinh
[i
]);
973 copy_mat(ekind
->tcstat
[i
].ekinf
,ekinstate
->ekinf
[i
]);
974 copy_mat(ekind
->tcstat
[i
].ekinh_old
,ekinstate
->ekinh_old
[i
]);
975 ekinstate
->ekinscalef_nhc
[i
] = ekind
->tcstat
[i
].ekinscalef_nhc
;
976 ekinstate
->ekinscaleh_nhc
[i
] = ekind
->tcstat
[i
].ekinscaleh_nhc
;
977 ekinstate
->vscale_nhc
[i
] = ekind
->tcstat
[i
].vscale_nhc
;
980 copy_mat(ekind
->ekin
,ekinstate
->ekin_total
);
981 ekinstate
->dekindl
= ekind
->dekindl
;
982 ekinstate
->mvcos
= ekind
->cosacc
.mvcos
;
986 void restore_ekinstate_from_state(t_commrec
*cr
,
987 gmx_ekindata_t
*ekind
,ekinstate_t
*ekinstate
)
993 for(i
=0;i
<ekinstate
->ekin_n
;i
++)
995 copy_mat(ekinstate
->ekinh
[i
],ekind
->tcstat
[i
].ekinh
);
996 copy_mat(ekinstate
->ekinf
[i
],ekind
->tcstat
[i
].ekinf
);
997 copy_mat(ekinstate
->ekinh_old
[i
],ekind
->tcstat
[i
].ekinh_old
);
998 ekind
->tcstat
[i
].ekinscalef_nhc
= ekinstate
->ekinscalef_nhc
[i
];
999 ekind
->tcstat
[i
].ekinscaleh_nhc
= ekinstate
->ekinscaleh_nhc
[i
];
1000 ekind
->tcstat
[i
].vscale_nhc
= ekinstate
->vscale_nhc
[i
];
1003 copy_mat(ekinstate
->ekin_total
,ekind
->ekin
);
1005 ekind
->dekindl
= ekinstate
->dekindl
;
1006 ekind
->cosacc
.mvcos
= ekinstate
->mvcos
;
1007 n
= ekinstate
->ekin_n
;
1012 gmx_bcast(sizeof(n
),&n
,cr
);
1015 gmx_bcast(DIM
*DIM
*sizeof(ekind
->tcstat
[i
].ekinh
[0][0]),
1016 ekind
->tcstat
[i
].ekinh
[0],cr
);
1017 gmx_bcast(DIM
*DIM
*sizeof(ekind
->tcstat
[i
].ekinf
[0][0]),
1018 ekind
->tcstat
[i
].ekinf
[0],cr
);
1019 gmx_bcast(DIM
*DIM
*sizeof(ekind
->tcstat
[i
].ekinh_old
[0][0]),
1020 ekind
->tcstat
[i
].ekinh_old
[0],cr
);
1022 gmx_bcast(sizeof(ekind
->tcstat
[i
].ekinscalef_nhc
),
1023 &(ekind
->tcstat
[i
].ekinscalef_nhc
),cr
);
1024 gmx_bcast(sizeof(ekind
->tcstat
[i
].ekinscaleh_nhc
),
1025 &(ekind
->tcstat
[i
].ekinscaleh_nhc
),cr
);
1026 gmx_bcast(sizeof(ekind
->tcstat
[i
].vscale_nhc
),
1027 &(ekind
->tcstat
[i
].vscale_nhc
),cr
);
1029 gmx_bcast(DIM
*DIM
*sizeof(ekind
->ekin
[0][0]),
1032 gmx_bcast(sizeof(ekind
->dekindl
),&ekind
->dekindl
,cr
);
1033 gmx_bcast(sizeof(ekind
->cosacc
.mvcos
),&ekind
->cosacc
.mvcos
,cr
);
1037 void set_deform_reference_box(gmx_update_t upd
,gmx_large_int_t step
,matrix box
)
1039 upd
->deformref_step
= step
;
1040 copy_mat(box
,upd
->deformref_box
);
1043 static void deform(gmx_update_t upd
,
1044 int start
,int homenr
,rvec x
[],matrix box
,matrix
*scale_tot
,
1045 const t_inputrec
*ir
,gmx_large_int_t step
)
1047 matrix bnew
,invbox
,mu
;
1051 elapsed_time
= (step
+ 1 - upd
->deformref_step
)*ir
->delta_t
;
1053 for(i
=0; i
<DIM
; i
++)
1055 for(j
=0; j
<DIM
; j
++)
1057 if (ir
->deform
[i
][j
] != 0)
1060 upd
->deformref_box
[i
][j
] + elapsed_time
*ir
->deform
[i
][j
];
1064 /* We correct the off-diagonal elements,
1065 * which can grow indefinitely during shearing,
1066 * so the shifts do not get messed up.
1068 for(i
=1; i
<DIM
; i
++)
1070 for(j
=i
-1; j
>=0; j
--)
1072 while (bnew
[i
][j
] - box
[i
][j
] > 0.5*bnew
[j
][j
])
1074 rvec_dec(bnew
[i
],bnew
[j
]);
1076 while (bnew
[i
][j
] - box
[i
][j
] < -0.5*bnew
[j
][j
])
1078 rvec_inc(bnew
[i
],bnew
[j
]);
1082 m_inv_ur0(box
,invbox
);
1084 mmul_ur0(box
,invbox
,mu
);
1086 for(i
=start
; i
<start
+homenr
; i
++)
1088 x
[i
][XX
] = mu
[XX
][XX
]*x
[i
][XX
]+mu
[YY
][XX
]*x
[i
][YY
]+mu
[ZZ
][XX
]*x
[i
][ZZ
];
1089 x
[i
][YY
] = mu
[YY
][YY
]*x
[i
][YY
]+mu
[ZZ
][YY
]*x
[i
][ZZ
];
1090 x
[i
][ZZ
] = mu
[ZZ
][ZZ
]*x
[i
][ZZ
];
1094 /* The transposes of the scaling matrices are stored,
1095 * so we need to do matrix multiplication in the inverse order.
1097 mmul_ur0(*scale_tot
,mu
,*scale_tot
);
1101 static void combine_forces(int nstlist
,
1102 gmx_constr_t constr
,
1103 t_inputrec
*ir
,t_mdatoms
*md
,t_idef
*idef
,
1104 t_commrec
*cr
,gmx_large_int_t step
,t_state
*state
,
1105 int start
,int nrend
,
1106 rvec f
[],rvec f_lr
[],
1111 /* f contains the short-range forces + the long range forces
1112 * which are stored separately in f_lr.
1115 if (constr
!= NULL
&& !(ir
->eConstrAlg
== econtSHAKE
&& ir
->epc
== epcNO
))
1117 /* We need to constrain the LR forces separately,
1118 * because due to the different pre-factor for the SR and LR
1119 * forces in the update algorithm, we can not determine
1120 * the constraint force for the coordinate constraining.
1121 * Constrain only the additional LR part of the force.
1123 /* MRS -- need to make sure this works with trotter integration -- the constraint calls may not be right.*/
1124 constrain(NULL
,FALSE
,FALSE
,constr
,idef
,ir
,NULL
,cr
,step
,0,md
,
1125 state
->x
,f_lr
,f_lr
,state
->box
,state
->lambda
,NULL
,
1126 NULL
,NULL
,nrnb
,econqForce
,ir
->epc
==epcMTTK
,state
->veta
,state
->veta
);
1129 /* Add nstlist-1 times the LR force to the sum of both forces
1130 * and store the result in forces_lr.
1133 for(i
=start
; i
<nrend
; i
++)
1135 for(d
=0; d
<DIM
; d
++)
1137 f_lr
[i
][d
] = f
[i
][d
] + nm1
*f_lr
[i
][d
];
1142 void update_tcouple(FILE *fplog
,
1143 gmx_large_int_t step
,
1144 t_inputrec
*inputrec
,
1146 gmx_ekindata_t
*ekind
,
1147 gmx_wallcycle_t wcycle
,
1153 gmx_bool bTCouple
=FALSE
;
1155 int i
,start
,end
,homenr
;
1157 /* if using vv, we do this elsewhere in the code */
1158 if (inputrec
->etc
!= etcNO
&&
1159 !(IR_NVT_TROTTER(inputrec
) || IR_NPT_TROTTER(inputrec
)))
1161 /* We should only couple after a step where energies were determined */
1162 bTCouple
= (inputrec
->nsttcouple
== 1 ||
1163 do_per_step(step
+inputrec
->nsttcouple
-1,
1164 inputrec
->nsttcouple
));
1169 dttc
= inputrec
->nsttcouple
*inputrec
->delta_t
;
1171 switch (inputrec
->etc
)
1176 berendsen_tcoupl(inputrec
,ekind
,dttc
);
1179 nosehoover_tcoupl(&(inputrec
->opts
),ekind
,dttc
,
1180 state
->nosehoover_xi
,state
->nosehoover_vxi
,MassQ
);
1183 vrescale_tcoupl(inputrec
,ekind
,dttc
,
1184 state
->therm_integral
,upd
->sd
->gaussrand
);
1187 /* rescale in place here */
1188 if (EI_VV(inputrec
->eI
))
1190 rescale_velocities(ekind
,md
,md
->start
,md
->start
+md
->homenr
,state
->v
);
1195 /* Set the T scaling lambda to 1 to have no scaling */
1196 for(i
=0; (i
<inputrec
->opts
.ngtc
); i
++)
1198 ekind
->tcstat
[i
].lambda
= 1.0;
1203 void update_pcouple(FILE *fplog
,
1204 gmx_large_int_t step
,
1205 t_inputrec
*inputrec
,
1209 gmx_wallcycle_t wcycle
,
1213 gmx_bool bPCouple
=FALSE
;
1217 /* if using vv, we do this elsewhere in the code */
1218 if (inputrec
->epc
!= epcNO
&&
1219 !(IR_NVT_TROTTER(inputrec
) || IR_NPT_TROTTER(inputrec
)))
1221 /* We should only couple after a step where energies were determined */
1222 bPCouple
= (inputrec
->nstpcouple
== 1 ||
1223 do_per_step(step
+inputrec
->nstpcouple
-1,
1224 inputrec
->nstpcouple
));
1227 clear_mat(pcoupl_mu
);
1228 for(i
=0; i
<DIM
; i
++)
1230 pcoupl_mu
[i
][i
] = 1.0;
1237 dtpc
= inputrec
->nstpcouple
*inputrec
->delta_t
;
1239 switch (inputrec
->epc
)
1241 /* We can always pcoupl, even if we did not sum the energies
1242 * the previous step, since state->pres_prev is only updated
1243 * when the energies have been summed.
1247 case (epcBERENDSEN
):
1250 berendsen_pcoupl(fplog
,step
,inputrec
,dtpc
,state
->pres_prev
,state
->box
,
1254 case (epcPARRINELLORAHMAN
):
1255 parrinellorahman_pcoupl(fplog
,step
,inputrec
,dtpc
,state
->pres_prev
,
1256 state
->box
,state
->box_rel
,state
->boxv
,
1257 M
,pcoupl_mu
,bInitStep
);
1265 static rvec
*get_xprime(const t_state
*state
,gmx_update_t upd
)
1267 if (state
->nalloc
> upd
->xp_nalloc
)
1269 upd
->xp_nalloc
= state
->nalloc
;
1270 srenew(upd
->xp
,upd
->xp_nalloc
);
1276 void update_constraints(FILE *fplog
,
1277 gmx_large_int_t step
,
1278 real
*dvdlambda
, /* FEP stuff */
1279 t_inputrec
*inputrec
, /* input record and box stuff */
1280 gmx_ekindata_t
*ekind
,
1284 rvec force
[], /* forces on home particles */
1287 tensor vir
, /* tensors for virial and ekin, needed for computing */
1290 gmx_wallcycle_t wcycle
,
1292 gmx_constr_t constr
,
1294 gmx_bool bFirstHalf
,
1298 gmx_bool bExtended
,bTrotter
,bLastStep
,bLog
=FALSE
,bEner
=FALSE
,bDoConstr
=FALSE
;
1301 int start
,homenr
,nrend
,i
,n
,m
,g
,d
;
1303 rvec
*vbuf
,*xprime
=NULL
;
1305 if (constr
) {bDoConstr
=TRUE
;}
1306 if (bFirstHalf
&& !EI_VV(inputrec
->eI
)) {bDoConstr
=FALSE
;}
1308 /* for now, SD update is here -- though it really seems like it
1309 should be reformulated as a velocity verlet method, since it has two parts */
1312 homenr
= md
->homenr
;
1313 nrend
= start
+homenr
;
1315 dt
= inputrec
->delta_t
;
1320 * APPLY CONSTRAINTS:
1323 * When doing PR pressure coupling we have to constrain the
1324 * bonds in each iteration. If we are only using Nose-Hoover tcoupling
1325 * it is enough to do this once though, since the relative velocities
1326 * after this will be normal to the bond vector
1331 /* clear out constraints before applying */
1332 clear_mat(vir_part
);
1334 xprime
= get_xprime(state
,upd
);
1336 bLastStep
= (step
== inputrec
->init_step
+inputrec
->nsteps
);
1337 bLog
= (do_per_step(step
,inputrec
->nstlog
) || bLastStep
|| (step
< 0));
1338 bEner
= (do_per_step(step
,inputrec
->nstenergy
) || bLastStep
);
1339 /* Constrain the coordinates xprime */
1340 wallcycle_start(wcycle
,ewcCONSTR
);
1341 if (EI_VV(inputrec
->eI
) && bFirstHalf
)
1343 constrain(NULL
,bLog
,bEner
,constr
,idef
,
1344 inputrec
,ekind
,cr
,step
,1,md
,
1345 state
->x
,state
->v
,state
->v
,
1346 state
->box
,state
->lambda
,dvdlambda
,
1347 NULL
,bCalcVir
? &vir_con
: NULL
,nrnb
,econqVeloc
,
1348 inputrec
->epc
==epcMTTK
,state
->veta
,vetanew
);
1352 constrain(NULL
,bLog
,bEner
,constr
,idef
,
1353 inputrec
,ekind
,cr
,step
,1,md
,
1354 state
->x
,xprime
,NULL
,
1355 state
->box
,state
->lambda
,dvdlambda
,
1356 state
->v
,bCalcVir
? &vir_con
: NULL
,nrnb
,econqCoord
,
1357 inputrec
->epc
==epcMTTK
,state
->veta
,state
->veta
);
1359 wallcycle_stop(wcycle
,ewcCONSTR
);
1363 dump_it_all(fplog
,"After Shake",
1364 state
->natoms
,state
->x
,xprime
,state
->v
,force
);
1368 if (inputrec
->eI
== eiSD2
)
1370 /* A correction factor eph is needed for the SD constraint force */
1371 /* Here we can, unfortunately, not have proper corrections
1372 * for different friction constants, so we use the first one.
1374 for(i
=0; i
<DIM
; i
++)
1376 for(m
=0; m
<DIM
; m
++)
1378 vir_part
[i
][m
] += upd
->sd
->sdc
[0].eph
*vir_con
[i
][m
];
1384 m_add(vir_part
,vir_con
,vir_part
);
1388 pr_rvecs(debug
,0,"constraint virial",vir_part
,DIM
);
1394 if ((inputrec
->eI
== eiSD2
) && !(bFirstHalf
))
1396 xprime
= get_xprime(state
,upd
);
1398 /* The second part of the SD integration */
1399 do_update_sd2(upd
->sd
,FALSE
,start
,homenr
,
1400 inputrec
->opts
.acc
,inputrec
->opts
.nFreeze
,
1401 md
->invmass
,md
->ptype
,
1402 md
->cFREEZE
,md
->cACC
,md
->cTC
,
1403 state
->x
,xprime
,state
->v
,force
,state
->sd_X
,
1404 inputrec
->opts
.ngtc
,inputrec
->opts
.tau_t
,
1405 inputrec
->opts
.ref_t
,FALSE
);
1406 inc_nrnb(nrnb
, eNR_UPDATE
, homenr
);
1410 /* Constrain the coordinates xprime */
1411 wallcycle_start(wcycle
,ewcCONSTR
);
1412 constrain(NULL
,bLog
,bEner
,constr
,idef
,
1413 inputrec
,NULL
,cr
,step
,1,md
,
1414 state
->x
,xprime
,NULL
,
1415 state
->box
,state
->lambda
,dvdlambda
,
1416 NULL
,NULL
,nrnb
,econqCoord
,FALSE
,0,0);
1417 wallcycle_stop(wcycle
,ewcCONSTR
);
1421 /* We must always unshift after updating coordinates; if we did not shake
1422 x was shifted in do_force */
1424 if (!(bFirstHalf
)) /* in the first half of vv, no shift. */
1426 if (graph
&& (graph
->nnodes
> 0))
1428 unshift_x(graph
,state
->box
,state
->x
,upd
->xp
);
1429 if (TRICLINIC(state
->box
))
1431 inc_nrnb(nrnb
,eNR_SHIFTX
,2*graph
->nnodes
);
1435 inc_nrnb(nrnb
,eNR_SHIFTX
,graph
->nnodes
);
1437 copy_rvecn(upd
->xp
,state
->x
,start
,graph
->start
);
1438 copy_rvecn(upd
->xp
,state
->x
,graph
->start
+graph
->nnodes
,nrend
);
1442 copy_rvecn(upd
->xp
,state
->x
,start
,nrend
);
1445 dump_it_all(fplog
,"After unshift",
1446 state
->natoms
,state
->x
,upd
->xp
,state
->v
,force
);
1448 /* ############# END the update of velocities and positions ######### */
1451 void update_box(FILE *fplog
,
1452 gmx_large_int_t step
,
1453 t_inputrec
*inputrec
, /* input record and box stuff */
1457 rvec force
[], /* forces on home particles */
1461 gmx_wallcycle_t wcycle
,
1464 gmx_bool bFirstHalf
)
1466 gmx_bool bExtended
,bTrotter
,bLastStep
,bLog
=FALSE
,bEner
=FALSE
;
1469 int start
,homenr
,nrend
,i
,n
,m
,g
;
1473 homenr
= md
->homenr
;
1474 nrend
= start
+homenr
;
1477 (inputrec
->etc
== etcNOSEHOOVER
) ||
1478 (inputrec
->epc
== epcPARRINELLORAHMAN
) ||
1479 (inputrec
->epc
== epcMTTK
);
1481 dt
= inputrec
->delta_t
;
1485 /* now update boxes */
1486 switch (inputrec
->epc
) {
1489 case (epcBERENDSEN
):
1490 berendsen_pscale(inputrec
,pcoupl_mu
,state
->box
,state
->box_rel
,
1491 start
,homenr
,state
->x
,md
->cFREEZE
,nrnb
);
1493 case (epcPARRINELLORAHMAN
):
1494 /* The box velocities were updated in do_pr_pcoupl in the update
1495 * iteration, but we dont change the box vectors until we get here
1496 * since we need to be able to shift/unshift above.
1502 state
->box
[i
][m
] += dt
*state
->boxv
[i
][m
];
1505 preserve_box_shape(inputrec
,state
->box_rel
,state
->box
);
1507 /* Scale the coordinates */
1508 for(n
=start
; (n
<start
+homenr
); n
++)
1510 tmvmul_ur0(pcoupl_mu
,state
->x
[n
],state
->x
[n
]);
1514 switch (inputrec
->epct
)
1516 case (epctISOTROPIC
):
1517 /* DIM * eta = ln V. so DIM*eta_new = DIM*eta_old + DIM*dt*veta =>
1518 ln V_new = ln V_old + 3*dt*veta => V_new = V_old*exp(3*dt*veta) =>
1519 Side length scales as exp(veta*dt) */
1521 msmul(state
->box
,exp(state
->veta
*dt
),state
->box
);
1523 /* Relate veta to boxv. veta = d(eta)/dT = (1/DIM)*1/V dV/dT.
1524 o If we assume isotropic scaling, and box length scaling
1525 factor L, then V = L^DIM (det(M)). So dV/dt = DIM
1526 L^(DIM-1) dL/dt det(M), and veta = (1/L) dL/dt. The
1527 determinant of B is L^DIM det(M), and the determinant
1528 of dB/dt is (dL/dT)^DIM det (M). veta will be
1529 (det(dB/dT)/det(B))^(1/3). Then since M =
1530 B_new*(vol_new)^(1/3), dB/dT_new = (veta_new)*B(new). */
1532 msmul(state
->box
,state
->veta
,state
->boxv
);
1542 if ((!(IR_NPT_TROTTER(inputrec
))) && scale_tot
)
1544 /* The transposes of the scaling matrices are stored,
1545 * therefore we need to reverse the order in the multiplication.
1547 mmul_ur0(*scale_tot
,pcoupl_mu
,*scale_tot
);
1550 if (DEFORM(*inputrec
))
1552 deform(upd
,start
,homenr
,state
->x
,state
->box
,scale_tot
,inputrec
,step
);
1555 dump_it_all(fplog
,"After update",
1556 state
->natoms
,state
->x
,upd
->xp
,state
->v
,force
);
1559 void update_coords(FILE *fplog
,
1560 gmx_large_int_t step
,
1561 t_inputrec
*inputrec
, /* input record and box stuff */
1564 rvec
*f
, /* forces on home particles */
1568 gmx_ekindata_t
*ekind
,
1570 gmx_wallcycle_t wcycle
,
1574 t_commrec
*cr
, /* these shouldn't be here -- need to think about it */
1576 gmx_constr_t constr
,
1579 gmx_bool bExtended
,bNH
,bPR
,bTrotter
,bLastStep
,bLog
=FALSE
,bEner
=FALSE
;
1581 real
*imass
,*imassin
;
1584 int start
,homenr
,nrend
,i
,j
,d
,n
,m
,g
;
1585 int blen0
,blen1
,iatom
,jatom
,nshake
,nsettle
,nconstr
,nexpand
;
1588 rvec
*vcom
,*xcom
,*vall
,*xall
,*xin
,*vin
,*forcein
,*fall
,*xpall
,*xprimein
,*xprime
;
1591 /* Running the velocity half does nothing except for velocity verlet */
1592 if ((UpdatePart
== etrtVELOCITY1
|| UpdatePart
== etrtVELOCITY2
) &&
1593 !EI_VV(inputrec
->eI
))
1595 gmx_incons("update_coords called for velocity without VV integrator");
1599 homenr
= md
->homenr
;
1600 nrend
= start
+homenr
;
1602 xprime
= get_xprime(state
,upd
);
1604 dt
= inputrec
->delta_t
;
1607 /* We need to update the NMR restraint history when time averaging is used */
1608 if (state
->flags
& (1<<estDISRE_RM3TAV
))
1610 update_disres_history(fcd
,&state
->hist
);
1612 if (state
->flags
& (1<<estORIRE_DTAV
))
1614 update_orires_history(fcd
,&state
->hist
);
1617 bNH
= inputrec
->etc
== etcNOSEHOOVER
;
1618 bPR
= ((inputrec
->epc
== epcPARRINELLORAHMAN
) || (inputrec
->epc
== epcMTTK
));
1620 bExtended
= bNH
|| bPR
;
1622 if (bDoLR
&& inputrec
->nstlist
> 1 && !EI_VV(inputrec
->eI
)) /* get this working with VV? */
1624 /* Store the total force + nstlist-1 times the LR force
1625 * in forces_lr, so it can be used in a normal update algorithm
1626 * to produce twin time stepping.
1628 /* is this correct in the new construction? MRS */
1629 combine_forces(inputrec
->nstlist
,constr
,inputrec
,md
,idef
,cr
,step
,state
,
1630 start
,nrend
,f
,f_lr
,nrnb
);
1638 /* ############# START The update of velocities and positions ######### */
1640 dump_it_all(fplog
,"Before update",
1641 state
->natoms
,state
->x
,xprime
,state
->v
,force
);
1643 switch (inputrec
->eI
) {
1645 if (ekind
->cosacc
.cos_accel
== 0) {
1646 /* use normal version of update */
1647 do_update_md(start
,nrend
,dt
,
1648 ekind
->tcstat
,ekind
->grpstat
,state
->nosehoover_vxi
,
1649 inputrec
->opts
.acc
,inputrec
->opts
.nFreeze
,md
->invmass
,md
->ptype
,
1650 md
->cFREEZE
,md
->cACC
,md
->cTC
,
1651 state
->x
,xprime
,state
->v
,force
,M
,
1656 do_update_visc(start
,nrend
,dt
,
1657 ekind
->tcstat
,md
->invmass
,state
->nosehoover_vxi
,
1658 md
->ptype
,md
->cTC
,state
->x
,xprime
,state
->v
,force
,M
,
1659 state
->box
,ekind
->cosacc
.cos_accel
,ekind
->cosacc
.vcos
,bNH
,bPR
);
1663 do_update_sd1(upd
->sd
,start
,homenr
,dt
,
1664 inputrec
->opts
.acc
,inputrec
->opts
.nFreeze
,
1665 md
->invmass
,md
->ptype
,
1666 md
->cFREEZE
,md
->cACC
,md
->cTC
,
1667 state
->x
,xprime
,state
->v
,force
,state
->sd_X
,
1668 inputrec
->opts
.ngtc
,inputrec
->opts
.tau_t
,inputrec
->opts
.ref_t
);
1671 /* The SD update is done in 2 parts, because an extra constraint step
1674 do_update_sd2(upd
->sd
,bInitStep
,start
,homenr
,
1675 inputrec
->opts
.acc
,inputrec
->opts
.nFreeze
,
1676 md
->invmass
,md
->ptype
,
1677 md
->cFREEZE
,md
->cACC
,md
->cTC
,
1678 state
->x
,xprime
,state
->v
,force
,state
->sd_X
,
1679 inputrec
->opts
.ngtc
,inputrec
->opts
.tau_t
,inputrec
->opts
.ref_t
,
1683 do_update_bd(start
,nrend
,dt
,
1684 inputrec
->opts
.nFreeze
,md
->invmass
,md
->ptype
,
1685 md
->cFREEZE
,md
->cTC
,
1686 state
->x
,xprime
,state
->v
,force
,
1688 inputrec
->opts
.ngtc
,inputrec
->opts
.tau_t
,inputrec
->opts
.ref_t
,
1689 upd
->sd
->bd_rf
,upd
->sd
->gaussrand
);
1693 alpha
= 1.0 + DIM
/((double)inputrec
->opts
.nrdf
[0]); /* assuming barostat coupled to group 0. */
1694 switch (UpdatePart
) {
1697 do_update_vv_vel(start
,nrend
,dt
,
1698 ekind
->tcstat
,ekind
->grpstat
,
1699 inputrec
->opts
.acc
,inputrec
->opts
.nFreeze
,
1700 md
->invmass
,md
->ptype
,
1701 md
->cFREEZE
,md
->cACC
,
1703 bExtended
,state
->veta
,alpha
);
1706 do_update_vv_pos(start
,nrend
,dt
,
1707 ekind
->tcstat
,ekind
->grpstat
,
1708 inputrec
->opts
.acc
,inputrec
->opts
.nFreeze
,
1709 md
->invmass
,md
->ptype
,md
->cFREEZE
,
1710 state
->x
,xprime
,state
->v
,force
,
1711 bExtended
,state
->veta
,alpha
);
1716 gmx_fatal(FARGS
,"Don't know how to update coordinates");
1722 void correct_ekin(FILE *log
,int start
,int end
,rvec v
[],rvec vcm
,real mass
[],
1723 real tmass
,tensor ekin
)
1726 * This is a debugging routine. It should not be called for production code
1728 * The kinetic energy should calculated according to:
1729 * Ekin = 1/2 m (v-vcm)^2
1730 * However the correction is not always applied, since vcm may not be
1731 * known in time and we compute
1732 * Ekin' = 1/2 m v^2 instead
1733 * This can be corrected afterwards by computing
1734 * Ekin = Ekin' + 1/2 m ( -2 v vcm + vcm^2)
1736 * Ekin = Ekin' - m v vcm + 1/2 m vcm^2
1743 /* Local particles */
1746 /* Processor dependent part. */
1748 for(i
=start
; (i
<end
); i
++)
1752 for(j
=0; (j
<DIM
); j
++)
1758 svmul(1/tmass
,vcm
,vcm
);
1759 svmul(0.5,vcm
,hvcm
);
1761 for(j
=0; (j
<DIM
); j
++)
1763 for(k
=0; (k
<DIM
); k
++)
1765 dekin
[j
][k
] += vcm
[k
]*(tm
*hvcm
[j
]-mv
[j
]);
1768 pr_rvecs(log
,0,"dekin",dekin
,DIM
);
1769 pr_rvecs(log
,0," ekin", ekin
,DIM
);
1770 fprintf(log
,"dekin = %g, ekin = %g vcm = (%8.4f %8.4f %8.4f)\n",
1771 trace(dekin
),trace(ekin
),vcm
[XX
],vcm
[YY
],vcm
[ZZ
]);
1772 fprintf(log
,"mv = (%8.4f %8.4f %8.4f)\n",
1773 mv
[XX
],mv
[YY
],mv
[ZZ
]);