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