2 * This file is part of the GROMACS molecular simulation package.
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, 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.
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.
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.
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.
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.
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.
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"
78 #define NTROTTERPARTS 3
80 /* Suzuki-Yoshida Constants, for n=3 and n=5, for symplectic integration */
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 */
85 #define MAX_SUZUKI_YOSHIDA_NUM 5
86 #define SUZUKI_YOSHIDA_NUM 5
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 };
93 static const double* sy_const
[] = { nullptr, sy_const_1
, nullptr, sy_const_3
, nullptr, sy_const_5
};
96 void update_tcouple(int64_t step
,
97 const t_inputrec
* inputrec
,
99 gmx_ekindata_t
* ekind
,
100 const t_extmass
* MassQ
,
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");
110 bool doTemperatureCoupling
= false;
112 // For VV temperature coupling parameters are updated on the current
113 // step, for the others - one step before.
114 if (inputrec
->etc
== etcNO
)
116 doTemperatureCoupling
= false;
118 else if (EI_VV(inputrec
->eI
))
120 doTemperatureCoupling
= do_per_step(step
, inputrec
->nsttcouple
);
124 doTemperatureCoupling
= do_per_step(step
+ inputrec
->nsttcouple
- 1, inputrec
->nsttcouple
);
127 if (doTemperatureCoupling
)
129 real dttc
= inputrec
->nsttcouple
* inputrec
->delta_t
;
131 // TODO: berendsen_tcoupl(...), nosehoover_tcoupl(...) and vrescale_tcoupl(...) update
132 // temperature coupling parameters, which should be reflected in the name of these
134 switch (inputrec
->etc
)
138 berendsen_tcoupl(inputrec
, ekind
, dttc
, state
->therm_integral
);
141 nosehoover_tcoupl(&(inputrec
->opts
), ekind
, dttc
, state
->nosehoover_xi
.data(),
142 state
->nosehoover_vxi
.data(), MassQ
);
145 vrescale_tcoupl(inputrec
, step
, ekind
, dttc
, state
->therm_integral
.data());
148 /* rescale in place here */
149 if (EI_VV(inputrec
->eI
))
151 rescale_velocities(ekind
, md
, 0, md
->homenr
, state
->v
.rvec_array());
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
++)
160 ekind
->tcstat
[i
].lambda
= 1.0;
165 void update_pcouple_before_coordinates(FILE* fplog
,
167 const t_inputrec
* inputrec
,
169 matrix parrinellorahmanMu
,
173 /* Berendsen P-coupling is completely handled after the coordinate update.
174 * Trotter P-coupling is handled by separate calls to trotter_update().
176 if (inputrec
->epc
== epcPARRINELLORAHMAN
177 && do_per_step(step
+ inputrec
->nstpcouple
- 1, inputrec
->nstpcouple
))
179 real dtpc
= inputrec
->nstpcouple
* inputrec
->delta_t
;
181 parrinellorahman_pcoupl(fplog
, step
, inputrec
, dtpc
, state
->pres_prev
, state
->box
,
182 state
->box_rel
, state
->boxv
, M
, parrinellorahmanMu
, bInitStep
);
186 void update_pcouple_after_coordinates(FILE* fplog
,
188 const t_inputrec
* inputrec
,
190 const matrix pressure
,
191 const matrix forceVirial
,
192 const matrix constraintVirial
,
193 matrix pressureCouplingMu
,
196 gmx::BoxDeformation
* boxDeformation
,
197 const bool scaleCoordinates
)
200 int homenr
= md
->homenr
;
202 /* Cast to real for faster code, no loss in precision (see comment above) */
203 real dt
= inputrec
->delta_t
;
206 /* now update boxes */
207 switch (inputrec
->epc
)
211 if (do_per_step(step
, inputrec
->nstpcouple
))
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
);
221 if (do_per_step(step
, inputrec
->nstpcouple
))
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
);
231 case (epcPARRINELLORAHMAN
):
232 if (do_per_step(step
+ inputrec
->nstpcouple
- 1, inputrec
->nstpcouple
))
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.
238 real dtpc
= inputrec
->nstpcouple
* dt
;
239 for (int i
= 0; i
< DIM
; i
++)
241 for (int m
= 0; m
<= i
; m
++)
243 state
->box
[i
][m
] += dtpc
* state
->boxv
[i
][m
];
246 preserve_box_shape(inputrec
, state
->box_rel
, state
->box
);
248 /* Scale the coordinates */
249 if (scaleCoordinates
)
251 auto x
= state
->x
.rvec_array();
252 for (int n
= start
; n
< start
+ homenr
; n
++)
254 tmvmul_ur0(pressureCouplingMu
, x
[n
], x
[n
]);
260 switch (inputrec
->epct
)
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) */
267 msmul(state
->box
, std::exp(state
->veta
* dt
), state
->box
);
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). */
278 msmul(state
->box
, state
->veta
, state
->boxv
);
288 auto localX
= makeArrayRef(state
->x
).subArray(start
, homenr
);
289 boxDeformation
->apply(localX
, state
->box
, step
);
293 extern gmx_bool
update_randomize_velocities(const t_inputrec
* ir
,
297 gmx::ArrayRef
<gmx::RVec
> v
,
298 const gmx::Update
* upd
,
299 const gmx::Constraints
* constr
)
302 real rate
= (ir
->delta_t
) / ir
->opts
.tau_t
[0];
304 if (ir
->etc
== etcANDERSEN
&& constr
!= nullptr)
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. */
313 "Normal Andersen is currently not supported with constraints, use massive "
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
)))
321 andersen_tcoupl(ir
, step
, cr
, md
, v
, rate
, upd
->getAndersenRandomizeGroup(),
322 upd
->getBoltzmanFactor());
329 static const double sy_const[MAX_SUZUKI_YOSHIDA_NUM+1][MAX_SUZUKI_YOSHIDA_NUM+1] = {
333 {0.828981543588751,-0.657963087177502,0.828981543588751},
335 {0.2967324292201065,0.2967324292201065,-0.186929716880426,0.2967324292201065,0.2967324292201065}
338 /* these integration routines are only referenced inside this file */
339 static void NHC_trotter(const t_grpopts
* opts
,
341 const gmx_ekindata_t
* ekind
,
347 const t_extmass
* MassQ
,
348 gmx_bool bEkinAveVel
)
351 /* general routine for both barostat and thermostat nose hoover chains */
354 double Ekin
, Efac
, reft
, kT
, nd
;
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
;
364 mstepsi
= mstepsj
= ns
;
366 /* if scalefac is NULL, we are doing the NHC of the barostat */
369 if (scalefac
== nullptr)
374 for (i
= 0; i
< nvar
; i
++)
377 /* make it easier to iterate by selecting
378 out the sub-array that corresponds to this T group */
382 gmx::ArrayRef
<const double> iQinv
;
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
;
392 iQinv
= gmx::arrayRefFromArray(&MassQ
->Qinv
[i
* nh
], nh
);
393 const t_grp_tcstat
* tcstat
= &ekind
->tcstat
[i
];
395 reft
= std::max
<real
>(0, opts
->ref_t
[i
]);
398 Ekin
= 2 * trace(tcstat
->ekinf
) * tcstat
->ekinscalef_nhc
;
402 Ekin
= 2 * trace(tcstat
->ekinh
) * tcstat
->ekinscaleh_nhc
;
407 for (mi
= 0; mi
< mstepsi
; mi
++)
409 for (mj
= 0; mj
< mstepsj
; mj
++)
411 /* weighting for this step using Suzuki-Yoshida integration - fixed at 5 */
412 dt
= sy_const
[ns
][mj
] * dtfull
/ mstepsi
;
414 /* compute the thermal forces */
415 GQ
[0] = iQinv
[0] * (Ekin
- nd
* kT
);
417 for (j
= 0; j
< nh
- 1; j
++)
419 if (iQinv
[j
+ 1] > 0)
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
);
431 ivxi
[nh
- 1] += 0.25 * dt
* GQ
[nh
- 1];
432 for (j
= nh
- 1; j
> 0; j
--)
434 Efac
= exp(-0.125 * dt
* ivxi
[j
]);
435 ivxi
[j
- 1] = Efac
* (ivxi
[j
- 1] * Efac
+ 0.25 * dt
* GQ
[j
- 1]);
438 Efac
= exp(-0.5 * dt
* ivxi
[0]);
447 Ekin
*= (Efac
* Efac
);
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 . . . */
453 GQ
[0] = iQinv
[0] * (Ekin
- nd
* kT
);
455 /* update thermostat positions */
456 for (j
= 0; j
< nh
; j
++)
458 ixi
[j
] += 0.5 * dt
* ivxi
[j
];
461 for (j
= 0; j
< nh
- 1; j
++)
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)
467 GQ
[j
+ 1] = iQinv
[j
+ 1] * ((gmx::square(ivxi
[j
]) / iQinv
[j
]) - kT
);
474 ivxi
[nh
- 1] += 0.25 * dt
* GQ
[nh
- 1];
481 static void boxv_trotter(const t_inputrec
* ir
,
485 const gmx_ekindata_t
* ekind
,
488 const t_extmass
* MassQ
)
495 tensor ekinmod
, localpres
;
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}
501 if (ir
->epct
== epctSEMIISOTROPIC
)
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. */
514 /* we take the partial pressure tensors, modify the
515 kinetic energy tensor, and recovert to pressure */
517 if (ir
->opts
.nrdf
[0] == 0)
519 gmx_fatal(FARGS
, "Barostat is coupled to a T-group with no degrees of freedom\n");
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? */
528 pscal
= calc_pres(ir
->pbcType
, nwall
, box
, ekinmod
, vir
, localpres
) + pcorr
;
531 GW
= (vol
* (MassQ
->Winv
/ PRESFAC
)) * (DIM
* pscal
- trace(ir
->ref_p
)); /* W is in ps^2 * bar * nm^3 */
533 *veta
+= 0.5 * dt
* GW
;
537 * This file implements temperature and pressure coupling algorithms:
538 * For now only the Weak coupling and the modified weak coupling.
540 * Furthermore computation of pressure and temperature is done here
544 real
calc_pres(PbcType pbcType
, int nwall
, const matrix box
, const tensor ekin
, const tensor vir
, tensor pres
)
549 if (pbcType
== PbcType::No
|| (pbcType
== PbcType::XY
&& nwall
!= 2))
555 /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
556 * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
560 fac
= PRESFAC
* 2.0 / det(box
);
561 for (n
= 0; (n
< DIM
); n
++)
563 for (m
= 0; (m
< DIM
); m
++)
565 pres
[n
][m
] = (ekin
[n
][m
] - vir
[n
][m
]) * fac
;
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
);
577 return trace(pres
) / DIM
;
580 real
calc_temp(real ekin
, real nrdf
)
584 return (2.0 * ekin
) / (nrdf
* BOLTZ
);
592 /*! \brief Sets 1/mass for Parrinello-Rahman in wInv; NOTE: PRESFAC is not included, so not in GROMACS units! */
593 static void calcParrinelloRahmanInvMass(const t_inputrec
* ir
, const matrix box
, tensor wInv
)
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
]);
601 for (int d
= 0; d
< DIM
; d
++)
603 for (int n
= 0; n
< DIM
; n
++)
605 wInv
[d
][n
] = (4 * M_PI
* M_PI
* ir
->compress
[d
][n
]) / (3 * ir
->tau_p
* ir
->tau_p
* maxBoxLength
);
610 void parrinellorahman_pcoupl(FILE* fplog
,
612 const t_inputrec
* ir
,
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.
628 * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
631 * M_nk = (h') * (h' * h + h' h) * h
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:
639 * M_gmx = M_nk' = b * (b * b' + b * b') * b'
641 * where b is the gromacs box matrix.
642 * Our box accelerations are given by
644 * b = vol/W inv(box') * (P-ref_P) (=h')
647 real vol
= box
[XX
][XX
] * box
[YY
][YY
] * box
[ZZ
][ZZ
];
648 real atot
, arel
, change
;
649 tensor invbox
, pdiff
, t1
, t2
;
651 gmx::invertBoxMatrix(box
, invbox
);
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.
660 calcParrinelloRahmanInvMass(ir
, box
, winv
);
662 m_sub(pres
, ir
->ref_p
, pdiff
);
664 if (ir
->epct
== epctSURFACETENSION
)
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.
670 real xy_pressure
= 0.5 * (pres
[XX
][XX
] + pres
[YY
][YY
]);
671 for (int d
= 0; d
< ZZ
; d
++)
673 pdiff
[d
][d
] = (xy_pressure
- (pres
[ZZ
][ZZ
] - ir
->ref_p
[d
][d
] / box
[d
][d
]));
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.
681 for (int d
= 0; d
< DIM
; d
++)
683 for (int n
= 0; n
< d
; n
++)
685 t1
[d
][n
] += t1
[n
][d
];
692 case epctANISOTROPIC
:
693 for (int d
= 0; d
< DIM
; d
++)
695 for (int n
= 0; n
<= d
; n
++)
697 t1
[d
][n
] *= winv
[d
][n
] * vol
;
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
708 for (int d
= 0; d
< DIM
; d
++)
710 for (int n
= 0; n
<= d
; n
++)
712 t1
[d
][n
] = winv
[0][0] * vol
* arel
* box
[d
][n
];
716 case epctSEMIISOTROPIC
:
717 case epctSURFACETENSION
:
718 /* Note the correction to pdiff above for surftens. coupling */
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
++)
727 for (int n
= 0; n
<= d
; n
++)
729 t1
[d
][n
] = winv
[d
][n
] * vol
* arel
* box
[d
][n
];
732 for (int n
= 0; n
< DIM
; n
++)
734 t1
[ZZ
][n
] *= winv
[ZZ
][n
] * vol
;
739 "Parrinello-Rahman pressure coupling type %s "
740 "not supported yet\n",
741 EPCOUPLTYPETYPE(ir
->epct
));
745 for (int d
= 0; d
< DIM
; d
++)
747 for (int n
= 0; n
<= d
; n
++)
749 boxv
[d
][n
] += dt
* t1
[d
][n
];
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.
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
762 change
= std::fabs(dt
* boxv
[d
][n
] / box
[d
][d
]);
764 if (change
> maxchange
)
771 if (maxchange
> 0.01 && 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
));
784 preserve_box_shape(ir
, box_rel
, boxv
);
786 mtmul(boxv
, box
, t1
); /* t1=boxv * b' */
787 mmul(invbox
, t1
, t2
);
788 mtmul(t2
, invbox
, M
);
790 /* Determine the scaling matrix mu for the coordinates */
791 for (int d
= 0; d
< DIM
; d
++)
793 for (int n
= 0; n
<= d
; n
++)
795 t1
[d
][n
] = box
[d
][n
] + dt
* boxv
[d
][n
];
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
);
803 void berendsen_pcoupl(FILE* fplog
,
805 const t_inputrec
* ir
,
809 const matrix force_vir
,
810 const matrix constraint_vir
,
812 double* baros_integral
)
814 real scalar_pressure
, xy_pressure
, p_corr_z
;
818 * Calculate the scaling matrix mu
822 for (int d
= 0; d
< DIM
; d
++)
824 scalar_pressure
+= pres
[d
][d
] / DIM
;
827 xy_pressure
+= pres
[d
][d
] / (DIM
- 1);
830 /* Pressure is now in bar, everywhere. */
831 #define factor(d, m) (ir->compress[d][m] * dt / ir->tau_p)
833 /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
834 * necessary for triclinic scaling
840 for (int d
= 0; d
< DIM
; d
++)
842 mu
[d
][d
] = 1.0 - factor(d
, d
) * (ir
->ref_p
[d
][d
] - scalar_pressure
) / DIM
;
845 case epctSEMIISOTROPIC
:
846 for (int d
= 0; d
< ZZ
; d
++)
848 mu
[d
][d
] = 1.0 - factor(d
, d
) * (ir
->ref_p
[d
][d
] - xy_pressure
) / DIM
;
850 mu
[ZZ
][ZZ
] = 1.0 - factor(ZZ
, ZZ
) * (ir
->ref_p
[ZZ
][ZZ
] - pres
[ZZ
][ZZ
]) / DIM
;
852 case epctANISOTROPIC
:
853 for (int d
= 0; d
< DIM
; d
++)
855 for (int n
= 0; n
< DIM
; n
++)
857 mu
[d
][n
] = (d
== n
? 1.0 : 0.0) - factor(d
, n
) * (ir
->ref_p
[d
][n
] - pres
[d
][n
]) / DIM
;
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
)
866 p_corr_z
= dt
/ ir
->tau_p
* (ir
->ref_p
[ZZ
][ZZ
] - pres
[ZZ
][ZZ
]);
870 /* when the compressibity is zero, set the pressure correction *
871 * in the z-direction to zero to get the correct surface tension */
874 mu
[ZZ
][ZZ
] = 1.0 - ir
->compress
[ZZ
][ZZ
] * p_corr_z
;
875 for (int d
= 0; d
< DIM
- 1; d
++)
879 * (ir
->ref_p
[d
][d
] / (mu
[ZZ
][ZZ
] * box
[ZZ
][ZZ
])
880 - (pres
[ZZ
][ZZ
] + p_corr_z
- xy_pressure
))
885 gmx_fatal(FARGS
, "Berendsen pressure coupling type %s not supported yet\n",
886 EPCOUPLTYPETYPE(ir
->epct
));
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.
892 mu
[YY
][XX
] += mu
[XX
][YY
];
893 mu
[ZZ
][XX
] += mu
[XX
][ZZ
];
894 mu
[ZZ
][YY
] += mu
[YY
][ZZ
];
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.
908 for (int d
= 0; d
< DIM
; d
++)
910 for (int n
= 0; n
<= d
; n
++)
913 2 * (mu
[d
][n
] - (n
== d
? 1 : 0)) * (force_vir
[d
][n
] + constraint_vir
[d
][n
]);
919 pr_rvecs(debug
, 0, "PC: pres ", pres
, 3);
920 pr_rvecs(debug
, 0, "PC: mu ", mu
, 3);
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)
928 "\nStep %s Warning: pressure scaling more than 1%%, "
930 gmx_step_str(step
, buf2
), mu
[XX
][XX
], mu
[YY
][YY
], mu
[ZZ
][ZZ
]);
933 fprintf(fplog
, "%s", buf
);
935 fprintf(stderr
, "%s", buf
);
939 void crescale_pcoupl(FILE* fplog
,
941 const t_inputrec
* ir
,
945 const matrix force_vir
,
946 const matrix constraint_vir
,
948 double* baros_integral
)
951 * Calculate the scaling matrix mu
953 real scalar_pressure
= 0;
954 real xy_pressure
= 0;
955 for (int d
= 0; d
< DIM
; d
++)
957 scalar_pressure
+= pres
[d
][d
] / DIM
;
960 xy_pressure
+= pres
[d
][d
] / (DIM
- 1);
965 gmx::ThreeFry2x64
<64> rng(ir
->ld_seed
, gmx::RandomDomain::Barostat
);
966 gmx::NormalDistribution
<real
> normalDist
;
967 rng
.restart(step
, 0);
969 for (int d
= 0; d
< DIM
; d
++)
975 real kt
= ir
->opts
.ref_t
[0] * BOLTZ
;
984 gauss
= normalDist(rng
);
985 for (int d
= 0; d
< DIM
; d
++)
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
)
993 case epctSEMIISOTROPIC
:
994 gauss
= normalDist(rng
);
995 gauss2
= normalDist(rng
);
996 for (int d
= 0; d
< ZZ
; d
++)
998 const real compressibilityFactor
= ir
->compress
[d
][d
] * dt
/ ir
->tau_p
;
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
);
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
);
1011 case epctSURFACETENSION
:
1012 gauss
= normalDist(rng
);
1013 gauss2
= normalDist(rng
);
1014 for (int d
= 0; d
< ZZ
; d
++)
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
);
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
);
1032 gmx_fatal(FARGS
, "C-rescale pressure coupling type %s not supported yet\n",
1033 EPCOUPLTYPETYPE(ir
->epct
));
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.
1039 mu
[YY
][XX
] += mu
[XX
][YY
];
1040 mu
[ZZ
][XX
] += mu
[XX
][ZZ
];
1041 mu
[ZZ
][YY
] += mu
[YY
][ZZ
];
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.
1055 for (int d
= 0; d
< DIM
; d
++)
1057 for (int n
= 0; n
<= d
; n
++)
1060 2 * (mu
[d
][n
] - (n
== d
? 1 : 0)) * (force_vir
[d
][n
] + constraint_vir
[d
][n
]);
1066 pr_rvecs(debug
, 0, "PC: pres ", pres
, 3);
1067 pr_rvecs(debug
, 0, "PC: mu ", mu
, 3);
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)
1076 "\nStep %s Warning: pressure scaling more than 1%%, "
1078 gmx_step_str(step
, buf2
), mu
[XX
][XX
], mu
[YY
][YY
], mu
[ZZ
][ZZ
]);
1081 fprintf(fplog
, "%s", buf
);
1083 fprintf(stderr
, "%s", buf
);
1087 void crescale_pscale(const t_inputrec
* ir
,
1095 const unsigned short cFREEZE
[],
1097 const bool scaleCoordinates
)
1099 ivec
* nFreeze
= ir
->opts
.nFreeze
;
1100 int nthreads gmx_unused
;
1103 #ifndef __clang_analyzer__
1104 nthreads
= gmx_omp_nthreads_get(emntUpdate
);
1107 gmx::invertBoxMatrix(mu
, inv_mu
);
1109 /* Scale the positions and the velocities */
1110 if (scaleCoordinates
)
1112 #pragma omp parallel for num_threads(nthreads) schedule(static)
1113 for (int n
= start
; n
< start
+ nr_atoms
; n
++)
1115 // Trivial OpenMP region that does not throw
1118 if (cFREEZE
== nullptr)
1127 if (!nFreeze
[g
][XX
])
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
];
1133 if (!nFreeze
[g
][YY
])
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
];
1138 if (!nFreeze
[g
][ZZ
])
1140 x
[n
][ZZ
] = mu
[ZZ
][ZZ
] * x
[n
][ZZ
];
1141 v
[n
][ZZ
] = inv_mu
[ZZ
][ZZ
] * v
[n
][ZZ
];
1145 /* compute final boxlengths */
1146 for (int d
= 0; d
< DIM
; d
++)
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
];
1153 preserve_box_shape(ir
, box_rel
, box
);
1155 /* (un)shifting should NOT be done after this,
1156 * since the box vectors might have changed
1158 inc_nrnb(nrnb
, eNR_PCOUPL
, nr_atoms
);
1161 void berendsen_pscale(const t_inputrec
* ir
,
1168 const unsigned short cFREEZE
[],
1170 const bool scaleCoordinates
)
1172 ivec
* nFreeze
= ir
->opts
.nFreeze
;
1174 int nthreads gmx_unused
;
1176 #ifndef __clang_analyzer__
1177 nthreads
= gmx_omp_nthreads_get(emntUpdate
);
1180 /* Scale the positions */
1181 if (scaleCoordinates
)
1183 #pragma omp parallel for num_threads(nthreads) schedule(static)
1184 for (int n
= start
; n
< start
+ nr_atoms
; n
++)
1186 // Trivial OpenMP region that does not throw
1189 if (cFREEZE
== nullptr)
1198 if (!nFreeze
[g
][XX
])
1200 x
[n
][XX
] = mu
[XX
][XX
] * x
[n
][XX
] + mu
[YY
][XX
] * x
[n
][YY
] + mu
[ZZ
][XX
] * x
[n
][ZZ
];
1202 if (!nFreeze
[g
][YY
])
1204 x
[n
][YY
] = mu
[YY
][YY
] * x
[n
][YY
] + mu
[ZZ
][YY
] * x
[n
][ZZ
];
1206 if (!nFreeze
[g
][ZZ
])
1208 x
[n
][ZZ
] = mu
[ZZ
][ZZ
] * x
[n
][ZZ
];
1212 /* compute final boxlengths */
1213 for (d
= 0; d
< DIM
; d
++)
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
];
1220 preserve_box_shape(ir
, box_rel
, box
);
1222 /* (un)shifting should NOT be done after this,
1223 * since the box vectors might have changed
1225 inc_nrnb(nrnb
, eNR_PCOUPL
, nr_atoms
);
1228 void berendsen_tcoupl(const t_inputrec
* ir
, gmx_ekindata_t
* ekind
, real dt
, std::vector
<double>& therm_integral
)
1230 const t_grpopts
* opts
= &ir
->opts
;
1232 for (int i
= 0; (i
< opts
->ngtc
); i
++)
1238 Ek
= trace(ekind
->tcstat
[i
].ekinf
);
1239 T
= ekind
->tcstat
[i
].T
;
1243 Ek
= trace(ekind
->tcstat
[i
].ekinh
);
1244 T
= ekind
->tcstat
[i
].Th
;
1247 if ((opts
->tau_t
[i
] > 0) && (T
> 0.0))
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);
1255 ekind
->tcstat
[i
].lambda
= 1.0;
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
;
1263 fprintf(debug
, "TC: group %d: T: %g, Lambda: %g\n", i
, T
, ekind
->tcstat
[i
].lambda
);
1268 void andersen_tcoupl(const t_inputrec
* ir
,
1270 const t_commrec
* cr
,
1271 const t_mdatoms
* md
,
1272 gmx::ArrayRef
<gmx::RVec
> v
,
1274 const std::vector
<bool>& randomize
,
1275 gmx::ArrayRef
<const real
> boltzfac
)
1277 const int* gatindex
= (DOMAINDECOMP(cr
) ? cr
->dd
->globalAtomIndices
.data() : nullptr);
1280 gmx::ThreeFry2x64
<0> rng(ir
->andersen_seed
, gmx::RandomDomain::Thermostat
);
1281 gmx::UniformRealDistribution
<real
> uniformDist
;
1282 gmx::TabulatedNormalDistribution
<real
, 14> normalDist
;
1284 /* randomize the velocities of the selected particles */
1286 for (i
= 0; i
< md
->homenr
; i
++) /* now loop over the list of atoms */
1288 int ng
= gatindex
? gatindex
[i
] : i
;
1289 gmx_bool bRandomize
;
1291 rng
.restart(step
, ng
);
1295 gc
= md
->cTC
[i
]; /* assign the atom to a temperature group if there are more than one */
1299 if (ir
->etc
== etcANDERSENMASSIVE
)
1301 /* Randomize particle always */
1306 /* Randomize particle probabilistically */
1307 uniformDist
.reset();
1308 bRandomize
= uniformDist(rng
) < rate
;
1315 scal
= std::sqrt(boltzfac
[gc
] * md
->invmass
[i
]);
1319 for (d
= 0; d
< DIM
; d
++)
1321 v
[i
][d
] = scal
* normalDist(rng
);
1329 void nosehoover_tcoupl(const t_grpopts
* opts
,
1330 const gmx_ekindata_t
* ekind
,
1334 const t_extmass
* MassQ
)
1339 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
1341 for (i
= 0; (i
< opts
->ngtc
); i
++)
1343 reft
= std::max
<real
>(0, opts
->ref_t
[i
]);
1345 vxi
[i
] += dt
* MassQ
->Qinv
[i
] * (ekind
->tcstat
[i
].Th
- reft
);
1346 xi
[i
] += dt
* (oldvxi
+ vxi
[i
]) * 0.5;
1350 void trotter_update(const t_inputrec
* ir
,
1352 gmx_ekindata_t
* ekind
,
1353 const gmx_enerdata_t
* enerd
,
1356 const t_mdatoms
* md
,
1357 const t_extmass
* MassQ
,
1358 gmx::ArrayRef
<std::vector
<int>> trotter_seqlist
,
1362 int n
, i
, d
, ngtc
, gc
= 0, t
;
1363 t_grp_tcstat
* tcstat
;
1364 const t_grpopts
* opts
;
1367 double * scalefac
, dtc
;
1368 rvec sumv
= { 0, 0, 0 };
1371 if (trotter_seqno
<= ettTSEQ2
)
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*/
1382 bCouple
= (ir
->nsttcouple
== 1 || do_per_step(step_eff
+ ir
->nsttcouple
, ir
->nsttcouple
));
1384 const gmx::ArrayRef
<const int> trotter_seq
= trotter_seqlist
[trotter_seqno
];
1386 if ((trotter_seq
[0] == etrtSKIPALL
) || (!bCouple
))
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 */
1394 snew(scalefac
, opts
->ngtc
);
1395 for (i
= 0; i
< ngtc
; i
++)
1399 /* execute the series of trotter updates specified in the trotterpart array */
1401 for (i
= 0; i
< NTROTTERPARTS
; i
++)
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
))
1414 auto v
= makeArrayRef(state
->v
);
1415 switch (trotter_seq
[i
])
1419 boxv_trotter(ir
, &(state
->veta
), dt
, state
->box
, ekind
, vir
,
1420 enerd
->term
[F_PDISPCORR
], MassQ
);
1424 NHC_trotter(opts
, state
->nnhpres
, ekind
, dt
, state
->nhpres_xi
.data(),
1425 state
->nhpres_vxi
.data(), nullptr, &(state
->veta
), MassQ
, FALSE
);
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. */
1435 for (t
= 0; t
< ngtc
; t
++)
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
]);
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? */
1445 /* modify the velocities as well */
1446 for (n
= 0; n
< md
->homenr
; n
++)
1448 if (md
->cTC
) /* does this conditional need to be here? is this always true?*/
1452 for (d
= 0; d
< DIM
; d
++)
1454 v
[n
][d
] *= scalefac
[gc
];
1459 for (d
= 0; d
< DIM
; d
++)
1461 sumv
[d
] += (v
[n
][d
]) / md
->invmass
[n
];
1469 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
1474 extern void init_npt_masses(const t_inputrec
* ir
, t_state
* state
, t_extmass
* MassQ
, gmx_bool bInit
)
1476 int n
, i
, j
, d
, ngtc
, nh
;
1477 const t_grpopts
* opts
;
1478 real reft
, kT
, ndj
, nd
;
1480 opts
= &(ir
->opts
); /* just for ease of referencing */
1481 ngtc
= ir
->opts
.ngtc
;
1482 nh
= state
->nhchainlength
;
1488 MassQ
->Qinv
.resize(ngtc
);
1490 for (i
= 0; (i
< ngtc
); i
++)
1492 if ((opts
->tau_t
[i
] > 0) && (opts
->ref_t
[i
] > 0))
1494 MassQ
->Qinv
[i
] = 1.0 / (gmx::square(opts
->tau_t
[i
] / M_2PI
) * opts
->ref_t
[i
]);
1498 MassQ
->Qinv
[i
] = 0.0;
1502 else if (EI_VV(ir
->eI
))
1504 /* Set pressure variables */
1508 if (state
->vol0
== 0)
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. */
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
++)
1525 for (n
= 0; n
< DIM
; n
++)
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.*/
1533 /* Allocate space for thermostat variables */
1536 MassQ
->Qinv
.resize(ngtc
* nh
);
1539 /* now, set temperature variables */
1540 for (i
= 0; i
< ngtc
; i
++)
1542 if (opts
->tau_t
[i
] > 0 && opts
->ref_t
[i
] > 0 && opts
->nrdf
[i
] > 0)
1544 reft
= std::max
<real
>(0, opts
->ref_t
[i
]);
1547 for (j
= 0; j
< nh
; j
++)
1557 MassQ
->Qinv
[i
* nh
+ j
] = 1.0 / (gmx::square(opts
->tau_t
[i
] / M_2PI
) * ndj
* kT
);
1562 for (j
= 0; j
< nh
; j
++)
1564 MassQ
->Qinv
[i
* nh
+ j
] = 0.0;
1571 std::array
<std::vector
<int>, ettTSEQMAX
>
1572 init_npt_vars(const t_inputrec
* ir
, t_state
* state
, t_extmass
* MassQ
, gmx_bool bTrotter
)
1574 int i
, j
, nnhpres
, nh
;
1575 const t_grpopts
* opts
;
1576 real bmass
, qmass
, reft
, kT
;
1578 opts
= &(ir
->opts
); /* just for ease of referencing */
1579 nnhpres
= state
->nnhpres
;
1580 nh
= state
->nhchainlength
;
1582 if (EI_VV(ir
->eI
) && (ir
->epc
== epcMTTK
) && (ir
->etc
!= etcNOSEHOOVER
))
1585 "Cannot do MTTK pressure coupling without Nose-Hoover temperature control");
1588 init_npt_masses(ir
, state
, MassQ
, TRUE
);
1590 /* first, initialize clear all the trotter calls */
1591 std::array
<std::vector
<int>, ettTSEQMAX
> trotter_seq
;
1592 for (i
= 0; i
< ettTSEQMAX
; i
++)
1594 trotter_seq
[i
].resize(NTROTTERPARTS
, etrtNONE
);
1595 trotter_seq
[i
][0] = etrtSKIPALL
;
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
1607 /* compute the kinetic energy by using the half step velocities or
1608 * the kinetic energies, depending on the order of the trotter calls */
1612 if (inputrecNptTrotter(ir
))
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
;
1619 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1621 /* The first half trotter update */
1622 trotter_seq
[2][0] = etrtBAROV
;
1623 trotter_seq
[2][1] = etrtNHC
;
1624 trotter_seq
[2][2] = etrtBARONHC
;
1626 /* The second half trotter update */
1627 trotter_seq
[3][0] = etrtBARONHC
;
1628 trotter_seq
[3][1] = etrtNHC
;
1629 trotter_seq
[3][2] = etrtBAROV
;
1631 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1633 else if (inputrecNvtTrotter(ir
))
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
;
1640 else if (inputrecNphTrotter(ir
))
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
;
1647 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1649 /* The first half trotter update */
1650 trotter_seq
[2][0] = etrtBAROV
;
1651 trotter_seq
[2][1] = etrtBARONHC
;
1653 /* The second half trotter update */
1654 trotter_seq
[3][0] = etrtBARONHC
;
1655 trotter_seq
[3][1] = etrtBAROV
;
1657 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1660 else if (ir
->eI
== eiVVAK
)
1662 if (inputrecNptTrotter(ir
))
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
;
1669 /* The first half trotter update, part 1 -- double update, because it commutes */
1670 trotter_seq
[1][0] = etrtNHC
;
1672 /* The first half trotter update, part 2 */
1673 trotter_seq
[2][0] = etrtBAROV
;
1674 trotter_seq
[2][1] = etrtBARONHC
;
1676 /* The second half trotter update, part 1 */
1677 trotter_seq
[3][0] = etrtBARONHC
;
1678 trotter_seq
[3][1] = etrtBAROV
;
1680 /* The second half trotter update */
1681 trotter_seq
[4][0] = etrtNHC
;
1683 else if (inputrecNvtTrotter(ir
))
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
;
1690 else if (inputrecNphTrotter(ir
))
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
;
1697 /* The first half trotter update, part 1 -- leave zero */
1698 trotter_seq
[1][0] = etrtNHC
;
1700 /* The first half trotter update, part 2 */
1701 trotter_seq
[2][0] = etrtBAROV
;
1702 trotter_seq
[2][1] = etrtBARONHC
;
1704 /* The second half trotter update, part 1 */
1705 trotter_seq
[3][0] = etrtBARONHC
;
1706 trotter_seq
[3][1] = etrtBAROV
;
1708 /* The second half trotter update -- blank for now */
1715 default: bmass
= DIM
* DIM
; /* recommended mass parameters for isotropic barostat */
1718 MassQ
->QPinv
.resize(nnhpres
* opts
->nhchainlength
);
1720 /* barostat temperature */
1721 if ((ir
->tau_p
> 0) && (opts
->ref_t
[0] > 0))
1723 reft
= std::max
<real
>(0, opts
->ref_t
[0]);
1725 for (i
= 0; i
< nnhpres
; i
++)
1727 for (j
= 0; j
< nh
; j
++)
1737 MassQ
->QPinv
[i
* opts
->nhchainlength
+ j
] =
1738 1.0 / (gmx::square(opts
->tau_t
[0] / M_2PI
) * qmass
* kT
);
1744 for (i
= 0; i
< nnhpres
; i
++)
1746 for (j
= 0; j
< nh
; j
++)
1748 MassQ
->QPinv
[i
* nh
+ j
] = 0.0;
1755 static real
energyNoseHoover(const t_inputrec
* ir
, const t_state
* state
, const t_extmass
* MassQ
)
1759 int nh
= state
->nhchainlength
;
1761 for (int i
= 0; i
< ir
->opts
.ngtc
; i
++)
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
]);
1767 int nd
= static_cast<int>(ir
->opts
.nrdf
[i
]);
1768 real reft
= std::max
<real
>(ir
->opts
.ref_t
[i
], 0);
1769 real kT
= BOLTZ
* reft
;
1773 if (inputrecNvtTrotter(ir
) || inputrecNptTrotter(ir
))
1775 /* contribution from the thermal momenta of the NH chain */
1776 for (int j
= 0; j
< nh
; j
++)
1780 energy
+= 0.5 * gmx::square(ivxi
[j
]) / iQinv
[j
];
1781 /* contribution from the thermal variable of the NH chain */
1791 energy
+= ndj
* ixi
[j
] * kT
;
1795 else /* Other non Trotter temperature NH control -- no chains yet. */
1797 energy
+= 0.5 * BOLTZ
* nd
* gmx::square(ivxi
[0]) / iQinv
[0];
1798 energy
+= nd
* ixi
[0] * kT
;
1806 /* Returns the energy from the barostat thermostat chain */
1807 static real
energyPressureMTTK(const t_inputrec
* ir
, const t_state
* state
, const t_extmass
* MassQ
)
1811 int nh
= state
->nhchainlength
;
1813 for (int i
= 0; i
< state
->nnhpres
; i
++)
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
;
1819 for (int j
= 0; j
< nh
; j
++)
1821 double iQinv
= MassQ
->QPinv
[i
* nh
+ j
];
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
;
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
]);
1839 /* Returns the energy accumulated by the V-rescale or Berendsen thermostat */
1840 static real
energyVrescale(const t_inputrec
* ir
, const t_state
* state
)
1843 for (int i
= 0; i
< ir
->opts
.ngtc
; i
++)
1845 energy
+= state
->therm_integral
[i
];
1851 real
NPT_energy(const t_inputrec
* ir
, const t_state
* state
, const t_extmass
* MassQ
)
1855 if (ir
->epc
!= epcNO
)
1857 /* Compute the contribution of the pressure to the conserved quantity*/
1859 real vol
= det(state
->box
);
1863 case epcPARRINELLORAHMAN
:
1865 /* contribution from the pressure momenta */
1867 calcParrinelloRahmanInvMass(ir
, state
->box
, invMass
);
1868 for (int d
= 0; d
< DIM
; d
++)
1870 for (int n
= 0; n
<= d
; n
++)
1872 if (invMass
[d
][n
] > 0)
1874 energyNPT
+= 0.5 * gmx::square(state
->boxv
[d
][n
]) / (invMass
[d
][n
] * PRESFAC
);
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().
1886 energyNPT
+= vol
* trace(ir
->ref_p
) / (DIM
* PRESFAC
);
1890 /* contribution from the pressure momenta */
1891 energyNPT
+= 0.5 * gmx::square(state
->veta
) / MassQ
->Winv
;
1893 /* contribution from the PV term */
1894 energyNPT
+= vol
* trace(ir
->ref_p
) / (DIM
* PRESFAC
);
1896 if (ir
->epc
== epcMTTK
)
1898 /* contribution from the MTTK chain */
1899 energyNPT
+= energyPressureMTTK(ir
, state
, MassQ
);
1903 case epcCRESCALE
: energyNPT
+= state
->baros_integral
; break;
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().");
1917 case etcBERENDSEN
: energyNPT
+= energyVrescale(ir
, state
); break;
1918 case etcNOSEHOOVER
: energyNPT
+= energyNoseHoover(ir
, state
, MassQ
); break;
1920 case etcANDERSENMASSIVE
:
1921 // Not supported, excluded in integratorHasConservedEnergyQuantity()
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().");
1935 static real
vrescale_sumnoises(real nn
, gmx::ThreeFry2x64
<>* rng
, gmx::NormalDistribution
<real
>* normalDist
)
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).
1942 const real ndeg_tol
= 0.0001;
1944 gmx::GammaDistribution
<real
> gammaDist(0.5 * nn
, 1.0);
1946 if (nn
< 2 + ndeg_tol
)
1951 nn_int
= gmx::roundToInt(nn
);
1953 if (nn
- nn_int
< -ndeg_tol
|| nn
- nn_int
> ndeg_tol
)
1956 "The v-rescale thermostat was called with a group with #DOF=%f, but for "
1957 "#DOF<3 only integer #DOF are supported",
1962 for (i
= 0; i
< nn_int
; i
++)
1964 gauss
= (*normalDist
)(*rng
);
1970 /* Use a gamma distribution for any real nn > 2 */
1971 r
= 2.0 * gammaDist(*rng
);
1977 real
vrescale_resamplekin(real kk
, real sigma
, real ndeg
, real taut
, int64_t step
, int64_t seed
)
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'
1987 real factor
, rr
, ekin_new
;
1988 gmx::ThreeFry2x64
<64> rng(seed
, gmx::RandomDomain::Thermostat
);
1989 gmx::NormalDistribution
<real
> normalDist
;
1993 factor
= exp(-1.0 / taut
);
2000 rng
.restart(step
, 0);
2002 rr
= normalDist(rng
);
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
);
2012 void vrescale_tcoupl(const t_inputrec
* ir
, int64_t step
, gmx_ekindata_t
* ekind
, real dt
, double therm_integral
[])
2014 const t_grpopts
* opts
;
2016 real Ek
, Ek_ref1
, Ek_ref
, Ek_new
;
2020 for (i
= 0; (i
< opts
->ngtc
); i
++)
2024 Ek
= trace(ekind
->tcstat
[i
].ekinf
);
2028 Ek
= trace(ekind
->tcstat
[i
].ekinh
);
2031 if (opts
->tau_t
[i
] >= 0 && opts
->nrdf
[i
] > 0 && Ek
> 0)
2033 Ek_ref1
= 0.5 * opts
->ref_t
[i
] * BOLTZ
;
2034 Ek_ref
= Ek_ref1
* opts
->nrdf
[i
];
2036 Ek_new
= vrescale_resamplekin(Ek
, Ek_ref
, opts
->nrdf
[i
], opts
->tau_t
[i
] / dt
, step
, ir
->ld_seed
);
2038 /* Analytically Ek_new>=0, but we check for rounding errors */
2041 ekind
->tcstat
[i
].lambda
= 0.0;
2045 ekind
->tcstat
[i
].lambda
= std::sqrt(Ek_new
/ Ek
);
2048 therm_integral
[i
] -= Ek_new
- Ek
;
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
);
2058 ekind
->tcstat
[i
].lambda
= 1.0;
2063 void rescale_velocities(const gmx_ekindata_t
* ekind
, const t_mdatoms
* mdatoms
, int start
, int end
, rvec v
[])
2065 unsigned short *cACC
, *cTC
;
2072 gmx::ArrayRef
<const t_grp_tcstat
> tcstat
= ekind
->tcstat
;
2076 gmx::ArrayRef
<const t_grp_acc
> gstat
= ekind
->grpstat
;
2077 cACC
= mdatoms
->cACC
;
2081 for (n
= start
; n
< end
; n
++)
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
++)
2096 v
[n
][d
] = gstat
[ga
].u
[d
] + lg
* vrel
[d
];
2103 for (n
= start
; n
< end
; n
++)
2109 lg
= tcstat
[gt
].lambda
;
2110 for (d
= 0; d
< DIM
; d
++)
2118 //! Check whether we're doing simulated annealing
2119 bool doSimulatedAnnealing(const t_inputrec
* ir
)
2121 for (int i
= 0; i
< ir
->opts
.ngtc
; i
++)
2123 /* set bSimAnn if any group is being annealed */
2124 if (ir
->opts
.annealing
[i
] != eannNO
)
2132 // TODO If we keep simulated annealing, make a proper module that
2133 // does not rely on changing inputrec.
2134 bool initSimulatedAnnealing(t_inputrec
* ir
, gmx::Update
* upd
)
2136 bool doSimAnnealing
= doSimulatedAnnealing(ir
);
2139 update_annealing_target_temp(ir
, ir
->init_t
, upd
);
2141 return doSimAnnealing
;
2144 /* set target temperatures if we are annealing */
2145 void update_annealing_target_temp(t_inputrec
* ir
, real t
, gmx::Update
* upd
)
2147 int i
, j
, n
, npoints
;
2148 real pert
, thist
= 0, x
;
2150 for (i
= 0; i
< ir
->opts
.ngtc
; i
++)
2152 npoints
= ir
->opts
.anneal_npoints
[i
];
2153 switch (ir
->opts
.annealing
[i
])
2155 case eannNO
: continue;
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)
2167 case eannSINGLE
: thist
= t
; break;
2169 gmx_fatal(FARGS
, "Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)",
2170 i
, ir
->opts
.ngtc
, npoints
);
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 */
2176 while ((j
< npoints
- 1) && (thist
> (ir
->opts
.anneal_time
[i
][j
+ 1])))
2180 if (j
< npoints
- 1)
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.
2186 if ((ir
->opts
.anneal_time
[i
][j
+ 1] - ir
->opts
.anneal_time
[i
][j
]) < GMX_REAL_EPS
* 100)
2188 ir
->opts
.ref_t
[i
] = ir
->opts
.anneal_temp
[i
][j
+ 1];
2192 x
= ((thist
- ir
->opts
.anneal_time
[i
][j
])
2193 / (ir
->opts
.anneal_time
[i
][j
+ 1] - ir
->opts
.anneal_time
[i
][j
]));
2195 x
* ir
->opts
.anneal_temp
[i
][j
+ 1] + (1 - x
) * ir
->opts
.anneal_temp
[i
][j
];
2200 ir
->opts
.ref_t
[i
] = ir
->opts
.anneal_temp
[i
][npoints
- 1];
2204 upd
->update_temperature_constants(*ir
);
2207 void pleaseCiteCouplingAlgorithms(FILE* fplog
, const t_inputrec
& ir
)
2209 if (EI_DYNAMICS(ir
.eI
))
2211 if (ir
.etc
== etcBERENDSEN
)
2213 please_cite(fplog
, "Berendsen84a");
2215 if (ir
.etc
== etcVRESCALE
)
2217 please_cite(fplog
, "Bussi2007a");
2219 if (ir
.epc
== epcCRESCALE
)
2221 please_cite(fplog
, "Bernetti2020");
2223 // TODO this is actually an integrator, not a coupling algorithm
2226 please_cite(fplog
, "Goga2012");