?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fixing handling of perturbation mass changes.
[gromacs.git] / src / mdlib / update.c
blob4c32f9661855df6c92560148cb0b1a93216fe588
1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012,2013, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
11 *
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
16 *
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
21 *
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 *
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
34 *
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
37 */
38 #ifdef HAVE_CONFIG_H
39 #include <config.h>
40 #endif
41
42
43 #include <stdio.h>
44 #include <math.h>
45
46 #include "types/commrec.h"
47 #include "sysstuff.h"
48 #include "smalloc.h"
49 #include "typedefs.h"
50 #include "nrnb.h"
51 #include "physics.h"
52 #include "macros.h"
53 #include "vec.h"
54 #include "main.h"
55 #include "confio.h"
56 #include "update.h"
57 #include "gmx_random.h"
58 #include "futil.h"
59 #include "mshift.h"
60 #include "tgroup.h"
61 #include "force.h"
62 #include "names.h"
63 #include "txtdump.h"
64 #include "mdrun.h"
65 #include "copyrite.h"
66 #include "constr.h"
67 #include "edsam.h"
68 #include "pull.h"
69 #include "disre.h"
70 #include "orires.h"
71 #include "gmx_wallcycle.h"
72 #include "gmx_omp_nthreads.h"
73 #include "gmx_omp.h"
74
75 /*For debugging, start at v(-dt/2) for velolcity verlet -- uncomment next line */
76 /*#define STARTFROMDT2*/
77
78 typedef struct {
79 double gdt;
80 double eph;
81 double emh;
82 double em;
83 double b;
84 double c;
85 double d;
86 } gmx_sd_const_t;
87
88 typedef struct {
89 real V;
90 real X;
91 real Yv;
92 real Yx;
93 } gmx_sd_sigma_t;
94
95 typedef struct {
96 /* The random state for ngaussrand threads.
97 * Normal thermostats need just 1 random number generator,
98 * but SD and BD with OpenMP parallelization need 1 for each thread.
99 */
100 int ngaussrand;
101 gmx_rng_t *gaussrand;
102 /* BD stuff */
103 real *bd_rf;
104 /* SD stuff */
105 gmx_sd_const_t *sdc;
106 gmx_sd_sigma_t *sdsig;
107 rvec *sd_V;
108 int sd_V_nalloc;
109 /* andersen temperature control stuff */
110 gmx_bool *randomize_group;
111 real *boltzfac;
112 } gmx_stochd_t;
113
114 typedef struct gmx_update
115 {
116 gmx_stochd_t *sd;
117 /* xprime for constraint algorithms */
118 rvec *xp;
119 int xp_nalloc;
120
121 /* variable size arrays for andersen */
122 gmx_bool *randatom;
123 int *randatom_list;
124 gmx_bool randatom_list_init;
125
126 /* Variables for the deform algorithm */
127 gmx_large_int_t deformref_step;
128 matrix deformref_box;
129 } t_gmx_update;
130
131
132 static void do_update_md(int start, int nrend, double dt,
133 t_grp_tcstat *tcstat,
134 double nh_vxi[],
135 gmx_bool bNEMD, t_grp_acc *gstat, rvec accel[],
136 ivec nFreeze[],
137 real invmass[],
138 unsigned short ptype[], unsigned short cFREEZE[],
139 unsigned short cACC[], unsigned short cTC[],
140 rvec x[], rvec xprime[], rvec v[],
141 rvec f[], matrix M,
142 gmx_bool bNH, gmx_bool bPR)
143 {
144 double imass, w_dt;
145 int gf = 0, ga = 0, gt = 0;
146 rvec vrel;
147 real vn, vv, va, vb, vnrel;
148 real lg, vxi = 0, u;
149 int n, d;
150
151 if (bNH || bPR)
152 {
153 /* Update with coupling to extended ensembles, used for
154 * Nose-Hoover and Parrinello-Rahman coupling
155 * Nose-Hoover uses the reversible leap-frog integrator from
156 * Holian et al. Phys Rev E 52(3) : 2338, 1995
157 */
158 for (n = start; n < nrend; n++)
159 {
160 imass = invmass[n];
161 if (cFREEZE)
162 {
163 gf = cFREEZE[n];
164 }
165 if (cACC)
166 {
167 ga = cACC[n];
168 }
169 if (cTC)
170 {
171 gt = cTC[n];
172 }
173 lg = tcstat[gt].lambda;
174 if (bNH)
175 {
176 vxi = nh_vxi[gt];
177 }
178 rvec_sub(v[n], gstat[ga].u, vrel);
179
180 for (d = 0; d < DIM; d++)
181 {
182 if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
183 {
184 vnrel = (lg*vrel[d] + dt*(imass*f[n][d] - 0.5*vxi*vrel[d]
185 - iprod(M[d], vrel)))/(1 + 0.5*vxi*dt);
186 /* do not scale the mean velocities u */
187 vn = gstat[ga].u[d] + accel[ga][d]*dt + vnrel;
188 v[n][d] = vn;
189 xprime[n][d] = x[n][d]+vn*dt;
190 }
191 else
192 {
193 v[n][d] = 0.0;
194 xprime[n][d] = x[n][d];
195 }
196 }
197 }
198 }
199 else if (cFREEZE != NULL ||
200 nFreeze[0][XX] || nFreeze[0][YY] || nFreeze[0][ZZ] ||
201 bNEMD)
202 {
203 /* Update with Berendsen/v-rescale coupling and freeze or NEMD */
204 for (n = start; n < nrend; n++)
205 {
206 w_dt = invmass[n]*dt;
207 if (cFREEZE)
208 {
209 gf = cFREEZE[n];
210 }
211 if (cACC)
212 {
213 ga = cACC[n];
214 }
215 if (cTC)
216 {
217 gt = cTC[n];
218 }
219 lg = tcstat[gt].lambda;
220
221 for (d = 0; d < DIM; d++)
222 {
223 vn = v[n][d];
224 if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
225 {
226 vv = lg*vn + f[n][d]*w_dt;
227
228 /* do not scale the mean velocities u */
229 u = gstat[ga].u[d];
230 va = vv + accel[ga][d]*dt;
231 vb = va + (1.0-lg)*u;
232 v[n][d] = vb;
233 xprime[n][d] = x[n][d]+vb*dt;
234 }
235 else
236 {
237 v[n][d] = 0.0;
238 xprime[n][d] = x[n][d];
239 }
240 }
241 }
242 }
243 else
244 {
245 /* Plain update with Berendsen/v-rescale coupling */
246 for (n = start; n < nrend; n++)
247 {
248 if ((ptype[n] != eptVSite) && (ptype[n] != eptShell))
249 {
250 w_dt = invmass[n]*dt;
251 if (cTC)
252 {
253 gt = cTC[n];
254 }
255 lg = tcstat[gt].lambda;
256
257 for (d = 0; d < DIM; d++)
258 {
259 vn = lg*v[n][d] + f[n][d]*w_dt;
260 v[n][d] = vn;
261 xprime[n][d] = x[n][d] + vn*dt;
262 }
263 }
264 else
265 {
266 for (d = 0; d < DIM; d++)
267 {
268 v[n][d] = 0.0;
269 xprime[n][d] = x[n][d];
270 }
271 }
272 }
273 }
274 }
275
276 static void do_update_vv_vel(int start, int nrend, double dt,
277 t_grp_tcstat *tcstat, t_grp_acc *gstat,
278 rvec accel[], ivec nFreeze[], real invmass[],
279 unsigned short ptype[], unsigned short cFREEZE[],
280 unsigned short cACC[], rvec v[], rvec f[],
281 gmx_bool bExtended, real veta, real alpha)
282 {
283 double imass, w_dt;
284 int gf = 0, ga = 0;
285 rvec vrel;
286 real u, vn, vv, va, vb, vnrel;
287 int n, d;
288 double g, mv1, mv2;
289
290 if (bExtended)
291 {
292 g = 0.25*dt*veta*alpha;
293 mv1 = exp(-g);
294 mv2 = series_sinhx(g);
295 }
296 else
297 {
298 mv1 = 1.0;
299 mv2 = 1.0;
300 }
301 for (n = start; n < nrend; n++)
302 {
303 w_dt = invmass[n]*dt;
304 if (cFREEZE)
305 {
306 gf = cFREEZE[n];
307 }
308 if (cACC)
309 {
310 ga = cACC[n];
311 }
312
313 for (d = 0; d < DIM; d++)
314 {
315 if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
316 {
317 v[n][d] = mv1*(mv1*v[n][d] + 0.5*(w_dt*mv2*f[n][d]))+0.5*accel[ga][d]*dt;
318 }
319 else
320 {
321 v[n][d] = 0.0;
322 }
323 }
324 }
325 } /* do_update_vv_vel */
326
327 static void do_update_vv_pos(int start, int nrend, double dt,
328 t_grp_tcstat *tcstat, t_grp_acc *gstat,
329 rvec accel[], ivec nFreeze[], real invmass[],
330 unsigned short ptype[], unsigned short cFREEZE[],
331 rvec x[], rvec xprime[], rvec v[],
332 rvec f[], gmx_bool bExtended, real veta, real alpha)
333 {
334 double imass, w_dt;
335 int gf = 0;
336 int n, d;
337 double g, mr1, mr2;
338
339 /* Would it make more sense if Parrinello-Rahman was put here? */
340 if (bExtended)
341 {
342 g = 0.5*dt*veta;
343 mr1 = exp(g);
344 mr2 = series_sinhx(g);
345 }
346 else
347 {
348 mr1 = 1.0;
349 mr2 = 1.0;
350 }
351
352 for (n = start; n < nrend; n++)
353 {
354
355 if (cFREEZE)
356 {
357 gf = cFREEZE[n];
358 }
359
360 for (d = 0; d < DIM; d++)
361 {
362 if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
363 {
364 xprime[n][d] = mr1*(mr1*x[n][d]+mr2*dt*v[n][d]);
365 }
366 else
367 {
368 xprime[n][d] = x[n][d];
369 }
370 }
371 }
372 } /* do_update_vv_pos */
373
374 static void do_update_visc(int start, int nrend, double dt,
375 t_grp_tcstat *tcstat,
376 double nh_vxi[],
377 real invmass[],
378 unsigned short ptype[], unsigned short cTC[],
379 rvec x[], rvec xprime[], rvec v[],
380 rvec f[], matrix M, matrix box, real
381 cos_accel, real vcos,
382 gmx_bool bNH, gmx_bool bPR)
383 {
384 double imass, w_dt;
385 int gt = 0;
386 real vn, vc;
387 real lg, vxi = 0, vv;
388 real fac, cosz;
389 rvec vrel;
390 int n, d;
391
392 fac = 2*M_PI/(box[ZZ][ZZ]);
393
394 if (bNH || bPR)
395 {
396 /* Update with coupling to extended ensembles, used for
397 * Nose-Hoover and Parrinello-Rahman coupling
398 */
399 for (n = start; n < nrend; n++)
400 {
401 imass = invmass[n];
402 if (cTC)
403 {
404 gt = cTC[n];
405 }
406 lg = tcstat[gt].lambda;
407 cosz = cos(fac*x[n][ZZ]);
408
409 copy_rvec(v[n], vrel);
410
411 vc = cosz*vcos;
412 vrel[XX] -= vc;
413 if (bNH)
414 {
415 vxi = nh_vxi[gt];
416 }
417 for (d = 0; d < DIM; d++)
418 {
419 vn = v[n][d];
420
421 if ((ptype[n] != eptVSite) && (ptype[n] != eptShell))
422 {
423 vn = (lg*vrel[d] + dt*(imass*f[n][d] - 0.5*vxi*vrel[d]
424 - iprod(M[d], vrel)))/(1 + 0.5*vxi*dt);
425 if (d == XX)
426 {
427 vn += vc + dt*cosz*cos_accel;
428 }
429 v[n][d] = vn;
430 xprime[n][d] = x[n][d]+vn*dt;
431 }
432 else
433 {
434 xprime[n][d] = x[n][d];
435 }
436 }
437 }
438 }
439 else
440 {
441 /* Classic version of update, used with berendsen coupling */
442 for (n = start; n < nrend; n++)
443 {
444 w_dt = invmass[n]*dt;
445 if (cTC)
446 {
447 gt = cTC[n];
448 }
449 lg = tcstat[gt].lambda;
450 cosz = cos(fac*x[n][ZZ]);
451
452 for (d = 0; d < DIM; d++)
453 {
454 vn = v[n][d];
455
456 if ((ptype[n] != eptVSite) && (ptype[n] != eptShell))
457 {
458 if (d == XX)
459 {
460 vc = cosz*vcos;
461 /* Do not scale the cosine velocity profile */
462 vv = vc + lg*(vn - vc + f[n][d]*w_dt);
463 /* Add the cosine accelaration profile */
464 vv += dt*cosz*cos_accel;
465 }
466 else
467 {
468 vv = lg*(vn + f[n][d]*w_dt);
469 }
470 v[n][d] = vv;
471 xprime[n][d] = x[n][d]+vv*dt;
472 }
473 else
474 {
475 v[n][d] = 0.0;
476 xprime[n][d] = x[n][d];
477 }
478 }
479 }
480 }
481 }
482
483 /* Allocates and initializes sd->gaussrand[i] for i=1, i<sd->ngaussrand,
484 * Using seeds generated from sd->gaussrand[0].
485 */
486 static void init_multiple_gaussrand(gmx_stochd_t *sd)
487 {
488 int ngr, i;
489 unsigned int *seed;
490
491 ngr = sd->ngaussrand;
492 snew(seed, ngr);
493
494 for (i = 1; i < ngr; i++)
495 {
496 seed[i] = gmx_rng_uniform_uint32(sd->gaussrand[0]);
497 }
498
499 #pragma omp parallel num_threads(ngr)
500 {
501 int th;
502
503 th = gmx_omp_get_thread_num();
504 if (th > 0)
505 {
506 /* Initialize on each thread to have thread-local memory alloced */
507 sd->gaussrand[th] = gmx_rng_init(seed[th]);
508 }
509 }
510
511 sfree(seed);
512 }
513
514 static gmx_stochd_t *init_stochd(FILE *fplog, t_inputrec *ir, int nthreads)
515 {
516 gmx_stochd_t *sd;
517 gmx_sd_const_t *sdc;
518 int ngtc, n, th;
519 real y;
520
521 snew(sd, 1);
522
523 /* Initiate random number generator for langevin type dynamics,
524 * for BD, SD or velocity rescaling temperature coupling.
525 */
526 if (ir->eI == eiBD || EI_SD(ir->eI))
527 {
528 sd->ngaussrand = nthreads;
529 }
530 else
531 {
532 sd->ngaussrand = 1;
533 }
534 snew(sd->gaussrand, sd->ngaussrand);
535
536 /* Initialize the first random generator */
537 sd->gaussrand[0] = gmx_rng_init(ir->ld_seed);
538
539 if (sd->ngaussrand > 1)
540 {
541 /* Initialize the rest of the random number generators,
542 * using the first one to generate seeds.
543 */
544 init_multiple_gaussrand(sd);
545 }
546
547 ngtc = ir->opts.ngtc;
548
549 if (ir->eI == eiBD)
550 {
551 snew(sd->bd_rf, ngtc);
552 }
553 else if (EI_SD(ir->eI))
554 {
555 snew(sd->sdc, ngtc);
556 snew(sd->sdsig, ngtc);
557
558 sdc = sd->sdc;
559 for (n = 0; n < ngtc; n++)
560 {
561 if (ir->opts.tau_t[n] > 0)
562 {
563 sdc[n].gdt = ir->delta_t/ir->opts.tau_t[n];
564 sdc[n].eph = exp(sdc[n].gdt/2);
565 sdc[n].emh = exp(-sdc[n].gdt/2);
566 sdc[n].em = exp(-sdc[n].gdt);
567 }
568 else
569 {
570 /* No friction and noise on this group */
571 sdc[n].gdt = 0;
572 sdc[n].eph = 1;
573 sdc[n].emh = 1;
574 sdc[n].em = 1;
575 }
576 if (sdc[n].gdt >= 0.05)
577 {
578 sdc[n].b = sdc[n].gdt*(sdc[n].eph*sdc[n].eph - 1)
579 - 4*(sdc[n].eph - 1)*(sdc[n].eph - 1);
580 sdc[n].c = sdc[n].gdt - 3 + 4*sdc[n].emh - sdc[n].em;
581 sdc[n].d = 2 - sdc[n].eph - sdc[n].emh;
582 }
583 else
584 {
585 y = sdc[n].gdt/2;
586 /* Seventh order expansions for small y */
587 sdc[n].b = y*y*y*y*(1/3.0+y*(1/3.0+y*(17/90.0+y*7/9.0)));
588 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))));
589 sdc[n].d = y*y*(-1+y*y*(-1/12.0-y*y/360.0));
590 }
591 if (debug)
592 {
593 fprintf(debug, "SD const tc-grp %d: b %g c %g d %g\n",
594 n, sdc[n].b, sdc[n].c, sdc[n].d);
595 }
596 }
597 }
598 else if (ETC_ANDERSEN(ir->etc))
599 {
600 int ngtc;
601 t_grpopts *opts;
602 real reft;
603
604 opts = &ir->opts;
605 ngtc = opts->ngtc;
606
607 snew(sd->randomize_group, ngtc);
608 snew(sd->boltzfac, ngtc);
609
610 /* for now, assume that all groups, if randomized, are randomized at the same rate, i.e. tau_t is the same. */
611 /* since constraint groups don't necessarily match up with temperature groups! This is checked in readir.c */
612
613 for (n = 0; n < ngtc; n++)
614 {
615 reft = max(0.0, opts->ref_t[n]);
616 if ((opts->tau_t[n] > 0) && (reft > 0)) /* tau_t or ref_t = 0 means that no randomization is done */
617 {
618 sd->randomize_group[n] = TRUE;
619 sd->boltzfac[n] = BOLTZ*opts->ref_t[n];
620 }
621 else
622 {
623 sd->randomize_group[n] = FALSE;
624 }
625 }
626 }
627 return sd;
628 }
629
630 void get_stochd_state(gmx_update_t upd, t_state *state)
631 {
632 /* Note that we only get the state of the first random generator,
633 * even if there are multiple. This avoids repetition.
634 */
635 gmx_rng_get_state(upd->sd->gaussrand[0], state->ld_rng, state->ld_rngi);
636 }
637
638 void set_stochd_state(gmx_update_t upd, t_state *state)
639 {
640 gmx_stochd_t *sd;
641 int i;
642
643 sd = upd->sd;
644
645 gmx_rng_set_state(sd->gaussrand[0], state->ld_rng, state->ld_rngi[0]);
646
647 if (sd->ngaussrand > 1)
648 {
649 /* We only end up here with SD or BD with OpenMP.
650 * Destroy and reinitialize the rest of the random number generators,
651 * using seeds generated from the first one.
652 * Although this doesn't recover the previous state,
653 * it at least avoids repetition, which is most important.
654 * Exaclty restoring states with all MPI+OpenMP setups is difficult
655 * and as the integrator is random to start with, doesn't gain us much.
656 */
657 for (i = 1; i < sd->ngaussrand; i++)
658 {
659 gmx_rng_destroy(sd->gaussrand[i]);
660 }
661
662 init_multiple_gaussrand(sd);
663 }
664 }
665
666 gmx_update_t init_update(FILE *fplog, t_inputrec *ir)
667 {
668 t_gmx_update *upd;
669
670 snew(upd, 1);
671
672 if (ir->eI == eiBD || EI_SD(ir->eI) || ir->etc == etcVRESCALE || ETC_ANDERSEN(ir->etc))
673 {
674 upd->sd = init_stochd(fplog, ir, gmx_omp_nthreads_get(emntUpdate));
675 }
676
677 upd->xp = NULL;
678 upd->xp_nalloc = 0;
679 upd->randatom = NULL;
680 upd->randatom_list = NULL;
681 upd->randatom_list_init = FALSE; /* we have not yet cleared the data structure at this point */
682
683 return upd;
684 }
685
686 static void do_update_sd1(gmx_stochd_t *sd,
687 gmx_rng_t gaussrand,
688 int start, int nrend, double dt,
689 rvec accel[], ivec nFreeze[],
690 real invmass[], unsigned short ptype[],
691 unsigned short cFREEZE[], unsigned short cACC[],
692 unsigned short cTC[],
693 rvec x[], rvec xprime[], rvec v[], rvec f[],
694 rvec sd_X[],
695 int ngtc, real tau_t[], real ref_t[])
696 {
697 gmx_sd_const_t *sdc;
698 gmx_sd_sigma_t *sig;
699 real kT;
700 int gf = 0, ga = 0, gt = 0;
701 real ism, sd_V;
702 int n, d;
703
704 sdc = sd->sdc;
705 sig = sd->sdsig;
706
707 for (n = 0; n < ngtc; n++)
708 {
709 kT = BOLTZ*ref_t[n];
710 /* The mass is encounted for later, since this differs per atom */
711 sig[n].V = sqrt(kT*(1 - sdc[n].em*sdc[n].em));
712 }
713
714 for (n = start; n < nrend; n++)
715 {
716 ism = sqrt(invmass[n]);
717 if (cFREEZE)
718 {
719 gf = cFREEZE[n];
720 }
721 if (cACC)
722 {
723 ga = cACC[n];
724 }
725 if (cTC)
726 {
727 gt = cTC[n];
728 }
729
730 for (d = 0; d < DIM; d++)
731 {
732 if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
733 {
734 sd_V = ism*sig[gt].V*gmx_rng_gaussian_table(gaussrand);
735
736 v[n][d] = v[n][d]*sdc[gt].em
737 + (invmass[n]*f[n][d] + accel[ga][d])*tau_t[gt]*(1 - sdc[gt].em)
738 + sd_V;
739
740 xprime[n][d] = x[n][d] + v[n][d]*dt;
741 }
742 else
743 {
744 v[n][d] = 0.0;
745 xprime[n][d] = x[n][d];
746 }
747 }
748 }
749 }
750
751 static void check_sd2_work_data_allocation(gmx_stochd_t *sd, int nrend)
752 {
753 if (nrend > sd->sd_V_nalloc)
754 {
755 sd->sd_V_nalloc = over_alloc_dd(nrend);
756 srenew(sd->sd_V, sd->sd_V_nalloc);
757 }
758 }
759
760 static void do_update_sd2(gmx_stochd_t *sd,
761 gmx_rng_t gaussrand,
762 gmx_bool bInitStep,
763 int start, int nrend,
764 rvec accel[], ivec nFreeze[],
765 real invmass[], unsigned short ptype[],
766 unsigned short cFREEZE[], unsigned short cACC[],
767 unsigned short cTC[],
768 rvec x[], rvec xprime[], rvec v[], rvec f[],
769 rvec sd_X[],
770 int ngtc, real tau_t[], real ref_t[],
771 gmx_bool bFirstHalf)
772 {
773 gmx_sd_const_t *sdc;
774 gmx_sd_sigma_t *sig;
775 /* The random part of the velocity update, generated in the first
776 * half of the update, needs to be remembered for the second half.
777 */
778 rvec *sd_V;
779 real kT;
780 int gf = 0, ga = 0, gt = 0;
781 real vn = 0, Vmh, Xmh;
782 real ism;
783 int n, d;
784
785 sdc = sd->sdc;
786 sig = sd->sdsig;
787 sd_V = sd->sd_V;
788
789 if (bFirstHalf)
790 {
791 for (n = 0; n < ngtc; n++)
792 {
793 kT = BOLTZ*ref_t[n];
794 /* The mass is encounted for later, since this differs per atom */
795 sig[n].V = sqrt(kT*(1-sdc[n].em));
796 sig[n].X = sqrt(kT*sqr(tau_t[n])*sdc[n].c);
797 sig[n].Yv = sqrt(kT*sdc[n].b/sdc[n].c);
798 sig[n].Yx = sqrt(kT*sqr(tau_t[n])*sdc[n].b/(1-sdc[n].em));
799 }
800 }
801
802 for (n = start; n < nrend; n++)
803 {
804 ism = sqrt(invmass[n]);
805 if (cFREEZE)
806 {
807 gf = cFREEZE[n];
808 }
809 if (cACC)
810 {
811 ga = cACC[n];
812 }
813 if (cTC)
814 {
815 gt = cTC[n];
816 }
817
818 for (d = 0; d < DIM; d++)
819 {
820 if (bFirstHalf)
821 {
822 vn = v[n][d];
823 }
824 if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
825 {
826 if (bFirstHalf)
827 {
828 if (bInitStep)
829 {
830 sd_X[n][d] = ism*sig[gt].X*gmx_rng_gaussian_table(gaussrand);
831 }
832 Vmh = sd_X[n][d]*sdc[gt].d/(tau_t[gt]*sdc[gt].c)
833 + ism*sig[gt].Yv*gmx_rng_gaussian_table(gaussrand);
834 sd_V[n][d] = ism*sig[gt].V*gmx_rng_gaussian_table(gaussrand);
835
836 v[n][d] = vn*sdc[gt].em
837 + (invmass[n]*f[n][d] + accel[ga][d])*tau_t[gt]*(1 - sdc[gt].em)
838 + sd_V[n][d] - sdc[gt].em*Vmh;
839
840 xprime[n][d] = x[n][d] + v[n][d]*tau_t[gt]*(sdc[gt].eph - sdc[gt].emh);
841 }
842 else
843 {
844
845 /* Correct the velocities for the constraints.
846 * This operation introduces some inaccuracy,
847 * since the velocity is determined from differences in coordinates.
848 */
849 v[n][d] =
850 (xprime[n][d] - x[n][d])/(tau_t[gt]*(sdc[gt].eph - sdc[gt].emh));
851
852 Xmh = sd_V[n][d]*tau_t[gt]*sdc[gt].d/(sdc[gt].em-1)
853 + ism*sig[gt].Yx*gmx_rng_gaussian_table(gaussrand);
854 sd_X[n][d] = ism*sig[gt].X*gmx_rng_gaussian_table(gaussrand);
855
856 xprime[n][d] += sd_X[n][d] - Xmh;
857
858 }
859 }
860 else
861 {
862 if (bFirstHalf)
863 {
864 v[n][d] = 0.0;
865 xprime[n][d] = x[n][d];
866 }
867 }
868 }
869 }
870 }
871
872 static void do_update_bd(int start, int nrend, double dt,
873 ivec nFreeze[],
874 real invmass[], unsigned short ptype[],
875 unsigned short cFREEZE[], unsigned short cTC[],
876 rvec x[], rvec xprime[], rvec v[],
877 rvec f[], real friction_coefficient,
878 int ngtc, real tau_t[], real ref_t[],
879 real *rf, gmx_rng_t gaussrand)
880 {
881 /* note -- these appear to be full step velocities . . . */
882 int gf = 0, gt = 0;
883 real vn;
884 real invfr = 0;
885 int n, d;
886
887 if (friction_coefficient != 0)
888 {
889 invfr = 1.0/friction_coefficient;
890 for (n = 0; n < ngtc; n++)
891 {
892 rf[n] = sqrt(2.0*BOLTZ*ref_t[n]/(friction_coefficient*dt));
893 }
894 }
895 else
896 {
897 for (n = 0; n < ngtc; n++)
898 {
899 rf[n] = sqrt(2.0*BOLTZ*ref_t[n]);
900 }
901 }
902 for (n = start; (n < nrend); n++)
903 {
904 if (cFREEZE)
905 {
906 gf = cFREEZE[n];
907 }
908 if (cTC)
909 {
910 gt = cTC[n];
911 }
912 for (d = 0; (d < DIM); d++)
913 {
914 if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
915 {
916 if (friction_coefficient != 0)
917 {
918 vn = invfr*f[n][d] + rf[gt]*gmx_rng_gaussian_table(gaussrand);
919 }
920 else
921 {
922 /* NOTE: invmass = 2/(mass*friction_constant*dt) */
923 vn = 0.5*invmass[n]*f[n][d]*dt
924 + sqrt(0.5*invmass[n])*rf[gt]*gmx_rng_gaussian_table(gaussrand);
925 }
926
927 v[n][d] = vn;
928 xprime[n][d] = x[n][d]+vn*dt;
929 }
930 else
931 {
932 v[n][d] = 0.0;
933 xprime[n][d] = x[n][d];
934 }
935 }
936 }
937 }
938
939 static void dump_it_all(FILE *fp, const char *title,
940 int natoms, rvec x[], rvec xp[], rvec v[], rvec f[])
941 {
942 #ifdef DEBUG
943 if (fp)
944 {
945 fprintf(fp, "%s\n", title);
946 pr_rvecs(fp, 0, "x", x, natoms);
947 pr_rvecs(fp, 0, "xp", xp, natoms);
948 pr_rvecs(fp, 0, "v", v, natoms);
949 pr_rvecs(fp, 0, "f", f, natoms);
950 }
951 #endif
952 }
953
954 static void calc_ke_part_normal(rvec v[], t_grpopts *opts, t_mdatoms *md,
955 gmx_ekindata_t *ekind, t_nrnb *nrnb, gmx_bool bEkinAveVel,
956 gmx_bool bSaveEkinOld)
957 {
958 int g;
959 t_grp_tcstat *tcstat = ekind->tcstat;
960 t_grp_acc *grpstat = ekind->grpstat;
961 int nthread, thread;
962
963 /* three main: VV with AveVel, vv with AveEkin, leap with AveEkin. Leap with AveVel is also
964 an option, but not supported now. Additionally, if we are doing iterations.
965 bEkinAveVel: If TRUE, we sum into ekin, if FALSE, into ekinh.
966 bSavEkinOld: If TRUE (in the case of iteration = bIterate is TRUE), we don't copy over the ekinh_old.
967 If FALSE, we overrwrite it.
968 */
969
970 /* group velocities are calculated in update_ekindata and
971 * accumulated in acumulate_groups.
972 * Now the partial global and groups ekin.
973 */
974 for (g = 0; (g < opts->ngtc); g++)
975 {
976
977 if (!bSaveEkinOld)
978 {
979 copy_mat(tcstat[g].ekinh, tcstat[g].ekinh_old);
980 }
981 if (bEkinAveVel)
982 {
983 clear_mat(tcstat[g].ekinf);
984 }
985 else
986 {
987 clear_mat(tcstat[g].ekinh);
988 }
989 if (bEkinAveVel)
990 {
991 tcstat[g].ekinscalef_nhc = 1.0; /* need to clear this -- logic is complicated! */
992 }
993 }
994 ekind->dekindl_old = ekind->dekindl;
995
996 nthread = gmx_omp_nthreads_get(emntUpdate);
997
998 #pragma omp parallel for num_threads(nthread) schedule(static)
999 for (thread = 0; thread < nthread; thread++)
1000 {
1001 int start_t, end_t, n;
1002 int ga, gt;
1003 rvec v_corrt;
1004 real hm;
1005 int d, m;
1006 matrix *ekin_sum;
1007 real *dekindl_sum;
1008
1009 start_t = md->start + ((thread+0)*md->homenr)/nthread;
1010 end_t = md->start + ((thread+1)*md->homenr)/nthread;
1011
1012 ekin_sum = ekind->ekin_work[thread];
1013 dekindl_sum = ekind->dekindl_work[thread];
1014
1015 for (gt = 0; gt < opts->ngtc; gt++)
1016 {
1017 clear_mat(ekin_sum[gt]);
1018 }
1019 *dekindl_sum = 0.0;
1020
1021 ga = 0;
1022 gt = 0;
1023 for (n = start_t; n < end_t; n++)
1024 {
1025 if (md->cACC)
1026 {
1027 ga = md->cACC[n];
1028 }
1029 if (md->cTC)
1030 {
1031 gt = md->cTC[n];
1032 }
1033 hm = 0.5*md->massT[n];
1034
1035 for (d = 0; (d < DIM); d++)
1036 {
1037 v_corrt[d] = v[n][d] - grpstat[ga].u[d];
1038 }
1039 for (d = 0; (d < DIM); d++)
1040 {
1041 for (m = 0; (m < DIM); m++)
1042 {
1043 /* if we're computing a full step velocity, v_corrt[d] has v(t). Otherwise, v(t+dt/2) */
1044 ekin_sum[gt][m][d] += hm*v_corrt[m]*v_corrt[d];
1045 }
1046 }
1047 if (md->nMassPerturbed && md->bPerturbed[n])
1048 {
1049 *dekindl_sum +=
1050 0.5*(md->massB[n] - md->massA[n])*iprod(v_corrt, v_corrt);
1051 }
1052 }
1053 }
1054
1055 ekind->dekindl = 0;
1056 for (thread = 0; thread < nthread; thread++)
1057 {
1058 for (g = 0; g < opts->ngtc; g++)
1059 {
1060 if (bEkinAveVel)
1061 {
1062 m_add(tcstat[g].ekinf, ekind->ekin_work[thread][g],
1063 tcstat[g].ekinf);
1064 }
1065 else
1066 {
1067 m_add(tcstat[g].ekinh, ekind->ekin_work[thread][g],
1068 tcstat[g].ekinh);
1069 }
1070 }
1071
1072 ekind->dekindl += *ekind->dekindl_work[thread];
1073 }
1074
1075 inc_nrnb(nrnb, eNR_EKIN, md->homenr);
1076 }
1077
1078 static void calc_ke_part_visc(matrix box, rvec x[], rvec v[],
1079 t_grpopts *opts, t_mdatoms *md,
1080 gmx_ekindata_t *ekind,
1081 t_nrnb *nrnb, gmx_bool bEkinAveVel, gmx_bool bSaveEkinOld)
1082 {
1083 int start = md->start, homenr = md->homenr;
1084 int g, d, n, m, gt = 0;
1085 rvec v_corrt;
1086 real hm;
1087 t_grp_tcstat *tcstat = ekind->tcstat;
1088 t_cos_acc *cosacc = &(ekind->cosacc);
1089 real dekindl;
1090 real fac, cosz;
1091 double mvcos;
1092
1093 for (g = 0; g < opts->ngtc; g++)
1094 {
1095 copy_mat(ekind->tcstat[g].ekinh, ekind->tcstat[g].ekinh_old);
1096 clear_mat(ekind->tcstat[g].ekinh);
1097 }
1098 ekind->dekindl_old = ekind->dekindl;
1099
1100 fac = 2*M_PI/box[ZZ][ZZ];
1101 mvcos = 0;
1102 dekindl = 0;
1103 for (n = start; n < start+homenr; n++)
1104 {
1105 if (md->cTC)
1106 {
1107 gt = md->cTC[n];
1108 }
1109 hm = 0.5*md->massT[n];
1110
1111 /* Note that the times of x and v differ by half a step */
1112 /* MRS -- would have to be changed for VV */
1113 cosz = cos(fac*x[n][ZZ]);
1114 /* Calculate the amplitude of the new velocity profile */
1115 mvcos += 2*cosz*md->massT[n]*v[n][XX];
1116
1117 copy_rvec(v[n], v_corrt);
1118 /* Subtract the profile for the kinetic energy */
1119 v_corrt[XX] -= cosz*cosacc->vcos;
1120 for (d = 0; (d < DIM); d++)
1121 {
1122 for (m = 0; (m < DIM); m++)
1123 {
1124 /* if we're computing a full step velocity, v_corrt[d] has v(t). Otherwise, v(t+dt/2) */
1125 if (bEkinAveVel)
1126 {
1127 tcstat[gt].ekinf[m][d] += hm*v_corrt[m]*v_corrt[d];
1128 }
1129 else
1130 {
1131 tcstat[gt].ekinh[m][d] += hm*v_corrt[m]*v_corrt[d];
1132 }
1133 }
1134 }
1135 if (md->nPerturbed && md->bPerturbed[n])
1136 {
1137 dekindl += 0.5*(md->massB[n] - md->massA[n])*iprod(v_corrt, v_corrt);
1138 }
1139 }
1140 ekind->dekindl = dekindl;
1141 cosacc->mvcos = mvcos;
1142
1143 inc_nrnb(nrnb, eNR_EKIN, homenr);
1144 }
1145
1146 void calc_ke_part(t_state *state, t_grpopts *opts, t_mdatoms *md,
1147 gmx_ekindata_t *ekind, t_nrnb *nrnb, gmx_bool bEkinAveVel, gmx_bool bSaveEkinOld)
1148 {
1149 if (ekind->cosacc.cos_accel == 0)
1150 {
1151 calc_ke_part_normal(state->v, opts, md, ekind, nrnb, bEkinAveVel, bSaveEkinOld);
1152 }
1153 else
1154 {
1155 calc_ke_part_visc(state->box, state->x, state->v, opts, md, ekind, nrnb, bEkinAveVel, bSaveEkinOld);
1156 }
1157 }
1158
1159 extern void init_ekinstate(ekinstate_t *ekinstate, const t_inputrec *ir)
1160 {
1161 ekinstate->ekin_n = ir->opts.ngtc;
1162 snew(ekinstate->ekinh, ekinstate->ekin_n);
1163 snew(ekinstate->ekinf, ekinstate->ekin_n);
1164 snew(ekinstate->ekinh_old, ekinstate->ekin_n);
1165 snew(ekinstate->ekinscalef_nhc, ekinstate->ekin_n);
1166 snew(ekinstate->ekinscaleh_nhc, ekinstate->ekin_n);
1167 snew(ekinstate->vscale_nhc, ekinstate->ekin_n);
1168 ekinstate->dekindl = 0;
1169 ekinstate->mvcos = 0;
1170 }
1171
1172 void update_ekinstate(ekinstate_t *ekinstate, gmx_ekindata_t *ekind)
1173 {
1174 int i;
1175
1176 for (i = 0; i < ekinstate->ekin_n; i++)
1177 {
1178 copy_mat(ekind->tcstat[i].ekinh, ekinstate->ekinh[i]);
1179 copy_mat(ekind->tcstat[i].ekinf, ekinstate->ekinf[i]);
1180 copy_mat(ekind->tcstat[i].ekinh_old, ekinstate->ekinh_old[i]);
1181 ekinstate->ekinscalef_nhc[i] = ekind->tcstat[i].ekinscalef_nhc;
1182 ekinstate->ekinscaleh_nhc[i] = ekind->tcstat[i].ekinscaleh_nhc;
1183 ekinstate->vscale_nhc[i] = ekind->tcstat[i].vscale_nhc;
1184 }
1185
1186 copy_mat(ekind->ekin, ekinstate->ekin_total);
1187 ekinstate->dekindl = ekind->dekindl;
1188 ekinstate->mvcos = ekind->cosacc.mvcos;
1189
1190 }
1191
1192 void restore_ekinstate_from_state(t_commrec *cr,
1193 gmx_ekindata_t *ekind, ekinstate_t *ekinstate)
1194 {
1195 int i, n;
1196
1197 if (MASTER(cr))
1198 {
1199 for (i = 0; i < ekinstate->ekin_n; i++)
1200 {
1201 copy_mat(ekinstate->ekinh[i], ekind->tcstat[i].ekinh);
1202 copy_mat(ekinstate->ekinf[i], ekind->tcstat[i].ekinf);
1203 copy_mat(ekinstate->ekinh_old[i], ekind->tcstat[i].ekinh_old);
1204 ekind->tcstat[i].ekinscalef_nhc = ekinstate->ekinscalef_nhc[i];
1205 ekind->tcstat[i].ekinscaleh_nhc = ekinstate->ekinscaleh_nhc[i];
1206 ekind->tcstat[i].vscale_nhc = ekinstate->vscale_nhc[i];
1207 }
1208
1209 copy_mat(ekinstate->ekin_total, ekind->ekin);
1210
1211 ekind->dekindl = ekinstate->dekindl;
1212 ekind->cosacc.mvcos = ekinstate->mvcos;
1213 n = ekinstate->ekin_n;
1214 }
1215
1216 if (PAR(cr))
1217 {
1218 gmx_bcast(sizeof(n), &n, cr);
1219 for (i = 0; i < n; i++)
1220 {
1221 gmx_bcast(DIM*DIM*sizeof(ekind->tcstat[i].ekinh[0][0]),
1222 ekind->tcstat[i].ekinh[0], cr);
1223 gmx_bcast(DIM*DIM*sizeof(ekind->tcstat[i].ekinf[0][0]),
1224 ekind->tcstat[i].ekinf[0], cr);
1225 gmx_bcast(DIM*DIM*sizeof(ekind->tcstat[i].ekinh_old[0][0]),
1226 ekind->tcstat[i].ekinh_old[0], cr);
1227
1228 gmx_bcast(sizeof(ekind->tcstat[i].ekinscalef_nhc),
1229 &(ekind->tcstat[i].ekinscalef_nhc), cr);
1230 gmx_bcast(sizeof(ekind->tcstat[i].ekinscaleh_nhc),
1231 &(ekind->tcstat[i].ekinscaleh_nhc), cr);
1232 gmx_bcast(sizeof(ekind->tcstat[i].vscale_nhc),
1233 &(ekind->tcstat[i].vscale_nhc), cr);
1234 }
1235 gmx_bcast(DIM*DIM*sizeof(ekind->ekin[0][0]),
1236 ekind->ekin[0], cr);
1237
1238 gmx_bcast(sizeof(ekind->dekindl), &ekind->dekindl, cr);
1239 gmx_bcast(sizeof(ekind->cosacc.mvcos), &ekind->cosacc.mvcos, cr);
1240 }
1241 }
1242
1243 void set_deform_reference_box(gmx_update_t upd, gmx_large_int_t step, matrix box)
1244 {
1245 upd->deformref_step = step;
1246 copy_mat(box, upd->deformref_box);
1247 }
1248
1249 static void deform(gmx_update_t upd,
1250 int start, int homenr, rvec x[], matrix box, matrix *scale_tot,
1251 const t_inputrec *ir, gmx_large_int_t step)
1252 {
1253 matrix bnew, invbox, mu;
1254 real elapsed_time;
1255 int i, j;
1256
1257 elapsed_time = (step + 1 - upd->deformref_step)*ir->delta_t;
1258 copy_mat(box, bnew);
1259 for (i = 0; i < DIM; i++)
1260 {
1261 for (j = 0; j < DIM; j++)
1262 {
1263 if (ir->deform[i][j] != 0)
1264 {
1265 bnew[i][j] =
1266 upd->deformref_box[i][j] + elapsed_time*ir->deform[i][j];
1267 }
1268 }
1269 }
1270 /* We correct the off-diagonal elements,
1271 * which can grow indefinitely during shearing,
1272 * so the shifts do not get messed up.
1273 */
1274 for (i = 1; i < DIM; i++)
1275 {
1276 for (j = i-1; j >= 0; j--)
1277 {
1278 while (bnew[i][j] - box[i][j] > 0.5*bnew[j][j])
1279 {
1280 rvec_dec(bnew[i], bnew[j]);
1281 }
1282 while (bnew[i][j] - box[i][j] < -0.5*bnew[j][j])
1283 {
1284 rvec_inc(bnew[i], bnew[j]);
1285 }
1286 }
1287 }
1288 m_inv_ur0(box, invbox);
1289 copy_mat(bnew, box);
1290 mmul_ur0(box, invbox, mu);
1291
1292 for (i = start; i < start+homenr; i++)
1293 {
1294 x[i][XX] = mu[XX][XX]*x[i][XX]+mu[YY][XX]*x[i][YY]+mu[ZZ][XX]*x[i][ZZ];
1295 x[i][YY] = mu[YY][YY]*x[i][YY]+mu[ZZ][YY]*x[i][ZZ];
1296 x[i][ZZ] = mu[ZZ][ZZ]*x[i][ZZ];
1297 }
1298 if (*scale_tot)
1299 {
1300 /* The transposes of the scaling matrices are stored,
1301 * so we need to do matrix multiplication in the inverse order.
1302 */
1303 mmul_ur0(*scale_tot, mu, *scale_tot);
1304 }
1305 }
1306
1307 static void combine_forces(int nstcalclr,
1308 gmx_constr_t constr,
1309 t_inputrec *ir, t_mdatoms *md, t_idef *idef,
1310 t_commrec *cr,
1311 gmx_large_int_t step,
1312 t_state *state, gmx_bool bMolPBC,
1313 int start, int nrend,
1314 rvec f[], rvec f_lr[],
1315 t_nrnb *nrnb)
1316 {
1317 int i, d, nm1;
1318
1319 /* f contains the short-range forces + the long range forces
1320 * which are stored separately in f_lr.
1321 */
1322
1323 if (constr != NULL && !(ir->eConstrAlg == econtSHAKE && ir->epc == epcNO))
1324 {
1325 /* We need to constrain the LR forces separately,
1326 * because due to the different pre-factor for the SR and LR
1327 * forces in the update algorithm, we can not determine
1328 * the constraint force for the coordinate constraining.
1329 * Constrain only the additional LR part of the force.
1330 */
1331 /* MRS -- need to make sure this works with trotter integration -- the constraint calls may not be right.*/
1332 constrain(NULL, FALSE, FALSE, constr, idef, ir, NULL, cr, step, 0, md,
1333 state->x, f_lr, f_lr, bMolPBC, state->box, state->lambda[efptBONDED], NULL,
1334 NULL, NULL, nrnb, econqForce, ir->epc == epcMTTK, state->veta, state->veta);
1335 }
1336
1337 /* Add nstcalclr-1 times the LR force to the sum of both forces
1338 * and store the result in forces_lr.
1339 */
1340 nm1 = nstcalclr - 1;
1341 for (i = start; i < nrend; i++)
1342 {
1343 for (d = 0; d < DIM; d++)
1344 {
1345 f_lr[i][d] = f[i][d] + nm1*f_lr[i][d];
1346 }
1347 }
1348 }
1349
1350 void update_tcouple(FILE *fplog,
1351 gmx_large_int_t step,
1352 t_inputrec *inputrec,
1353 t_state *state,
1354 gmx_ekindata_t *ekind,
1355 gmx_wallcycle_t wcycle,
1356 gmx_update_t upd,
1357 t_extmass *MassQ,
1358 t_mdatoms *md)
1359
1360 {
1361 gmx_bool bTCouple = FALSE;
1362 real dttc;
1363 int i, start, end, homenr, offset;
1364
1365 /* if using vv with trotter decomposition methods, we do this elsewhere in the code */
1366 if (inputrec->etc != etcNO &&
1367 !(IR_NVT_TROTTER(inputrec) || IR_NPT_TROTTER(inputrec) || IR_NPH_TROTTER(inputrec)))
1368 {
1369 /* We should only couple after a step where energies were determined (for leapfrog versions)
1370 or the step energies are determined, for velocity verlet versions */
1371
1372 if (EI_VV(inputrec->eI))
1373 {
1374 offset = 0;
1375 }
1376 else
1377 {
1378 offset = 1;
1379 }
1380 bTCouple = (inputrec->nsttcouple == 1 ||
1381 do_per_step(step+inputrec->nsttcouple-offset,
1382 inputrec->nsttcouple));
1383 }
1384
1385 if (bTCouple)
1386 {
1387 dttc = inputrec->nsttcouple*inputrec->delta_t;
1388
1389 switch (inputrec->etc)
1390 {
1391 case etcNO:
1392 break;
1393 case etcBERENDSEN:
1394 berendsen_tcoupl(inputrec, ekind, dttc);
1395 break;
1396 case etcNOSEHOOVER:
1397 nosehoover_tcoupl(&(inputrec->opts), ekind, dttc,
1398 state->nosehoover_xi, state->nosehoover_vxi, MassQ);
1399 break;
1400 case etcVRESCALE:
1401 vrescale_tcoupl(inputrec, ekind, dttc,
1402 state->therm_integral, upd->sd->gaussrand[0]);
1403 break;
1404 }
1405 /* rescale in place here */
1406 if (EI_VV(inputrec->eI))
1407 {
1408 rescale_velocities(ekind, md, md->start, md->start+md->homenr, state->v);
1409 }
1410 }
1411 else
1412 {
1413 /* Set the T scaling lambda to 1 to have no scaling */
1414 for (i = 0; (i < inputrec->opts.ngtc); i++)
1415 {
1416 ekind->tcstat[i].lambda = 1.0;
1417 }
1418 }
1419 }
1420
1421 void update_pcouple(FILE *fplog,
1422 gmx_large_int_t step,
1423 t_inputrec *inputrec,
1424 t_state *state,
1425 matrix pcoupl_mu,
1426 matrix M,
1427 gmx_wallcycle_t wcycle,
1428 gmx_update_t upd,
1429 gmx_bool bInitStep)
1430 {
1431 gmx_bool bPCouple = FALSE;
1432 real dtpc = 0;
1433 int i;
1434
1435 /* if using Trotter pressure, we do this in coupling.c, so we leave it false. */
1436 if (inputrec->epc != epcNO && (!(IR_NPT_TROTTER(inputrec) || IR_NPH_TROTTER(inputrec))))
1437 {
1438 /* We should only couple after a step where energies were determined */
1439 bPCouple = (inputrec->nstpcouple == 1 ||
1440 do_per_step(step+inputrec->nstpcouple-1,
1441 inputrec->nstpcouple));
1442 }
1443
1444 clear_mat(pcoupl_mu);
1445 for (i = 0; i < DIM; i++)
1446 {
1447 pcoupl_mu[i][i] = 1.0;
1448 }
1449
1450 clear_mat(M);
1451
1452 if (bPCouple)
1453 {
1454 dtpc = inputrec->nstpcouple*inputrec->delta_t;
1455
1456 switch (inputrec->epc)
1457 {
1458 /* We can always pcoupl, even if we did not sum the energies
1459 * the previous step, since state->pres_prev is only updated
1460 * when the energies have been summed.
1461 */
1462 case (epcNO):
1463 break;
1464 case (epcBERENDSEN):
1465 if (!bInitStep)
1466 {
1467 berendsen_pcoupl(fplog, step, inputrec, dtpc, state->pres_prev, state->box,
1468 pcoupl_mu);
1469 }
1470 break;
1471 case (epcPARRINELLORAHMAN):
1472 parrinellorahman_pcoupl(fplog, step, inputrec, dtpc, state->pres_prev,
1473 state->box, state->box_rel, state->boxv,
1474 M, pcoupl_mu, bInitStep);
1475 break;
1476 default:
1477 break;
1478 }
1479 }
1480 }
1481
1482 static rvec *get_xprime(const t_state *state, gmx_update_t upd)
1483 {
1484 if (state->nalloc > upd->xp_nalloc)
1485 {
1486 upd->xp_nalloc = state->nalloc;
1487 srenew(upd->xp, upd->xp_nalloc);
1488 }
1489
1490 return upd->xp;
1491 }
1492
1493 void update_constraints(FILE *fplog,
1494 gmx_large_int_t step,
1495 real *dvdlambda, /* the contribution to be added to the bonded interactions */
1496 t_inputrec *inputrec, /* input record and box stuff */
1497 gmx_ekindata_t *ekind,
1498 t_mdatoms *md,
1499 t_state *state,
1500 gmx_bool bMolPBC,
1501 t_graph *graph,
1502 rvec force[], /* forces on home particles */
1503 t_idef *idef,
1504 tensor vir_part,
1505 tensor vir, /* tensors for virial and ekin, needed for computing */
1506 t_commrec *cr,
1507 t_nrnb *nrnb,
1508 gmx_wallcycle_t wcycle,
1509 gmx_update_t upd,
1510 gmx_constr_t constr,
1511 gmx_bool bInitStep,
1512 gmx_bool bFirstHalf,
1513 gmx_bool bCalcVir,
1514 real vetanew)
1515 {
1516 gmx_bool bExtended, bLastStep, bLog = FALSE, bEner = FALSE, bDoConstr = FALSE;
1517 double dt;
1518 real dt_1;
1519 int start, homenr, nrend, i, n, m, g, d;
1520 tensor vir_con;
1521 rvec *vbuf, *xprime = NULL;
1522 int nth, th;
1523
1524 if (constr)
1525 {
1526 bDoConstr = TRUE;
1527 }
1528 if (bFirstHalf && !EI_VV(inputrec->eI))
1529 {
1530 bDoConstr = FALSE;
1531 }
1532
1533 /* for now, SD update is here -- though it really seems like it
1534 should be reformulated as a velocity verlet method, since it has two parts */
1535
1536 start = md->start;
1537 homenr = md->homenr;
1538 nrend = start+homenr;
1539
1540 dt = inputrec->delta_t;
1541 dt_1 = 1.0/dt;
1542
1543 /*
1544 * Steps (7C, 8C)
1545 * APPLY CONSTRAINTS:
1546 * BLOCK SHAKE
1547
1548 * When doing PR pressure coupling we have to constrain the
1549 * bonds in each iteration. If we are only using Nose-Hoover tcoupling
1550 * it is enough to do this once though, since the relative velocities
1551 * after this will be normal to the bond vector
1552 */
1553
1554 if (bDoConstr)
1555 {
1556 /* clear out constraints before applying */
1557 clear_mat(vir_part);
1558
1559 xprime = get_xprime(state, upd);
1560
1561 bLastStep = (step == inputrec->init_step+inputrec->nsteps);
1562 bLog = (do_per_step(step, inputrec->nstlog) || bLastStep || (step < 0));
1563 bEner = (do_per_step(step, inputrec->nstenergy) || bLastStep);
1564 /* Constrain the coordinates xprime */
1565 wallcycle_start(wcycle, ewcCONSTR);
1566 if (EI_VV(inputrec->eI) && bFirstHalf)
1567 {
1568 constrain(NULL, bLog, bEner, constr, idef,
1569 inputrec, ekind, cr, step, 1, md,
1570 state->x, state->v, state->v,
1571 bMolPBC, state->box,
1572 state->lambda[efptBONDED], dvdlambda,
1573 NULL, bCalcVir ? &vir_con : NULL, nrnb, econqVeloc,
1574 inputrec->epc == epcMTTK, state->veta, vetanew);
1575 }
1576 else
1577 {
1578 constrain(NULL, bLog, bEner, constr, idef,
1579 inputrec, ekind, cr, step, 1, md,
1580 state->x, xprime, NULL,
1581 bMolPBC, state->box,
1582 state->lambda[efptBONDED], dvdlambda,
1583 state->v, bCalcVir ? &vir_con : NULL, nrnb, econqCoord,
1584 inputrec->epc == epcMTTK, state->veta, state->veta);
1585 }
1586 wallcycle_stop(wcycle, ewcCONSTR);
1587
1588 where();
1589
1590 dump_it_all(fplog, "After Shake",
1591 state->natoms, state->x, xprime, state->v, force);
1592
1593 if (bCalcVir)
1594 {
1595 if (inputrec->eI == eiSD2)
1596 {
1597 /* A correction factor eph is needed for the SD constraint force */
1598 /* Here we can, unfortunately, not have proper corrections
1599 * for different friction constants, so we use the first one.
1600 */
1601 for (i = 0; i < DIM; i++)
1602 {
1603 for (m = 0; m < DIM; m++)
1604 {
1605 vir_part[i][m] += upd->sd->sdc[0].eph*vir_con[i][m];
1606 }
1607 }
1608 }
1609 else
1610 {
1611 m_add(vir_part, vir_con, vir_part);
1612 }
1613 if (debug)
1614 {
1615 pr_rvecs(debug, 0, "constraint virial", vir_part, DIM);
1616 }
1617 }
1618 }
1619
1620 where();
1621 if ((inputrec->eI == eiSD2) && !(bFirstHalf))
1622 {
1623 xprime = get_xprime(state, upd);
1624
1625 nth = gmx_omp_nthreads_get(emntUpdate);
1626
1627 #pragma omp parallel for num_threads(nth) schedule(static)
1628 for (th = 0; th < nth; th++)
1629 {
1630 int start_th, end_th;
1631
1632 start_th = start + ((nrend-start)* th )/nth;
1633 end_th = start + ((nrend-start)*(th+1))/nth;
1634
1635 /* The second part of the SD integration */
1636 do_update_sd2(upd->sd, upd->sd->gaussrand[th],
1637 FALSE, start_th, end_th,
1638 inputrec->opts.acc, inputrec->opts.nFreeze,
1639 md->invmass, md->ptype,
1640 md->cFREEZE, md->cACC, md->cTC,
1641 state->x, xprime, state->v, force, state->sd_X,
1642 inputrec->opts.ngtc, inputrec->opts.tau_t,
1643 inputrec->opts.ref_t, FALSE);
1644 }
1645 inc_nrnb(nrnb, eNR_UPDATE, homenr);
1646
1647 if (bDoConstr)
1648 {
1649 /* Constrain the coordinates xprime */
1650 wallcycle_start(wcycle, ewcCONSTR);
1651 constrain(NULL, bLog, bEner, constr, idef,
1652 inputrec, NULL, cr, step, 1, md,
1653 state->x, xprime, NULL,
1654 bMolPBC, state->box,
1655 state->lambda[efptBONDED], dvdlambda,
1656 NULL, NULL, nrnb, econqCoord, FALSE, 0, 0);
1657 wallcycle_stop(wcycle, ewcCONSTR);
1658 }
1659 }
1660
1661 /* We must always unshift after updating coordinates; if we did not shake
1662 x was shifted in do_force */
1663
1664 if (!(bFirstHalf)) /* in the first half of vv, no shift. */
1665 {
1666 if (graph && (graph->nnodes > 0))
1667 {
1668 unshift_x(graph, state->box, state->x, upd->xp);
1669 if (TRICLINIC(state->box))
1670 {
1671 inc_nrnb(nrnb, eNR_SHIFTX, 2*graph->nnodes);
1672 }
1673 else
1674 {
1675 inc_nrnb(nrnb, eNR_SHIFTX, graph->nnodes);
1676 }
1677 }
1678 else
1679 {
1680 #pragma omp parallel for num_threads(gmx_omp_nthreads_get(emntUpdate)) schedule(static)
1681 for (i = start; i < nrend; i++)
1682 {
1683 copy_rvec(upd->xp[i], state->x[i]);
1684 }
1685 }
1686
1687 dump_it_all(fplog, "After unshift",
1688 state->natoms, state->x, upd->xp, state->v, force);
1689 }
1690 /* ############# END the update of velocities and positions ######### */
1691 }
1692
1693 void update_box(FILE *fplog,
1694 gmx_large_int_t step,
1695 t_inputrec *inputrec, /* input record and box stuff */
1696 t_mdatoms *md,
1697 t_state *state,
1698 t_graph *graph,
1699 rvec force[], /* forces on home particles */
1700 matrix *scale_tot,
1701 matrix pcoupl_mu,
1702 t_nrnb *nrnb,
1703 gmx_wallcycle_t wcycle,
1704 gmx_update_t upd,
1705 gmx_bool bInitStep,
1706 gmx_bool bFirstHalf)
1707 {
1708 gmx_bool bExtended, bLastStep, bLog = FALSE, bEner = FALSE;
1709 double dt;
1710 real dt_1;
1711 int start, homenr, nrend, i, n, m, g;
1712 tensor vir_con;
1713
1714 start = md->start;
1715 homenr = md->homenr;
1716 nrend = start+homenr;
1717
1718 bExtended =
1719 (inputrec->etc == etcNOSEHOOVER) ||
1720 (inputrec->epc == epcPARRINELLORAHMAN) ||
1721 (inputrec->epc == epcMTTK);
1722
1723 dt = inputrec->delta_t;
1724
1725 where();
1726
1727 /* now update boxes */
1728 switch (inputrec->epc)
1729 {
1730 case (epcNO):
1731 break;
1732 case (epcBERENDSEN):
1733 berendsen_pscale(inputrec, pcoupl_mu, state->box, state->box_rel,
1734 start, homenr, state->x, md->cFREEZE, nrnb);
1735 break;
1736 case (epcPARRINELLORAHMAN):
1737 /* The box velocities were updated in do_pr_pcoupl in the update
1738 * iteration, but we dont change the box vectors until we get here
1739 * since we need to be able to shift/unshift above.
1740 */
1741 for (i = 0; i < DIM; i++)
1742 {
1743 for (m = 0; m <= i; m++)
1744 {
1745 state->box[i][m] += dt*state->boxv[i][m];
1746 }
1747 }
1748 preserve_box_shape(inputrec, state->box_rel, state->box);
1749
1750 /* Scale the coordinates */
1751 for (n = start; (n < start+homenr); n++)
1752 {
1753 tmvmul_ur0(pcoupl_mu, state->x[n], state->x[n]);
1754 }
1755 break;
1756 case (epcMTTK):
1757 switch (inputrec->epct)
1758 {
1759 case (epctISOTROPIC):
1760 /* DIM * eta = ln V. so DIM*eta_new = DIM*eta_old + DIM*dt*veta =>
1761 ln V_new = ln V_old + 3*dt*veta => V_new = V_old*exp(3*dt*veta) =>
1762 Side length scales as exp(veta*dt) */
1763
1764 msmul(state->box, exp(state->veta*dt), state->box);
1765
1766 /* Relate veta to boxv. veta = d(eta)/dT = (1/DIM)*1/V dV/dT.
1767 o If we assume isotropic scaling, and box length scaling
1768 factor L, then V = L^DIM (det(M)). So dV/dt = DIM
1769 L^(DIM-1) dL/dt det(M), and veta = (1/L) dL/dt. The
1770 determinant of B is L^DIM det(M), and the determinant
1771 of dB/dt is (dL/dT)^DIM det (M). veta will be
1772 (det(dB/dT)/det(B))^(1/3). Then since M =
1773 B_new*(vol_new)^(1/3), dB/dT_new = (veta_new)*B(new). */
1774
1775 msmul(state->box, state->veta, state->boxv);
1776 break;
1777 default:
1778 break;
1779 }
1780 break;
1781 default:
1782 break;
1783 }
1784
1785 if ((!(IR_NPT_TROTTER(inputrec) || IR_NPH_TROTTER(inputrec))) && scale_tot)
1786 {
1787 /* The transposes of the scaling matrices are stored,
1788 * therefore we need to reverse the order in the multiplication.
1789 */
1790 mmul_ur0(*scale_tot, pcoupl_mu, *scale_tot);
1791 }
1792
1793 if (DEFORM(*inputrec))
1794 {
1795 deform(upd, start, homenr, state->x, state->box, scale_tot, inputrec, step);
1796 }
1797 where();
1798 dump_it_all(fplog, "After update",
1799 state->natoms, state->x, upd->xp, state->v, force);
1800 }
1801
1802 void update_coords(FILE *fplog,
1803 gmx_large_int_t step,
1804 t_inputrec *inputrec, /* input record and box stuff */
1805 t_mdatoms *md,
1806 t_state *state,
1807 gmx_bool bMolPBC,
1808 rvec *f, /* forces on home particles */
1809 gmx_bool bDoLR,
1810 rvec *f_lr,
1811 t_fcdata *fcd,
1812 gmx_ekindata_t *ekind,
1813 matrix M,
1814 gmx_wallcycle_t wcycle,
1815 gmx_update_t upd,
1816 gmx_bool bInitStep,
1817 int UpdatePart,
1818 t_commrec *cr, /* these shouldn't be here -- need to think about it */
1819 t_nrnb *nrnb,
1820 gmx_constr_t constr,
1821 t_idef *idef)
1822 {
1823 gmx_bool bNH, bPR, bLastStep, bLog = FALSE, bEner = FALSE;
1824 double dt, alpha;
1825 real *imass, *imassin;
1826 rvec *force;
1827 real dt_1;
1828 int start, homenr, nrend, i, j, d, n, m, g;
1829 int blen0, blen1, iatom, jatom, nshake, nsettle, nconstr, nexpand;
1830 int *icom = NULL;
1831 tensor vir_con;
1832 rvec *vcom, *xcom, *vall, *xall, *xin, *vin, *forcein, *fall, *xpall, *xprimein, *xprime;
1833 int nth, th;
1834
1835 /* Running the velocity half does nothing except for velocity verlet */
1836 if ((UpdatePart == etrtVELOCITY1 || UpdatePart == etrtVELOCITY2) &&
1837 !EI_VV(inputrec->eI))
1838 {
1839 gmx_incons("update_coords called for velocity without VV integrator");
1840 }
1841
1842 start = md->start;
1843 homenr = md->homenr;
1844 nrend = start+homenr;
1845
1846 xprime = get_xprime(state, upd);
1847
1848 dt = inputrec->delta_t;
1849 dt_1 = 1.0/dt;
1850
1851 /* We need to update the NMR restraint history when time averaging is used */
1852 if (state->flags & (1<<estDISRE_RM3TAV))
1853 {
1854 update_disres_history(fcd, &state->hist);
1855 }
1856 if (state->flags & (1<<estORIRE_DTAV))
1857 {
1858 update_orires_history(fcd, &state->hist);
1859 }
1860
1861
1862 bNH = inputrec->etc == etcNOSEHOOVER;
1863 bPR = ((inputrec->epc == epcPARRINELLORAHMAN) || (inputrec->epc == epcMTTK));
1864
1865 if (bDoLR && inputrec->nstcalclr > 1 && !EI_VV(inputrec->eI)) /* get this working with VV? */
1866 {
1867 /* Store the total force + nstcalclr-1 times the LR force
1868 * in forces_lr, so it can be used in a normal update algorithm
1869 * to produce twin time stepping.
1870 */
1871 /* is this correct in the new construction? MRS */
1872 combine_forces(inputrec->nstcalclr, constr, inputrec, md, idef, cr,
1873 step, state, bMolPBC,
1874 start, nrend, f, f_lr, nrnb);
1875 force = f_lr;
1876 }
1877 else
1878 {
1879 force = f;
1880 }
1881
1882 /* ############# START The update of velocities and positions ######### */
1883 where();
1884 dump_it_all(fplog, "Before update",
1885 state->natoms, state->x, xprime, state->v, force);
1886
1887 if (EI_RANDOM(inputrec->eI))
1888 {
1889 /* We still need to take care of generating random seeds properly
1890 * when multi-threading.
1891 */
1892 nth = 1;
1893 }
1894 else
1895 {
1896 nth = gmx_omp_nthreads_get(emntUpdate);
1897 }
1898
1899 if (inputrec->eI == eiSD2)
1900 {
1901 check_sd2_work_data_allocation(upd->sd, nrend);
1902 }
1903
1904 #pragma omp parallel for num_threads(nth) schedule(static) private(alpha)
1905 for (th = 0; th < nth; th++)
1906 {
1907 int start_th, end_th;
1908
1909 start_th = start + ((nrend-start)* th )/nth;
1910 end_th = start + ((nrend-start)*(th+1))/nth;
1911
1912 switch (inputrec->eI)
1913 {
1914 case (eiMD):
1915 if (ekind->cosacc.cos_accel == 0)
1916 {
1917 do_update_md(start_th, end_th, dt,
1918 ekind->tcstat, state->nosehoover_vxi,
1919 ekind->bNEMD, ekind->grpstat, inputrec->opts.acc,
1920 inputrec->opts.nFreeze,
1921 md->invmass, md->ptype,
1922 md->cFREEZE, md->cACC, md->cTC,
1923 state->x, xprime, state->v, force, M,
1924 bNH, bPR);
1925 }
1926 else
1927 {
1928 do_update_visc(start_th, end_th, dt,
1929 ekind->tcstat, state->nosehoover_vxi,
1930 md->invmass, md->ptype,
1931 md->cTC, state->x, xprime, state->v, force, M,
1932 state->box,
1933 ekind->cosacc.cos_accel,
1934 ekind->cosacc.vcos,
1935 bNH, bPR);
1936 }
1937 break;
1938 case (eiSD1):
1939 do_update_sd1(upd->sd, upd->sd->gaussrand[th],
1940 start_th, end_th, dt,
1941 inputrec->opts.acc, inputrec->opts.nFreeze,
1942 md->invmass, md->ptype,
1943 md->cFREEZE, md->cACC, md->cTC,
1944 state->x, xprime, state->v, force, state->sd_X,
1945 inputrec->opts.ngtc, inputrec->opts.tau_t, inputrec->opts.ref_t);
1946 break;
1947 case (eiSD2):
1948 /* The SD update is done in 2 parts, because an extra constraint step
1949 * is needed
1950 */
1951 do_update_sd2(upd->sd, upd->sd->gaussrand[th],
1952 bInitStep, start_th, end_th,
1953 inputrec->opts.acc, inputrec->opts.nFreeze,
1954 md->invmass, md->ptype,
1955 md->cFREEZE, md->cACC, md->cTC,
1956 state->x, xprime, state->v, force, state->sd_X,
1957 inputrec->opts.ngtc, inputrec->opts.tau_t, inputrec->opts.ref_t,
1958 TRUE);
1959 break;
1960 case (eiBD):
1961 do_update_bd(start_th, end_th, dt,
1962 inputrec->opts.nFreeze, md->invmass, md->ptype,
1963 md->cFREEZE, md->cTC,
1964 state->x, xprime, state->v, force,
1965 inputrec->bd_fric,
1966 inputrec->opts.ngtc, inputrec->opts.tau_t, inputrec->opts.ref_t,
1967 upd->sd->bd_rf, upd->sd->gaussrand[th]);
1968 break;
1969 case (eiVV):
1970 case (eiVVAK):
1971 alpha = 1.0 + DIM/((double)inputrec->opts.nrdf[0]); /* assuming barostat coupled to group 0. */
1972 switch (UpdatePart)
1973 {
1974 case etrtVELOCITY1:
1975 case etrtVELOCITY2:
1976 do_update_vv_vel(start_th, end_th, dt,
1977 ekind->tcstat, ekind->grpstat,
1978 inputrec->opts.acc, inputrec->opts.nFreeze,
1979 md->invmass, md->ptype,
1980 md->cFREEZE, md->cACC,
1981 state->v, force,
1982 (bNH || bPR), state->veta, alpha);
1983 break;
1984 case etrtPOSITION:
1985 do_update_vv_pos(start_th, end_th, dt,
1986 ekind->tcstat, ekind->grpstat,
1987 inputrec->opts.acc, inputrec->opts.nFreeze,
1988 md->invmass, md->ptype, md->cFREEZE,
1989 state->x, xprime, state->v, force,
1990 (bNH || bPR), state->veta, alpha);
1991 break;
1992 }
1993 break;
1994 default:
1995 gmx_fatal(FARGS, "Don't know how to update coordinates");
1996 break;
1997 }
1998 }
1999
2000 }
2001
2002
2003 void correct_ekin(FILE *log, int start, int end, rvec v[], rvec vcm, real mass[],
2004 real tmass, tensor ekin)
2005 {
2006 /*
2007 * This is a debugging routine. It should not be called for production code
2008 *
2009 * The kinetic energy should calculated according to:
2010 * Ekin = 1/2 m (v-vcm)^2
2011 * However the correction is not always applied, since vcm may not be
2012 * known in time and we compute
2013 * Ekin' = 1/2 m v^2 instead
2014 * This can be corrected afterwards by computing
2015 * Ekin = Ekin' + 1/2 m ( -2 v vcm + vcm^2)
2016 * or in hsorthand:
2017 * Ekin = Ekin' - m v vcm + 1/2 m vcm^2
2018 */
2019 int i, j, k;
2020 real m, tm;
2021 rvec hvcm, mv;
2022 tensor dekin;
2023
2024 /* Local particles */
2025 clear_rvec(mv);
2026
2027 /* Processor dependent part. */
2028 tm = 0;
2029 for (i = start; (i < end); i++)
2030 {
2031 m = mass[i];
2032 tm += m;
2033 for (j = 0; (j < DIM); j++)
2034 {
2035 mv[j] += m*v[i][j];
2036 }
2037 }
2038 /* Shortcut */
2039 svmul(1/tmass, vcm, vcm);
2040 svmul(0.5, vcm, hvcm);
2041 clear_mat(dekin);
2042 for (j = 0; (j < DIM); j++)
2043 {
2044 for (k = 0; (k < DIM); k++)
2045 {
2046 dekin[j][k] += vcm[k]*(tm*hvcm[j]-mv[j]);
2047 }
2048 }
2049 pr_rvecs(log, 0, "dekin", dekin, DIM);
2050 pr_rvecs(log, 0, " ekin", ekin, DIM);
2051 fprintf(log, "dekin = %g, ekin = %g vcm = (%8.4f %8.4f %8.4f)\n",
2052 trace(dekin), trace(ekin), vcm[XX], vcm[YY], vcm[ZZ]);
2053 fprintf(log, "mv = (%8.4f %8.4f %8.4f)\n",
2054 mv[XX], mv[YY], mv[ZZ]);
2055 }
2056
2057 extern gmx_bool update_randomize_velocities(t_inputrec *ir, gmx_large_int_t step, t_mdatoms *md, t_state *state, gmx_update_t upd, t_idef *idef, gmx_constr_t constr)
2058 {
2059
2060 int i;
2061 real rate = (ir->delta_t)/ir->opts.tau_t[0];
2062 /* proceed with andersen if 1) it's fixed probability per
2063 particle andersen or 2) it's massive andersen and it's tau_t/dt */
2064 if ((ir->etc == etcANDERSEN) || do_per_step(step, (int)(1.0/rate)))
2065 {
2066 srenew(upd->randatom, state->nalloc);
2067 srenew(upd->randatom_list, state->nalloc);
2068 if (upd->randatom_list_init == FALSE)
2069 {
2070 for (i = 0; i < state->nalloc; i++)
2071 {
2072 upd->randatom[i] = FALSE;
2073 upd->randatom_list[i] = 0;
2074 }
2075 upd->randatom_list_init = TRUE;
2076 }
2077 andersen_tcoupl(ir, md, state, upd->sd->gaussrand[0], rate,
2078 (ir->etc == etcANDERSEN) ? idef : NULL,
2079 constr ? get_nblocks(constr) : 0,
2080 constr ? get_sblock(constr) : NULL,
2081 upd->randatom, upd->randatom_list,
2082 upd->sd->randomize_group, upd->sd->boltzfac);
2083 return TRUE;
2084 }
2085 return FALSE;
2086 }
The ultimate molecular dynamics simulation package