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,2018,2019, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
44 #include "gromacs/domdec/domdec_struct.h"
45 #include "gromacs/gmxlib/nrnb.h"
46 #include "gromacs/math/functions.h"
47 #include "gromacs/math/invertmatrix.h"
48 #include "gromacs/math/units.h"
49 #include "gromacs/math/vec.h"
50 #include "gromacs/math/vecdump.h"
51 #include "gromacs/mdlib/expanded.h"
52 #include "gromacs/mdlib/gmx_omp_nthreads.h"
53 #include "gromacs/mdlib/stat.h"
54 #include "gromacs/mdlib/update.h"
55 #include "gromacs/mdtypes/commrec.h"
56 #include "gromacs/mdtypes/enerdata.h"
57 #include "gromacs/mdtypes/group.h"
58 #include "gromacs/mdtypes/inputrec.h"
59 #include "gromacs/mdtypes/md_enums.h"
60 #include "gromacs/mdtypes/mdatom.h"
61 #include "gromacs/mdtypes/state.h"
62 #include "gromacs/pbcutil/boxutilities.h"
63 #include "gromacs/pbcutil/pbc.h"
64 #include "gromacs/random/gammadistribution.h"
65 #include "gromacs/random/normaldistribution.h"
66 #include "gromacs/random/tabulatednormaldistribution.h"
67 #include "gromacs/random/threefry.h"
68 #include "gromacs/random/uniformrealdistribution.h"
69 #include "gromacs/utility/cstringutil.h"
70 #include "gromacs/utility/fatalerror.h"
71 #include "gromacs/utility/pleasecite.h"
72 #include "gromacs/utility/smalloc.h"
74 #define NTROTTERPARTS 3
76 /* Suzuki-Yoshida Constants, for n=3 and n=5, for symplectic integration */
78 /* for n=3, w0 = w2 = 1/(2-2^-(1/3)), w1 = 1-2*w0 */
79 /* for n=5, w0 = w1 = w3 = w4 = 1/(4-4^-(1/3)), w1 = 1-4*w0 */
81 #define MAX_SUZUKI_YOSHIDA_NUM 5
82 #define SUZUKI_YOSHIDA_NUM 5
84 static const double sy_const_1
[] = { 1. };
85 static const double sy_const_3
[] = { 0.828981543588751, -0.657963087177502, 0.828981543588751 };
86 static const double sy_const_5
[] = { 0.2967324292201065, 0.2967324292201065, -0.186929716880426, 0.2967324292201065, 0.2967324292201065 };
88 static const double* sy_const
[] = {
98 static const double sy_const[MAX_SUZUKI_YOSHIDA_NUM+1][MAX_SUZUKI_YOSHIDA_NUM+1] = {
102 {0.828981543588751,-0.657963087177502,0.828981543588751},
104 {0.2967324292201065,0.2967324292201065,-0.186929716880426,0.2967324292201065,0.2967324292201065}
107 /* these integration routines are only referenced inside this file */
108 static void NHC_trotter(const t_grpopts
*opts
, int nvar
, const gmx_ekindata_t
*ekind
, real dtfull
,
109 double xi
[], double vxi
[], double scalefac
[], real
*veta
, const t_extmass
*MassQ
, gmx_bool bEkinAveVel
)
112 /* general routine for both barostat and thermostat nose hoover chains */
115 double Ekin
, Efac
, reft
, kT
, nd
;
120 int mstepsi
, mstepsj
;
121 int ns
= SUZUKI_YOSHIDA_NUM
; /* set the degree of integration in the types/state.h file */
122 int nh
= opts
->nhchainlength
;
125 mstepsi
= mstepsj
= ns
;
127 /* if scalefac is NULL, we are doing the NHC of the barostat */
130 if (scalefac
== nullptr)
135 for (i
= 0; i
< nvar
; i
++)
138 /* make it easier to iterate by selecting
139 out the sub-array that corresponds to this T group */
143 gmx::ArrayRef
<const double> iQinv
;
146 iQinv
= gmx::arrayRefFromArray(&MassQ
->QPinv
[i
*nh
], nh
);
147 nd
= 1.0; /* THIS WILL CHANGE IF NOT ISOTROPIC */
148 reft
= std::max
<real
>(0, opts
->ref_t
[0]);
149 Ekin
= gmx::square(*veta
)/MassQ
->Winv
;
153 iQinv
= gmx::arrayRefFromArray(&MassQ
->Qinv
[i
*nh
], nh
);
154 const t_grp_tcstat
*tcstat
= &ekind
->tcstat
[i
];
156 reft
= std::max
<real
>(0, opts
->ref_t
[i
]);
159 Ekin
= 2*trace(tcstat
->ekinf
)*tcstat
->ekinscalef_nhc
;
163 Ekin
= 2*trace(tcstat
->ekinh
)*tcstat
->ekinscaleh_nhc
;
168 for (mi
= 0; mi
< mstepsi
; mi
++)
170 for (mj
= 0; mj
< mstepsj
; mj
++)
172 /* weighting for this step using Suzuki-Yoshida integration - fixed at 5 */
173 dt
= sy_const
[ns
][mj
] * dtfull
/ mstepsi
;
175 /* compute the thermal forces */
176 GQ
[0] = iQinv
[0]*(Ekin
- nd
*kT
);
178 for (j
= 0; j
< nh
-1; j
++)
182 /* we actually don't need to update here if we save the
183 state of the GQ, but it's easier to just recompute*/
184 GQ
[j
+1] = iQinv
[j
+1]*((gmx::square(ivxi
[j
])/iQinv
[j
])-kT
);
192 ivxi
[nh
-1] += 0.25*dt
*GQ
[nh
-1];
193 for (j
= nh
-1; j
> 0; j
--)
195 Efac
= exp(-0.125*dt
*ivxi
[j
]);
196 ivxi
[j
-1] = Efac
*(ivxi
[j
-1]*Efac
+ 0.25*dt
*GQ
[j
-1]);
199 Efac
= exp(-0.5*dt
*ivxi
[0]);
210 /* Issue - if the KE is an average of the last and the current temperatures, then we might not be
211 able to scale the kinetic energy directly with this factor. Might take more bookkeeping -- have to
212 think about this a bit more . . . */
214 GQ
[0] = iQinv
[0]*(Ekin
- nd
*kT
);
216 /* update thermostat positions */
217 for (j
= 0; j
< nh
; j
++)
219 ixi
[j
] += 0.5*dt
*ivxi
[j
];
222 for (j
= 0; j
< nh
-1; j
++)
224 Efac
= exp(-0.125*dt
*ivxi
[j
+1]);
225 ivxi
[j
] = Efac
*(ivxi
[j
]*Efac
+ 0.25*dt
*GQ
[j
]);
228 GQ
[j
+1] = iQinv
[j
+1]*((gmx::square(ivxi
[j
])/iQinv
[j
])-kT
);
235 ivxi
[nh
-1] += 0.25*dt
*GQ
[nh
-1];
242 static void boxv_trotter(const t_inputrec
*ir
, real
*veta
, real dt
, const tensor box
,
243 const gmx_ekindata_t
*ekind
, const tensor vir
, real pcorr
, const t_extmass
*MassQ
)
250 tensor ekinmod
, localpres
;
252 /* The heat bath is coupled to a separate barostat, the last temperature group. In the
253 2006 Tuckerman et al paper., the order is iL_{T_baro} iL {T_part}
256 if (ir
->epct
== epctSEMIISOTROPIC
)
265 /* eta is in pure units. veta is in units of ps^-1. GW is in
266 units of ps^-2. However, eta has a reference of 1 nm^3, so care must be
267 taken to use only RATIOS of eta in updating the volume. */
269 /* we take the partial pressure tensors, modify the
270 kinetic energy tensor, and recovert to pressure */
272 if (ir
->opts
.nrdf
[0] == 0)
274 gmx_fatal(FARGS
, "Barostat is coupled to a T-group with no degrees of freedom\n");
276 /* alpha factor for phase space volume, then multiply by the ekin scaling factor. */
277 alpha
= 1.0 + DIM
/(static_cast<double>(ir
->opts
.nrdf
[0]));
278 alpha
*= ekind
->tcstat
[0].ekinscalef_nhc
;
279 msmul(ekind
->ekin
, alpha
, ekinmod
);
280 /* for now, we use Elr = 0, because if you want to get it right, you
281 really should be using PME. Maybe print a warning? */
283 pscal
= calc_pres(ir
->ePBC
, nwall
, box
, ekinmod
, vir
, localpres
)+pcorr
;
286 GW
= (vol
*(MassQ
->Winv
/PRESFAC
))*(DIM
*pscal
- trace(ir
->ref_p
)); /* W is in ps^2 * bar * nm^3 */
292 * This file implements temperature and pressure coupling algorithms:
293 * For now only the Weak coupling and the modified weak coupling.
295 * Furthermore computation of pressure and temperature is done here
299 real
calc_pres(int ePBC
, int nwall
, const matrix box
, const tensor ekin
, const tensor vir
,
305 if (ePBC
== epbcNONE
|| (ePBC
== epbcXY
&& nwall
!= 2))
311 /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
312 * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
316 fac
= PRESFAC
*2.0/det(box
);
317 for (n
= 0; (n
< DIM
); n
++)
319 for (m
= 0; (m
< DIM
); m
++)
321 pres
[n
][m
] = (ekin
[n
][m
] - vir
[n
][m
])*fac
;
327 pr_rvecs(debug
, 0, "PC: pres", pres
, DIM
);
328 pr_rvecs(debug
, 0, "PC: ekin", ekin
, DIM
);
329 pr_rvecs(debug
, 0, "PC: vir ", vir
, DIM
);
330 pr_rvecs(debug
, 0, "PC: box ", box
, DIM
);
333 return trace(pres
)/DIM
;
336 real
calc_temp(real ekin
, real nrdf
)
340 return (2.0*ekin
)/(nrdf
*BOLTZ
);
348 /*! \brief Sets 1/mass for Parrinello-Rahman in wInv; NOTE: PRESFAC is not included, so not in GROMACS units! */
349 static void calcParrinelloRahmanInvMass(const t_inputrec
*ir
, const matrix box
,
354 /* TODO: See if we can make the mass independent of the box size */
355 maxBoxLength
= std::max(box
[XX
][XX
], box
[YY
][YY
]);
356 maxBoxLength
= std::max(maxBoxLength
, box
[ZZ
][ZZ
]);
358 for (int d
= 0; d
< DIM
; d
++)
360 for (int n
= 0; n
< DIM
; n
++)
362 wInv
[d
][n
] = (4*M_PI
*M_PI
*ir
->compress
[d
][n
])/(3*ir
->tau_p
*ir
->tau_p
*maxBoxLength
);
367 void parrinellorahman_pcoupl(FILE *fplog
, int64_t step
,
368 const t_inputrec
*ir
, real dt
, const tensor pres
,
369 const tensor box
, tensor box_rel
, tensor boxv
,
370 tensor M
, matrix mu
, gmx_bool bFirstStep
)
372 /* This doesn't do any coordinate updating. It just
373 * integrates the box vector equations from the calculated
374 * acceleration due to pressure difference. We also compute
375 * the tensor M which is used in update to couple the particle
376 * coordinates to the box vectors.
378 * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
381 * M_nk = (h') * (h' * h + h' h) * h
383 * with the dots denoting time derivatives and h is the transformation from
384 * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
385 * This also goes for the pressure and M tensors - they are transposed relative
386 * to ours. Our equation thus becomes:
389 * M_gmx = M_nk' = b * (b * b' + b * b') * b'
391 * where b is the gromacs box matrix.
392 * Our box accelerations are given by
394 * b = vol/W inv(box') * (P-ref_P) (=h')
397 real vol
= box
[XX
][XX
]*box
[YY
][YY
]*box
[ZZ
][ZZ
];
398 real atot
, arel
, change
;
399 tensor invbox
, pdiff
, t1
, t2
;
401 gmx::invertBoxMatrix(box
, invbox
);
405 /* Note that PRESFAC does not occur here.
406 * The pressure and compressibility always occur as a product,
407 * therefore the pressure unit drops out.
410 calcParrinelloRahmanInvMass(ir
, box
, winv
);
412 m_sub(pres
, ir
->ref_p
, pdiff
);
414 if (ir
->epct
== epctSURFACETENSION
)
416 /* Unlike Berendsen coupling it might not be trivial to include a z
417 * pressure correction here? On the other hand we don't scale the
418 * box momentarily, but change accelerations, so it might not be crucial.
420 real xy_pressure
= 0.5*(pres
[XX
][XX
]+pres
[YY
][YY
]);
421 for (int d
= 0; d
< ZZ
; d
++)
423 pdiff
[d
][d
] = (xy_pressure
-(pres
[ZZ
][ZZ
]-ir
->ref_p
[d
][d
]/box
[d
][d
]));
427 tmmul(invbox
, pdiff
, t1
);
428 /* Move the off-diagonal elements of the 'force' to one side to ensure
429 * that we obey the box constraints.
431 for (int d
= 0; d
< DIM
; d
++)
433 for (int n
= 0; n
< d
; n
++)
435 t1
[d
][n
] += t1
[n
][d
];
442 case epctANISOTROPIC
:
443 for (int d
= 0; d
< DIM
; d
++)
445 for (int n
= 0; n
<= d
; n
++)
447 t1
[d
][n
] *= winv
[d
][n
]*vol
;
452 /* calculate total volume acceleration */
453 atot
= box
[XX
][XX
]*box
[YY
][YY
]*t1
[ZZ
][ZZ
]+
454 box
[XX
][XX
]*t1
[YY
][YY
]*box
[ZZ
][ZZ
]+
455 t1
[XX
][XX
]*box
[YY
][YY
]*box
[ZZ
][ZZ
];
457 /* set all RELATIVE box accelerations equal, and maintain total V
459 for (int d
= 0; d
< DIM
; d
++)
461 for (int n
= 0; n
<= d
; n
++)
463 t1
[d
][n
] = winv
[0][0]*vol
*arel
*box
[d
][n
];
467 case epctSEMIISOTROPIC
:
468 case epctSURFACETENSION
:
469 /* Note the correction to pdiff above for surftens. coupling */
471 /* calculate total XY volume acceleration */
472 atot
= box
[XX
][XX
]*t1
[YY
][YY
]+t1
[XX
][XX
]*box
[YY
][YY
];
473 arel
= atot
/(2*box
[XX
][XX
]*box
[YY
][YY
]);
474 /* set RELATIVE XY box accelerations equal, and maintain total V
475 * change speed. Dont change the third box vector accelerations */
476 for (int d
= 0; d
< ZZ
; d
++)
478 for (int n
= 0; n
<= d
; n
++)
480 t1
[d
][n
] = winv
[d
][n
]*vol
*arel
*box
[d
][n
];
483 for (int n
= 0; n
< DIM
; n
++)
485 t1
[ZZ
][n
] *= winv
[ZZ
][n
]*vol
;
489 gmx_fatal(FARGS
, "Parrinello-Rahman pressure coupling type %s "
490 "not supported yet\n", EPCOUPLTYPETYPE(ir
->epct
));
494 for (int d
= 0; d
< DIM
; d
++)
496 for (int n
= 0; n
<= d
; n
++)
498 boxv
[d
][n
] += dt
*t1
[d
][n
];
500 /* We do NOT update the box vectors themselves here, since
501 * we need them for shifting later. It is instead done last
502 * in the update() routine.
505 /* Calculate the change relative to diagonal elements-
506 since it's perfectly ok for the off-diagonal ones to
507 be zero it doesn't make sense to check the change relative
511 change
= std::fabs(dt
*boxv
[d
][n
]/box
[d
][d
]);
513 if (change
> maxchange
)
520 if (maxchange
> 0.01 && fplog
)
524 "\nStep %s Warning: Pressure scaling more than 1%%. "
525 "This may mean your system\n is not yet equilibrated. "
526 "Use of Parrinello-Rahman pressure coupling during\n"
527 "equilibration can lead to simulation instability, "
528 "and is discouraged.\n",
529 gmx_step_str(step
, buf
));
533 preserve_box_shape(ir
, box_rel
, boxv
);
535 mtmul(boxv
, box
, t1
); /* t1=boxv * b' */
536 mmul(invbox
, t1
, t2
);
537 mtmul(t2
, invbox
, M
);
539 /* Determine the scaling matrix mu for the coordinates */
540 for (int d
= 0; d
< DIM
; d
++)
542 for (int n
= 0; n
<= d
; n
++)
544 t1
[d
][n
] = box
[d
][n
] + dt
*boxv
[d
][n
];
547 preserve_box_shape(ir
, box_rel
, t1
);
548 /* t1 is the box at t+dt, determine mu as the relative change */
549 mmul_ur0(invbox
, t1
, mu
);
552 void berendsen_pcoupl(FILE *fplog
, int64_t step
,
553 const t_inputrec
*ir
, real dt
,
554 const tensor pres
, const matrix box
,
555 const matrix force_vir
, const matrix constraint_vir
,
556 matrix mu
, double *baros_integral
)
559 real scalar_pressure
, xy_pressure
, p_corr_z
;
563 * Calculate the scaling matrix mu
567 for (d
= 0; d
< DIM
; d
++)
569 scalar_pressure
+= pres
[d
][d
]/DIM
;
572 xy_pressure
+= pres
[d
][d
]/(DIM
-1);
575 /* Pressure is now in bar, everywhere. */
576 #define factor(d, m) (ir->compress[d][m]*dt/ir->tau_p)
578 /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
579 * necessary for triclinic scaling
585 for (d
= 0; d
< DIM
; d
++)
587 mu
[d
][d
] = 1.0 - factor(d
, d
)*(ir
->ref_p
[d
][d
] - scalar_pressure
) /DIM
;
590 case epctSEMIISOTROPIC
:
591 for (d
= 0; d
< ZZ
; d
++)
593 mu
[d
][d
] = 1.0 - factor(d
, d
)*(ir
->ref_p
[d
][d
]-xy_pressure
)/DIM
;
596 1.0 - factor(ZZ
, ZZ
)*(ir
->ref_p
[ZZ
][ZZ
] - pres
[ZZ
][ZZ
])/DIM
;
598 case epctANISOTROPIC
:
599 for (d
= 0; d
< DIM
; d
++)
601 for (n
= 0; n
< DIM
; n
++)
603 mu
[d
][n
] = (d
== n
? 1.0 : 0.0)
604 -factor(d
, n
)*(ir
->ref_p
[d
][n
] - pres
[d
][n
])/DIM
;
608 case epctSURFACETENSION
:
609 /* ir->ref_p[0/1] is the reference surface-tension times *
610 * the number of surfaces */
611 if (ir
->compress
[ZZ
][ZZ
] != 0.0F
)
613 p_corr_z
= dt
/ir
->tau_p
*(ir
->ref_p
[ZZ
][ZZ
] - pres
[ZZ
][ZZ
]);
617 /* when the compressibity is zero, set the pressure correction *
618 * in the z-direction to zero to get the correct surface tension */
621 mu
[ZZ
][ZZ
] = 1.0 - ir
->compress
[ZZ
][ZZ
]*p_corr_z
;
622 for (d
= 0; d
< DIM
-1; d
++)
624 mu
[d
][d
] = 1.0 + factor(d
, d
)*(ir
->ref_p
[d
][d
]/(mu
[ZZ
][ZZ
]*box
[ZZ
][ZZ
])
625 - (pres
[ZZ
][ZZ
]+p_corr_z
- xy_pressure
))/(DIM
-1);
629 gmx_fatal(FARGS
, "Berendsen pressure coupling type %s not supported yet\n",
630 EPCOUPLTYPETYPE(ir
->epct
));
632 /* To fullfill the orientation restrictions on triclinic boxes
633 * we will set mu_yx, mu_zx and mu_zy to 0 and correct
634 * the other elements of mu to first order.
636 mu
[YY
][XX
] += mu
[XX
][YY
];
637 mu
[ZZ
][XX
] += mu
[XX
][ZZ
];
638 mu
[ZZ
][YY
] += mu
[YY
][ZZ
];
643 /* Keep track of the work the barostat applies on the system.
644 * Without constraints force_vir tells us how Epot changes when scaling.
645 * With constraints constraint_vir gives us the constraint contribution
646 * to both Epot and Ekin. Although we are not scaling velocities, scaling
647 * the coordinates leads to scaling of distances involved in constraints.
648 * This in turn changes the angular momentum (even if the constrained
649 * distances are corrected at the next step). The kinetic component
650 * of the constraint virial captures the angular momentum change.
652 for (int d
= 0; d
< DIM
; d
++)
654 for (int n
= 0; n
<= d
; n
++)
656 *baros_integral
-= 2*(mu
[d
][n
] - (n
== d
? 1 : 0))*(force_vir
[d
][n
] + constraint_vir
[d
][n
]);
662 pr_rvecs(debug
, 0, "PC: pres ", pres
, 3);
663 pr_rvecs(debug
, 0, "PC: mu ", mu
, 3);
666 if (mu
[XX
][XX
] < 0.99 || mu
[XX
][XX
] > 1.01 ||
667 mu
[YY
][YY
] < 0.99 || mu
[YY
][YY
] > 1.01 ||
668 mu
[ZZ
][ZZ
] < 0.99 || mu
[ZZ
][ZZ
] > 1.01)
671 sprintf(buf
, "\nStep %s Warning: pressure scaling more than 1%%, "
673 gmx_step_str(step
, buf2
), mu
[XX
][XX
], mu
[YY
][YY
], mu
[ZZ
][ZZ
]);
676 fprintf(fplog
, "%s", buf
);
678 fprintf(stderr
, "%s", buf
);
682 void berendsen_pscale(const t_inputrec
*ir
, const matrix mu
,
683 matrix box
, matrix box_rel
,
684 int start
, int nr_atoms
,
685 rvec x
[], const unsigned short cFREEZE
[],
688 ivec
*nFreeze
= ir
->opts
.nFreeze
;
690 int nthreads gmx_unused
;
692 #ifndef __clang_analyzer__
693 nthreads
= gmx_omp_nthreads_get(emntUpdate
);
696 /* Scale the positions */
697 #pragma omp parallel for num_threads(nthreads) schedule(static)
698 for (int n
= start
; n
< start
+nr_atoms
; n
++)
700 // Trivial OpenMP region that does not throw
703 if (cFREEZE
== nullptr)
714 x
[n
][XX
] = mu
[XX
][XX
]*x
[n
][XX
]+mu
[YY
][XX
]*x
[n
][YY
]+mu
[ZZ
][XX
]*x
[n
][ZZ
];
718 x
[n
][YY
] = mu
[YY
][YY
]*x
[n
][YY
]+mu
[ZZ
][YY
]*x
[n
][ZZ
];
722 x
[n
][ZZ
] = mu
[ZZ
][ZZ
]*x
[n
][ZZ
];
725 /* compute final boxlengths */
726 for (d
= 0; d
< DIM
; d
++)
728 box
[d
][XX
] = mu
[XX
][XX
]*box
[d
][XX
]+mu
[YY
][XX
]*box
[d
][YY
]+mu
[ZZ
][XX
]*box
[d
][ZZ
];
729 box
[d
][YY
] = mu
[YY
][YY
]*box
[d
][YY
]+mu
[ZZ
][YY
]*box
[d
][ZZ
];
730 box
[d
][ZZ
] = mu
[ZZ
][ZZ
]*box
[d
][ZZ
];
733 preserve_box_shape(ir
, box_rel
, box
);
735 /* (un)shifting should NOT be done after this,
736 * since the box vectors might have changed
738 inc_nrnb(nrnb
, eNR_PCOUPL
, nr_atoms
);
741 void berendsen_tcoupl(const t_inputrec
*ir
, gmx_ekindata_t
*ekind
, real dt
,
742 std::vector
<double> &therm_integral
)
744 const t_grpopts
*opts
= &ir
->opts
;
746 for (int i
= 0; (i
< opts
->ngtc
); i
++)
752 Ek
= trace(ekind
->tcstat
[i
].ekinf
);
753 T
= ekind
->tcstat
[i
].T
;
757 Ek
= trace(ekind
->tcstat
[i
].ekinh
);
758 T
= ekind
->tcstat
[i
].Th
;
761 if ((opts
->tau_t
[i
] > 0) && (T
> 0.0))
763 real reft
= std::max
<real
>(0, opts
->ref_t
[i
]);
764 real lll
= std::sqrt(1.0 + (dt
/opts
->tau_t
[i
])*(reft
/T
-1.0));
765 ekind
->tcstat
[i
].lambda
= std::max
<real
>(std::min
<real
>(lll
, 1.25), 0.8);
769 ekind
->tcstat
[i
].lambda
= 1.0;
772 /* Keep track of the amount of energy we are adding to the system */
773 therm_integral
[i
] -= (gmx::square(ekind
->tcstat
[i
].lambda
) - 1)*Ek
;
777 fprintf(debug
, "TC: group %d: T: %g, Lambda: %g\n",
778 i
, T
, ekind
->tcstat
[i
].lambda
);
783 void andersen_tcoupl(const t_inputrec
*ir
, int64_t step
,
784 const t_commrec
*cr
, const t_mdatoms
*md
,
785 gmx::ArrayRef
<gmx::RVec
> v
,
786 real rate
, const std::vector
<bool> &randomize
,
787 gmx::ArrayRef
<const real
> boltzfac
)
789 const int *gatindex
= (DOMAINDECOMP(cr
) ? cr
->dd
->globalAtomIndices
.data() : nullptr);
792 gmx::ThreeFry2x64
<0> rng(ir
->andersen_seed
, gmx::RandomDomain::Thermostat
);
793 gmx::UniformRealDistribution
<real
> uniformDist
;
794 gmx::TabulatedNormalDistribution
<real
, 14> normalDist
;
796 /* randomize the velocities of the selected particles */
798 for (i
= 0; i
< md
->homenr
; i
++) /* now loop over the list of atoms */
800 int ng
= gatindex
? gatindex
[i
] : i
;
803 rng
.restart(step
, ng
);
807 gc
= md
->cTC
[i
]; /* assign the atom to a temperature group if there are more than one */
811 if (ir
->etc
== etcANDERSENMASSIVE
)
813 /* Randomize particle always */
818 /* Randomize particle probabilistically */
820 bRandomize
= uniformDist(rng
) < rate
;
827 scal
= std::sqrt(boltzfac
[gc
]*md
->invmass
[i
]);
831 for (d
= 0; d
< DIM
; d
++)
833 v
[i
][d
] = scal
*normalDist(rng
);
841 void nosehoover_tcoupl(const t_grpopts
*opts
, const gmx_ekindata_t
*ekind
, real dt
,
842 double xi
[], double vxi
[], const t_extmass
*MassQ
)
847 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
849 for (i
= 0; (i
< opts
->ngtc
); i
++)
851 reft
= std::max
<real
>(0, opts
->ref_t
[i
]);
853 vxi
[i
] += dt
*MassQ
->Qinv
[i
]*(ekind
->tcstat
[i
].Th
- reft
);
854 xi
[i
] += dt
*(oldvxi
+ vxi
[i
])*0.5;
858 void trotter_update(const t_inputrec
*ir
, int64_t step
, gmx_ekindata_t
*ekind
,
859 const gmx_enerdata_t
*enerd
, t_state
*state
,
860 const tensor vir
, const t_mdatoms
*md
,
861 const t_extmass
*MassQ
, gmx::ArrayRef
< std::vector
< int>> trotter_seqlist
,
865 int n
, i
, d
, ngtc
, gc
= 0, t
;
866 t_grp_tcstat
*tcstat
;
867 const t_grpopts
*opts
;
870 double *scalefac
, dtc
;
871 rvec sumv
= {0, 0, 0};
874 if (trotter_seqno
<= ettTSEQ2
)
876 step_eff
= step
-1; /* the velocity verlet calls are actually out of order -- the first half step
877 is actually the last half step from the previous step. Thus the first half step
878 actually corresponds to the n-1 step*/
886 bCouple
= (ir
->nsttcouple
== 1 ||
887 do_per_step(step_eff
+ir
->nsttcouple
, ir
->nsttcouple
));
889 const gmx::ArrayRef
<const int> trotter_seq
= trotter_seqlist
[trotter_seqno
];
891 if ((trotter_seq
[0] == etrtSKIPALL
) || (!bCouple
))
895 dtc
= ir
->nsttcouple
*ir
->delta_t
; /* This is OK for NPT, because nsttcouple == nstpcouple is enforcesd */
896 opts
= &(ir
->opts
); /* just for ease of referencing */
899 snew(scalefac
, opts
->ngtc
);
900 for (i
= 0; i
< ngtc
; i
++)
904 /* execute the series of trotter updates specified in the trotterpart array */
906 for (i
= 0; i
< NTROTTERPARTS
; i
++)
908 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
909 if ((trotter_seq
[i
] == etrtBAROV2
) || (trotter_seq
[i
] == etrtBARONHC2
) || (trotter_seq
[i
] == etrtNHC2
))
918 auto v
= makeArrayRef(state
->v
);
919 switch (trotter_seq
[i
])
923 boxv_trotter(ir
, &(state
->veta
), dt
, state
->box
, ekind
, vir
,
924 enerd
->term
[F_PDISPCORR
], MassQ
);
928 NHC_trotter(opts
, state
->nnhpres
, ekind
, dt
, state
->nhpres_xi
.data(),
929 state
->nhpres_vxi
.data(), nullptr, &(state
->veta
), MassQ
, FALSE
);
933 NHC_trotter(opts
, opts
->ngtc
, ekind
, dt
, state
->nosehoover_xi
.data(),
934 state
->nosehoover_vxi
.data(), scalefac
, nullptr, MassQ
, (ir
->eI
== eiVV
));
935 /* need to rescale the kinetic energies and velocities here. Could
936 scale the velocities later, but we need them scaled in order to
937 produce the correct outputs, so we'll scale them here. */
939 for (t
= 0; t
< ngtc
; t
++)
941 tcstat
= &ekind
->tcstat
[t
];
942 tcstat
->vscale_nhc
= scalefac
[t
];
943 tcstat
->ekinscaleh_nhc
*= (scalefac
[t
]*scalefac
[t
]);
944 tcstat
->ekinscalef_nhc
*= (scalefac
[t
]*scalefac
[t
]);
946 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
947 /* but do we actually need the total? */
949 /* modify the velocities as well */
950 for (n
= 0; n
< md
->homenr
; n
++)
952 if (md
->cTC
) /* does this conditional need to be here? is this always true?*/
956 for (d
= 0; d
< DIM
; d
++)
958 v
[n
][d
] *= scalefac
[gc
];
963 for (d
= 0; d
< DIM
; d
++)
965 sumv
[d
] += (v
[n
][d
])/md
->invmass
[n
];
974 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
979 extern void init_npt_masses(const t_inputrec
*ir
, t_state
*state
, t_extmass
*MassQ
, gmx_bool bInit
)
981 int n
, i
, j
, d
, ngtc
, nh
;
982 const t_grpopts
*opts
;
983 real reft
, kT
, ndj
, nd
;
985 opts
= &(ir
->opts
); /* just for ease of referencing */
986 ngtc
= ir
->opts
.ngtc
;
987 nh
= state
->nhchainlength
;
993 MassQ
->Qinv
.resize(ngtc
);
995 for (i
= 0; (i
< ngtc
); i
++)
997 if ((opts
->tau_t
[i
] > 0) && (opts
->ref_t
[i
] > 0))
999 MassQ
->Qinv
[i
] = 1.0/(gmx::square(opts
->tau_t
[i
]/M_2PI
)*opts
->ref_t
[i
]);
1003 MassQ
->Qinv
[i
] = 0.0;
1007 else if (EI_VV(ir
->eI
))
1009 /* Set pressure variables */
1013 if (state
->vol0
== 0)
1015 state
->vol0
= det(state
->box
);
1016 /* because we start by defining a fixed
1017 compressibility, we need the volume at this
1018 compressibility to solve the problem. */
1022 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
1023 /* Consider evaluating eventually if this the right mass to use. All are correct, some might be more stable */
1024 MassQ
->Winv
= (PRESFAC
*trace(ir
->compress
)*BOLTZ
*opts
->ref_t
[0])/(DIM
*state
->vol0
*gmx::square(ir
->tau_p
/M_2PI
));
1025 /* An alternate mass definition, from Tuckerman et al. */
1026 /* MassQ->Winv = 1.0/(gmx::square(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
1027 for (d
= 0; d
< DIM
; d
++)
1029 for (n
= 0; n
< DIM
; n
++)
1031 MassQ
->Winvm
[d
][n
] = PRESFAC
*ir
->compress
[d
][n
]/(state
->vol0
*gmx::square(ir
->tau_p
/M_2PI
));
1032 /* not clear this is correct yet for the anisotropic case. Will need to reevaluate
1033 before using MTTK for anisotropic states.*/
1036 /* Allocate space for thermostat variables */
1039 MassQ
->Qinv
.resize(ngtc
* nh
);
1042 /* now, set temperature variables */
1043 for (i
= 0; i
< ngtc
; i
++)
1045 if (opts
->tau_t
[i
] > 0 && opts
->ref_t
[i
] > 0 && opts
->nrdf
[i
] > 0)
1047 reft
= std::max
<real
>(0, opts
->ref_t
[i
]);
1050 for (j
= 0; j
< nh
; j
++)
1060 MassQ
->Qinv
[i
*nh
+j
] = 1.0/(gmx::square(opts
->tau_t
[i
]/M_2PI
)*ndj
*kT
);
1065 for (j
= 0; j
< nh
; j
++)
1067 MassQ
->Qinv
[i
*nh
+j
] = 0.0;
1074 std::array
< std::vector
< int>, ettTSEQMAX
> init_npt_vars(const t_inputrec
*ir
, t_state
*state
,
1075 t_extmass
*MassQ
, gmx_bool bTrotter
)
1077 int i
, j
, nnhpres
, nh
;
1078 const t_grpopts
*opts
;
1079 real bmass
, qmass
, reft
, kT
;
1081 opts
= &(ir
->opts
); /* just for ease of referencing */
1082 nnhpres
= state
->nnhpres
;
1083 nh
= state
->nhchainlength
;
1085 if (EI_VV(ir
->eI
) && (ir
->epc
== epcMTTK
) && (ir
->etc
!= etcNOSEHOOVER
))
1087 gmx_fatal(FARGS
, "Cannot do MTTK pressure coupling without Nose-Hoover temperature control");
1090 init_npt_masses(ir
, state
, MassQ
, TRUE
);
1092 /* first, initialize clear all the trotter calls */
1093 std::array
< std::vector
< int>, ettTSEQMAX
> trotter_seq
;
1094 for (i
= 0; i
< ettTSEQMAX
; i
++)
1096 trotter_seq
[i
].resize(NTROTTERPARTS
, etrtNONE
);
1097 trotter_seq
[i
][0] = etrtSKIPALL
;
1102 /* no trotter calls, so we never use the values in the array.
1103 * We access them (so we need to define them, but ignore
1109 /* compute the kinetic energy by using the half step velocities or
1110 * the kinetic energies, depending on the order of the trotter calls */
1114 if (inputrecNptTrotter(ir
))
1116 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1117 We start with the initial one. */
1118 /* first, a round that estimates veta. */
1119 trotter_seq
[0][0] = etrtBAROV
;
1121 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1123 /* The first half trotter update */
1124 trotter_seq
[2][0] = etrtBAROV
;
1125 trotter_seq
[2][1] = etrtNHC
;
1126 trotter_seq
[2][2] = etrtBARONHC
;
1128 /* The second half trotter update */
1129 trotter_seq
[3][0] = etrtBARONHC
;
1130 trotter_seq
[3][1] = etrtNHC
;
1131 trotter_seq
[3][2] = etrtBAROV
;
1133 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1136 else if (inputrecNvtTrotter(ir
))
1138 /* This is the easy version - there are only two calls, both the same.
1139 Otherwise, even easier -- no calls */
1140 trotter_seq
[2][0] = etrtNHC
;
1141 trotter_seq
[3][0] = etrtNHC
;
1143 else if (inputrecNphTrotter(ir
))
1145 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1146 We start with the initial one. */
1147 /* first, a round that estimates veta. */
1148 trotter_seq
[0][0] = etrtBAROV
;
1150 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1152 /* The first half trotter update */
1153 trotter_seq
[2][0] = etrtBAROV
;
1154 trotter_seq
[2][1] = etrtBARONHC
;
1156 /* The second half trotter update */
1157 trotter_seq
[3][0] = etrtBARONHC
;
1158 trotter_seq
[3][1] = etrtBAROV
;
1160 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1163 else if (ir
->eI
== eiVVAK
)
1165 if (inputrecNptTrotter(ir
))
1167 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1168 We start with the initial one. */
1169 /* first, a round that estimates veta. */
1170 trotter_seq
[0][0] = etrtBAROV
;
1172 /* The first half trotter update, part 1 -- double update, because it commutes */
1173 trotter_seq
[1][0] = etrtNHC
;
1175 /* The first half trotter update, part 2 */
1176 trotter_seq
[2][0] = etrtBAROV
;
1177 trotter_seq
[2][1] = etrtBARONHC
;
1179 /* The second half trotter update, part 1 */
1180 trotter_seq
[3][0] = etrtBARONHC
;
1181 trotter_seq
[3][1] = etrtBAROV
;
1183 /* The second half trotter update */
1184 trotter_seq
[4][0] = etrtNHC
;
1186 else if (inputrecNvtTrotter(ir
))
1188 /* This is the easy version - there is only one call, both the same.
1189 Otherwise, even easier -- no calls */
1190 trotter_seq
[1][0] = etrtNHC
;
1191 trotter_seq
[4][0] = etrtNHC
;
1193 else if (inputrecNphTrotter(ir
))
1195 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1196 We start with the initial one. */
1197 /* first, a round that estimates veta. */
1198 trotter_seq
[0][0] = etrtBAROV
;
1200 /* The first half trotter update, part 1 -- leave zero */
1201 trotter_seq
[1][0] = etrtNHC
;
1203 /* The first half trotter update, part 2 */
1204 trotter_seq
[2][0] = etrtBAROV
;
1205 trotter_seq
[2][1] = etrtBARONHC
;
1207 /* The second half trotter update, part 1 */
1208 trotter_seq
[3][0] = etrtBARONHC
;
1209 trotter_seq
[3][1] = etrtBAROV
;
1211 /* The second half trotter update -- blank for now */
1219 bmass
= DIM
*DIM
; /* recommended mass parameters for isotropic barostat */
1222 MassQ
->QPinv
.resize(nnhpres
*opts
->nhchainlength
);
1224 /* barostat temperature */
1225 if ((ir
->tau_p
> 0) && (opts
->ref_t
[0] > 0))
1227 reft
= std::max
<real
>(0, opts
->ref_t
[0]);
1229 for (i
= 0; i
< nnhpres
; i
++)
1231 for (j
= 0; j
< nh
; j
++)
1241 MassQ
->QPinv
[i
*opts
->nhchainlength
+j
] = 1.0/(gmx::square(opts
->tau_t
[0]/M_2PI
)*qmass
*kT
);
1247 for (i
= 0; i
< nnhpres
; i
++)
1249 for (j
= 0; j
< nh
; j
++)
1251 MassQ
->QPinv
[i
*nh
+j
] = 0.0;
1258 static real
energyNoseHoover(const t_inputrec
*ir
, const t_state
*state
, const t_extmass
*MassQ
)
1262 int nh
= state
->nhchainlength
;
1264 for (int i
= 0; i
< ir
->opts
.ngtc
; i
++)
1266 const double *ixi
= &state
->nosehoover_xi
[i
*nh
];
1267 const double *ivxi
= &state
->nosehoover_vxi
[i
*nh
];
1268 const double *iQinv
= &(MassQ
->Qinv
[i
*nh
]);
1270 int nd
= static_cast<int>(ir
->opts
.nrdf
[i
]);
1271 real reft
= std::max
<real
>(ir
->opts
.ref_t
[i
], 0);
1272 real kT
= BOLTZ
* reft
;
1276 if (inputrecNvtTrotter(ir
))
1278 /* contribution from the thermal momenta of the NH chain */
1279 for (int j
= 0; j
< nh
; j
++)
1283 energy
+= 0.5*gmx::square(ivxi
[j
])/iQinv
[j
];
1284 /* contribution from the thermal variable of the NH chain */
1294 energy
+= ndj
*ixi
[j
]*kT
;
1298 else /* Other non Trotter temperature NH control -- no chains yet. */
1300 energy
+= 0.5*BOLTZ
*nd
*gmx::square(ivxi
[0])/iQinv
[0];
1301 energy
+= nd
*ixi
[0]*kT
;
1309 /* Returns the energy from the barostat thermostat chain */
1310 static real
energyPressureMTTK(const t_inputrec
*ir
, const t_state
*state
, const t_extmass
*MassQ
)
1314 int nh
= state
->nhchainlength
;
1316 for (int i
= 0; i
< state
->nnhpres
; i
++)
1318 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1319 real reft
= std::max
<real
>(ir
->opts
.ref_t
[0], 0.0); /* using 'System' temperature */
1320 real kT
= BOLTZ
* reft
;
1322 for (int j
= 0; j
< nh
; j
++)
1324 double iQinv
= MassQ
->QPinv
[i
*nh
+ j
];
1327 energy
+= 0.5*gmx::square(state
->nhpres_vxi
[i
*nh
+ j
]/iQinv
);
1328 /* contribution from the thermal variable of the NH chain */
1329 energy
+= state
->nhpres_xi
[i
*nh
+ j
]*kT
;
1333 fprintf(debug
, "P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f", i
, j
, state
->nhpres_vxi
[i
*nh
+ j
], state
->nhpres_xi
[i
*nh
+ j
]);
1341 /* Returns the energy accumulated by the V-rescale or Berendsen thermostat */
1342 static real
energyVrescale(const t_inputrec
*ir
, const t_state
*state
)
1345 for (int i
= 0; i
< ir
->opts
.ngtc
; i
++)
1347 energy
+= state
->therm_integral
[i
];
1353 real
NPT_energy(const t_inputrec
*ir
, const t_state
*state
, const t_extmass
*MassQ
)
1357 if (ir
->epc
!= epcNO
)
1359 /* Compute the contribution of the pressure to the conserved quantity*/
1361 real vol
= det(state
->box
);
1365 case epcPARRINELLORAHMAN
:
1367 /* contribution from the pressure momenta */
1369 calcParrinelloRahmanInvMass(ir
, state
->box
, invMass
);
1370 for (int d
= 0; d
< DIM
; d
++)
1372 for (int n
= 0; n
<= d
; n
++)
1374 if (invMass
[d
][n
] > 0)
1376 energyNPT
+= 0.5*gmx::square(state
->boxv
[d
][n
])/(invMass
[d
][n
]*PRESFAC
);
1381 /* Contribution from the PV term.
1382 * Not that with non-zero off-diagonal reference pressures,
1383 * i.e. applied shear stresses, there are additional terms.
1384 * We don't support this here, since that requires keeping
1385 * track of unwrapped box diagonal elements. This case is
1386 * excluded in integratorHasConservedEnergyQuantity().
1388 energyNPT
+= vol
*trace(ir
->ref_p
)/(DIM
*PRESFAC
);
1392 /* contribution from the pressure momenta */
1393 energyNPT
+= 0.5*gmx::square(state
->veta
)/MassQ
->Winv
;
1395 /* contribution from the PV term */
1396 energyNPT
+= vol
*trace(ir
->ref_p
)/(DIM
*PRESFAC
);
1398 if (ir
->epc
== epcMTTK
)
1400 /* contribution from the MTTK chain */
1401 energyNPT
+= energyPressureMTTK(ir
, state
, MassQ
);
1405 energyNPT
+= state
->baros_integral
;
1408 GMX_RELEASE_ASSERT(false, "Conserved energy quantity for pressure coupling is not handled. A case should be added with either the conserved quantity added or nothing added and an exclusion added to integratorHasConservedEnergyQuantity().");
1418 energyNPT
+= energyVrescale(ir
, state
);
1421 energyNPT
+= energyNoseHoover(ir
, state
, MassQ
);
1424 case etcANDERSENMASSIVE
:
1425 // Not supported, excluded in integratorHasConservedEnergyQuantity()
1428 GMX_RELEASE_ASSERT(false, "Conserved energy quantity for temperature coupling is not handled. A case should be added with either the conserved quantity added or nothing added and an exclusion added to integratorHasConservedEnergyQuantity().");
1435 static real
vrescale_sumnoises(real nn
,
1436 gmx::ThreeFry2x64
<> *rng
,
1437 gmx::NormalDistribution
<real
> *normalDist
)
1440 * Returns the sum of nn independent gaussian noises squared
1441 * (i.e. equivalent to summing the square of the return values
1442 * of nn calls to a normal distribution).
1444 const real ndeg_tol
= 0.0001;
1446 gmx::GammaDistribution
<real
> gammaDist(0.5*nn
, 1.0);
1448 if (nn
< 2 + ndeg_tol
)
1453 nn_int
= gmx::roundToInt(nn
);
1455 if (nn
- nn_int
< -ndeg_tol
|| nn
- nn_int
> ndeg_tol
)
1457 gmx_fatal(FARGS
, "The v-rescale thermostat was called with a group with #DOF=%f, but for #DOF<3 only integer #DOF are supported", nn
+ 1);
1461 for (i
= 0; i
< nn_int
; i
++)
1463 gauss
= (*normalDist
)(*rng
);
1469 /* Use a gamma distribution for any real nn > 2 */
1470 r
= 2.0*gammaDist(*rng
);
1476 static real
vrescale_resamplekin(real kk
, real sigma
, real ndeg
, real taut
,
1477 int64_t step
, int64_t seed
)
1480 * Generates a new value for the kinetic energy,
1481 * according to Bussi et al JCP (2007), Eq. (A7)
1482 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1483 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1484 * ndeg: number of degrees of freedom of the atoms to be thermalized
1485 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1487 real factor
, rr
, ekin_new
;
1488 gmx::ThreeFry2x64
<64> rng(seed
, gmx::RandomDomain::Thermostat
);
1489 gmx::NormalDistribution
<real
> normalDist
;
1493 factor
= exp(-1.0/taut
);
1500 rng
.restart(step
, 0);
1502 rr
= normalDist(rng
);
1506 (1.0 - factor
)*(sigma
*(vrescale_sumnoises(ndeg
-1, &rng
, &normalDist
) + rr
*rr
)/ndeg
- kk
) +
1507 2.0*rr
*std::sqrt(kk
*sigma
/ndeg
*(1.0 - factor
)*factor
);
1512 void vrescale_tcoupl(const t_inputrec
*ir
, int64_t step
,
1513 gmx_ekindata_t
*ekind
, real dt
,
1514 double therm_integral
[])
1516 const t_grpopts
*opts
;
1518 real Ek
, Ek_ref1
, Ek_ref
, Ek_new
;
1522 for (i
= 0; (i
< opts
->ngtc
); i
++)
1526 Ek
= trace(ekind
->tcstat
[i
].ekinf
);
1530 Ek
= trace(ekind
->tcstat
[i
].ekinh
);
1533 if (opts
->tau_t
[i
] >= 0 && opts
->nrdf
[i
] > 0 && Ek
> 0)
1535 Ek_ref1
= 0.5*opts
->ref_t
[i
]*BOLTZ
;
1536 Ek_ref
= Ek_ref1
*opts
->nrdf
[i
];
1538 Ek_new
= vrescale_resamplekin(Ek
, Ek_ref
, opts
->nrdf
[i
],
1542 /* Analytically Ek_new>=0, but we check for rounding errors */
1545 ekind
->tcstat
[i
].lambda
= 0.0;
1549 ekind
->tcstat
[i
].lambda
= std::sqrt(Ek_new
/Ek
);
1552 therm_integral
[i
] -= Ek_new
- Ek
;
1556 fprintf(debug
, "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
1557 i
, Ek_ref
, Ek
, Ek_new
, ekind
->tcstat
[i
].lambda
);
1562 ekind
->tcstat
[i
].lambda
= 1.0;
1567 void rescale_velocities(const gmx_ekindata_t
*ekind
, const t_mdatoms
*mdatoms
,
1568 int start
, int end
, rvec v
[])
1570 unsigned short *cACC
, *cTC
;
1577 gmx::ArrayRef
<const t_grp_tcstat
> tcstat
= ekind
->tcstat
;
1581 gmx::ArrayRef
<const t_grp_acc
> gstat
= ekind
->grpstat
;
1582 cACC
= mdatoms
->cACC
;
1586 for (n
= start
; n
< end
; n
++)
1596 /* Only scale the velocity component relative to the COM velocity */
1597 rvec_sub(v
[n
], gstat
[ga
].u
, vrel
);
1598 lg
= tcstat
[gt
].lambda
;
1599 for (d
= 0; d
< DIM
; d
++)
1601 v
[n
][d
] = gstat
[ga
].u
[d
] + lg
*vrel
[d
];
1608 for (n
= start
; n
< end
; n
++)
1614 lg
= tcstat
[gt
].lambda
;
1615 for (d
= 0; d
< DIM
; d
++)
1623 // TODO If we keep simulated annealing, make a proper module that
1624 // does not rely on changing inputrec.
1625 bool initSimulatedAnnealing(t_inputrec
*ir
,
1628 bool doSimulatedAnnealing
= false;
1629 for (int i
= 0; i
< ir
->opts
.ngtc
; i
++)
1631 /* set bSimAnn if any group is being annealed */
1632 if (ir
->opts
.annealing
[i
] != eannNO
)
1634 doSimulatedAnnealing
= true;
1637 if (doSimulatedAnnealing
)
1639 update_annealing_target_temp(ir
, ir
->init_t
, upd
);
1641 return doSimulatedAnnealing
;
1644 /* set target temperatures if we are annealing */
1645 void update_annealing_target_temp(t_inputrec
*ir
, real t
, gmx::Update
*upd
)
1647 int i
, j
, n
, npoints
;
1648 real pert
, thist
= 0, x
;
1650 for (i
= 0; i
< ir
->opts
.ngtc
; i
++)
1652 npoints
= ir
->opts
.anneal_npoints
[i
];
1653 switch (ir
->opts
.annealing
[i
])
1658 /* calculate time modulo the period */
1659 pert
= ir
->opts
.anneal_time
[i
][npoints
-1];
1660 n
= static_cast<int>(t
/ pert
);
1661 thist
= t
- n
*pert
; /* modulo time */
1662 /* Make sure rounding didn't get us outside the interval */
1663 if (std::fabs(thist
-pert
) < GMX_REAL_EPS
*100)
1672 gmx_fatal(FARGS
, "Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)", i
, ir
->opts
.ngtc
, npoints
);
1674 /* We are doing annealing for this group if we got here,
1675 * and we have the (relative) time as thist.
1676 * calculate target temp */
1678 while ((j
< npoints
-1) && (thist
> (ir
->opts
.anneal_time
[i
][j
+1])))
1684 /* Found our position between points j and j+1.
1685 * Interpolate: x is the amount from j+1, (1-x) from point j
1686 * First treat possible jumps in temperature as a special case.
1688 if ((ir
->opts
.anneal_time
[i
][j
+1]-ir
->opts
.anneal_time
[i
][j
]) < GMX_REAL_EPS
*100)
1690 ir
->opts
.ref_t
[i
] = ir
->opts
.anneal_temp
[i
][j
+1];
1694 x
= ((thist
-ir
->opts
.anneal_time
[i
][j
])/
1695 (ir
->opts
.anneal_time
[i
][j
+1]-ir
->opts
.anneal_time
[i
][j
]));
1696 ir
->opts
.ref_t
[i
] = x
*ir
->opts
.anneal_temp
[i
][j
+1]+(1-x
)*ir
->opts
.anneal_temp
[i
][j
];
1701 ir
->opts
.ref_t
[i
] = ir
->opts
.anneal_temp
[i
][npoints
-1];
1705 update_temperature_constants(upd
->sd(), ir
);
1708 void pleaseCiteCouplingAlgorithms(FILE *fplog
,
1709 const t_inputrec
&ir
)
1711 if (EI_DYNAMICS(ir
.eI
))
1713 if (ir
.etc
== etcBERENDSEN
)
1715 please_cite(fplog
, "Berendsen84a");
1717 if (ir
.etc
== etcVRESCALE
)
1719 please_cite(fplog
, "Bussi2007a");
1721 // TODO this is actually an integrator, not a coupling algorithm
1724 please_cite(fplog
, "Goga2012");