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, 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.
43 #include "gromacs/domdec/domdec_struct.h"
44 #include "gromacs/gmxlib/nrnb.h"
45 #include "gromacs/math/functions.h"
46 #include "gromacs/math/invertmatrix.h"
47 #include "gromacs/math/units.h"
48 #include "gromacs/math/vec.h"
49 #include "gromacs/math/vecdump.h"
50 #include "gromacs/mdlib/expanded.h"
51 #include "gromacs/mdlib/gmx_omp_nthreads.h"
52 #include "gromacs/mdlib/sim_util.h"
53 #include "gromacs/mdlib/update.h"
54 #include "gromacs/mdtypes/commrec.h"
55 #include "gromacs/mdtypes/group.h"
56 #include "gromacs/mdtypes/inputrec.h"
57 #include "gromacs/mdtypes/md_enums.h"
58 #include "gromacs/mdtypes/state.h"
59 #include "gromacs/pbcutil/boxutilities.h"
60 #include "gromacs/pbcutil/pbc.h"
61 #include "gromacs/random/gammadistribution.h"
62 #include "gromacs/random/normaldistribution.h"
63 #include "gromacs/random/tabulatednormaldistribution.h"
64 #include "gromacs/random/threefry.h"
65 #include "gromacs/random/uniformrealdistribution.h"
66 #include "gromacs/trajectory/energy.h"
67 #include "gromacs/utility/cstringutil.h"
68 #include "gromacs/utility/fatalerror.h"
69 #include "gromacs/utility/smalloc.h"
71 #define NTROTTERPARTS 3
73 /* Suzuki-Yoshida Constants, for n=3 and n=5, for symplectic integration */
75 /* for n=3, w0 = w2 = 1/(2-2^-(1/3)), w1 = 1-2*w0 */
76 /* for n=5, w0 = w1 = w3 = w4 = 1/(4-4^-(1/3)), w1 = 1-4*w0 */
78 #define MAX_SUZUKI_YOSHIDA_NUM 5
79 #define SUZUKI_YOSHIDA_NUM 5
81 static const double sy_const_1
[] = { 1. };
82 static const double sy_const_3
[] = { 0.828981543588751, -0.657963087177502, 0.828981543588751 };
83 static const double sy_const_5
[] = { 0.2967324292201065, 0.2967324292201065, -0.186929716880426, 0.2967324292201065, 0.2967324292201065 };
85 static const double* sy_const
[] = {
95 static const double sy_const[MAX_SUZUKI_YOSHIDA_NUM+1][MAX_SUZUKI_YOSHIDA_NUM+1] = {
99 {0.828981543588751,-0.657963087177502,0.828981543588751},
101 {0.2967324292201065,0.2967324292201065,-0.186929716880426,0.2967324292201065,0.2967324292201065}
104 /* these integration routines are only referenced inside this file */
105 static void NHC_trotter(t_grpopts
*opts
, int nvar
, gmx_ekindata_t
*ekind
, real dtfull
,
106 double xi
[], double vxi
[], double scalefac
[], real
*veta
, t_extmass
*MassQ
, gmx_bool bEkinAveVel
)
109 /* general routine for both barostat and thermostat nose hoover chains */
112 double Ekin
, Efac
, reft
, kT
, nd
;
114 t_grp_tcstat
*tcstat
;
119 int mstepsi
, mstepsj
;
120 int ns
= SUZUKI_YOSHIDA_NUM
; /* set the degree of integration in the types/state.h file */
121 int nh
= opts
->nhchainlength
;
124 mstepsi
= mstepsj
= ns
;
126 /* if scalefac is NULL, we are doing the NHC of the barostat */
129 if (scalefac
== nullptr)
134 for (i
= 0; i
< nvar
; i
++)
137 /* make it easier to iterate by selecting
138 out the sub-array that corresponds to this T group */
144 iQinv
= &(MassQ
->QPinv
[i
*nh
]);
145 nd
= 1.0; /* THIS WILL CHANGE IF NOT ISOTROPIC */
146 reft
= std::max
<real
>(0, opts
->ref_t
[0]);
147 Ekin
= gmx::square(*veta
)/MassQ
->Winv
;
151 iQinv
= &(MassQ
->Qinv
[i
*nh
]);
152 tcstat
= &ekind
->tcstat
[i
];
154 reft
= std::max
<real
>(0, opts
->ref_t
[i
]);
157 Ekin
= 2*trace(tcstat
->ekinf
)*tcstat
->ekinscalef_nhc
;
161 Ekin
= 2*trace(tcstat
->ekinh
)*tcstat
->ekinscaleh_nhc
;
166 for (mi
= 0; mi
< mstepsi
; mi
++)
168 for (mj
= 0; mj
< mstepsj
; mj
++)
170 /* weighting for this step using Suzuki-Yoshida integration - fixed at 5 */
171 dt
= sy_const
[ns
][mj
] * dtfull
/ mstepsi
;
173 /* compute the thermal forces */
174 GQ
[0] = iQinv
[0]*(Ekin
- nd
*kT
);
176 for (j
= 0; j
< nh
-1; j
++)
180 /* we actually don't need to update here if we save the
181 state of the GQ, but it's easier to just recompute*/
182 GQ
[j
+1] = iQinv
[j
+1]*((gmx::square(ivxi
[j
])/iQinv
[j
])-kT
);
190 ivxi
[nh
-1] += 0.25*dt
*GQ
[nh
-1];
191 for (j
= nh
-1; j
> 0; j
--)
193 Efac
= exp(-0.125*dt
*ivxi
[j
]);
194 ivxi
[j
-1] = Efac
*(ivxi
[j
-1]*Efac
+ 0.25*dt
*GQ
[j
-1]);
197 Efac
= exp(-0.5*dt
*ivxi
[0]);
208 /* Issue - if the KE is an average of the last and the current temperatures, then we might not be
209 able to scale the kinetic energy directly with this factor. Might take more bookkeeping -- have to
210 think about this a bit more . . . */
212 GQ
[0] = iQinv
[0]*(Ekin
- nd
*kT
);
214 /* update thermostat positions */
215 for (j
= 0; j
< nh
; j
++)
217 ixi
[j
] += 0.5*dt
*ivxi
[j
];
220 for (j
= 0; j
< nh
-1; j
++)
222 Efac
= exp(-0.125*dt
*ivxi
[j
+1]);
223 ivxi
[j
] = Efac
*(ivxi
[j
]*Efac
+ 0.25*dt
*GQ
[j
]);
226 GQ
[j
+1] = iQinv
[j
+1]*((gmx::square(ivxi
[j
])/iQinv
[j
])-kT
);
233 ivxi
[nh
-1] += 0.25*dt
*GQ
[nh
-1];
240 static void boxv_trotter(t_inputrec
*ir
, real
*veta
, real dt
, tensor box
,
241 gmx_ekindata_t
*ekind
, tensor vir
, real pcorr
, t_extmass
*MassQ
)
248 tensor ekinmod
, localpres
;
250 /* The heat bath is coupled to a separate barostat, the last temperature group. In the
251 2006 Tuckerman et al paper., the order is iL_{T_baro} iL {T_part}
254 if (ir
->epct
== epctSEMIISOTROPIC
)
263 /* eta is in pure units. veta is in units of ps^-1. GW is in
264 units of ps^-2. However, eta has a reference of 1 nm^3, so care must be
265 taken to use only RATIOS of eta in updating the volume. */
267 /* we take the partial pressure tensors, modify the
268 kinetic energy tensor, and recovert to pressure */
270 if (ir
->opts
.nrdf
[0] == 0)
272 gmx_fatal(FARGS
, "Barostat is coupled to a T-group with no degrees of freedom\n");
274 /* alpha factor for phase space volume, then multiply by the ekin scaling factor. */
275 alpha
= 1.0 + DIM
/((double)ir
->opts
.nrdf
[0]);
276 alpha
*= ekind
->tcstat
[0].ekinscalef_nhc
;
277 msmul(ekind
->ekin
, alpha
, ekinmod
);
278 /* for now, we use Elr = 0, because if you want to get it right, you
279 really should be using PME. Maybe print a warning? */
281 pscal
= calc_pres(ir
->ePBC
, nwall
, box
, ekinmod
, vir
, localpres
)+pcorr
;
284 GW
= (vol
*(MassQ
->Winv
/PRESFAC
))*(DIM
*pscal
- trace(ir
->ref_p
)); /* W is in ps^2 * bar * nm^3 */
290 * This file implements temperature and pressure coupling algorithms:
291 * For now only the Weak coupling and the modified weak coupling.
293 * Furthermore computation of pressure and temperature is done here
297 real
calc_pres(int ePBC
, int nwall
, matrix box
, tensor ekin
, tensor vir
,
303 if (ePBC
== epbcNONE
|| (ePBC
== epbcXY
&& nwall
!= 2))
309 /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
310 * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
314 fac
= PRESFAC
*2.0/det(box
);
315 for (n
= 0; (n
< DIM
); n
++)
317 for (m
= 0; (m
< DIM
); m
++)
319 pres
[n
][m
] = (ekin
[n
][m
] - vir
[n
][m
])*fac
;
325 pr_rvecs(debug
, 0, "PC: pres", pres
, DIM
);
326 pr_rvecs(debug
, 0, "PC: ekin", ekin
, DIM
);
327 pr_rvecs(debug
, 0, "PC: vir ", vir
, DIM
);
328 pr_rvecs(debug
, 0, "PC: box ", box
, DIM
);
331 return trace(pres
)/DIM
;
334 real
calc_temp(real ekin
, real nrdf
)
338 return (2.0*ekin
)/(nrdf
*BOLTZ
);
346 /*! \brief Sets 1/mass for Parrinello-Rahman in wInv; NOTE: PRESFAC is not included, so not in GROMACS units! */
347 static void calcParrinelloRahmanInvMass(const t_inputrec
*ir
, const matrix box
,
352 /* TODO: See if we can make the mass independent of the box size */
353 maxBoxLength
= std::max(box
[XX
][XX
], box
[YY
][YY
]);
354 maxBoxLength
= std::max(maxBoxLength
, box
[ZZ
][ZZ
]);
356 for (int d
= 0; d
< DIM
; d
++)
358 for (int n
= 0; n
< DIM
; n
++)
360 wInv
[d
][n
] = (4*M_PI
*M_PI
*ir
->compress
[d
][n
])/(3*ir
->tau_p
*ir
->tau_p
*maxBoxLength
);
365 void parrinellorahman_pcoupl(FILE *fplog
, gmx_int64_t step
,
366 const t_inputrec
*ir
, real dt
, const tensor pres
,
367 tensor box
, tensor box_rel
, tensor boxv
,
368 tensor M
, matrix mu
, gmx_bool bFirstStep
)
370 /* This doesn't do any coordinate updating. It just
371 * integrates the box vector equations from the calculated
372 * acceleration due to pressure difference. We also compute
373 * the tensor M which is used in update to couple the particle
374 * coordinates to the box vectors.
376 * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
379 * M_nk = (h') * (h' * h + h' h) * h
381 * with the dots denoting time derivatives and h is the transformation from
382 * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
383 * This also goes for the pressure and M tensors - they are transposed relative
384 * to ours. Our equation thus becomes:
387 * M_gmx = M_nk' = b * (b * b' + b * b') * b'
389 * where b is the gromacs box matrix.
390 * Our box accelerations are given by
392 * b = vol/W inv(box') * (P-ref_P) (=h')
395 real vol
= box
[XX
][XX
]*box
[YY
][YY
]*box
[ZZ
][ZZ
];
396 real atot
, arel
, change
;
397 tensor invbox
, pdiff
, t1
, t2
;
399 gmx::invertBoxMatrix(box
, invbox
);
403 /* Note that PRESFAC does not occur here.
404 * The pressure and compressibility always occur as a product,
405 * therefore the pressure unit drops out.
408 calcParrinelloRahmanInvMass(ir
, box
, winv
);
410 m_sub(pres
, ir
->ref_p
, pdiff
);
412 if (ir
->epct
== epctSURFACETENSION
)
414 /* Unlike Berendsen coupling it might not be trivial to include a z
415 * pressure correction here? On the other hand we don't scale the
416 * box momentarily, but change accelerations, so it might not be crucial.
418 real xy_pressure
= 0.5*(pres
[XX
][XX
]+pres
[YY
][YY
]);
419 for (int d
= 0; d
< ZZ
; d
++)
421 pdiff
[d
][d
] = (xy_pressure
-(pres
[ZZ
][ZZ
]-ir
->ref_p
[d
][d
]/box
[d
][d
]));
425 tmmul(invbox
, pdiff
, t1
);
426 /* Move the off-diagonal elements of the 'force' to one side to ensure
427 * that we obey the box constraints.
429 for (int d
= 0; d
< DIM
; d
++)
431 for (int n
= 0; n
< d
; n
++)
433 t1
[d
][n
] += t1
[n
][d
];
440 case epctANISOTROPIC
:
441 for (int d
= 0; d
< DIM
; d
++)
443 for (int n
= 0; n
<= d
; n
++)
445 t1
[d
][n
] *= winv
[d
][n
]*vol
;
450 /* calculate total volume acceleration */
451 atot
= box
[XX
][XX
]*box
[YY
][YY
]*t1
[ZZ
][ZZ
]+
452 box
[XX
][XX
]*t1
[YY
][YY
]*box
[ZZ
][ZZ
]+
453 t1
[XX
][XX
]*box
[YY
][YY
]*box
[ZZ
][ZZ
];
455 /* set all RELATIVE box accelerations equal, and maintain total V
457 for (int d
= 0; d
< DIM
; d
++)
459 for (int n
= 0; n
<= d
; n
++)
461 t1
[d
][n
] = winv
[0][0]*vol
*arel
*box
[d
][n
];
465 case epctSEMIISOTROPIC
:
466 case epctSURFACETENSION
:
467 /* Note the correction to pdiff above for surftens. coupling */
469 /* calculate total XY volume acceleration */
470 atot
= box
[XX
][XX
]*t1
[YY
][YY
]+t1
[XX
][XX
]*box
[YY
][YY
];
471 arel
= atot
/(2*box
[XX
][XX
]*box
[YY
][YY
]);
472 /* set RELATIVE XY box accelerations equal, and maintain total V
473 * change speed. Dont change the third box vector accelerations */
474 for (int d
= 0; d
< ZZ
; d
++)
476 for (int n
= 0; n
<= d
; n
++)
478 t1
[d
][n
] = winv
[d
][n
]*vol
*arel
*box
[d
][n
];
481 for (int n
= 0; n
< DIM
; n
++)
483 t1
[ZZ
][n
] *= winv
[ZZ
][n
]*vol
;
487 gmx_fatal(FARGS
, "Parrinello-Rahman pressure coupling type %s "
488 "not supported yet\n", EPCOUPLTYPETYPE(ir
->epct
));
493 for (int d
= 0; d
< DIM
; d
++)
495 for (int n
= 0; n
<= d
; n
++)
497 boxv
[d
][n
] += dt
*t1
[d
][n
];
499 /* We do NOT update the box vectors themselves here, since
500 * we need them for shifting later. It is instead done last
501 * in the update() routine.
504 /* Calculate the change relative to diagonal elements-
505 since it's perfectly ok for the off-diagonal ones to
506 be zero it doesn't make sense to check the change relative
510 change
= fabs(dt
*boxv
[d
][n
]/box
[d
][d
]);
512 if (change
> maxchange
)
519 if (maxchange
> 0.01 && fplog
)
523 "\nStep %s Warning: Pressure scaling more than 1%%. "
524 "This may mean your system\n is not yet equilibrated. "
525 "Use of Parrinello-Rahman pressure coupling during\n"
526 "equilibration can lead to simulation instability, "
527 "and is discouraged.\n",
528 gmx_step_str(step
, buf
));
532 preserve_box_shape(ir
, box_rel
, boxv
);
534 mtmul(boxv
, box
, t1
); /* t1=boxv * b' */
535 mmul(invbox
, t1
, t2
);
536 mtmul(t2
, invbox
, M
);
538 /* Determine the scaling matrix mu for the coordinates */
539 for (int d
= 0; d
< DIM
; d
++)
541 for (int n
= 0; n
<= d
; n
++)
543 t1
[d
][n
] = box
[d
][n
] + dt
*boxv
[d
][n
];
546 preserve_box_shape(ir
, box_rel
, t1
);
547 /* t1 is the box at t+dt, determine mu as the relative change */
548 mmul_ur0(invbox
, t1
, mu
);
551 void berendsen_pcoupl(FILE *fplog
, gmx_int64_t step
,
552 const t_inputrec
*ir
, real dt
,
553 const tensor pres
, const matrix box
,
554 const matrix force_vir
, const matrix constraint_vir
,
555 matrix mu
, double *baros_integral
)
558 real scalar_pressure
, xy_pressure
, p_corr_z
;
562 * Calculate the scaling matrix mu
566 for (d
= 0; d
< DIM
; d
++)
568 scalar_pressure
+= pres
[d
][d
]/DIM
;
571 xy_pressure
+= pres
[d
][d
]/(DIM
-1);
574 /* Pressure is now in bar, everywhere. */
575 #define factor(d, m) (ir->compress[d][m]*dt/ir->tau_p)
577 /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
578 * necessary for triclinic scaling
584 for (d
= 0; d
< DIM
; d
++)
586 mu
[d
][d
] = 1.0 - factor(d
, d
)*(ir
->ref_p
[d
][d
] - scalar_pressure
) /DIM
;
589 case epctSEMIISOTROPIC
:
590 for (d
= 0; d
< ZZ
; d
++)
592 mu
[d
][d
] = 1.0 - factor(d
, d
)*(ir
->ref_p
[d
][d
]-xy_pressure
)/DIM
;
595 1.0 - factor(ZZ
, ZZ
)*(ir
->ref_p
[ZZ
][ZZ
] - pres
[ZZ
][ZZ
])/DIM
;
597 case epctANISOTROPIC
:
598 for (d
= 0; d
< DIM
; d
++)
600 for (n
= 0; n
< DIM
; n
++)
602 mu
[d
][n
] = (d
== n
? 1.0 : 0.0)
603 -factor(d
, n
)*(ir
->ref_p
[d
][n
] - pres
[d
][n
])/DIM
;
607 case epctSURFACETENSION
:
608 /* ir->ref_p[0/1] is the reference surface-tension times *
609 * the number of surfaces */
610 if (ir
->compress
[ZZ
][ZZ
])
612 p_corr_z
= dt
/ir
->tau_p
*(ir
->ref_p
[ZZ
][ZZ
] - pres
[ZZ
][ZZ
]);
616 /* when the compressibity is zero, set the pressure correction *
617 * in the z-direction to zero to get the correct surface tension */
620 mu
[ZZ
][ZZ
] = 1.0 - ir
->compress
[ZZ
][ZZ
]*p_corr_z
;
621 for (d
= 0; d
< DIM
-1; d
++)
623 mu
[d
][d
] = 1.0 + factor(d
, d
)*(ir
->ref_p
[d
][d
]/(mu
[ZZ
][ZZ
]*box
[ZZ
][ZZ
])
624 - (pres
[ZZ
][ZZ
]+p_corr_z
- xy_pressure
))/(DIM
-1);
628 gmx_fatal(FARGS
, "Berendsen pressure coupling type %s not supported yet\n",
629 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 // cppcheck-suppress unreadVariable
694 nthreads
= gmx_omp_nthreads_get(emntUpdate
);
697 /* Scale the positions */
698 #pragma omp parallel for num_threads(nthreads) schedule(static)
699 for (n
= start
; n
< start
+nr_atoms
; n
++)
701 // Trivial OpenMP region that does not throw
704 if (cFREEZE
== nullptr)
715 x
[n
][XX
] = mu
[XX
][XX
]*x
[n
][XX
]+mu
[YY
][XX
]*x
[n
][YY
]+mu
[ZZ
][XX
]*x
[n
][ZZ
];
719 x
[n
][YY
] = mu
[YY
][YY
]*x
[n
][YY
]+mu
[ZZ
][YY
]*x
[n
][ZZ
];
723 x
[n
][ZZ
] = mu
[ZZ
][ZZ
]*x
[n
][ZZ
];
726 /* compute final boxlengths */
727 for (d
= 0; d
< DIM
; d
++)
729 box
[d
][XX
] = mu
[XX
][XX
]*box
[d
][XX
]+mu
[YY
][XX
]*box
[d
][YY
]+mu
[ZZ
][XX
]*box
[d
][ZZ
];
730 box
[d
][YY
] = mu
[YY
][YY
]*box
[d
][YY
]+mu
[ZZ
][YY
]*box
[d
][ZZ
];
731 box
[d
][ZZ
] = mu
[ZZ
][ZZ
]*box
[d
][ZZ
];
734 preserve_box_shape(ir
, box_rel
, box
);
736 /* (un)shifting should NOT be done after this,
737 * since the box vectors might have changed
739 inc_nrnb(nrnb
, eNR_PCOUPL
, nr_atoms
);
742 void berendsen_tcoupl(const t_inputrec
*ir
, const gmx_ekindata_t
*ekind
, real dt
,
743 std::vector
<double> &therm_integral
)
745 const t_grpopts
*opts
= &ir
->opts
;
747 for (int i
= 0; (i
< opts
->ngtc
); i
++)
753 Ek
= trace(ekind
->tcstat
[i
].ekinf
);
754 T
= ekind
->tcstat
[i
].T
;
758 Ek
= trace(ekind
->tcstat
[i
].ekinh
);
759 T
= ekind
->tcstat
[i
].Th
;
762 if ((opts
->tau_t
[i
] > 0) && (T
> 0.0))
764 real reft
= std::max
<real
>(0, opts
->ref_t
[i
]);
765 real lll
= std::sqrt(1.0 + (dt
/opts
->tau_t
[i
])*(reft
/T
-1.0));
766 ekind
->tcstat
[i
].lambda
= std::max
<real
>(std::min
<real
>(lll
, 1.25), 0.8);
770 ekind
->tcstat
[i
].lambda
= 1.0;
773 /* Keep track of the amount of energy we are adding to the system */
774 therm_integral
[i
] -= (gmx::square(ekind
->tcstat
[i
].lambda
) - 1)*Ek
;
778 fprintf(debug
, "TC: group %d: T: %g, Lambda: %g\n",
779 i
, T
, ekind
->tcstat
[i
].lambda
);
784 void andersen_tcoupl(t_inputrec
*ir
, gmx_int64_t step
,
785 const t_commrec
*cr
, const t_mdatoms
*md
, t_state
*state
, real rate
, const gmx_bool
*randomize
, const real
*boltzfac
)
787 const int *gatindex
= (DOMAINDECOMP(cr
) ? cr
->dd
->gatindex
: nullptr);
790 gmx::ThreeFry2x64
<0> rng(ir
->andersen_seed
, gmx::RandomDomain::Thermostat
);
791 gmx::UniformRealDistribution
<real
> uniformDist
;
792 gmx::TabulatedNormalDistribution
<real
, 14> normalDist
;
794 /* randomize the velocities of the selected particles */
796 for (i
= 0; i
< md
->homenr
; i
++) /* now loop over the list of atoms */
798 int ng
= gatindex
? gatindex
[i
] : i
;
801 rng
.restart(step
, ng
);
805 gc
= md
->cTC
[i
]; /* assign the atom to a temperature group if there are more than one */
809 if (ir
->etc
== etcANDERSENMASSIVE
)
811 /* Randomize particle always */
816 /* Randomize particle probabilistically */
818 bRandomize
= uniformDist(rng
) < rate
;
825 scal
= std::sqrt(boltzfac
[gc
]*md
->invmass
[i
]);
829 for (d
= 0; d
< DIM
; d
++)
831 state
->v
[i
][d
] = scal
*normalDist(rng
);
839 void nosehoover_tcoupl(t_grpopts
*opts
, gmx_ekindata_t
*ekind
, real dt
,
840 double xi
[], double vxi
[], t_extmass
*MassQ
)
845 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
847 for (i
= 0; (i
< opts
->ngtc
); i
++)
849 reft
= std::max
<real
>(0, opts
->ref_t
[i
]);
851 vxi
[i
] += dt
*MassQ
->Qinv
[i
]*(ekind
->tcstat
[i
].Th
- reft
);
852 xi
[i
] += dt
*(oldvxi
+ vxi
[i
])*0.5;
856 void trotter_update(t_inputrec
*ir
, gmx_int64_t step
, gmx_ekindata_t
*ekind
,
857 gmx_enerdata_t
*enerd
, t_state
*state
,
858 tensor vir
, t_mdatoms
*md
,
859 t_extmass
*MassQ
, int **trotter_seqlist
, int trotter_seqno
)
862 int n
, i
, d
, ngtc
, gc
= 0, t
;
863 t_grp_tcstat
*tcstat
;
865 gmx_int64_t step_eff
;
867 double *scalefac
, dtc
;
869 rvec sumv
= {0, 0, 0};
872 if (trotter_seqno
<= ettTSEQ2
)
874 step_eff
= step
-1; /* the velocity verlet calls are actually out of order -- the first half step
875 is actually the last half step from the previous step. Thus the first half step
876 actually corresponds to the n-1 step*/
884 bCouple
= (ir
->nsttcouple
== 1 ||
885 do_per_step(step_eff
+ir
->nsttcouple
, ir
->nsttcouple
));
887 trotter_seq
= trotter_seqlist
[trotter_seqno
];
889 if ((trotter_seq
[0] == etrtSKIPALL
) || (!bCouple
))
893 dtc
= ir
->nsttcouple
*ir
->delta_t
; /* This is OK for NPT, because nsttcouple == nstpcouple is enforcesd */
894 opts
= &(ir
->opts
); /* just for ease of referencing */
897 snew(scalefac
, opts
->ngtc
);
898 for (i
= 0; i
< ngtc
; i
++)
902 /* execute the series of trotter updates specified in the trotterpart array */
904 for (i
= 0; i
< NTROTTERPARTS
; i
++)
906 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
907 if ((trotter_seq
[i
] == etrtBAROV2
) || (trotter_seq
[i
] == etrtBARONHC2
) || (trotter_seq
[i
] == etrtNHC2
))
916 switch (trotter_seq
[i
])
920 boxv_trotter(ir
, &(state
->veta
), dt
, state
->box
, ekind
, vir
,
921 enerd
->term
[F_PDISPCORR
], MassQ
);
925 NHC_trotter(opts
, state
->nnhpres
, ekind
, dt
, state
->nhpres_xi
.data(),
926 state
->nhpres_vxi
.data(), nullptr, &(state
->veta
), MassQ
, FALSE
);
930 NHC_trotter(opts
, opts
->ngtc
, ekind
, dt
, state
->nosehoover_xi
.data(),
931 state
->nosehoover_vxi
.data(), scalefac
, nullptr, MassQ
, (ir
->eI
== eiVV
));
932 /* need to rescale the kinetic energies and velocities here. Could
933 scale the velocities later, but we need them scaled in order to
934 produce the correct outputs, so we'll scale them here. */
936 for (t
= 0; t
< ngtc
; t
++)
938 tcstat
= &ekind
->tcstat
[t
];
939 tcstat
->vscale_nhc
= scalefac
[t
];
940 tcstat
->ekinscaleh_nhc
*= (scalefac
[t
]*scalefac
[t
]);
941 tcstat
->ekinscalef_nhc
*= (scalefac
[t
]*scalefac
[t
]);
943 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
944 /* but do we actually need the total? */
946 /* modify the velocities as well */
947 for (n
= 0; n
< md
->homenr
; n
++)
949 if (md
->cTC
) /* does this conditional need to be here? is this always true?*/
953 for (d
= 0; d
< DIM
; d
++)
955 state
->v
[n
][d
] *= scalefac
[gc
];
960 for (d
= 0; d
< DIM
; d
++)
962 sumv
[d
] += (state
->v
[n
][d
])/md
->invmass
[n
];
971 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
977 for (d
= 0; d
< DIM
; d
++)
979 consk
[d
] = sumv
[d
]*exp((1 + 1.0/opts
->nrdf
[0])*((1.0/DIM
)*log(det(state
->box
)/state
->vol0
)) + state
->nosehoover_xi
[0]);
981 fprintf(debug
, "Conserved kappa: %15.8f %15.8f %15.8f\n", consk
[0], consk
[1], consk
[2]);
989 extern void init_npt_masses(t_inputrec
*ir
, t_state
*state
, t_extmass
*MassQ
, gmx_bool bInit
)
991 int n
, i
, j
, d
, ngtc
, nh
;
993 real reft
, kT
, ndj
, nd
;
995 opts
= &(ir
->opts
); /* just for ease of referencing */
996 ngtc
= ir
->opts
.ngtc
;
997 nh
= state
->nhchainlength
;
1003 snew(MassQ
->Qinv
, ngtc
);
1005 for (i
= 0; (i
< ngtc
); i
++)
1007 if ((opts
->tau_t
[i
] > 0) && (opts
->ref_t
[i
] > 0))
1009 MassQ
->Qinv
[i
] = 1.0/(gmx::square(opts
->tau_t
[i
]/M_2PI
)*opts
->ref_t
[i
]);
1013 MassQ
->Qinv
[i
] = 0.0;
1017 else if (EI_VV(ir
->eI
))
1019 /* Set pressure variables */
1023 if (state
->vol0
== 0)
1025 state
->vol0
= det(state
->box
);
1026 /* because we start by defining a fixed
1027 compressibility, we need the volume at this
1028 compressibility to solve the problem. */
1032 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
1033 /* Consider evaluating eventually if this the right mass to use. All are correct, some might be more stable */
1034 MassQ
->Winv
= (PRESFAC
*trace(ir
->compress
)*BOLTZ
*opts
->ref_t
[0])/(DIM
*state
->vol0
*gmx::square(ir
->tau_p
/M_2PI
));
1035 /* An alternate mass definition, from Tuckerman et al. */
1036 /* MassQ->Winv = 1.0/(gmx::square(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
1037 for (d
= 0; d
< DIM
; d
++)
1039 for (n
= 0; n
< DIM
; n
++)
1041 MassQ
->Winvm
[d
][n
] = PRESFAC
*ir
->compress
[d
][n
]/(state
->vol0
*gmx::square(ir
->tau_p
/M_2PI
));
1042 /* not clear this is correct yet for the anisotropic case. Will need to reevaluate
1043 before using MTTK for anisotropic states.*/
1046 /* Allocate space for thermostat variables */
1049 snew(MassQ
->Qinv
, ngtc
*nh
);
1052 /* now, set temperature variables */
1053 for (i
= 0; i
< ngtc
; i
++)
1055 if (opts
->tau_t
[i
] > 0 && opts
->ref_t
[i
] > 0 && opts
->nrdf
[i
] > 0)
1057 reft
= std::max
<real
>(0, opts
->ref_t
[i
]);
1060 for (j
= 0; j
< nh
; j
++)
1070 MassQ
->Qinv
[i
*nh
+j
] = 1.0/(gmx::square(opts
->tau_t
[i
]/M_2PI
)*ndj
*kT
);
1075 for (j
= 0; j
< nh
; j
++)
1077 MassQ
->Qinv
[i
*nh
+j
] = 0.0;
1084 int **init_npt_vars(t_inputrec
*ir
, t_state
*state
, t_extmass
*MassQ
, gmx_bool bTrotter
)
1086 int i
, j
, nnhpres
, nh
;
1088 real bmass
, qmass
, reft
, kT
;
1091 opts
= &(ir
->opts
); /* just for ease of referencing */
1092 nnhpres
= state
->nnhpres
;
1093 nh
= state
->nhchainlength
;
1095 if (EI_VV(ir
->eI
) && (ir
->epc
== epcMTTK
) && (ir
->etc
!= etcNOSEHOOVER
))
1097 gmx_fatal(FARGS
, "Cannot do MTTK pressure coupling without Nose-Hoover temperature control");
1100 init_npt_masses(ir
, state
, MassQ
, TRUE
);
1102 /* first, initialize clear all the trotter calls */
1103 snew(trotter_seq
, ettTSEQMAX
);
1104 for (i
= 0; i
< ettTSEQMAX
; i
++)
1106 snew(trotter_seq
[i
], NTROTTERPARTS
);
1107 for (j
= 0; j
< NTROTTERPARTS
; j
++)
1109 trotter_seq
[i
][j
] = etrtNONE
;
1111 trotter_seq
[i
][0] = etrtSKIPALL
;
1116 /* no trotter calls, so we never use the values in the array.
1117 * We access them (so we need to define them, but ignore
1123 /* compute the kinetic energy by using the half step velocities or
1124 * the kinetic energies, depending on the order of the trotter calls */
1128 if (inputrecNptTrotter(ir
))
1130 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1131 We start with the initial one. */
1132 /* first, a round that estimates veta. */
1133 trotter_seq
[0][0] = etrtBAROV
;
1135 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1137 /* The first half trotter update */
1138 trotter_seq
[2][0] = etrtBAROV
;
1139 trotter_seq
[2][1] = etrtNHC
;
1140 trotter_seq
[2][2] = etrtBARONHC
;
1142 /* The second half trotter update */
1143 trotter_seq
[3][0] = etrtBARONHC
;
1144 trotter_seq
[3][1] = etrtNHC
;
1145 trotter_seq
[3][2] = etrtBAROV
;
1147 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1150 else if (inputrecNvtTrotter(ir
))
1152 /* This is the easy version - there are only two calls, both the same.
1153 Otherwise, even easier -- no calls */
1154 trotter_seq
[2][0] = etrtNHC
;
1155 trotter_seq
[3][0] = etrtNHC
;
1157 else if (inputrecNphTrotter(ir
))
1159 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1160 We start with the initial one. */
1161 /* first, a round that estimates veta. */
1162 trotter_seq
[0][0] = etrtBAROV
;
1164 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1166 /* The first half trotter update */
1167 trotter_seq
[2][0] = etrtBAROV
;
1168 trotter_seq
[2][1] = etrtBARONHC
;
1170 /* The second half trotter update */
1171 trotter_seq
[3][0] = etrtBARONHC
;
1172 trotter_seq
[3][1] = etrtBAROV
;
1174 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1177 else if (ir
->eI
== eiVVAK
)
1179 if (inputrecNptTrotter(ir
))
1181 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1182 We start with the initial one. */
1183 /* first, a round that estimates veta. */
1184 trotter_seq
[0][0] = etrtBAROV
;
1186 /* The first half trotter update, part 1 -- double update, because it commutes */
1187 trotter_seq
[1][0] = etrtNHC
;
1189 /* The first half trotter update, part 2 */
1190 trotter_seq
[2][0] = etrtBAROV
;
1191 trotter_seq
[2][1] = etrtBARONHC
;
1193 /* The second half trotter update, part 1 */
1194 trotter_seq
[3][0] = etrtBARONHC
;
1195 trotter_seq
[3][1] = etrtBAROV
;
1197 /* The second half trotter update */
1198 trotter_seq
[4][0] = etrtNHC
;
1200 else if (inputrecNvtTrotter(ir
))
1202 /* This is the easy version - there is only one call, both the same.
1203 Otherwise, even easier -- no calls */
1204 trotter_seq
[1][0] = etrtNHC
;
1205 trotter_seq
[4][0] = etrtNHC
;
1207 else if (inputrecNphTrotter(ir
))
1209 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1210 We start with the initial one. */
1211 /* first, a round that estimates veta. */
1212 trotter_seq
[0][0] = etrtBAROV
;
1214 /* The first half trotter update, part 1 -- leave zero */
1215 trotter_seq
[1][0] = etrtNHC
;
1217 /* The first half trotter update, part 2 */
1218 trotter_seq
[2][0] = etrtBAROV
;
1219 trotter_seq
[2][1] = etrtBARONHC
;
1221 /* The second half trotter update, part 1 */
1222 trotter_seq
[3][0] = etrtBARONHC
;
1223 trotter_seq
[3][1] = etrtBAROV
;
1225 /* The second half trotter update -- blank for now */
1233 bmass
= DIM
*DIM
; /* recommended mass parameters for isotropic barostat */
1236 snew(MassQ
->QPinv
, nnhpres
*opts
->nhchainlength
);
1238 /* barostat temperature */
1239 if ((ir
->tau_p
> 0) && (opts
->ref_t
[0] > 0))
1241 reft
= std::max
<real
>(0, opts
->ref_t
[0]);
1243 for (i
= 0; i
< nnhpres
; i
++)
1245 for (j
= 0; j
< nh
; j
++)
1255 MassQ
->QPinv
[i
*opts
->nhchainlength
+j
] = 1.0/(gmx::square(opts
->tau_t
[0]/M_2PI
)*qmass
*kT
);
1261 for (i
= 0; i
< nnhpres
; i
++)
1263 for (j
= 0; j
< nh
; j
++)
1265 MassQ
->QPinv
[i
*nh
+j
] = 0.0;
1272 static real
energyNoseHoover(const t_inputrec
*ir
, const t_state
*state
, const t_extmass
*MassQ
)
1276 int nh
= state
->nhchainlength
;
1278 for (int i
= 0; i
< ir
->opts
.ngtc
; i
++)
1280 const double *ixi
= &state
->nosehoover_xi
[i
*nh
];
1281 const double *ivxi
= &state
->nosehoover_vxi
[i
*nh
];
1282 const double *iQinv
= &(MassQ
->Qinv
[i
*nh
]);
1284 int nd
= ir
->opts
.nrdf
[i
];
1285 real reft
= std::max
<real
>(ir
->opts
.ref_t
[i
], 0);
1286 real kT
= BOLTZ
* reft
;
1290 if (inputrecNvtTrotter(ir
))
1292 /* contribution from the thermal momenta of the NH chain */
1293 for (int j
= 0; j
< nh
; j
++)
1297 energy
+= 0.5*gmx::square(ivxi
[j
])/iQinv
[j
];
1298 /* contribution from the thermal variable of the NH chain */
1308 energy
+= ndj
*ixi
[j
]*kT
;
1312 else /* Other non Trotter temperature NH control -- no chains yet. */
1314 energy
+= 0.5*BOLTZ
*nd
*gmx::square(ivxi
[0])/iQinv
[0];
1315 energy
+= nd
*ixi
[0]*kT
;
1323 /* Returns the energy from the barostat thermostat chain */
1324 static real
energyPressureMTTK(const t_inputrec
*ir
, const t_state
*state
, const t_extmass
*MassQ
)
1328 int nh
= state
->nhchainlength
;
1330 for (int i
= 0; i
< state
->nnhpres
; i
++)
1332 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1333 real reft
= std::max
<real
>(ir
->opts
.ref_t
[0], 0.0); /* using 'System' temperature */
1334 real kT
= BOLTZ
* reft
;
1336 for (int j
= 0; j
< nh
; j
++)
1338 double iQinv
= MassQ
->QPinv
[i
*nh
+ j
];
1341 energy
+= 0.5*gmx::square(state
->nhpres_vxi
[i
*nh
+ j
]/iQinv
);
1342 /* contribution from the thermal variable of the NH chain */
1343 energy
+= state
->nhpres_xi
[i
*nh
+ j
]*kT
;
1347 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
]);
1355 /* Returns the energy accumulated by the V-rescale or Berendsen thermostat */
1356 static real
energyVrescale(const t_inputrec
*ir
, const t_state
*state
)
1359 for (int i
= 0; i
< ir
->opts
.ngtc
; i
++)
1361 energy
+= state
->therm_integral
[i
];
1367 real
NPT_energy(const t_inputrec
*ir
, const t_state
*state
, const t_extmass
*MassQ
)
1371 if (ir
->epc
!= epcNO
)
1373 /* Compute the contribution of the pressure to the conserved quantity*/
1375 real vol
= det(state
->box
);
1379 case epcPARRINELLORAHMAN
:
1381 /* contribution from the pressure momenta */
1383 calcParrinelloRahmanInvMass(ir
, state
->box
, invMass
);
1384 for (int d
= 0; d
< DIM
; d
++)
1386 for (int n
= 0; n
<= d
; n
++)
1388 if (invMass
[d
][n
] > 0)
1390 energyNPT
+= 0.5*gmx::square(state
->boxv
[d
][n
])/(invMass
[d
][n
]*PRESFAC
);
1395 /* Contribution from the PV term.
1396 * Not that with non-zero off-diagonal reference pressures,
1397 * i.e. applied shear stresses, there are additional terms.
1398 * We don't support this here, since that requires keeping
1399 * track of unwrapped box diagonal elements. This case is
1400 * excluded in integratorHasConservedEnergyQuantity().
1402 energyNPT
+= vol
*trace(ir
->ref_p
)/(DIM
*PRESFAC
);
1406 /* contribution from the pressure momenta */
1407 energyNPT
+= 0.5*gmx::square(state
->veta
)/MassQ
->Winv
;
1409 /* contribution from the PV term */
1410 energyNPT
+= vol
*trace(ir
->ref_p
)/(DIM
*PRESFAC
);
1412 if (ir
->epc
== epcMTTK
)
1414 /* contribution from the MTTK chain */
1415 energyNPT
+= energyPressureMTTK(ir
, state
, MassQ
);
1419 energyNPT
+= state
->baros_integral
;
1422 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().");
1432 energyNPT
+= energyVrescale(ir
, state
);
1435 energyNPT
+= energyNoseHoover(ir
, state
, MassQ
);
1438 case etcANDERSENMASSIVE
:
1439 // Not supported, excluded in integratorHasConservedEnergyQuantity()
1442 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().");
1449 static real
vrescale_sumnoises(real nn
,
1450 gmx::ThreeFry2x64
<> *rng
,
1451 gmx::NormalDistribution
<real
> *normalDist
)
1454 * Returns the sum of nn independent gaussian noises squared
1455 * (i.e. equivalent to summing the square of the return values
1456 * of nn calls to a normal distribution).
1458 const real ndeg_tol
= 0.0001;
1460 gmx::GammaDistribution
<real
> gammaDist(0.5*nn
, 1.0);
1462 if (nn
< 2 + ndeg_tol
)
1467 nn_int
= (int)(nn
+ 0.5);
1469 if (nn
- nn_int
< -ndeg_tol
|| nn
- nn_int
> ndeg_tol
)
1471 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);
1475 for (i
= 0; i
< nn_int
; i
++)
1477 gauss
= (*normalDist
)(*rng
);
1483 /* Use a gamma distribution for any real nn > 2 */
1484 r
= 2.0*gammaDist(*rng
);
1490 static real
vrescale_resamplekin(real kk
, real sigma
, real ndeg
, real taut
,
1491 gmx_int64_t step
, gmx_int64_t seed
)
1494 * Generates a new value for the kinetic energy,
1495 * according to Bussi et al JCP (2007), Eq. (A7)
1496 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1497 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1498 * ndeg: number of degrees of freedom of the atoms to be thermalized
1499 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1501 real factor
, rr
, ekin_new
;
1502 gmx::ThreeFry2x64
<64> rng(seed
, gmx::RandomDomain::Thermostat
);
1503 gmx::NormalDistribution
<real
> normalDist
;
1507 factor
= exp(-1.0/taut
);
1514 rng
.restart(step
, 0);
1516 rr
= normalDist(rng
);
1520 (1.0 - factor
)*(sigma
*(vrescale_sumnoises(ndeg
-1, &rng
, &normalDist
) + rr
*rr
)/ndeg
- kk
) +
1521 2.0*rr
*std::sqrt(kk
*sigma
/ndeg
*(1.0 - factor
)*factor
);
1526 void vrescale_tcoupl(t_inputrec
*ir
, gmx_int64_t step
,
1527 gmx_ekindata_t
*ekind
, real dt
,
1528 double therm_integral
[])
1532 real Ek
, Ek_ref1
, Ek_ref
, Ek_new
;
1536 for (i
= 0; (i
< opts
->ngtc
); i
++)
1540 Ek
= trace(ekind
->tcstat
[i
].ekinf
);
1544 Ek
= trace(ekind
->tcstat
[i
].ekinh
);
1547 if (opts
->tau_t
[i
] >= 0 && opts
->nrdf
[i
] > 0 && Ek
> 0)
1549 Ek_ref1
= 0.5*opts
->ref_t
[i
]*BOLTZ
;
1550 Ek_ref
= Ek_ref1
*opts
->nrdf
[i
];
1552 Ek_new
= vrescale_resamplekin(Ek
, Ek_ref
, opts
->nrdf
[i
],
1556 /* Analytically Ek_new>=0, but we check for rounding errors */
1559 ekind
->tcstat
[i
].lambda
= 0.0;
1563 ekind
->tcstat
[i
].lambda
= std::sqrt(Ek_new
/Ek
);
1566 therm_integral
[i
] -= Ek_new
- Ek
;
1570 fprintf(debug
, "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
1571 i
, Ek_ref
, Ek
, Ek_new
, ekind
->tcstat
[i
].lambda
);
1576 ekind
->tcstat
[i
].lambda
= 1.0;
1581 void rescale_velocities(gmx_ekindata_t
*ekind
, t_mdatoms
*mdatoms
,
1582 int start
, int end
, rvec v
[])
1585 t_grp_tcstat
*tcstat
;
1586 unsigned short *cACC
, *cTC
;
1591 tcstat
= ekind
->tcstat
;
1596 gstat
= ekind
->grpstat
;
1597 cACC
= mdatoms
->cACC
;
1601 for (n
= start
; n
< end
; n
++)
1611 /* Only scale the velocity component relative to the COM velocity */
1612 rvec_sub(v
[n
], gstat
[ga
].u
, vrel
);
1613 lg
= tcstat
[gt
].lambda
;
1614 for (d
= 0; d
< DIM
; d
++)
1616 v
[n
][d
] = gstat
[ga
].u
[d
] + lg
*vrel
[d
];
1623 for (n
= start
; n
< end
; n
++)
1629 lg
= tcstat
[gt
].lambda
;
1630 for (d
= 0; d
< DIM
; d
++)
1639 /* set target temperatures if we are annealing */
1640 void update_annealing_target_temp(t_inputrec
*ir
, real t
, gmx_update_t
*upd
)
1642 int i
, j
, n
, npoints
;
1643 real pert
, thist
= 0, x
;
1645 for (i
= 0; i
< ir
->opts
.ngtc
; i
++)
1647 npoints
= ir
->opts
.anneal_npoints
[i
];
1648 switch (ir
->opts
.annealing
[i
])
1653 /* calculate time modulo the period */
1654 pert
= ir
->opts
.anneal_time
[i
][npoints
-1];
1655 n
= static_cast<int>(t
/ pert
);
1656 thist
= t
- n
*pert
; /* modulo time */
1657 /* Make sure rounding didn't get us outside the interval */
1658 if (fabs(thist
-pert
) < GMX_REAL_EPS
*100)
1667 gmx_fatal(FARGS
, "Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)", i
, ir
->opts
.ngtc
, npoints
);
1669 /* We are doing annealing for this group if we got here,
1670 * and we have the (relative) time as thist.
1671 * calculate target temp */
1673 while ((j
< npoints
-1) && (thist
> (ir
->opts
.anneal_time
[i
][j
+1])))
1679 /* Found our position between points j and j+1.
1680 * Interpolate: x is the amount from j+1, (1-x) from point j
1681 * First treat possible jumps in temperature as a special case.
1683 if ((ir
->opts
.anneal_time
[i
][j
+1]-ir
->opts
.anneal_time
[i
][j
]) < GMX_REAL_EPS
*100)
1685 ir
->opts
.ref_t
[i
] = ir
->opts
.anneal_temp
[i
][j
+1];
1689 x
= ((thist
-ir
->opts
.anneal_time
[i
][j
])/
1690 (ir
->opts
.anneal_time
[i
][j
+1]-ir
->opts
.anneal_time
[i
][j
]));
1691 ir
->opts
.ref_t
[i
] = x
*ir
->opts
.anneal_temp
[i
][j
+1]+(1-x
)*ir
->opts
.anneal_temp
[i
][j
];
1696 ir
->opts
.ref_t
[i
] = ir
->opts
.anneal_temp
[i
][npoints
-1];
1700 update_temperature_constants(upd
, ir
);