src/gromacs/mdlib/coupling.cpp

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