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 * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source 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 #include "gmxpre.h"
39
40 #include "coupling.h"
41
42 #include <cassert>
43 #include <cmath>
44
45 #include <algorithm>
46
47 #include "gromacs/domdec/domdec_struct.h"
48 #include "gromacs/gmxlib/nrnb.h"
49 #include "gromacs/math/functions.h"
50 #include "gromacs/math/invertmatrix.h"
51 #include "gromacs/math/units.h"
52 #include "gromacs/math/vec.h"
53 #include "gromacs/math/vecdump.h"
54 #include "gromacs/mdlib/boxdeformation.h"
55 #include "gromacs/mdlib/expanded.h"
56 #include "gromacs/mdlib/gmx_omp_nthreads.h"
57 #include "gromacs/mdlib/stat.h"
58 #include "gromacs/mdlib/update.h"
59 #include "gromacs/mdtypes/commrec.h"
60 #include "gromacs/mdtypes/enerdata.h"
61 #include "gromacs/mdtypes/group.h"
62 #include "gromacs/mdtypes/inputrec.h"
63 #include "gromacs/mdtypes/md_enums.h"
64 #include "gromacs/mdtypes/mdatom.h"
65 #include "gromacs/mdtypes/state.h"
66 #include "gromacs/pbcutil/boxutilities.h"
67 #include "gromacs/pbcutil/pbc.h"
68 #include "gromacs/random/gammadistribution.h"
69 #include "gromacs/random/normaldistribution.h"
70 #include "gromacs/random/tabulatednormaldistribution.h"
71 #include "gromacs/random/threefry.h"
72 #include "gromacs/random/uniformrealdistribution.h"
73 #include "gromacs/utility/cstringutil.h"
74 #include "gromacs/utility/fatalerror.h"
75 #include "gromacs/utility/pleasecite.h"
76 #include "gromacs/utility/smalloc.h"
77
78 #define NTROTTERPARTS 3
79
80 /* Suzuki-Yoshida Constants, for n=3 and n=5, for symplectic integration */
81 /* for n=1, w0 = 1 */
82 /* for n=3, w0 = w2 = 1/(2-2^-(1/3)), w1 = 1-2*w0 */
83 /* for n=5, w0 = w1 = w3 = w4 = 1/(4-4^-(1/3)), w1 = 1-4*w0 */
84
85 #define MAX_SUZUKI_YOSHIDA_NUM 5
86 #define SUZUKI_YOSHIDA_NUM 5
87
88 static const double sy_const_1[] = { 1. };
89 static const double sy_const_3[] = { 0.828981543588751, -0.657963087177502, 0.828981543588751 };
90 static const double sy_const_5[] = { 0.2967324292201065, 0.2967324292201065, -0.186929716880426,
91 0.2967324292201065, 0.2967324292201065 };
92
93 static const double* sy_const[] = { nullptr, sy_const_1, nullptr, sy_const_3, nullptr, sy_const_5 };
94
95
update_tcouple(int64_t step,const t_inputrec * inputrec,t_state * state,gmx_ekindata_t * ekind,const t_extmass * MassQ,const t_mdatoms * md)96 void update_tcouple(int64_t step,
97 const t_inputrec* inputrec,
98 t_state* state,
99 gmx_ekindata_t* ekind,
100 const t_extmass* MassQ,
101 const t_mdatoms* md)
102
103 {
104 // This condition was explicitly checked in previous version, but should have never been satisfied
105 GMX_ASSERT(!(EI_VV(inputrec->eI)
106 && (inputrecNvtTrotter(inputrec) || inputrecNptTrotter(inputrec)
107 || inputrecNphTrotter(inputrec))),
108 "Temperature coupling was requested with velocity verlet and trotter");
109
110 bool doTemperatureCoupling = false;
111
112 // For VV temperature coupling parameters are updated on the current
113 // step, for the others - one step before.
114 if (inputrec->etc == etcNO)
115 {
116 doTemperatureCoupling = false;
117 }
118 else if (EI_VV(inputrec->eI))
119 {
120 doTemperatureCoupling = do_per_step(step, inputrec->nsttcouple);
121 }
122 else
123 {
124 doTemperatureCoupling = do_per_step(step + inputrec->nsttcouple - 1, inputrec->nsttcouple);
125 }
126
127 if (doTemperatureCoupling)
128 {
129 real dttc = inputrec->nsttcouple * inputrec->delta_t;
130
131 // TODO: berendsen_tcoupl(...), nosehoover_tcoupl(...) and vrescale_tcoupl(...) update
132 // temperature coupling parameters, which should be reflected in the name of these
133 // subroutines
134 switch (inputrec->etc)
135 {
136 case etcNO: break;
137 case etcBERENDSEN:
138 berendsen_tcoupl(inputrec, ekind, dttc, state->therm_integral);
139 break;
140 case etcNOSEHOOVER:
141 nosehoover_tcoupl(&(inputrec->opts), ekind, dttc, state->nosehoover_xi.data(),
142 state->nosehoover_vxi.data(), MassQ);
143 break;
144 case etcVRESCALE:
145 vrescale_tcoupl(inputrec, step, ekind, dttc, state->therm_integral.data());
146 break;
147 }
148 /* rescale in place here */
149 if (EI_VV(inputrec->eI))
150 {
151 rescale_velocities(ekind, md, 0, md->homenr, state->v.rvec_array());
152 }
153 }
154 else
155 {
156 // Set the T scaling lambda to 1 to have no scaling
157 // TODO: Do we have to do it on every non-t-couple step?
158 for (int i = 0; (i < inputrec->opts.ngtc); i++)
159 {
160 ekind->tcstat[i].lambda = 1.0;
161 }
162 }
163 }
164
update_pcouple_before_coordinates(FILE * fplog,int64_t step,const t_inputrec * inputrec,t_state * state,matrix parrinellorahmanMu,matrix M,gmx_bool bInitStep)165 void update_pcouple_before_coordinates(FILE* fplog,
166 int64_t step,
167 const t_inputrec* inputrec,
168 t_state* state,
169 matrix parrinellorahmanMu,
170 matrix M,
171 gmx_bool bInitStep)
172 {
173 /* Berendsen P-coupling is completely handled after the coordinate update.
174 * Trotter P-coupling is handled by separate calls to trotter_update().
175 */
176 if (inputrec->epc == epcPARRINELLORAHMAN
177 && do_per_step(step + inputrec->nstpcouple - 1, inputrec->nstpcouple))
178 {
179 real dtpc = inputrec->nstpcouple * inputrec->delta_t;
180
181 parrinellorahman_pcoupl(fplog, step, inputrec, dtpc, state->pres_prev, state->box,
182 state->box_rel, state->boxv, M, parrinellorahmanMu, bInitStep);
183 }
184 }
185
update_pcouple_after_coordinates(FILE * fplog,int64_t step,const t_inputrec * inputrec,const t_mdatoms * md,const matrix pressure,const matrix forceVirial,const matrix constraintVirial,matrix pressureCouplingMu,t_state * state,t_nrnb * nrnb,gmx::BoxDeformation * boxDeformation,const bool scaleCoordinates)186 void update_pcouple_after_coordinates(FILE* fplog,
187 int64_t step,
188 const t_inputrec* inputrec,
189 const t_mdatoms* md,
190 const matrix pressure,
191 const matrix forceVirial,
192 const matrix constraintVirial,
193 matrix pressureCouplingMu,
194 t_state* state,
195 t_nrnb* nrnb,
196 gmx::BoxDeformation* boxDeformation,
197 const bool scaleCoordinates)
198 {
199 int start = 0;
200 int homenr = md->homenr;
201
202 /* Cast to real for faster code, no loss in precision (see comment above) */
203 real dt = inputrec->delta_t;
204
205
206 /* now update boxes */
207 switch (inputrec->epc)
208 {
209 case (epcNO): break;
210 case (epcBERENDSEN):
211 if (do_per_step(step, inputrec->nstpcouple))
212 {
213 real dtpc = inputrec->nstpcouple * dt;
214 berendsen_pcoupl(fplog, step, inputrec, dtpc, pressure, state->box, forceVirial,
215 constraintVirial, pressureCouplingMu, &state->baros_integral);
216 berendsen_pscale(inputrec, pressureCouplingMu, state->box, state->box_rel, start,
217 homenr, state->x.rvec_array(), md->cFREEZE, nrnb, scaleCoordinates);
218 }
219 break;
220 case (epcCRESCALE):
221 if (do_per_step(step, inputrec->nstpcouple))
222 {
223 real dtpc = inputrec->nstpcouple * dt;
224 crescale_pcoupl(fplog, step, inputrec, dtpc, pressure, state->box, forceVirial,
225 constraintVirial, pressureCouplingMu, &state->baros_integral);
226 crescale_pscale(inputrec, pressureCouplingMu, state->box, state->box_rel, start,
227 homenr, state->x.rvec_array(), state->v.rvec_array(), md->cFREEZE,
228 nrnb, scaleCoordinates);
229 }
230 break;
231 case (epcPARRINELLORAHMAN):
232 if (do_per_step(step + inputrec->nstpcouple - 1, inputrec->nstpcouple))
233 {
234 /* The box velocities were updated in do_pr_pcoupl,
235 * but we dont change the box vectors until we get here
236 * since we need to be able to shift/unshift above.
237 */
238 real dtpc = inputrec->nstpcouple * dt;
239 for (int i = 0; i < DIM; i++)
240 {
241 for (int m = 0; m <= i; m++)
242 {
243 state->box[i][m] += dtpc * state->boxv[i][m];
244 }
245 }
246 preserve_box_shape(inputrec, state->box_rel, state->box);
247
248 /* Scale the coordinates */
249 if (scaleCoordinates)
250 {
251 auto x = state->x.rvec_array();
252 for (int n = start; n < start + homenr; n++)
253 {
254 tmvmul_ur0(pressureCouplingMu, x[n], x[n]);
255 }
256 }
257 }
258 break;
259 case (epcMTTK):
260 switch (inputrec->epct)
261 {
262 case (epctISOTROPIC):
263 /* DIM * eta = ln V. so DIM*eta_new = DIM*eta_old + DIM*dt*veta =>
264 ln V_new = ln V_old + 3*dt*veta => V_new = V_old*exp(3*dt*veta) =>
265 Side length scales as exp(veta*dt) */
266
267 msmul(state->box, std::exp(state->veta * dt), state->box);
268
269 /* Relate veta to boxv. veta = d(eta)/dT = (1/DIM)*1/V dV/dT.
270 o If we assume isotropic scaling, and box length scaling
271 factor L, then V = L^DIM (det(M)). So dV/dt = DIM
272 L^(DIM-1) dL/dt det(M), and veta = (1/L) dL/dt. The
273 determinant of B is L^DIM det(M), and the determinant
274 of dB/dt is (dL/dT)^DIM det (M). veta will be
275 (det(dB/dT)/det(B))^(1/3). Then since M =
276 B_new*(vol_new)^(1/3), dB/dT_new = (veta_new)*B(new). */
277
278 msmul(state->box, state->veta, state->boxv);
279 break;
280 default: break;
281 }
282 break;
283 default: break;
284 }
285
286 if (boxDeformation)
287 {
288 auto localX = makeArrayRef(state->x).subArray(start, homenr);
289 boxDeformation->apply(localX, state->box, step);
290 }
291 }
292
update_randomize_velocities(const t_inputrec * ir,int64_t step,const t_commrec * cr,const t_mdatoms * md,gmx::ArrayRef<gmx::RVec> v,const gmx::Update * upd,const gmx::Constraints * constr)293 extern gmx_bool update_randomize_velocities(const t_inputrec* ir,
294 int64_t step,
295 const t_commrec* cr,
296 const t_mdatoms* md,
297 gmx::ArrayRef<gmx::RVec> v,
298 const gmx::Update* upd,
299 const gmx::Constraints* constr)
300 {
301
302 real rate = (ir->delta_t) / ir->opts.tau_t[0];
303
304 if (ir->etc == etcANDERSEN && constr != nullptr)
305 {
306 /* Currently, Andersen thermostat does not support constrained
307 systems. Functionality exists in the andersen_tcoupl
308 function in GROMACS 4.5.7 to allow this combination. That
309 code could be ported to the current random-number
310 generation approach, but has not yet been done because of
311 lack of time and resources. */
312 gmx_fatal(FARGS,
313 "Normal Andersen is currently not supported with constraints, use massive "
314 "Andersen instead");
315 }
316
317 /* proceed with andersen if 1) it's fixed probability per
318 particle andersen or 2) it's massive andersen and it's tau_t/dt */
319 if ((ir->etc == etcANDERSEN) || do_per_step(step, gmx::roundToInt(1.0 / rate)))
320 {
321 andersen_tcoupl(ir, step, cr, md, v, rate, upd->getAndersenRandomizeGroup(),
322 upd->getBoltzmanFactor());
323 return TRUE;
324 }
325 return FALSE;
326 }
327
328 /*
329 static const double sy_const[MAX_SUZUKI_YOSHIDA_NUM+1][MAX_SUZUKI_YOSHIDA_NUM+1] = {
330 {},
331 {1},
332 {},
333 {0.828981543588751,-0.657963087177502,0.828981543588751},
334 {},
335 {0.2967324292201065,0.2967324292201065,-0.186929716880426,0.2967324292201065,0.2967324292201065}
336 };*/
337
338 /* these integration routines are only referenced inside this file */
NHC_trotter(const t_grpopts * opts,int nvar,const gmx_ekindata_t * ekind,real dtfull,double xi[],double vxi[],double scalefac[],real * veta,const t_extmass * MassQ,gmx_bool bEkinAveVel)339 static void NHC_trotter(const t_grpopts* opts,
340 int nvar,
341 const gmx_ekindata_t* ekind,
342 real dtfull,
343 double xi[],
344 double vxi[],
345 double scalefac[],
346 real* veta,
347 const t_extmass* MassQ,
348 gmx_bool bEkinAveVel)
349
350 {
351 /* general routine for both barostat and thermostat nose hoover chains */
352
353 int i, j, mi, mj;
354 double Ekin, Efac, reft, kT, nd;
355 double dt;
356 double * ivxi, *ixi;
357 double* GQ;
358 gmx_bool bBarostat;
359 int mstepsi, mstepsj;
360 int ns = SUZUKI_YOSHIDA_NUM; /* set the degree of integration in the types/state.h file */
361 int nh = opts->nhchainlength;
362
363 snew(GQ, nh);
364 mstepsi = mstepsj = ns;
365
366 /* if scalefac is NULL, we are doing the NHC of the barostat */
367
368 bBarostat = FALSE;
369 if (scalefac == nullptr)
370 {
371 bBarostat = TRUE;
372 }
373
374 for (i = 0; i < nvar; i++)
375 {
376
377 /* make it easier to iterate by selecting
378 out the sub-array that corresponds to this T group */
379
380 ivxi = &vxi[i * nh];
381 ixi = &xi[i * nh];
382 gmx::ArrayRef<const double> iQinv;
383 if (bBarostat)
384 {
385 iQinv = gmx::arrayRefFromArray(&MassQ->QPinv[i * nh], nh);
386 nd = 1.0; /* THIS WILL CHANGE IF NOT ISOTROPIC */
387 reft = std::max<real>(0, opts->ref_t[0]);
388 Ekin = gmx::square(*veta) / MassQ->Winv;
389 }
390 else
391 {
392 iQinv = gmx::arrayRefFromArray(&MassQ->Qinv[i * nh], nh);
393 const t_grp_tcstat* tcstat = &ekind->tcstat[i];
394 nd = opts->nrdf[i];
395 reft = std::max<real>(0, opts->ref_t[i]);
396 if (bEkinAveVel)
397 {
398 Ekin = 2 * trace(tcstat->ekinf) * tcstat->ekinscalef_nhc;
399 }
400 else
401 {
402 Ekin = 2 * trace(tcstat->ekinh) * tcstat->ekinscaleh_nhc;
403 }
404 }
405 kT = BOLTZ * reft;
406
407 for (mi = 0; mi < mstepsi; mi++)
408 {
409 for (mj = 0; mj < mstepsj; mj++)
410 {
411 /* weighting for this step using Suzuki-Yoshida integration - fixed at 5 */
412 dt = sy_const[ns][mj] * dtfull / mstepsi;
413
414 /* compute the thermal forces */
415 GQ[0] = iQinv[0] * (Ekin - nd * kT);
416
417 for (j = 0; j < nh - 1; j++)
418 {
419 if (iQinv[j + 1] > 0)
420 {
421 /* we actually don't need to update here if we save the
422 state of the GQ, but it's easier to just recompute*/
423 GQ[j + 1] = iQinv[j + 1] * ((gmx::square(ivxi[j]) / iQinv[j]) - kT);
424 }
425 else
426 {
427 GQ[j + 1] = 0;
428 }
429 }
430
431 ivxi[nh - 1] += 0.25 * dt * GQ[nh - 1];
432 for (j = nh - 1; j > 0; j--)
433 {
434 Efac = exp(-0.125 * dt * ivxi[j]);
435 ivxi[j - 1] = Efac * (ivxi[j - 1] * Efac + 0.25 * dt * GQ[j - 1]);
436 }
437
438 Efac = exp(-0.5 * dt * ivxi[0]);
439 if (bBarostat)
440 {
441 *veta *= Efac;
442 }
443 else
444 {
445 scalefac[i] *= Efac;
446 }
447 Ekin *= (Efac * Efac);
448
449 /* Issue - if the KE is an average of the last and the current temperatures, then we
450 might not be able to scale the kinetic energy directly with this factor. Might
451 take more bookkeeping -- have to think about this a bit more . . . */
452
453 GQ[0] = iQinv[0] * (Ekin - nd * kT);
454
455 /* update thermostat positions */
456 for (j = 0; j < nh; j++)
457 {
458 ixi[j] += 0.5 * dt * ivxi[j];
459 }
460
461 for (j = 0; j < nh - 1; j++)
462 {
463 Efac = exp(-0.125 * dt * ivxi[j + 1]);
464 ivxi[j] = Efac * (ivxi[j] * Efac + 0.25 * dt * GQ[j]);
465 if (iQinv[j + 1] > 0)
466 {
467 GQ[j + 1] = iQinv[j + 1] * ((gmx::square(ivxi[j]) / iQinv[j]) - kT);
468 }
469 else
470 {
471 GQ[j + 1] = 0;
472 }
473 }
474 ivxi[nh - 1] += 0.25 * dt * GQ[nh - 1];
475 }
476 }
477 }
478 sfree(GQ);
479 }
480
boxv_trotter(const t_inputrec * ir,real * veta,real dt,const tensor box,const gmx_ekindata_t * ekind,const tensor vir,real pcorr,const t_extmass * MassQ)481 static void boxv_trotter(const t_inputrec* ir,
482 real* veta,
483 real dt,
484 const tensor box,
485 const gmx_ekindata_t* ekind,
486 const tensor vir,
487 real pcorr,
488 const t_extmass* MassQ)
489 {
490
491 real pscal;
492 double alpha;
493 int nwall;
494 real GW, vol;
495 tensor ekinmod, localpres;
496
497 /* The heat bath is coupled to a separate barostat, the last temperature group. In the
498 2006 Tuckerman et al paper., the order is iL_{T_baro} iL {T_part}
499 */
500
501 if (ir->epct == epctSEMIISOTROPIC)
502 {
503 nwall = 2;
504 }
505 else
506 {
507 nwall = 3;
508 }
509
510 /* eta is in pure units. veta is in units of ps^-1. GW is in
511 units of ps^-2. However, eta has a reference of 1 nm^3, so care must be
512 taken to use only RATIOS of eta in updating the volume. */
513
514 /* we take the partial pressure tensors, modify the
515 kinetic energy tensor, and recovert to pressure */
516
517 if (ir->opts.nrdf[0] == 0)
518 {
519 gmx_fatal(FARGS, "Barostat is coupled to a T-group with no degrees of freedom\n");
520 }
521 /* alpha factor for phase space volume, then multiply by the ekin scaling factor. */
522 alpha = 1.0 + DIM / (static_cast<double>(ir->opts.nrdf[0]));
523 alpha *= ekind->tcstat[0].ekinscalef_nhc;
524 msmul(ekind->ekin, alpha, ekinmod);
525 /* for now, we use Elr = 0, because if you want to get it right, you
526 really should be using PME. Maybe print a warning? */
527
528 pscal = calc_pres(ir->pbcType, nwall, box, ekinmod, vir, localpres) + pcorr;
529
530 vol = det(box);
531 GW = (vol * (MassQ->Winv / PRESFAC)) * (DIM * pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
532
533 *veta += 0.5 * dt * GW;
534 }
535
536 /*
537 * This file implements temperature and pressure coupling algorithms:
538 * For now only the Weak coupling and the modified weak coupling.
539 *
540 * Furthermore computation of pressure and temperature is done here
541 *
542 */
543
calc_pres(PbcType pbcType,int nwall,const matrix box,const tensor ekin,const tensor vir,tensor pres)544 real calc_pres(PbcType pbcType, int nwall, const matrix box, const tensor ekin, const tensor vir, tensor pres)
545 {
546 int n, m;
547 real fac;
548
549 if (pbcType == PbcType::No || (pbcType == PbcType::XY && nwall != 2))
550 {
551 clear_mat(pres);
552 }
553 else
554 {
555 /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
556 * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
557 * het systeem...
558 */
559
560 fac = PRESFAC * 2.0 / det(box);
561 for (n = 0; (n < DIM); n++)
562 {
563 for (m = 0; (m < DIM); m++)
564 {
565 pres[n][m] = (ekin[n][m] - vir[n][m]) * fac;
566 }
567 }
568
569 if (debug)
570 {
571 pr_rvecs(debug, 0, "PC: pres", pres, DIM);
572 pr_rvecs(debug, 0, "PC: ekin", ekin, DIM);
573 pr_rvecs(debug, 0, "PC: vir ", vir, DIM);
574 pr_rvecs(debug, 0, "PC: box ", box, DIM);
575 }
576 }
577 return trace(pres) / DIM;
578 }
579
calc_temp(real ekin,real nrdf)580 real calc_temp(real ekin, real nrdf)
581 {
582 if (nrdf > 0)
583 {
584 return (2.0 * ekin) / (nrdf * BOLTZ);
585 }
586 else
587 {
588 return 0;
589 }
590 }
591
592 /*! \brief Sets 1/mass for Parrinello-Rahman in wInv; NOTE: PRESFAC is not included, so not in GROMACS units! */
calcParrinelloRahmanInvMass(const t_inputrec * ir,const matrix box,tensor wInv)593 static void calcParrinelloRahmanInvMass(const t_inputrec* ir, const matrix box, tensor wInv)
594 {
595 real maxBoxLength;
596
597 /* TODO: See if we can make the mass independent of the box size */
598 maxBoxLength = std::max(box[XX][XX], box[YY][YY]);
599 maxBoxLength = std::max(maxBoxLength, box[ZZ][ZZ]);
600
601 for (int d = 0; d < DIM; d++)
602 {
603 for (int n = 0; n < DIM; n++)
604 {
605 wInv[d][n] = (4 * M_PI * M_PI * ir->compress[d][n]) / (3 * ir->tau_p * ir->tau_p * maxBoxLength);
606 }
607 }
608 }
609
parrinellorahman_pcoupl(FILE * fplog,int64_t step,const t_inputrec * ir,real dt,const tensor pres,const tensor box,tensor box_rel,tensor boxv,tensor M,matrix mu,gmx_bool bFirstStep)610 void parrinellorahman_pcoupl(FILE* fplog,
611 int64_t step,
612 const t_inputrec* ir,
613 real dt,
614 const tensor pres,
615 const tensor box,
616 tensor box_rel,
617 tensor boxv,
618 tensor M,
619 matrix mu,
620 gmx_bool bFirstStep)
621 {
622 /* This doesn't do any coordinate updating. It just
623 * integrates the box vector equations from the calculated
624 * acceleration due to pressure difference. We also compute
625 * the tensor M which is used in update to couple the particle
626 * coordinates to the box vectors.
627 *
628 * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
629 * given as
630 * -1 . . -1
631 * M_nk = (h') * (h' * h + h' h) * h
632 *
633 * with the dots denoting time derivatives and h is the transformation from
634 * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
635 * This also goes for the pressure and M tensors - they are transposed relative
636 * to ours. Our equation thus becomes:
637 *
638 * -1 . . -1
639 * M_gmx = M_nk' = b * (b * b' + b * b') * b'
640 *
641 * where b is the gromacs box matrix.
642 * Our box accelerations are given by
643 * .. ..
644 * b = vol/W inv(box') * (P-ref_P) (=h')
645 */
646
647 real vol = box[XX][XX] * box[YY][YY] * box[ZZ][ZZ];
648 real atot, arel, change;
649 tensor invbox, pdiff, t1, t2;
650
651 gmx::invertBoxMatrix(box, invbox);
652
653 if (!bFirstStep)
654 {
655 /* Note that PRESFAC does not occur here.
656 * The pressure and compressibility always occur as a product,
657 * therefore the pressure unit drops out.
658 */
659 tensor winv;
660 calcParrinelloRahmanInvMass(ir, box, winv);
661
662 m_sub(pres, ir->ref_p, pdiff);
663
664 if (ir->epct == epctSURFACETENSION)
665 {
666 /* Unlike Berendsen coupling it might not be trivial to include a z
667 * pressure correction here? On the other hand we don't scale the
668 * box momentarily, but change accelerations, so it might not be crucial.
669 */
670 real xy_pressure = 0.5 * (pres[XX][XX] + pres[YY][YY]);
671 for (int d = 0; d < ZZ; d++)
672 {
673 pdiff[d][d] = (xy_pressure - (pres[ZZ][ZZ] - ir->ref_p[d][d] / box[d][d]));
674 }
675 }
676
677 tmmul(invbox, pdiff, t1);
678 /* Move the off-diagonal elements of the 'force' to one side to ensure
679 * that we obey the box constraints.
680 */
681 for (int d = 0; d < DIM; d++)
682 {
683 for (int n = 0; n < d; n++)
684 {
685 t1[d][n] += t1[n][d];
686 t1[n][d] = 0;
687 }
688 }
689
690 switch (ir->epct)
691 {
692 case epctANISOTROPIC:
693 for (int d = 0; d < DIM; d++)
694 {
695 for (int n = 0; n <= d; n++)
696 {
697 t1[d][n] *= winv[d][n] * vol;
698 }
699 }
700 break;
701 case epctISOTROPIC:
702 /* calculate total volume acceleration */
703 atot = box[XX][XX] * box[YY][YY] * t1[ZZ][ZZ] + box[XX][XX] * t1[YY][YY] * box[ZZ][ZZ]
704 + t1[XX][XX] * box[YY][YY] * box[ZZ][ZZ];
705 arel = atot / (3 * vol);
706 /* set all RELATIVE box accelerations equal, and maintain total V
707 * change speed */
708 for (int d = 0; d < DIM; d++)
709 {
710 for (int n = 0; n <= d; n++)
711 {
712 t1[d][n] = winv[0][0] * vol * arel * box[d][n];
713 }
714 }
715 break;
716 case epctSEMIISOTROPIC:
717 case epctSURFACETENSION:
718 /* Note the correction to pdiff above for surftens. coupling */
719
720 /* calculate total XY volume acceleration */
721 atot = box[XX][XX] * t1[YY][YY] + t1[XX][XX] * box[YY][YY];
722 arel = atot / (2 * box[XX][XX] * box[YY][YY]);
723 /* set RELATIVE XY box accelerations equal, and maintain total V
724 * change speed. Dont change the third box vector accelerations */
725 for (int d = 0; d < ZZ; d++)
726 {
727 for (int n = 0; n <= d; n++)
728 {
729 t1[d][n] = winv[d][n] * vol * arel * box[d][n];
730 }
731 }
732 for (int n = 0; n < DIM; n++)
733 {
734 t1[ZZ][n] *= winv[ZZ][n] * vol;
735 }
736 break;
737 default:
738 gmx_fatal(FARGS,
739 "Parrinello-Rahman pressure coupling type %s "
740 "not supported yet\n",
741 EPCOUPLTYPETYPE(ir->epct));
742 }
743
744 real maxchange = 0;
745 for (int d = 0; d < DIM; d++)
746 {
747 for (int n = 0; n <= d; n++)
748 {
749 boxv[d][n] += dt * t1[d][n];
750
751 /* We do NOT update the box vectors themselves here, since
752 * we need them for shifting later. It is instead done last
753 * in the update() routine.
754 */
755
756 /* Calculate the change relative to diagonal elements-
757 since it's perfectly ok for the off-diagonal ones to
758 be zero it doesn't make sense to check the change relative
759 to its current size.
760 */
761
762 change = std::fabs(dt * boxv[d][n] / box[d][d]);
763
764 if (change > maxchange)
765 {
766 maxchange = change;
767 }
768 }
769 }
770
771 if (maxchange > 0.01 && fplog)
772 {
773 char buf[22];
774 fprintf(fplog,
775 "\nStep %s Warning: Pressure scaling more than 1%%. "
776 "This may mean your system\n is not yet equilibrated. "
777 "Use of Parrinello-Rahman pressure coupling during\n"
778 "equilibration can lead to simulation instability, "
779 "and is discouraged.\n",
780 gmx_step_str(step, buf));
781 }
782 }
783
784 preserve_box_shape(ir, box_rel, boxv);
785
786 mtmul(boxv, box, t1); /* t1=boxv * b' */
787 mmul(invbox, t1, t2);
788 mtmul(t2, invbox, M);
789
790 /* Determine the scaling matrix mu for the coordinates */
791 for (int d = 0; d < DIM; d++)
792 {
793 for (int n = 0; n <= d; n++)
794 {
795 t1[d][n] = box[d][n] + dt * boxv[d][n];
796 }
797 }
798 preserve_box_shape(ir, box_rel, t1);
799 /* t1 is the box at t+dt, determine mu as the relative change */
800 mmul_ur0(invbox, t1, mu);
801 }
802
berendsen_pcoupl(FILE * fplog,int64_t step,const t_inputrec * ir,real dt,const tensor pres,const matrix box,const matrix force_vir,const matrix constraint_vir,matrix mu,double * baros_integral)803 void berendsen_pcoupl(FILE* fplog,
804 int64_t step,
805 const t_inputrec* ir,
806 real dt,
807 const tensor pres,
808 const matrix box,
809 const matrix force_vir,
810 const matrix constraint_vir,
811 matrix mu,
812 double* baros_integral)
813 {
814 real scalar_pressure, xy_pressure, p_corr_z;
815 char buf[STRLEN];
816
817 /*
818 * Calculate the scaling matrix mu
819 */
820 scalar_pressure = 0;
821 xy_pressure = 0;
822 for (int d = 0; d < DIM; d++)
823 {
824 scalar_pressure += pres[d][d] / DIM;
825 if (d != ZZ)
826 {
827 xy_pressure += pres[d][d] / (DIM - 1);
828 }
829 }
830 /* Pressure is now in bar, everywhere. */
831 #define factor(d, m) (ir->compress[d][m] * dt / ir->tau_p)
832
833 /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
834 * necessary for triclinic scaling
835 */
836 clear_mat(mu);
837 switch (ir->epct)
838 {
839 case epctISOTROPIC:
840 for (int d = 0; d < DIM; d++)
841 {
842 mu[d][d] = 1.0 - factor(d, d) * (ir->ref_p[d][d] - scalar_pressure) / DIM;
843 }
844 break;
845 case epctSEMIISOTROPIC:
846 for (int d = 0; d < ZZ; d++)
847 {
848 mu[d][d] = 1.0 - factor(d, d) * (ir->ref_p[d][d] - xy_pressure) / DIM;
849 }
850 mu[ZZ][ZZ] = 1.0 - factor(ZZ, ZZ) * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) / DIM;
851 break;
852 case epctANISOTROPIC:
853 for (int d = 0; d < DIM; d++)
854 {
855 for (int n = 0; n < DIM; n++)
856 {
857 mu[d][n] = (d == n ? 1.0 : 0.0) - factor(d, n) * (ir->ref_p[d][n] - pres[d][n]) / DIM;
858 }
859 }
860 break;
861 case epctSURFACETENSION:
862 /* ir->ref_p[0/1] is the reference surface-tension times *
863 * the number of surfaces */
864 if (ir->compress[ZZ][ZZ] != 0.0F)
865 {
866 p_corr_z = dt / ir->tau_p * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]);
867 }
868 else
869 {
870 /* when the compressibity is zero, set the pressure correction *
871 * in the z-direction to zero to get the correct surface tension */
872 p_corr_z = 0;
873 }
874 mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ] * p_corr_z;
875 for (int d = 0; d < DIM - 1; d++)
876 {
877 mu[d][d] = 1.0
878 + factor(d, d)
879 * (ir->ref_p[d][d] / (mu[ZZ][ZZ] * box[ZZ][ZZ])
880 - (pres[ZZ][ZZ] + p_corr_z - xy_pressure))
881 / (DIM - 1);
882 }
883 break;
884 default:
885 gmx_fatal(FARGS, "Berendsen pressure coupling type %s not supported yet\n",
886 EPCOUPLTYPETYPE(ir->epct));
887 }
888 /* To fullfill the orientation restrictions on triclinic boxes
889 * we will set mu_yx, mu_zx and mu_zy to 0 and correct
890 * the other elements of mu to first order.
891 */
892 mu[YY][XX] += mu[XX][YY];
893 mu[ZZ][XX] += mu[XX][ZZ];
894 mu[ZZ][YY] += mu[YY][ZZ];
895 mu[XX][YY] = 0;
896 mu[XX][ZZ] = 0;
897 mu[YY][ZZ] = 0;
898
899 /* Keep track of the work the barostat applies on the system.
900 * Without constraints force_vir tells us how Epot changes when scaling.
901 * With constraints constraint_vir gives us the constraint contribution
902 * to both Epot and Ekin. Although we are not scaling velocities, scaling
903 * the coordinates leads to scaling of distances involved in constraints.
904 * This in turn changes the angular momentum (even if the constrained
905 * distances are corrected at the next step). The kinetic component
906 * of the constraint virial captures the angular momentum change.
907 */
908 for (int d = 0; d < DIM; d++)
909 {
910 for (int n = 0; n <= d; n++)
911 {
912 *baros_integral -=
913 2 * (mu[d][n] - (n == d ? 1 : 0)) * (force_vir[d][n] + constraint_vir[d][n]);
914 }
915 }
916
917 if (debug)
918 {
919 pr_rvecs(debug, 0, "PC: pres ", pres, 3);
920 pr_rvecs(debug, 0, "PC: mu ", mu, 3);
921 }
922
923 if (mu[XX][XX] < 0.99 || mu[XX][XX] > 1.01 || mu[YY][YY] < 0.99 || mu[YY][YY] > 1.01
924 || mu[ZZ][ZZ] < 0.99 || mu[ZZ][ZZ] > 1.01)
925 {
926 char buf2[22];
927 sprintf(buf,
928 "\nStep %s Warning: pressure scaling more than 1%%, "
929 "mu: %g %g %g\n",
930 gmx_step_str(step, buf2), mu[XX][XX], mu[YY][YY], mu[ZZ][ZZ]);
931 if (fplog)
932 {
933 fprintf(fplog, "%s", buf);
934 }
935 fprintf(stderr, "%s", buf);
936 }
937 }
938
crescale_pcoupl(FILE * fplog,int64_t step,const t_inputrec * ir,real dt,const tensor pres,const matrix box,const matrix force_vir,const matrix constraint_vir,matrix mu,double * baros_integral)939 void crescale_pcoupl(FILE* fplog,
940 int64_t step,
941 const t_inputrec* ir,
942 real dt,
943 const tensor pres,
944 const matrix box,
945 const matrix force_vir,
946 const matrix constraint_vir,
947 matrix mu,
948 double* baros_integral)
949 {
950 /*
951 * Calculate the scaling matrix mu
952 */
953 real scalar_pressure = 0;
954 real xy_pressure = 0;
955 for (int d = 0; d < DIM; d++)
956 {
957 scalar_pressure += pres[d][d] / DIM;
958 if (d != ZZ)
959 {
960 xy_pressure += pres[d][d] / (DIM - 1);
961 }
962 }
963 clear_mat(mu);
964
965 gmx::ThreeFry2x64<64> rng(ir->ld_seed, gmx::RandomDomain::Barostat);
966 gmx::NormalDistribution<real> normalDist;
967 rng.restart(step, 0);
968 real vol = 1.0;
969 for (int d = 0; d < DIM; d++)
970 {
971 vol *= box[d][d];
972 }
973 real gauss;
974 real gauss2;
975 real kt = ir->opts.ref_t[0] * BOLTZ;
976 if (kt < 0.0)
977 {
978 kt = 0.0;
979 }
980
981 switch (ir->epct)
982 {
983 case epctISOTROPIC:
984 gauss = normalDist(rng);
985 for (int d = 0; d < DIM; d++)
986 {
987 const real compressibilityFactor = ir->compress[d][d] * dt / ir->tau_p;
988 mu[d][d] = std::exp(-compressibilityFactor * (ir->ref_p[d][d] - scalar_pressure) / DIM
989 + std::sqrt(2.0 * kt * compressibilityFactor * PRESFAC / vol)
990 * gauss / DIM);
991 }
992 break;
993 case epctSEMIISOTROPIC:
994 gauss = normalDist(rng);
995 gauss2 = normalDist(rng);
996 for (int d = 0; d < ZZ; d++)
997 {
998 const real compressibilityFactor = ir->compress[d][d] * dt / ir->tau_p;
999 mu[d][d] = std::exp(
1000 -compressibilityFactor * (ir->ref_p[d][d] - xy_pressure) / DIM
1001 + std::sqrt((DIM - 1) * 2.0 * kt * compressibilityFactor * PRESFAC / vol / DIM)
1002 / (DIM - 1) * gauss);
1003 }
1004 {
1005 const real compressibilityFactor = ir->compress[ZZ][ZZ] * dt / ir->tau_p;
1006 mu[ZZ][ZZ] = std::exp(
1007 -compressibilityFactor * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) / DIM
1008 + std::sqrt(2.0 * kt * compressibilityFactor * PRESFAC / vol / DIM) * gauss2);
1009 }
1010 break;
1011 case epctSURFACETENSION:
1012 gauss = normalDist(rng);
1013 gauss2 = normalDist(rng);
1014 for (int d = 0; d < ZZ; d++)
1015 {
1016 const real compressibilityFactor = ir->compress[d][d] * dt / ir->tau_p;
1017 /* Notice: we here use ref_p[ZZ][ZZ] as isotropic pressure and ir->ref_p[d][d] as surface tension */
1018 mu[d][d] = std::exp(
1019 -compressibilityFactor
1020 * (ir->ref_p[ZZ][ZZ] - ir->ref_p[d][d] / box[ZZ][ZZ] - xy_pressure) / DIM
1021 + std::sqrt(4.0 / 3.0 * kt * compressibilityFactor * PRESFAC / vol)
1022 / (DIM - 1) * gauss);
1023 }
1024 {
1025 const real compressibilityFactor = ir->compress[ZZ][ZZ] * dt / ir->tau_p;
1026 mu[ZZ][ZZ] = std::exp(
1027 -compressibilityFactor * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) / DIM
1028 + std::sqrt(2.0 / 3.0 * kt * compressibilityFactor * PRESFAC / vol) * gauss2);
1029 }
1030 break;
1031 default:
1032 gmx_fatal(FARGS, "C-rescale pressure coupling type %s not supported yet\n",
1033 EPCOUPLTYPETYPE(ir->epct));
1034 }
1035 /* To fullfill the orientation restrictions on triclinic boxes
1036 * we will set mu_yx, mu_zx and mu_zy to 0 and correct
1037 * the other elements of mu to first order.
1038 */
1039 mu[YY][XX] += mu[XX][YY];
1040 mu[ZZ][XX] += mu[XX][ZZ];
1041 mu[ZZ][YY] += mu[YY][ZZ];
1042 mu[XX][YY] = 0;
1043 mu[XX][ZZ] = 0;
1044 mu[YY][ZZ] = 0;
1045
1046 /* Keep track of the work the barostat applies on the system.
1047 * Without constraints force_vir tells us how Epot changes when scaling.
1048 * With constraints constraint_vir gives us the constraint contribution
1049 * to both Epot and Ekin. Although we are not scaling velocities, scaling
1050 * the coordinates leads to scaling of distances involved in constraints.
1051 * This in turn changes the angular momentum (even if the constrained
1052 * distances are corrected at the next step). The kinetic component
1053 * of the constraint virial captures the angular momentum change.
1054 */
1055 for (int d = 0; d < DIM; d++)
1056 {
1057 for (int n = 0; n <= d; n++)
1058 {
1059 *baros_integral -=
1060 2 * (mu[d][n] - (n == d ? 1 : 0)) * (force_vir[d][n] + constraint_vir[d][n]);
1061 }
1062 }
1063
1064 if (debug)
1065 {
1066 pr_rvecs(debug, 0, "PC: pres ", pres, 3);
1067 pr_rvecs(debug, 0, "PC: mu ", mu, 3);
1068 }
1069
1070 if (mu[XX][XX] < 0.99 || mu[XX][XX] > 1.01 || mu[YY][YY] < 0.99 || mu[YY][YY] > 1.01
1071 || mu[ZZ][ZZ] < 0.99 || mu[ZZ][ZZ] > 1.01)
1072 {
1073 char buf[STRLEN];
1074 char buf2[22];
1075 sprintf(buf,
1076 "\nStep %s Warning: pressure scaling more than 1%%, "
1077 "mu: %g %g %g\n",
1078 gmx_step_str(step, buf2), mu[XX][XX], mu[YY][YY], mu[ZZ][ZZ]);
1079 if (fplog)
1080 {
1081 fprintf(fplog, "%s", buf);
1082 }
1083 fprintf(stderr, "%s", buf);
1084 }
1085 }
1086
crescale_pscale(const t_inputrec * ir,const matrix mu,matrix box,matrix box_rel,int start,int nr_atoms,rvec x[],rvec v[],const unsigned short cFREEZE[],t_nrnb * nrnb,const bool scaleCoordinates)1087 void crescale_pscale(const t_inputrec* ir,
1088 const matrix mu,
1089 matrix box,
1090 matrix box_rel,
1091 int start,
1092 int nr_atoms,
1093 rvec x[],
1094 rvec v[],
1095 const unsigned short cFREEZE[],
1096 t_nrnb* nrnb,
1097 const bool scaleCoordinates)
1098 {
1099 ivec* nFreeze = ir->opts.nFreeze;
1100 int nthreads gmx_unused;
1101 matrix inv_mu;
1102
1103 #ifndef __clang_analyzer__
1104 nthreads = gmx_omp_nthreads_get(emntUpdate);
1105 #endif
1106
1107 gmx::invertBoxMatrix(mu, inv_mu);
1108
1109 /* Scale the positions and the velocities */
1110 if (scaleCoordinates)
1111 {
1112 #pragma omp parallel for num_threads(nthreads) schedule(static)
1113 for (int n = start; n < start + nr_atoms; n++)
1114 {
1115 // Trivial OpenMP region that does not throw
1116 int g;
1117
1118 if (cFREEZE == nullptr)
1119 {
1120 g = 0;
1121 }
1122 else
1123 {
1124 g = cFREEZE[n];
1125 }
1126
1127 if (!nFreeze[g][XX])
1128 {
1129 x[n][XX] = mu[XX][XX] * x[n][XX] + mu[YY][XX] * x[n][YY] + mu[ZZ][XX] * x[n][ZZ];
1130 v[n][XX] = inv_mu[XX][XX] * v[n][XX] + inv_mu[YY][XX] * v[n][YY]
1131 + inv_mu[ZZ][XX] * v[n][ZZ];
1132 }
1133 if (!nFreeze[g][YY])
1134 {
1135 x[n][YY] = mu[YY][YY] * x[n][YY] + mu[ZZ][YY] * x[n][ZZ];
1136 v[n][YY] = inv_mu[YY][YY] * v[n][YY] + inv_mu[ZZ][YY] * v[n][ZZ];
1137 }
1138 if (!nFreeze[g][ZZ])
1139 {
1140 x[n][ZZ] = mu[ZZ][ZZ] * x[n][ZZ];
1141 v[n][ZZ] = inv_mu[ZZ][ZZ] * v[n][ZZ];
1142 }
1143 }
1144 }
1145 /* compute final boxlengths */
1146 for (int d = 0; d < DIM; d++)
1147 {
1148 box[d][XX] = mu[XX][XX] * box[d][XX] + mu[YY][XX] * box[d][YY] + mu[ZZ][XX] * box[d][ZZ];
1149 box[d][YY] = mu[YY][YY] * box[d][YY] + mu[ZZ][YY] * box[d][ZZ];
1150 box[d][ZZ] = mu[ZZ][ZZ] * box[d][ZZ];
1151 }
1152
1153 preserve_box_shape(ir, box_rel, box);
1154
1155 /* (un)shifting should NOT be done after this,
1156 * since the box vectors might have changed
1157 */
1158 inc_nrnb(nrnb, eNR_PCOUPL, nr_atoms);
1159 }
1160
berendsen_pscale(const t_inputrec * ir,const matrix mu,matrix box,matrix box_rel,int start,int nr_atoms,rvec x[],const unsigned short cFREEZE[],t_nrnb * nrnb,const bool scaleCoordinates)1161 void berendsen_pscale(const t_inputrec* ir,
1162 const matrix mu,
1163 matrix box,
1164 matrix box_rel,
1165 int start,
1166 int nr_atoms,
1167 rvec x[],
1168 const unsigned short cFREEZE[],
1169 t_nrnb* nrnb,
1170 const bool scaleCoordinates)
1171 {
1172 ivec* nFreeze = ir->opts.nFreeze;
1173 int d;
1174 int nthreads gmx_unused;
1175
1176 #ifndef __clang_analyzer__
1177 nthreads = gmx_omp_nthreads_get(emntUpdate);
1178 #endif
1179
1180 /* Scale the positions */
1181 if (scaleCoordinates)
1182 {
1183 #pragma omp parallel for num_threads(nthreads) schedule(static)
1184 for (int n = start; n < start + nr_atoms; n++)
1185 {
1186 // Trivial OpenMP region that does not throw
1187 int g;
1188
1189 if (cFREEZE == nullptr)
1190 {
1191 g = 0;
1192 }
1193 else
1194 {
1195 g = cFREEZE[n];
1196 }
1197
1198 if (!nFreeze[g][XX])
1199 {
1200 x[n][XX] = mu[XX][XX] * x[n][XX] + mu[YY][XX] * x[n][YY] + mu[ZZ][XX] * x[n][ZZ];
1201 }
1202 if (!nFreeze[g][YY])
1203 {
1204 x[n][YY] = mu[YY][YY] * x[n][YY] + mu[ZZ][YY] * x[n][ZZ];
1205 }
1206 if (!nFreeze[g][ZZ])
1207 {
1208 x[n][ZZ] = mu[ZZ][ZZ] * x[n][ZZ];
1209 }
1210 }
1211 }
1212 /* compute final boxlengths */
1213 for (d = 0; d < DIM; d++)
1214 {
1215 box[d][XX] = mu[XX][XX] * box[d][XX] + mu[YY][XX] * box[d][YY] + mu[ZZ][XX] * box[d][ZZ];
1216 box[d][YY] = mu[YY][YY] * box[d][YY] + mu[ZZ][YY] * box[d][ZZ];
1217 box[d][ZZ] = mu[ZZ][ZZ] * box[d][ZZ];
1218 }
1219
1220 preserve_box_shape(ir, box_rel, box);
1221
1222 /* (un)shifting should NOT be done after this,
1223 * since the box vectors might have changed
1224 */
1225 inc_nrnb(nrnb, eNR_PCOUPL, nr_atoms);
1226 }
1227
berendsen_tcoupl(const t_inputrec * ir,gmx_ekindata_t * ekind,real dt,std::vector<double> & therm_integral)1228 void berendsen_tcoupl(const t_inputrec* ir, gmx_ekindata_t* ekind, real dt, std::vector<double>& therm_integral)
1229 {
1230 const t_grpopts* opts = &ir->opts;
1231
1232 for (int i = 0; (i < opts->ngtc); i++)
1233 {
1234 real Ek, T;
1235
1236 if (ir->eI == eiVV)
1237 {
1238 Ek = trace(ekind->tcstat[i].ekinf);
1239 T = ekind->tcstat[i].T;
1240 }
1241 else
1242 {
1243 Ek = trace(ekind->tcstat[i].ekinh);
1244 T = ekind->tcstat[i].Th;
1245 }
1246
1247 if ((opts->tau_t[i] > 0) && (T > 0.0))
1248 {
1249 real reft = std::max<real>(0, opts->ref_t[i]);
1250 real lll = std::sqrt(1.0 + (dt / opts->tau_t[i]) * (reft / T - 1.0));
1251 ekind->tcstat[i].lambda = std::max<real>(std::min<real>(lll, 1.25), 0.8);
1252 }
1253 else
1254 {
1255 ekind->tcstat[i].lambda = 1.0;
1256 }
1257
1258 /* Keep track of the amount of energy we are adding to the system */
1259 therm_integral[i] -= (gmx::square(ekind->tcstat[i].lambda) - 1) * Ek;
1260
1261 if (debug)
1262 {
1263 fprintf(debug, "TC: group %d: T: %g, Lambda: %g\n", i, T, ekind->tcstat[i].lambda);
1264 }
1265 }
1266 }
1267
andersen_tcoupl(const t_inputrec * ir,int64_t step,const t_commrec * cr,const t_mdatoms * md,gmx::ArrayRef<gmx::RVec> v,real rate,const std::vector<bool> & randomize,gmx::ArrayRef<const real> boltzfac)1268 void andersen_tcoupl(const t_inputrec* ir,
1269 int64_t step,
1270 const t_commrec* cr,
1271 const t_mdatoms* md,
1272 gmx::ArrayRef<gmx::RVec> v,
1273 real rate,
1274 const std::vector<bool>& randomize,
1275 gmx::ArrayRef<const real> boltzfac)
1276 {
1277 const int* gatindex = (DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
1278 int i;
1279 int gc = 0;
1280 gmx::ThreeFry2x64<0> rng(ir->andersen_seed, gmx::RandomDomain::Thermostat);
1281 gmx::UniformRealDistribution<real> uniformDist;
1282 gmx::TabulatedNormalDistribution<real, 14> normalDist;
1283
1284 /* randomize the velocities of the selected particles */
1285
1286 for (i = 0; i < md->homenr; i++) /* now loop over the list of atoms */
1287 {
1288 int ng = gatindex ? gatindex[i] : i;
1289 gmx_bool bRandomize;
1290
1291 rng.restart(step, ng);
1292
1293 if (md->cTC)
1294 {
1295 gc = md->cTC[i]; /* assign the atom to a temperature group if there are more than one */
1296 }
1297 if (randomize[gc])
1298 {
1299 if (ir->etc == etcANDERSENMASSIVE)
1300 {
1301 /* Randomize particle always */
1302 bRandomize = TRUE;
1303 }
1304 else
1305 {
1306 /* Randomize particle probabilistically */
1307 uniformDist.reset();
1308 bRandomize = uniformDist(rng) < rate;
1309 }
1310 if (bRandomize)
1311 {
1312 real scal;
1313 int d;
1314
1315 scal = std::sqrt(boltzfac[gc] * md->invmass[i]);
1316
1317 normalDist.reset();
1318
1319 for (d = 0; d < DIM; d++)
1320 {
1321 v[i][d] = scal * normalDist(rng);
1322 }
1323 }
1324 }
1325 }
1326 }
1327
1328
nosehoover_tcoupl(const t_grpopts * opts,const gmx_ekindata_t * ekind,real dt,double xi[],double vxi[],const t_extmass * MassQ)1329 void nosehoover_tcoupl(const t_grpopts* opts,
1330 const gmx_ekindata_t* ekind,
1331 real dt,
1332 double xi[],
1333 double vxi[],
1334 const t_extmass* MassQ)
1335 {
1336 int i;
1337 real reft, oldvxi;
1338
1339 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
1340
1341 for (i = 0; (i < opts->ngtc); i++)
1342 {
1343 reft = std::max<real>(0, opts->ref_t[i]);
1344 oldvxi = vxi[i];
1345 vxi[i] += dt * MassQ->Qinv[i] * (ekind->tcstat[i].Th - reft);
1346 xi[i] += dt * (oldvxi + vxi[i]) * 0.5;
1347 }
1348 }
1349
trotter_update(const t_inputrec * ir,int64_t step,gmx_ekindata_t * ekind,const gmx_enerdata_t * enerd,t_state * state,const tensor vir,const t_mdatoms * md,const t_extmass * MassQ,gmx::ArrayRef<std::vector<int>> trotter_seqlist,int trotter_seqno)1350 void trotter_update(const t_inputrec* ir,
1351 int64_t step,
1352 gmx_ekindata_t* ekind,
1353 const gmx_enerdata_t* enerd,
1354 t_state* state,
1355 const tensor vir,
1356 const t_mdatoms* md,
1357 const t_extmass* MassQ,
1358 gmx::ArrayRef<std::vector<int>> trotter_seqlist,
1359 int trotter_seqno)
1360 {
1361
1362 int n, i, d, ngtc, gc = 0, t;
1363 t_grp_tcstat* tcstat;
1364 const t_grpopts* opts;
1365 int64_t step_eff;
1366 real dt;
1367 double * scalefac, dtc;
1368 rvec sumv = { 0, 0, 0 };
1369 gmx_bool bCouple;
1370
1371 if (trotter_seqno <= ettTSEQ2)
1372 {
1373 step_eff = step - 1; /* the velocity verlet calls are actually out of order -- the first
1374 half step is actually the last half step from the previous step.
1375 Thus the first half step actually corresponds to the n-1 step*/
1376 }
1377 else
1378 {
1379 step_eff = step;
1380 }
1381
1382 bCouple = (ir->nsttcouple == 1 || do_per_step(step_eff + ir->nsttcouple, ir->nsttcouple));
1383
1384 const gmx::ArrayRef<const int> trotter_seq = trotter_seqlist[trotter_seqno];
1385
1386 if ((trotter_seq[0] == etrtSKIPALL) || (!bCouple))
1387 {
1388 return;
1389 }
1390 dtc = ir->nsttcouple * ir->delta_t; /* This is OK for NPT, because nsttcouple == nstpcouple is enforcesd */
1391 opts = &(ir->opts); /* just for ease of referencing */
1392 ngtc = opts->ngtc;
1393 assert(ngtc > 0);
1394 snew(scalefac, opts->ngtc);
1395 for (i = 0; i < ngtc; i++)
1396 {
1397 scalefac[i] = 1;
1398 }
1399 /* execute the series of trotter updates specified in the trotterpart array */
1400
1401 for (i = 0; i < NTROTTERPARTS; i++)
1402 {
1403 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
1404 if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2)
1405 || (trotter_seq[i] == etrtNHC2))
1406 {
1407 dt = 2 * dtc;
1408 }
1409 else
1410 {
1411 dt = dtc;
1412 }
1413
1414 auto v = makeArrayRef(state->v);
1415 switch (trotter_seq[i])
1416 {
1417 case etrtBAROV:
1418 case etrtBAROV2:
1419 boxv_trotter(ir, &(state->veta), dt, state->box, ekind, vir,
1420 enerd->term[F_PDISPCORR], MassQ);
1421 break;
1422 case etrtBARONHC:
1423 case etrtBARONHC2:
1424 NHC_trotter(opts, state->nnhpres, ekind, dt, state->nhpres_xi.data(),
1425 state->nhpres_vxi.data(), nullptr, &(state->veta), MassQ, FALSE);
1426 break;
1427 case etrtNHC:
1428 case etrtNHC2:
1429 NHC_trotter(opts, opts->ngtc, ekind, dt, state->nosehoover_xi.data(),
1430 state->nosehoover_vxi.data(), scalefac, nullptr, MassQ, (ir->eI == eiVV));
1431 /* need to rescale the kinetic energies and velocities here. Could
1432 scale the velocities later, but we need them scaled in order to
1433 produce the correct outputs, so we'll scale them here. */
1434
1435 for (t = 0; t < ngtc; t++)
1436 {
1437 tcstat = &ekind->tcstat[t];
1438 tcstat->vscale_nhc = scalefac[t];
1439 tcstat->ekinscaleh_nhc *= (scalefac[t] * scalefac[t]);
1440 tcstat->ekinscalef_nhc *= (scalefac[t] * scalefac[t]);
1441 }
1442 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
1443 /* but do we actually need the total? */
1444
1445 /* modify the velocities as well */
1446 for (n = 0; n < md->homenr; n++)
1447 {
1448 if (md->cTC) /* does this conditional need to be here? is this always true?*/
1449 {
1450 gc = md->cTC[n];
1451 }
1452 for (d = 0; d < DIM; d++)
1453 {
1454 v[n][d] *= scalefac[gc];
1455 }
1456
1457 if (debug)
1458 {
1459 for (d = 0; d < DIM; d++)
1460 {
1461 sumv[d] += (v[n][d]) / md->invmass[n];
1462 }
1463 }
1464 }
1465 break;
1466 default: break;
1467 }
1468 }
1469 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
1470 sfree(scalefac);
1471 }
1472
1473
init_npt_masses(const t_inputrec * ir,t_state * state,t_extmass * MassQ,gmx_bool bInit)1474 extern void init_npt_masses(const t_inputrec* ir, t_state* state, t_extmass* MassQ, gmx_bool bInit)
1475 {
1476 int n, i, j, d, ngtc, nh;
1477 const t_grpopts* opts;
1478 real reft, kT, ndj, nd;
1479
1480 opts = &(ir->opts); /* just for ease of referencing */
1481 ngtc = ir->opts.ngtc;
1482 nh = state->nhchainlength;
1483
1484 if (ir->eI == eiMD)
1485 {
1486 if (bInit)
1487 {
1488 MassQ->Qinv.resize(ngtc);
1489 }
1490 for (i = 0; (i < ngtc); i++)
1491 {
1492 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
1493 {
1494 MassQ->Qinv[i] = 1.0 / (gmx::square(opts->tau_t[i] / M_2PI) * opts->ref_t[i]);
1495 }
1496 else
1497 {
1498 MassQ->Qinv[i] = 0.0;
1499 }
1500 }
1501 }
1502 else if (EI_VV(ir->eI))
1503 {
1504 /* Set pressure variables */
1505
1506 if (bInit)
1507 {
1508 if (state->vol0 == 0)
1509 {
1510 state->vol0 = det(state->box);
1511 /* because we start by defining a fixed
1512 compressibility, we need the volume at this
1513 compressibility to solve the problem. */
1514 }
1515 }
1516
1517 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
1518 /* Consider evaluating eventually if this the right mass to use. All are correct, some might be more stable */
1519 MassQ->Winv = (PRESFAC * trace(ir->compress) * BOLTZ * opts->ref_t[0])
1520 / (DIM * state->vol0 * gmx::square(ir->tau_p / M_2PI));
1521 /* An alternate mass definition, from Tuckerman et al. */
1522 /* MassQ->Winv = 1.0/(gmx::square(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
1523 for (d = 0; d < DIM; d++)
1524 {
1525 for (n = 0; n < DIM; n++)
1526 {
1527 MassQ->Winvm[d][n] =
1528 PRESFAC * ir->compress[d][n] / (state->vol0 * gmx::square(ir->tau_p / M_2PI));
1529 /* not clear this is correct yet for the anisotropic case. Will need to reevaluate
1530 before using MTTK for anisotropic states.*/
1531 }
1532 }
1533 /* Allocate space for thermostat variables */
1534 if (bInit)
1535 {
1536 MassQ->Qinv.resize(ngtc * nh);
1537 }
1538
1539 /* now, set temperature variables */
1540 for (i = 0; i < ngtc; i++)
1541 {
1542 if (opts->tau_t[i] > 0 && opts->ref_t[i] > 0 && opts->nrdf[i] > 0)
1543 {
1544 reft = std::max<real>(0, opts->ref_t[i]);
1545 nd = opts->nrdf[i];
1546 kT = BOLTZ * reft;
1547 for (j = 0; j < nh; j++)
1548 {
1549 if (j == 0)
1550 {
1551 ndj = nd;
1552 }
1553 else
1554 {
1555 ndj = 1;
1556 }
1557 MassQ->Qinv[i * nh + j] = 1.0 / (gmx::square(opts->tau_t[i] / M_2PI) * ndj * kT);
1558 }
1559 }
1560 else
1561 {
1562 for (j = 0; j < nh; j++)
1563 {
1564 MassQ->Qinv[i * nh + j] = 0.0;
1565 }
1566 }
1567 }
1568 }
1569 }
1570
1571 std::array<std::vector<int>, ettTSEQMAX>
init_npt_vars(const t_inputrec * ir,t_state * state,t_extmass * MassQ,gmx_bool bTrotter)1572 init_npt_vars(const t_inputrec* ir, t_state* state, t_extmass* MassQ, gmx_bool bTrotter)
1573 {
1574 int i, j, nnhpres, nh;
1575 const t_grpopts* opts;
1576 real bmass, qmass, reft, kT;
1577
1578 opts = &(ir->opts); /* just for ease of referencing */
1579 nnhpres = state->nnhpres;
1580 nh = state->nhchainlength;
1581
1582 if (EI_VV(ir->eI) && (ir->epc == epcMTTK) && (ir->etc != etcNOSEHOOVER))
1583 {
1584 gmx_fatal(FARGS,
1585 "Cannot do MTTK pressure coupling without Nose-Hoover temperature control");
1586 }
1587
1588 init_npt_masses(ir, state, MassQ, TRUE);
1589
1590 /* first, initialize clear all the trotter calls */
1591 std::array<std::vector<int>, ettTSEQMAX> trotter_seq;
1592 for (i = 0; i < ettTSEQMAX; i++)
1593 {
1594 trotter_seq[i].resize(NTROTTERPARTS, etrtNONE);
1595 trotter_seq[i][0] = etrtSKIPALL;
1596 }
1597
1598 if (!bTrotter)
1599 {
1600 /* no trotter calls, so we never use the values in the array.
1601 * We access them (so we need to define them, but ignore
1602 * then.*/
1603
1604 return trotter_seq;
1605 }
1606
1607 /* compute the kinetic energy by using the half step velocities or
1608 * the kinetic energies, depending on the order of the trotter calls */
1609
1610 if (ir->eI == eiVV)
1611 {
1612 if (inputrecNptTrotter(ir))
1613 {
1614 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1615 We start with the initial one. */
1616 /* first, a round that estimates veta. */
1617 trotter_seq[0][0] = etrtBAROV;
1618
1619 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1620
1621 /* The first half trotter update */
1622 trotter_seq[2][0] = etrtBAROV;
1623 trotter_seq[2][1] = etrtNHC;
1624 trotter_seq[2][2] = etrtBARONHC;
1625
1626 /* The second half trotter update */
1627 trotter_seq[3][0] = etrtBARONHC;
1628 trotter_seq[3][1] = etrtNHC;
1629 trotter_seq[3][2] = etrtBAROV;
1630
1631 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1632 }
1633 else if (inputrecNvtTrotter(ir))
1634 {
1635 /* This is the easy version - there are only two calls, both the same.
1636 Otherwise, even easier -- no calls */
1637 trotter_seq[2][0] = etrtNHC;
1638 trotter_seq[3][0] = etrtNHC;
1639 }
1640 else if (inputrecNphTrotter(ir))
1641 {
1642 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1643 We start with the initial one. */
1644 /* first, a round that estimates veta. */
1645 trotter_seq[0][0] = etrtBAROV;
1646
1647 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1648
1649 /* The first half trotter update */
1650 trotter_seq[2][0] = etrtBAROV;
1651 trotter_seq[2][1] = etrtBARONHC;
1652
1653 /* The second half trotter update */
1654 trotter_seq[3][0] = etrtBARONHC;
1655 trotter_seq[3][1] = etrtBAROV;
1656
1657 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1658 }
1659 }
1660 else if (ir->eI == eiVVAK)
1661 {
1662 if (inputrecNptTrotter(ir))
1663 {
1664 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1665 We start with the initial one. */
1666 /* first, a round that estimates veta. */
1667 trotter_seq[0][0] = etrtBAROV;
1668
1669 /* The first half trotter update, part 1 -- double update, because it commutes */
1670 trotter_seq[1][0] = etrtNHC;
1671
1672 /* The first half trotter update, part 2 */
1673 trotter_seq[2][0] = etrtBAROV;
1674 trotter_seq[2][1] = etrtBARONHC;
1675
1676 /* The second half trotter update, part 1 */
1677 trotter_seq[3][0] = etrtBARONHC;
1678 trotter_seq[3][1] = etrtBAROV;
1679
1680 /* The second half trotter update */
1681 trotter_seq[4][0] = etrtNHC;
1682 }
1683 else if (inputrecNvtTrotter(ir))
1684 {
1685 /* This is the easy version - there is only one call, both the same.
1686 Otherwise, even easier -- no calls */
1687 trotter_seq[1][0] = etrtNHC;
1688 trotter_seq[4][0] = etrtNHC;
1689 }
1690 else if (inputrecNphTrotter(ir))
1691 {
1692 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1693 We start with the initial one. */
1694 /* first, a round that estimates veta. */
1695 trotter_seq[0][0] = etrtBAROV;
1696
1697 /* The first half trotter update, part 1 -- leave zero */
1698 trotter_seq[1][0] = etrtNHC;
1699
1700 /* The first half trotter update, part 2 */
1701 trotter_seq[2][0] = etrtBAROV;
1702 trotter_seq[2][1] = etrtBARONHC;
1703
1704 /* The second half trotter update, part 1 */
1705 trotter_seq[3][0] = etrtBARONHC;
1706 trotter_seq[3][1] = etrtBAROV;
1707
1708 /* The second half trotter update -- blank for now */
1709 }
1710 }
1711
1712 switch (ir->epct)
1713 {
1714 case epctISOTROPIC:
1715 default: bmass = DIM * DIM; /* recommended mass parameters for isotropic barostat */
1716 }
1717
1718 MassQ->QPinv.resize(nnhpres * opts->nhchainlength);
1719
1720 /* barostat temperature */
1721 if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
1722 {
1723 reft = std::max<real>(0, opts->ref_t[0]);
1724 kT = BOLTZ * reft;
1725 for (i = 0; i < nnhpres; i++)
1726 {
1727 for (j = 0; j < nh; j++)
1728 {
1729 if (j == 0)
1730 {
1731 qmass = bmass;
1732 }
1733 else
1734 {
1735 qmass = 1;
1736 }
1737 MassQ->QPinv[i * opts->nhchainlength + j] =
1738 1.0 / (gmx::square(opts->tau_t[0] / M_2PI) * qmass * kT);
1739 }
1740 }
1741 }
1742 else
1743 {
1744 for (i = 0; i < nnhpres; i++)
1745 {
1746 for (j = 0; j < nh; j++)
1747 {
1748 MassQ->QPinv[i * nh + j] = 0.0;
1749 }
1750 }
1751 }
1752 return trotter_seq;
1753 }
1754
energyNoseHoover(const t_inputrec * ir,const t_state * state,const t_extmass * MassQ)1755 static real energyNoseHoover(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ)
1756 {
1757 real energy = 0;
1758
1759 int nh = state->nhchainlength;
1760
1761 for (int i = 0; i < ir->opts.ngtc; i++)
1762 {
1763 const double* ixi = &state->nosehoover_xi[i * nh];
1764 const double* ivxi = &state->nosehoover_vxi[i * nh];
1765 const double* iQinv = &(MassQ->Qinv[i * nh]);
1766
1767 real nd = ir->opts.nrdf[i];
1768 real reft = std::max<real>(ir->opts.ref_t[i], 0);
1769 real kT = BOLTZ * reft;
1770
1771 if (nd > 0.0)
1772 {
1773 if (inputrecNvtTrotter(ir) || inputrecNptTrotter(ir))
1774 {
1775 /* contribution from the thermal momenta of the NH chain */
1776 for (int j = 0; j < nh; j++)
1777 {
1778 if (iQinv[j] > 0)
1779 {
1780 energy += 0.5 * gmx::square(ivxi[j]) / iQinv[j];
1781 /* contribution from the thermal variable of the NH chain */
1782 real ndj = 0;
1783 if (j == 0)
1784 {
1785 ndj = nd;
1786 }
1787 else
1788 {
1789 ndj = 1;
1790 }
1791 energy += ndj * ixi[j] * kT;
1792 }
1793 }
1794 }
1795 else /* Other non Trotter temperature NH control -- no chains yet. */
1796 {
1797 energy += 0.5 * BOLTZ * nd * gmx::square(ivxi[0]) / iQinv[0];
1798 energy += nd * ixi[0] * kT;
1799 }
1800 }
1801 }
1802
1803 return energy;
1804 }
1805
1806 /* Returns the energy from the barostat thermostat chain */
energyPressureMTTK(const t_inputrec * ir,const t_state * state,const t_extmass * MassQ)1807 static real energyPressureMTTK(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ)
1808 {
1809 real energy = 0;
1810
1811 int nh = state->nhchainlength;
1812
1813 for (int i = 0; i < state->nnhpres; i++)
1814 {
1815 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1816 real reft = std::max<real>(ir->opts.ref_t[0], 0.0); /* using 'System' temperature */
1817 real kT = BOLTZ * reft;
1818
1819 for (int j = 0; j < nh; j++)
1820 {
1821 double iQinv = MassQ->QPinv[i * nh + j];
1822 if (iQinv > 0)
1823 {
1824 energy += 0.5 * gmx::square(state->nhpres_vxi[i * nh + j]) / iQinv;
1825 /* contribution from the thermal variable of the NH chain */
1826 energy += state->nhpres_xi[i * nh + j] * kT;
1827 }
1828 if (debug)
1829 {
1830 fprintf(debug, "P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f", i, j,
1831 state->nhpres_vxi[i * nh + j], state->nhpres_xi[i * nh + j]);
1832 }
1833 }
1834 }
1835
1836 return energy;
1837 }
1838
1839 /* Returns the energy accumulated by the V-rescale or Berendsen thermostat */
energyVrescale(const t_inputrec * ir,const t_state * state)1840 static real energyVrescale(const t_inputrec* ir, const t_state* state)
1841 {
1842 real energy = 0;
1843 for (int i = 0; i < ir->opts.ngtc; i++)
1844 {
1845 energy += state->therm_integral[i];
1846 }
1847
1848 return energy;
1849 }
1850
NPT_energy(const t_inputrec * ir,const t_state * state,const t_extmass * MassQ)1851 real NPT_energy(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ)
1852 {
1853 real energyNPT = 0;
1854
1855 if (ir->epc != epcNO)
1856 {
1857 /* Compute the contribution of the pressure to the conserved quantity*/
1858
1859 real vol = det(state->box);
1860
1861 switch (ir->epc)
1862 {
1863 case epcPARRINELLORAHMAN:
1864 {
1865 /* contribution from the pressure momenta */
1866 tensor invMass;
1867 calcParrinelloRahmanInvMass(ir, state->box, invMass);
1868 for (int d = 0; d < DIM; d++)
1869 {
1870 for (int n = 0; n <= d; n++)
1871 {
1872 if (invMass[d][n] > 0)
1873 {
1874 energyNPT += 0.5 * gmx::square(state->boxv[d][n]) / (invMass[d][n] * PRESFAC);
1875 }
1876 }
1877 }
1878
1879 /* Contribution from the PV term.
1880 * Not that with non-zero off-diagonal reference pressures,
1881 * i.e. applied shear stresses, there are additional terms.
1882 * We don't support this here, since that requires keeping
1883 * track of unwrapped box diagonal elements. This case is
1884 * excluded in integratorHasConservedEnergyQuantity().
1885 */
1886 energyNPT += vol * trace(ir->ref_p) / (DIM * PRESFAC);
1887 break;
1888 }
1889 case epcMTTK:
1890 /* contribution from the pressure momenta */
1891 energyNPT += 0.5 * gmx::square(state->veta) / MassQ->Winv;
1892
1893 /* contribution from the PV term */
1894 energyNPT += vol * trace(ir->ref_p) / (DIM * PRESFAC);
1895
1896 if (ir->epc == epcMTTK)
1897 {
1898 /* contribution from the MTTK chain */
1899 energyNPT += energyPressureMTTK(ir, state, MassQ);
1900 }
1901 break;
1902 case epcBERENDSEN:
1903 case epcCRESCALE: energyNPT += state->baros_integral; break;
1904 default:
1905 GMX_RELEASE_ASSERT(
1906 false,
1907 "Conserved energy quantity for pressure coupling is not handled. A case "
1908 "should be added with either the conserved quantity added or nothing added "
1909 "and an exclusion added to integratorHasConservedEnergyQuantity().");
1910 }
1911 }
1912
1913 switch (ir->etc)
1914 {
1915 case etcNO: break;
1916 case etcVRESCALE:
1917 case etcBERENDSEN: energyNPT += energyVrescale(ir, state); break;
1918 case etcNOSEHOOVER: energyNPT += energyNoseHoover(ir, state, MassQ); break;
1919 case etcANDERSEN:
1920 case etcANDERSENMASSIVE:
1921 // Not supported, excluded in integratorHasConservedEnergyQuantity()
1922 break;
1923 default:
1924 GMX_RELEASE_ASSERT(
1925 false,
1926 "Conserved energy quantity for temperature coupling is not handled. A case "
1927 "should be added with either the conserved quantity added or nothing added and "
1928 "an exclusion added to integratorHasConservedEnergyQuantity().");
1929 }
1930
1931 return energyNPT;
1932 }
1933
1934
vrescale_sumnoises(real nn,gmx::ThreeFry2x64<> * rng,gmx::NormalDistribution<real> * normalDist)1935 static real vrescale_sumnoises(real nn, gmx::ThreeFry2x64<>* rng, gmx::NormalDistribution<real>* normalDist)
1936 {
1937 /*
1938 * Returns the sum of nn independent gaussian noises squared
1939 * (i.e. equivalent to summing the square of the return values
1940 * of nn calls to a normal distribution).
1941 */
1942 const real ndeg_tol = 0.0001;
1943 real r;
1944 gmx::GammaDistribution<real> gammaDist(0.5 * nn, 1.0);
1945
1946 if (nn < 2 + ndeg_tol)
1947 {
1948 int nn_int, i;
1949 real gauss;
1950
1951 nn_int = gmx::roundToInt(nn);
1952
1953 if (nn - nn_int < -ndeg_tol || nn - nn_int > ndeg_tol)
1954 {
1955 gmx_fatal(FARGS,
1956 "The v-rescale thermostat was called with a group with #DOF=%f, but for "
1957 "#DOF<3 only integer #DOF are supported",
1958 nn + 1);
1959 }
1960
1961 r = 0;
1962 for (i = 0; i < nn_int; i++)
1963 {
1964 gauss = (*normalDist)(*rng);
1965 r += gauss * gauss;
1966 }
1967 }
1968 else
1969 {
1970 /* Use a gamma distribution for any real nn > 2 */
1971 r = 2.0 * gammaDist(*rng);
1972 }
1973
1974 return r;
1975 }
1976
vrescale_resamplekin(real kk,real sigma,real ndeg,real taut,int64_t step,int64_t seed)1977 real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut, int64_t step, int64_t seed)
1978 {
1979 /*
1980 * Generates a new value for the kinetic energy,
1981 * according to Bussi et al JCP (2007), Eq. (A7)
1982 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1983 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1984 * ndeg: number of degrees of freedom of the atoms to be thermalized
1985 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1986 */
1987 real factor, rr, ekin_new;
1988 gmx::ThreeFry2x64<64> rng(seed, gmx::RandomDomain::Thermostat);
1989 gmx::NormalDistribution<real> normalDist;
1990
1991 if (taut > 0.1)
1992 {
1993 factor = exp(-1.0 / taut);
1994 }
1995 else
1996 {
1997 factor = 0.0;
1998 }
1999
2000 rng.restart(step, 0);
2001
2002 rr = normalDist(rng);
2003
2004 ekin_new = kk
2005 + (1.0 - factor)
2006 * (sigma * (vrescale_sumnoises(ndeg - 1, &rng, &normalDist) + rr * rr) / ndeg - kk)
2007 + 2.0 * rr * std::sqrt(kk * sigma / ndeg * (1.0 - factor) * factor);
2008
2009 return ekin_new;
2010 }
2011
vrescale_tcoupl(const t_inputrec * ir,int64_t step,gmx_ekindata_t * ekind,real dt,double therm_integral[])2012 void vrescale_tcoupl(const t_inputrec* ir, int64_t step, gmx_ekindata_t* ekind, real dt, double therm_integral[])
2013 {
2014 const t_grpopts* opts;
2015 int i;
2016 real Ek, Ek_ref1, Ek_ref, Ek_new;
2017
2018 opts = &ir->opts;
2019
2020 for (i = 0; (i < opts->ngtc); i++)
2021 {
2022 if (ir->eI == eiVV)
2023 {
2024 Ek = trace(ekind->tcstat[i].ekinf);
2025 }
2026 else
2027 {
2028 Ek = trace(ekind->tcstat[i].ekinh);
2029 }
2030
2031 if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0)
2032 {
2033 Ek_ref1 = 0.5 * opts->ref_t[i] * BOLTZ;
2034 Ek_ref = Ek_ref1 * opts->nrdf[i];
2035
2036 Ek_new = vrescale_resamplekin(Ek, Ek_ref, opts->nrdf[i], opts->tau_t[i] / dt, step, ir->ld_seed);
2037
2038 /* Analytically Ek_new>=0, but we check for rounding errors */
2039 if (Ek_new <= 0)
2040 {
2041 ekind->tcstat[i].lambda = 0.0;
2042 }
2043 else
2044 {
2045 ekind->tcstat[i].lambda = std::sqrt(Ek_new / Ek);
2046 }
2047
2048 therm_integral[i] -= Ek_new - Ek;
2049
2050 if (debug)
2051 {
2052 fprintf(debug, "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n", i, Ek_ref,
2053 Ek, Ek_new, ekind->tcstat[i].lambda);
2054 }
2055 }
2056 else
2057 {
2058 ekind->tcstat[i].lambda = 1.0;
2059 }
2060 }
2061 }
2062
rescale_velocities(const gmx_ekindata_t * ekind,const t_mdatoms * mdatoms,int start,int end,rvec v[])2063 void rescale_velocities(const gmx_ekindata_t* ekind, const t_mdatoms* mdatoms, int start, int end, rvec v[])
2064 {
2065 unsigned short *cACC, *cTC;
2066 int ga, gt, n, d;
2067 real lg;
2068 rvec vrel;
2069
2070 cTC = mdatoms->cTC;
2071
2072 gmx::ArrayRef<const t_grp_tcstat> tcstat = ekind->tcstat;
2073
2074 if (ekind->bNEMD)
2075 {
2076 gmx::ArrayRef<const t_grp_acc> gstat = ekind->grpstat;
2077 cACC = mdatoms->cACC;
2078
2079 ga = 0;
2080 gt = 0;
2081 for (n = start; n < end; n++)
2082 {
2083 if (cACC)
2084 {
2085 ga = cACC[n];
2086 }
2087 if (cTC)
2088 {
2089 gt = cTC[n];
2090 }
2091 /* Only scale the velocity component relative to the COM velocity */
2092 rvec_sub(v[n], gstat[ga].u, vrel);
2093 lg = tcstat[gt].lambda;
2094 for (d = 0; d < DIM; d++)
2095 {
2096 v[n][d] = gstat[ga].u[d] + lg * vrel[d];
2097 }
2098 }
2099 }
2100 else
2101 {
2102 gt = 0;
2103 for (n = start; n < end; n++)
2104 {
2105 if (cTC)
2106 {
2107 gt = cTC[n];
2108 }
2109 lg = tcstat[gt].lambda;
2110 for (d = 0; d < DIM; d++)
2111 {
2112 v[n][d] *= lg;
2113 }
2114 }
2115 }
2116 }
2117
2118 //! Check whether we're doing simulated annealing
doSimulatedAnnealing(const t_inputrec * ir)2119 bool doSimulatedAnnealing(const t_inputrec* ir)
2120 {
2121 for (int i = 0; i < ir->opts.ngtc; i++)
2122 {
2123 /* set bSimAnn if any group is being annealed */
2124 if (ir->opts.annealing[i] != eannNO)
2125 {
2126 return true;
2127 }
2128 }
2129 return false;
2130 }
2131
2132 // TODO If we keep simulated annealing, make a proper module that
2133 // does not rely on changing inputrec.
initSimulatedAnnealing(t_inputrec * ir,gmx::Update * upd)2134 bool initSimulatedAnnealing(t_inputrec* ir, gmx::Update* upd)
2135 {
2136 bool doSimAnnealing = doSimulatedAnnealing(ir);
2137 if (doSimAnnealing)
2138 {
2139 update_annealing_target_temp(ir, ir->init_t, upd);
2140 }
2141 return doSimAnnealing;
2142 }
2143
2144 /* set target temperatures if we are annealing */
update_annealing_target_temp(t_inputrec * ir,real t,gmx::Update * upd)2145 void update_annealing_target_temp(t_inputrec* ir, real t, gmx::Update* upd)
2146 {
2147 int i, j, n, npoints;
2148 real pert, thist = 0, x;
2149
2150 for (i = 0; i < ir->opts.ngtc; i++)
2151 {
2152 npoints = ir->opts.anneal_npoints[i];
2153 switch (ir->opts.annealing[i])
2154 {
2155 case eannNO: continue;
2156 case eannPERIODIC:
2157 /* calculate time modulo the period */
2158 pert = ir->opts.anneal_time[i][npoints - 1];
2159 n = static_cast<int>(t / pert);
2160 thist = t - n * pert; /* modulo time */
2161 /* Make sure rounding didn't get us outside the interval */
2162 if (std::fabs(thist - pert) < GMX_REAL_EPS * 100)
2163 {
2164 thist = 0;
2165 }
2166 break;
2167 case eannSINGLE: thist = t; break;
2168 default:
2169 gmx_fatal(FARGS, "Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)",
2170 i, ir->opts.ngtc, npoints);
2171 }
2172 /* We are doing annealing for this group if we got here,
2173 * and we have the (relative) time as thist.
2174 * calculate target temp */
2175 j = 0;
2176 while ((j < npoints - 1) && (thist > (ir->opts.anneal_time[i][j + 1])))
2177 {
2178 j++;
2179 }
2180 if (j < npoints - 1)
2181 {
2182 /* Found our position between points j and j+1.
2183 * Interpolate: x is the amount from j+1, (1-x) from point j
2184 * First treat possible jumps in temperature as a special case.
2185 */
2186 if ((ir->opts.anneal_time[i][j + 1] - ir->opts.anneal_time[i][j]) < GMX_REAL_EPS * 100)
2187 {
2188 ir->opts.ref_t[i] = ir->opts.anneal_temp[i][j + 1];
2189 }
2190 else
2191 {
2192 x = ((thist - ir->opts.anneal_time[i][j])
2193 / (ir->opts.anneal_time[i][j + 1] - ir->opts.anneal_time[i][j]));
2194 ir->opts.ref_t[i] =
2195 x * ir->opts.anneal_temp[i][j + 1] + (1 - x) * ir->opts.anneal_temp[i][j];
2196 }
2197 }
2198 else
2199 {
2200 ir->opts.ref_t[i] = ir->opts.anneal_temp[i][npoints - 1];
2201 }
2202 }
2203
2204 upd->update_temperature_constants(*ir);
2205 }
2206
pleaseCiteCouplingAlgorithms(FILE * fplog,const t_inputrec & ir)2207 void pleaseCiteCouplingAlgorithms(FILE* fplog, const t_inputrec& ir)
2208 {
2209 if (EI_DYNAMICS(ir.eI))
2210 {
2211 if (ir.etc == etcBERENDSEN)
2212 {
2213 please_cite(fplog, "Berendsen84a");
2214 }
2215 if (ir.etc == etcVRESCALE)
2216 {
2217 please_cite(fplog, "Bussi2007a");
2218 }
2219 if (ir.epc == epcCRESCALE)
2220 {
2221 please_cite(fplog, "Bernetti2020");
2222 }
2223 // TODO this is actually an integrator, not a coupling algorithm
2224 if (ir.eI == eiSD1)
2225 {
2226 please_cite(fplog, "Goga2012");
2227 }
2228 }
2229 }
2230