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 = 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 / PRESFAC)) * (DIM * pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
 576
 577     *veta += 0.5 * dt * GW;
 578 }
 579
 580 /*
 581  * This file implements temperature and pressure coupling algorithms:
 582  * For now only the Weak coupling and the modified weak coupling.
 583  *
 584  * Furthermore computation of pressure and temperature is done here
 585  *
 586  */
 587
 588 real calc_pres(PbcType pbcType, int nwall, const matrix box, const tensor ekin, const tensor vir, tensor pres)
 589 {
 590     int  n, m;
 591     real fac;
 592
 593     if (pbcType == PbcType::No || (pbcType == PbcType::XY && nwall != 2))
 594     {
 595         clear_mat(pres);
 596     }
 597     else
 598     {
 599         /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
 600          * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
 601          * het systeem...
 602          */
 603
 604         fac = PRESFAC * 2.0 / det(box);
 605         for (n = 0; (n < DIM); n++)
 606         {
 607             for (m = 0; (m < DIM); m++)
 608             {
 609                 pres[n][m] = (ekin[n][m] - vir[n][m]) * fac;
 610             }
 611         }
 612
 613         if (debug)
 614         {
 615             pr_rvecs(debug, 0, "PC: pres", pres, DIM);
 616             pr_rvecs(debug, 0, "PC: ekin", ekin, DIM);
 617             pr_rvecs(debug, 0, "PC: vir ", vir, DIM);
 618             pr_rvecs(debug, 0, "PC: box ", box, DIM);
 619         }
 620     }
 621     return trace(pres) / DIM;
 622 }
 623
 624 real calc_temp(real ekin, real nrdf)
 625 {
 626     if (nrdf > 0)
 627     {
 628         return (2.0 * ekin) / (nrdf * BOLTZ);
 629     }
 630     else
 631     {
 632         return 0;
 633     }
 634 }
 635
 636 /*! \brief Sets 1/mass for Parrinello-Rahman in wInv; NOTE: PRESFAC is not included, so not in GROMACS units! */
 637 static void calcParrinelloRahmanInvMass(const t_inputrec* ir, const matrix box, tensor wInv)
 638 {
 639     real maxBoxLength;
 640
 641     /* TODO: See if we can make the mass independent of the box size */
 642     maxBoxLength = std::max(box[XX][XX], box[YY][YY]);
 643     maxBoxLength = std::max(maxBoxLength, box[ZZ][ZZ]);
 644
 645     for (int d = 0; d < DIM; d++)
 646     {
 647         for (int n = 0; n < DIM; n++)
 648         {
 649             wInv[d][n] = (4 * M_PI * M_PI * ir->compress[d][n]) / (3 * ir->tau_p * ir->tau_p * maxBoxLength);
 650         }
 651     }
 652 }
 653
 654 void parrinellorahman_pcoupl(FILE*             fplog,
 655                              int64_t           step,
 656                              const t_inputrec* ir,
 657                              real              dt,
 658                              const tensor      pres,
 659                              const tensor      box,
 660                              tensor            box_rel,
 661                              tensor            boxv,
 662                              tensor            M,
 663                              matrix            mu,
 664                              gmx_bool          bFirstStep)
 665 {
 666     /* This doesn't do any coordinate updating. It just
 667      * integrates the box vector equations from the calculated
 668      * acceleration due to pressure difference. We also compute
 669      * the tensor M which is used in update to couple the particle
 670      * coordinates to the box vectors.
 671      *
 672      * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
 673      * given as
 674      *            -1    .           .     -1
 675      * M_nk = (h')   * (h' * h + h' h) * h
 676      *
 677      * with the dots denoting time derivatives and h is the transformation from
 678      * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
 679      * This also goes for the pressure and M tensors - they are transposed relative
 680      * to ours. Our equation thus becomes:
 681      *
 682      *                  -1       .    .           -1
 683      * M_gmx = M_nk' = b  * (b * b' + b * b') * b'
 684      *
 685      * where b is the gromacs box matrix.
 686      * Our box accelerations are given by
 687      *   ..                                    ..
 688      *   b = vol/W inv(box') * (P-ref_P)     (=h')
 689      */
 690
 691     real   vol = box[XX][XX] * box[YY][YY] * box[ZZ][ZZ];
 692     real   atot, arel, change;
 693     tensor invbox, pdiff, t1, t2;
 694
 695     gmx::invertBoxMatrix(box, invbox);
 696
 697     if (!bFirstStep)
 698     {
 699         /* Note that PRESFAC does not occur here.
 700          * The pressure and compressibility always occur as a product,
 701          * therefore the pressure unit drops out.
 702          */
 703         tensor winv;
 704         calcParrinelloRahmanInvMass(ir, box, winv);
 705
 706         m_sub(pres, ir->ref_p, pdiff);
 707
 708         if (ir->epct == PressureCouplingType::SurfaceTension)
 709         {
 710             /* Unlike Berendsen coupling it might not be trivial to include a z
 711              * pressure correction here? On the other hand we don't scale the
 712              * box momentarily, but change accelerations, so it might not be crucial.
 713              */
 714             real xy_pressure = 0.5 * (pres[XX][XX] + pres[YY][YY]);
 715             for (int d = 0; d < ZZ; d++)
 716             {
 717                 pdiff[d][d] = (xy_pressure - (pres[ZZ][ZZ] - ir->ref_p[d][d] / box[d][d]));
 718             }
 719         }
 720
 721         tmmul(invbox, pdiff, t1);
 722         /* Move the off-diagonal elements of the 'force' to one side to ensure
 723          * that we obey the box constraints.
 724          */
 725         for (int d = 0; d < DIM; d++)
 726         {
 727             for (int n = 0; n < d; n++)
 728             {
 729                 t1[d][n] += t1[n][d];
 730                 t1[n][d] = 0;
 731             }
 732         }
 733
 734         switch (ir->epct)
 735         {
 736             case PressureCouplingType::Anisotropic:
 737                 for (int d = 0; d < DIM; d++)
 738                 {
 739                     for (int n = 0; n <= d; n++)
 740                     {
 741                         t1[d][n] *= winv[d][n] * vol;
 742                     }
 743                 }
 744                 break;
 745             case PressureCouplingType::Isotropic:
 746                 /* calculate total volume acceleration */
 747                 atot = box[XX][XX] * box[YY][YY] * t1[ZZ][ZZ] + box[XX][XX] * t1[YY][YY] * box[ZZ][ZZ]
 748                        + t1[XX][XX] * box[YY][YY] * box[ZZ][ZZ];
 749                 arel = atot / (3 * vol);
 750                 /* set all RELATIVE box accelerations equal, and maintain total V
 751                  * change speed */
 752                 for (int d = 0; d < DIM; d++)
 753                 {
 754                     for (int n = 0; n <= d; n++)
 755                     {
 756                         t1[d][n] = winv[0][0] * vol * arel * box[d][n];
 757                     }
 758                 }
 759                 break;
 760             case PressureCouplingType::SemiIsotropic:
 761             case PressureCouplingType::SurfaceTension:
 762                 /* Note the correction to pdiff above for surftens. coupling  */
 763
 764                 /* calculate total XY volume acceleration */
 765                 atot = box[XX][XX] * t1[YY][YY] + t1[XX][XX] * box[YY][YY];
 766                 arel = atot / (2 * box[XX][XX] * box[YY][YY]);
 767                 /* set RELATIVE XY box accelerations equal, and maintain total V
 768                  * change speed. Dont change the third box vector accelerations */
 769                 for (int d = 0; d < ZZ; d++)
 770                 {
 771                     for (int n = 0; n <= d; n++)
 772                     {
 773                         t1[d][n] = winv[d][n] * vol * arel * box[d][n];
 774                     }
 775                 }
 776                 for (int n = 0; n < DIM; n++)
 777                 {
 778                     t1[ZZ][n] *= winv[ZZ][n] * vol;
 779                 }
 780                 break;
 781             default:
 782                 gmx_fatal(FARGS,
 783                           "Parrinello-Rahman pressure coupling type %s "
 784                           "not supported yet\n",
 785                           enumValueToString(ir->epct));
 786         }
 787
 788         real maxchange = 0;
 789         for (int d = 0; d < DIM; d++)
 790         {
 791             for (int n = 0; n <= d; n++)
 792             {
 793                 boxv[d][n] += dt * t1[d][n];
 794
 795                 /* We do NOT update the box vectors themselves here, since
 796                  * we need them for shifting later. It is instead done last
 797                  * in the update() routine.
 798                  */
 799
 800                 /* Calculate the change relative to diagonal elements-
 801                    since it's perfectly ok for the off-diagonal ones to
 802                    be zero it doesn't make sense to check the change relative
 803                    to its current size.
 804                  */
 805
 806                 change = std::fabs(dt * boxv[d][n] / box[d][d]);
 807
 808                 if (change > maxchange)
 809                 {
 810                     maxchange = change;
 811                 }
 812             }
 813         }
 814
 815         if (maxchange > 0.01 && fplog)
 816         {
 817             char buf[22];
 818             fprintf(fplog,
 819                     "\nStep %s  Warning: Pressure scaling more than 1%%. "
 820                     "This may mean your system\n is not yet equilibrated. "
 821                     "Use of Parrinello-Rahman pressure coupling during\n"
 822                     "equilibration can lead to simulation instability, "
 823                     "and is discouraged.\n",
 824                     gmx_step_str(step, buf));
 825         }
 826     }
 827
 828     preserve_box_shape(ir, box_rel, boxv);
 829
 830     mtmul(boxv, box, t1); /* t1=boxv * b' */
 831     mmul(invbox, t1, t2);
 832     mtmul(t2, invbox, M);
 833
 834     /* Determine the scaling matrix mu for the coordinates */
 835     for (int d = 0; d < DIM; d++)
 836     {
 837         for (int n = 0; n <= d; n++)
 838         {
 839             t1[d][n] = box[d][n] + dt * boxv[d][n];
 840         }
 841     }
 842     preserve_box_shape(ir, box_rel, t1);
 843     /* t1 is the box at t+dt, determine mu as the relative change */
 844     mmul_ur0(invbox, t1, mu);
 845 }
 846
 847 //! Return compressibility factor for entry (i,j) of Berendsen / C-rescale scaling matrix
 848 static inline real compressibilityFactor(int i, int j, const t_inputrec* ir, real dt)
 849 {
 850     return ir->compress[i][j] * dt / ir->tau_p;
 851 }
 852
 853 //! Details of Berendsen / C-rescale scaling matrix calculation
 854 template<PressureCoupling pressureCouplingType>
 855 static void calculateScalingMatrixImplDetail(const t_inputrec* ir,
 856                                              matrix            mu,
 857                                              real              dt,
 858                                              const matrix      pres,
 859                                              const matrix      box,
 860                                              real              scalar_pressure,
 861                                              real              xy_pressure,
 862                                              int64_t           step);
 863
 864 //! Calculate Berendsen / C-rescale scaling matrix
 865 template<PressureCoupling pressureCouplingType>
 866 static void calculateScalingMatrixImpl(const t_inputrec* ir,
 867                                        matrix            mu,
 868                                        real              dt,
 869                                        const matrix      pres,
 870                                        const matrix      box,
 871                                        int64_t           step)
 872 {
 873     real scalar_pressure = 0;
 874     real xy_pressure     = 0;
 875     for (int d = 0; d < DIM; d++)
 876     {
 877         scalar_pressure += pres[d][d] / DIM;
 878         if (d != ZZ)
 879         {
 880             xy_pressure += pres[d][d] / (DIM - 1);
 881         }
 882     }
 883     clear_mat(mu);
 884     calculateScalingMatrixImplDetail<pressureCouplingType>(
 885             ir, mu, dt, pres, box, scalar_pressure, xy_pressure, step);
 886 }
 887
 888 template<>
 889 void calculateScalingMatrixImplDetail<PressureCoupling::Berendsen>(const t_inputrec* ir,
 890                                                                    matrix            mu,
 891                                                                    real              dt,
 892                                                                    const matrix      pres,
 893                                                                    const matrix      box,
 894                                                                    real    scalar_pressure,
 895                                                                    real    xy_pressure,
 896                                                                    int64_t gmx_unused step)
 897 {
 898     real p_corr_z = 0;
 899     switch (ir->epct)
 900     {
 901         case PressureCouplingType::Isotropic:
 902             for (int d = 0; d < DIM; d++)
 903             {
 904                 mu[d][d] = 1.0 - compressibilityFactor(d, d, ir, dt) * (ir->ref_p[d][d] - scalar_pressure) / DIM;
 905             }
 906             break;
 907         case PressureCouplingType::SemiIsotropic:
 908             for (int d = 0; d < ZZ; d++)
 909             {
 910                 mu[d][d] = 1.0 - compressibilityFactor(d, d, ir, dt) * (ir->ref_p[d][d] - xy_pressure) / DIM;
 911             }
 912             mu[ZZ][ZZ] =
 913                     1.0 - compressibilityFactor(ZZ, ZZ, ir, dt) * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) / DIM;
 914             break;
 915         case PressureCouplingType::Anisotropic:
 916             for (int d = 0; d < DIM; d++)
 917             {
 918                 for (int n = 0; n < DIM; n++)
 919                 {
 920                     mu[d][n] = (d == n ? 1.0 : 0.0)
 921                                - compressibilityFactor(d, n, ir, dt) * (ir->ref_p[d][n] - pres[d][n]) / DIM;
 922                 }
 923             }
 924             break;
 925         case PressureCouplingType::SurfaceTension:
 926             /* ir->ref_p[0/1] is the reference surface-tension times *
 927              * the number of surfaces                                */
 928             if (ir->compress[ZZ][ZZ] != 0.0F)
 929             {
 930                 p_corr_z = dt / ir->tau_p * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]);
 931             }
 932             else
 933             {
 934                 /* when the compressibity is zero, set the pressure correction   *
 935                  * in the z-direction to zero to get the correct surface tension */
 936                 p_corr_z = 0;
 937             }
 938             mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ] * p_corr_z;
 939             for (int d = 0; d < DIM - 1; d++)
 940             {
 941                 mu[d][d] = 1.0
 942                            + compressibilityFactor(d, d, ir, dt)
 943                                      * (ir->ref_p[d][d] / (mu[ZZ][ZZ] * box[ZZ][ZZ])
 944                                         - (pres[ZZ][ZZ] + p_corr_z - xy_pressure))
 945                                      / (DIM - 1);
 946             }
 947             break;
 948         default:
 949             gmx_fatal(FARGS,
 950                       "Berendsen pressure coupling type %s not supported yet\n",
 951                       enumValueToString(ir->epct));
 952     }
 953 }
 954
 955 template<>
 956 void calculateScalingMatrixImplDetail<PressureCoupling::CRescale>(const t_inputrec* ir,
 957                                                                   matrix            mu,
 958                                                                   real              dt,
 959                                                                   const matrix      pres,
 960                                                                   const matrix      box,
 961                                                                   real              scalar_pressure,
 962                                                                   real              xy_pressure,
 963                                                                   int64_t           step)
 964 {
 965     gmx::ThreeFry2x64<64>         rng(ir->ld_seed, gmx::RandomDomain::Barostat);
 966     gmx::NormalDistribution<real> normalDist;
 967     rng.restart(step, 0);
 968     real vol = 1.0;
 969     for (int d = 0; d < DIM; d++)
 970     {
 971         vol *= box[d][d];
 972     }
 973     real gauss  = 0;
 974     real gauss2 = 0;
 975     real kt     = ir->opts.ref_t[0] * BOLTZ;
 976     if (kt < 0.0)
 977     {
 978         kt = 0.0;
 979     }
 980
 981     switch (ir->epct)
 982     {
 983         case PressureCouplingType::Isotropic:
 984             gauss = normalDist(rng);
 985             for (int d = 0; d < DIM; d++)
 986             {
 987                 const real factor = compressibilityFactor(d, d, ir, dt);
 988                 mu[d][d]          = std::exp(-factor * (ir->ref_p[d][d] - scalar_pressure) / DIM
 989                                     + std::sqrt(2.0 * kt * factor * PRESFAC / vol) * gauss / DIM);
 990             }
 991             break;
 992         case PressureCouplingType::SemiIsotropic:
 993             gauss  = normalDist(rng);
 994             gauss2 = normalDist(rng);
 995             for (int d = 0; d < ZZ; d++)
 996             {
 997                 const real factor = compressibilityFactor(d, d, ir, dt);
 998                 mu[d][d]          = std::exp(-factor * (ir->ref_p[d][d] - xy_pressure) / DIM
 999                                     + std::sqrt((DIM - 1) * 2.0 * kt * factor * PRESFAC / vol / DIM)
1000                                               / (DIM - 1) * gauss);
1001             }
1002             {
1003                 const real factor = compressibilityFactor(ZZ, ZZ, ir, dt);
1004                 mu[ZZ][ZZ]        = std::exp(-factor * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) / DIM
1005                                       + std::sqrt(2.0 * kt * factor * PRESFAC / vol / DIM) * gauss2);
1006             }
1007             break;
1008         case PressureCouplingType::SurfaceTension:
1009             gauss  = normalDist(rng);
1010             gauss2 = normalDist(rng);
1011             for (int d = 0; d < ZZ; d++)
1012             {
1013                 const real factor = compressibilityFactor(d, d, ir, dt);
1014                 /* Notice: we here use ref_p[ZZ][ZZ] as isotropic pressure and ir->ref_p[d][d] as surface tension */
1015                 mu[d][d] = std::exp(
1016                         -factor * (ir->ref_p[ZZ][ZZ] - ir->ref_p[d][d] / box[ZZ][ZZ] - xy_pressure) / DIM
1017                         + std::sqrt(4.0 / 3.0 * kt * factor * PRESFAC / vol) / (DIM - 1) * gauss);
1018             }
1019             {
1020                 const real factor = compressibilityFactor(ZZ, ZZ, ir, dt);
1021                 mu[ZZ][ZZ]        = std::exp(-factor * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) / DIM
1022                                       + std::sqrt(2.0 / 3.0 * kt * factor * PRESFAC / vol) * gauss2);
1023             }
1024             break;
1025         default:
1026             gmx_fatal(FARGS,
1027                       "C-rescale pressure coupling type %s not supported yet\n",
1028                       enumValueToString(ir->epct));
1029     }
1030 }
1031
1032 template<PressureCoupling pressureCouplingType>
1033 void pressureCouplingCalculateScalingMatrix(FILE*             fplog,
1034                                             int64_t           step,
1035                                             const t_inputrec* ir,
1036                                             real              dt,
1037                                             const tensor      pres,
1038                                             const matrix      box,
1039                                             const matrix      force_vir,
1040                                             const matrix      constraint_vir,
1041                                             matrix            mu,
1042                                             double*           baros_integral)
1043 {
1044     static_assert(pressureCouplingType == PressureCoupling::Berendsen
1045                           || pressureCouplingType == PressureCoupling::CRescale,
1046                   "pressureCouplingCalculateScalingMatrix is only implemented for Berendsen and "
1047                   "C-rescale pressure coupling");
1048
1049     calculateScalingMatrixImpl<pressureCouplingType>(ir, mu, dt, pres, box, step);
1050
1051     /* To fulfill the orientation restrictions on triclinic boxes
1052      * we will set mu_yx, mu_zx and mu_zy to 0 and correct
1053      * the other elements of mu to first order.
1054      */
1055     mu[YY][XX] += mu[XX][YY];
1056     mu[ZZ][XX] += mu[XX][ZZ];
1057     mu[ZZ][YY] += mu[YY][ZZ];
1058     mu[XX][YY] = 0;
1059     mu[XX][ZZ] = 0;
1060     mu[YY][ZZ] = 0;
1061
1062     /* Keep track of the work the barostat applies on the system.
1063      * Without constraints force_vir tells us how Epot changes when scaling.
1064      * With constraints constraint_vir gives us the constraint contribution
1065      * to both Epot and Ekin. Although we are not scaling velocities, scaling
1066      * the coordinates leads to scaling of distances involved in constraints.
1067      * This in turn changes the angular momentum (even if the constrained
1068      * distances are corrected at the next step). The kinetic component
1069      * of the constraint virial captures the angular momentum change.
1070      */
1071     for (int d = 0; d < DIM; d++)
1072     {
1073         for (int n = 0; n <= d; n++)
1074         {
1075             *baros_integral -=
1076                     2 * (mu[d][n] - (n == d ? 1 : 0)) * (force_vir[d][n] + constraint_vir[d][n]);
1077         }
1078     }
1079
1080     if (debug)
1081     {
1082         pr_rvecs(debug, 0, "PC: pres ", pres, 3);
1083         pr_rvecs(debug, 0, "PC: mu   ", mu, 3);
1084     }
1085
1086     if (mu[XX][XX] < 0.99 || mu[XX][XX] > 1.01 || mu[YY][YY] < 0.99 || mu[YY][YY] > 1.01
1087         || mu[ZZ][ZZ] < 0.99 || mu[ZZ][ZZ] > 1.01)
1088     {
1089         char buf[STRLEN];
1090         char buf2[22];
1091         sprintf(buf,
1092                 "\nStep %s  Warning: pressure scaling more than 1%%, "
1093                 "mu: %g %g %g\n",
1094                 gmx_step_str(step, buf2),
1095                 mu[XX][XX],
1096                 mu[YY][YY],
1097                 mu[ZZ][ZZ]);
1098         if (fplog)
1099         {
1100             fprintf(fplog, "%s", buf);
1101         }
1102         fprintf(stderr, "%s", buf);
1103     }
1104 }
1105
1106 template<PressureCoupling pressureCouplingType>
1107 void pressureCouplingScaleBoxAndCoordinates(const t_inputrec*    ir,
1108                                             const matrix         mu,
1109                                             matrix               box,
1110                                             matrix               box_rel,
1111                                             int                  start,
1112                                             int                  nr_atoms,
1113                                             rvec                 x[],
1114                                             rvec                 v[],
1115                                             const unsigned short cFREEZE[],
1116                                             t_nrnb*              nrnb,
1117                                             const bool           scaleCoordinates)
1118 {
1119     static_assert(pressureCouplingType == PressureCoupling::Berendsen
1120                           || pressureCouplingType == PressureCoupling::CRescale,
1121                   "pressureCouplingScaleBoxAndCoordinates is only implemented for Berendsen and "
1122                   "C-rescale pressure coupling");
1123
1124     ivec*  nFreeze = ir->opts.nFreeze;
1125     matrix inv_mu;
1126     if (pressureCouplingType == PressureCoupling::CRescale)
1127     {
1128         gmx::invertBoxMatrix(mu, inv_mu);
1129     }
1130
1131     /* Scale the positions and the velocities */
1132     if (scaleCoordinates)
1133     {
1134         const int gmx_unused numThreads = gmx_omp_nthreads_get(emntUpdate);
1135 #pragma omp parallel for num_threads(numThreads) schedule(static)
1136         for (int n = start; n < start + nr_atoms; n++)
1137         {
1138             // Trivial OpenMP region that does not throw
1139             int g = 0;
1140             if (cFREEZE != nullptr)
1141             {
1142                 g = cFREEZE[n];
1143             }
1144
1145             if (!nFreeze[g][XX])
1146             {
1147                 x[n][XX] = mu[XX][XX] * x[n][XX] + mu[YY][XX] * x[n][YY] + mu[ZZ][XX] * x[n][ZZ];
1148                 if (pressureCouplingType == PressureCoupling::CRescale)
1149                 {
1150                     v[n][XX] = inv_mu[XX][XX] * v[n][XX] + inv_mu[YY][XX] * v[n][YY]
1151                                + inv_mu[ZZ][XX] * v[n][ZZ];
1152                 }
1153             }
1154             if (!nFreeze[g][YY])
1155             {
1156                 x[n][YY] = mu[YY][YY] * x[n][YY] + mu[ZZ][YY] * x[n][ZZ];
1157                 if (pressureCouplingType == PressureCoupling::CRescale)
1158                 {
1159                     v[n][YY] = inv_mu[YY][YY] * v[n][YY] + inv_mu[ZZ][YY] * v[n][ZZ];
1160                 }
1161             }
1162             if (!nFreeze[g][ZZ])
1163             {
1164                 x[n][ZZ] = mu[ZZ][ZZ] * x[n][ZZ];
1165                 if (pressureCouplingType == PressureCoupling::CRescale)
1166                 {
1167                     v[n][ZZ] = inv_mu[ZZ][ZZ] * v[n][ZZ];
1168                 }
1169             }
1170         }
1171     }
1172     /* compute final boxlengths */
1173     for (int d = 0; d < DIM; d++)
1174     {
1175         box[d][XX] = mu[XX][XX] * box[d][XX] + mu[YY][XX] * box[d][YY] + mu[ZZ][XX] * box[d][ZZ];
1176         box[d][YY] = mu[YY][YY] * box[d][YY] + mu[ZZ][YY] * box[d][ZZ];
1177         box[d][ZZ] = mu[ZZ][ZZ] * box[d][ZZ];
1178     }
1179
1180     preserve_box_shape(ir, box_rel, box);
1181
1182     /* (un)shifting should NOT be done after this,
1183      * since the box vectors might have changed
1184      */
1185     inc_nrnb(nrnb, eNR_PCOUPL, nr_atoms);
1186 }
1187
1188 void berendsen_tcoupl(const t_inputrec* ir, gmx_ekindata_t* ekind, real dt, std::vector<double>& therm_integral)
1189 {
1190     const t_grpopts* opts = &ir->opts;
1191
1192     for (int i = 0; (i < opts->ngtc); i++)
1193     {
1194         real Ek, T;
1195
1196         if (ir->eI == IntegrationAlgorithm::VV)
1197         {
1198             Ek = trace(ekind->tcstat[i].ekinf);
1199             T  = ekind->tcstat[i].T;
1200         }
1201         else
1202         {
1203             Ek = trace(ekind->tcstat[i].ekinh);
1204             T  = ekind->tcstat[i].Th;
1205         }
1206
1207         if ((opts->tau_t[i] > 0) && (T > 0.0))
1208         {
1209             real reft               = std::max<real>(0, opts->ref_t[i]);
1210             real lll                = std::sqrt(1.0 + (dt / opts->tau_t[i]) * (reft / T - 1.0));
1211             ekind->tcstat[i].lambda = std::max<real>(std::min<real>(lll, 1.25), 0.8);
1212         }
1213         else
1214         {
1215             ekind->tcstat[i].lambda = 1.0;
1216         }
1217
1218         /* Keep track of the amount of energy we are adding to the system */
1219         therm_integral[i] -= (gmx::square(ekind->tcstat[i].lambda) - 1) * Ek;
1220
1221         if (debug)
1222         {
1223             fprintf(debug, "TC: group %d: T: %g, Lambda: %g\n", i, T, ekind->tcstat[i].lambda);
1224         }
1225     }
1226 }
1227
1228 void andersen_tcoupl(const t_inputrec*         ir,
1229                      int64_t                   step,
1230                      const t_commrec*          cr,
1231                      const t_mdatoms*          md,
1232                      gmx::ArrayRef<gmx::RVec>  v,
1233                      real                      rate,
1234                      const std::vector<bool>&  randomize,
1235                      gmx::ArrayRef<const real> boltzfac)
1236 {
1237     const int*           gatindex = (DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
1238     int                  i;
1239     int                  gc = 0;
1240     gmx::ThreeFry2x64<0> rng(ir->andersen_seed, gmx::RandomDomain::Thermostat);
1241     gmx::UniformRealDistribution<real>         uniformDist;
1242     gmx::TabulatedNormalDistribution<real, 14> normalDist;
1243
1244     /* randomize the velocities of the selected particles */
1245
1246     for (i = 0; i < md->homenr; i++) /* now loop over the list of atoms */
1247     {
1248         int      ng = gatindex ? gatindex[i] : i;
1249         gmx_bool bRandomize;
1250
1251         rng.restart(step, ng);
1252
1253         if (md->cTC)
1254         {
1255             gc = md->cTC[i]; /* assign the atom to a temperature group if there are more than one */
1256         }
1257         if (randomize[gc])
1258         {
1259             if (ir->etc == TemperatureCoupling::AndersenMassive)
1260             {
1261                 /* Randomize particle always */
1262                 bRandomize = TRUE;
1263             }
1264             else
1265             {
1266                 /* Randomize particle probabilistically */
1267                 uniformDist.reset();
1268                 bRandomize = uniformDist(rng) < rate;
1269             }
1270             if (bRandomize)
1271             {
1272                 real scal;
1273                 int  d;
1274
1275                 scal = std::sqrt(boltzfac[gc] * md->invmass[i]);
1276
1277                 normalDist.reset();
1278
1279                 for (d = 0; d < DIM; d++)
1280                 {
1281                     v[i][d] = scal * normalDist(rng);
1282                 }
1283             }
1284         }
1285     }
1286 }
1287
1288
1289 void nosehoover_tcoupl(const t_grpopts*      opts,
1290                        const gmx_ekindata_t* ekind,
1291                        real                  dt,
1292                        double                xi[],
1293                        double                vxi[],
1294                        const t_extmass*      MassQ)
1295 {
1296     int  i;
1297     real reft, oldvxi;
1298
1299     /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
1300
1301     for (i = 0; (i < opts->ngtc); i++)
1302     {
1303         reft   = std::max<real>(0, opts->ref_t[i]);
1304         oldvxi = vxi[i];
1305         vxi[i] += dt * MassQ->Qinv[i] * (ekind->tcstat[i].Th - reft);
1306         xi[i] += dt * (oldvxi + vxi[i]) * 0.5;
1307     }
1308 }
1309
1310 void trotter_update(const t_inputrec*               ir,
1311                     int64_t                         step,
1312                     gmx_ekindata_t*                 ekind,
1313                     const gmx_enerdata_t*           enerd,
1314                     t_state*                        state,
1315                     const tensor                    vir,
1316                     const t_mdatoms*                md,
1317                     const t_extmass*                MassQ,
1318                     gmx::ArrayRef<std::vector<int>> trotter_seqlist,
1319                     int                             trotter_seqno)
1320 {
1321
1322     int              n, i, d, ngtc, gc = 0, t;
1323     t_grp_tcstat*    tcstat;
1324     const t_grpopts* opts;
1325     int64_t          step_eff;
1326     real             dt;
1327     double *         scalefac, dtc;
1328     rvec             sumv = { 0, 0, 0 };
1329     gmx_bool         bCouple;
1330
1331     if (trotter_seqno <= ettTSEQ2)
1332     {
1333         step_eff = step - 1; /* the velocity verlet calls are actually out of order -- the first
1334                                 half step is actually the last half step from the previous step.
1335                                 Thus the first half step actually corresponds to the n-1 step*/
1336     }
1337     else
1338     {
1339         step_eff = step;
1340     }
1341
1342     bCouple = (ir->nsttcouple == 1 || do_per_step(step_eff + ir->nsttcouple, ir->nsttcouple));
1343
1344     const gmx::ArrayRef<const int> trotter_seq = trotter_seqlist[trotter_seqno];
1345
1346     if ((trotter_seq[0] == etrtSKIPALL) || (!bCouple))
1347     {
1348         return;
1349     }
1350     dtc  = ir->nsttcouple * ir->delta_t; /* This is OK for NPT, because nsttcouple == nstpcouple is enforcesd */
1351     opts = &(ir->opts);                  /* just for ease of referencing */
1352     ngtc = opts->ngtc;
1353     assert(ngtc > 0);
1354     snew(scalefac, opts->ngtc);
1355     for (i = 0; i < ngtc; i++)
1356     {
1357         scalefac[i] = 1;
1358     }
1359     /* execute the series of trotter updates specified in the trotterpart array */
1360
1361     for (i = 0; i < NTROTTERPARTS; i++)
1362     {
1363         /* allow for doubled intgrators by doubling dt instead of making 2 calls */
1364         if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2)
1365             || (trotter_seq[i] == etrtNHC2))
1366         {
1367             dt = 2 * dtc;
1368         }
1369         else
1370         {
1371             dt = dtc;
1372         }
1373
1374         auto v = makeArrayRef(state->v);
1375         switch (trotter_seq[i])
1376         {
1377             case etrtBAROV:
1378             case etrtBAROV2:
1379                 boxv_trotter(ir, &(state->veta), dt, state->box, ekind, vir, enerd->term[F_PDISPCORR], MassQ);
1380                 break;
1381             case etrtBARONHC:
1382             case etrtBARONHC2:
1383                 NHC_trotter(opts,
1384                             state->nnhpres,
1385                             ekind,
1386                             dt,
1387                             state->nhpres_xi.data(),
1388                             state->nhpres_vxi.data(),
1389                             nullptr,
1390                             &(state->veta),
1391                             MassQ,
1392                             FALSE);
1393                 break;
1394             case etrtNHC:
1395             case etrtNHC2:
1396                 NHC_trotter(opts,
1397                             opts->ngtc,
1398                             ekind,
1399                             dt,
1400                             state->nosehoover_xi.data(),
1401                             state->nosehoover_vxi.data(),
1402                             scalefac,
1403                             nullptr,
1404                             MassQ,
1405                             (ir->eI == IntegrationAlgorithm::VV));
1406                 /* need to rescale the kinetic energies and velocities here.  Could
1407                    scale the velocities later, but we need them scaled in order to
1408                    produce the correct outputs, so we'll scale them here. */
1409
1410                 for (t = 0; t < ngtc; t++)
1411                 {
1412                     tcstat             = &ekind->tcstat[t];
1413                     tcstat->vscale_nhc = scalefac[t];
1414                     tcstat->ekinscaleh_nhc *= (scalefac[t] * scalefac[t]);
1415                     tcstat->ekinscalef_nhc *= (scalefac[t] * scalefac[t]);
1416                 }
1417                 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
1418                 /* but do we actually need the total? */
1419
1420                 /* modify the velocities as well */
1421                 for (n = 0; n < md->homenr; n++)
1422                 {
1423                     if (md->cTC) /* does this conditional need to be here? is this always true?*/
1424                     {
1425                         gc = md->cTC[n];
1426                     }
1427                     for (d = 0; d < DIM; d++)
1428                     {
1429                         v[n][d] *= scalefac[gc];
1430                     }
1431
1432                     if (debug)
1433                     {
1434                         for (d = 0; d < DIM; d++)
1435                         {
1436                             sumv[d] += (v[n][d]) / md->invmass[n];
1437                         }
1438                     }
1439                 }
1440                 break;
1441             default: break;
1442         }
1443     }
1444     /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
1445     sfree(scalefac);
1446 }
1447
1448
1449 extern void init_npt_masses(const t_inputrec* ir, t_state* state, t_extmass* MassQ, gmx_bool bInit)
1450 {
1451     int              n, i, j, d, ngtc, nh;
1452     const t_grpopts* opts;
1453     real             reft, kT, ndj, nd;
1454
1455     opts = &(ir->opts); /* just for ease of referencing */
1456     ngtc = ir->opts.ngtc;
1457     nh   = state->nhchainlength;
1458
1459     if (ir->eI == IntegrationAlgorithm::MD)
1460     {
1461         if (bInit)
1462         {
1463             MassQ->Qinv.resize(ngtc);
1464         }
1465         for (i = 0; (i < ngtc); i++)
1466         {
1467             if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
1468             {
1469                 MassQ->Qinv[i] = 1.0 / (gmx::square(opts->tau_t[i] / M_2PI) * opts->ref_t[i]);
1470             }
1471             else
1472             {
1473                 MassQ->Qinv[i] = 0.0;
1474             }
1475         }
1476     }
1477     else if (EI_VV(ir->eI))
1478     {
1479         /* Set pressure variables */
1480
1481         if (bInit)
1482         {
1483             if (state->vol0 == 0)
1484             {
1485                 state->vol0 = det(state->box);
1486                 /* because we start by defining a fixed
1487                    compressibility, we need the volume at this
1488                    compressibility to solve the problem. */
1489             }
1490         }
1491
1492         /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol  */
1493         /* Consider evaluating eventually if this the right mass to use.  All are correct, some might be more stable  */
1494         MassQ->Winv = (PRESFAC * trace(ir->compress) * BOLTZ * opts->ref_t[0])
1495                       / (DIM * state->vol0 * gmx::square(ir->tau_p / M_2PI));
1496         /* An alternate mass definition, from Tuckerman et al. */
1497         /* MassQ->Winv = 1.0/(gmx::square(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
1498         for (d = 0; d < DIM; d++)
1499         {
1500             for (n = 0; n < DIM; n++)
1501             {
1502                 MassQ->Winvm[d][n] =
1503                         PRESFAC * ir->compress[d][n] / (state->vol0 * gmx::square(ir->tau_p / M_2PI));
1504                 /* not clear this is correct yet for the anisotropic case. Will need to reevaluate
1505                    before using MTTK for anisotropic states.*/
1506             }
1507         }
1508         /* Allocate space for thermostat variables */
1509         if (bInit)
1510         {
1511             MassQ->Qinv.resize(ngtc * nh);
1512         }
1513
1514         /* now, set temperature variables */
1515         for (i = 0; i < ngtc; i++)
1516         {
1517             if (opts->tau_t[i] > 0 && opts->ref_t[i] > 0 && opts->nrdf[i] > 0)
1518             {
1519                 reft = std::max<real>(0, opts->ref_t[i]);
1520                 nd   = opts->nrdf[i];
1521                 kT   = BOLTZ * reft;
1522                 for (j = 0; j < nh; j++)
1523                 {
1524                     if (j == 0)
1525                     {
1526                         ndj = nd;
1527                     }
1528                     else
1529                     {
1530                         ndj = 1;
1531                     }
1532                     MassQ->Qinv[i * nh + j] = 1.0 / (gmx::square(opts->tau_t[i] / M_2PI) * ndj * kT);
1533                 }
1534             }
1535             else
1536             {
1537                 for (j = 0; j < nh; j++)
1538                 {
1539                     MassQ->Qinv[i * nh + j] = 0.0;
1540                 }
1541             }
1542         }
1543     }
1544 }
1545
1546 std::array<std::vector<int>, ettTSEQMAX>
1547 init_npt_vars(const t_inputrec* ir, t_state* state, t_extmass* MassQ, gmx_bool bTrotter)
1548 {
1549     int              i, j, nnhpres, nh;
1550     const t_grpopts* opts;
1551     real             bmass, qmass, reft, kT;
1552
1553     opts    = &(ir->opts); /* just for ease of referencing */
1554     nnhpres = state->nnhpres;
1555     nh      = state->nhchainlength;
1556
1557     if (EI_VV(ir->eI) && (ir->epc == PressureCoupling::Mttk) && (ir->etc != TemperatureCoupling::NoseHoover))
1558     {
1559         gmx_fatal(FARGS, "Cannot do MTTK pressure coupling without Nose-Hoover temperature control");
1560     }
1561
1562     init_npt_masses(ir, state, MassQ, TRUE);
1563
1564     /* first, initialize clear all the trotter calls */
1565     std::array<std::vector<int>, ettTSEQMAX> trotter_seq;
1566     for (i = 0; i < ettTSEQMAX; i++)
1567     {
1568         trotter_seq[i].resize(NTROTTERPARTS, etrtNONE);
1569         trotter_seq[i][0] = etrtSKIPALL;
1570     }
1571
1572     if (!bTrotter)
1573     {
1574         /* no trotter calls, so we never use the values in the array.
1575          * We access them (so we need to define them, but ignore
1576          * then.*/
1577
1578         return trotter_seq;
1579     }
1580
1581     /* compute the kinetic energy by using the half step velocities or
1582      * the kinetic energies, depending on the order of the trotter calls */
1583
1584     if (ir->eI == IntegrationAlgorithm::VV)
1585     {
1586         if (inputrecNptTrotter(ir))
1587         {
1588             /* This is the complicated version - there are 4 possible calls, depending on ordering.
1589                We start with the initial one. */
1590             /* first, a round that estimates veta. */
1591             trotter_seq[0][0] = etrtBAROV;
1592
1593             /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1594
1595             /* The first half trotter update */
1596             trotter_seq[2][0] = etrtBAROV;
1597             trotter_seq[2][1] = etrtNHC;
1598             trotter_seq[2][2] = etrtBARONHC;
1599
1600             /* The second half trotter update */
1601             trotter_seq[3][0] = etrtBARONHC;
1602             trotter_seq[3][1] = etrtNHC;
1603             trotter_seq[3][2] = etrtBAROV;
1604
1605             /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1606         }
1607         else if (inputrecNvtTrotter(ir))
1608         {
1609             /* This is the easy version - there are only two calls, both the same.
1610                Otherwise, even easier -- no calls  */
1611             trotter_seq[2][0] = etrtNHC;
1612             trotter_seq[3][0] = etrtNHC;
1613         }
1614         else if (inputrecNphTrotter(ir))
1615         {
1616             /* This is the complicated version - there are 4 possible calls, depending on ordering.
1617                We start with the initial one. */
1618             /* first, a round that estimates veta. */
1619             trotter_seq[0][0] = etrtBAROV;
1620
1621             /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1622
1623             /* The first half trotter update */
1624             trotter_seq[2][0] = etrtBAROV;
1625             trotter_seq[2][1] = etrtBARONHC;
1626
1627             /* The second half trotter update */
1628             trotter_seq[3][0] = etrtBARONHC;
1629             trotter_seq[3][1] = etrtBAROV;
1630
1631             /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1632         }
1633     }
1634     else if (ir->eI == IntegrationAlgorithm::VVAK)
1635     {
1636         if (inputrecNptTrotter(ir))
1637         {
1638             /* This is the complicated version - there are 4 possible calls, depending on ordering.
1639                We start with the initial one. */
1640             /* first, a round that estimates veta. */
1641             trotter_seq[0][0] = etrtBAROV;
1642
1643             /* The first half trotter update, part 1 -- double update, because it commutes */
1644             trotter_seq[1][0] = etrtNHC;
1645
1646             /* The first half trotter update, part 2 */
1647             trotter_seq[2][0] = etrtBAROV;
1648             trotter_seq[2][1] = etrtBARONHC;
1649
1650             /* The second half trotter update, part 1 */
1651             trotter_seq[3][0] = etrtBARONHC;
1652             trotter_seq[3][1] = etrtBAROV;
1653
1654             /* The second half trotter update */
1655             trotter_seq[4][0] = etrtNHC;
1656         }
1657         else if (inputrecNvtTrotter(ir))
1658         {
1659             /* This is the easy version - there is only one call, both the same.
1660                Otherwise, even easier -- no calls  */
1661             trotter_seq[1][0] = etrtNHC;
1662             trotter_seq[4][0] = etrtNHC;
1663         }
1664         else if (inputrecNphTrotter(ir))
1665         {
1666             /* This is the complicated version - there are 4 possible calls, depending on ordering.
1667                We start with the initial one. */
1668             /* first, a round that estimates veta. */
1669             trotter_seq[0][0] = etrtBAROV;
1670
1671             /* The first half trotter update, part 1 -- leave zero */
1672             trotter_seq[1][0] = etrtNHC;
1673
1674             /* The first half trotter update, part 2 */
1675             trotter_seq[2][0] = etrtBAROV;
1676             trotter_seq[2][1] = etrtBARONHC;
1677
1678             /* The second half trotter update, part 1 */
1679             trotter_seq[3][0] = etrtBARONHC;
1680             trotter_seq[3][1] = etrtBAROV;
1681
1682             /* The second half trotter update -- blank for now */
1683         }
1684     }
1685
1686     switch (ir->epct)
1687     {
1688         case PressureCouplingType::Isotropic:
1689         default: bmass = DIM * DIM; /* recommended mass parameters for isotropic barostat */
1690     }
1691
1692     MassQ->QPinv.resize(nnhpres * opts->nhchainlength);
1693
1694     /* barostat temperature */
1695     if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
1696     {
1697         reft = std::max<real>(0, opts->ref_t[0]);
1698         kT   = BOLTZ * reft;
1699         for (i = 0; i < nnhpres; i++)
1700         {
1701             for (j = 0; j < nh; j++)
1702             {
1703                 if (j == 0)
1704                 {
1705                     qmass = bmass;
1706                 }
1707                 else
1708                 {
1709                     qmass = 1;
1710                 }
1711                 MassQ->QPinv[i * opts->nhchainlength + j] =
1712                         1.0 / (gmx::square(opts->tau_t[0] / M_2PI) * qmass * kT);
1713             }
1714         }
1715     }
1716     else
1717     {
1718         for (i = 0; i < nnhpres; i++)
1719         {
1720             for (j = 0; j < nh; j++)
1721             {
1722                 MassQ->QPinv[i * nh + j] = 0.0;
1723             }
1724         }
1725     }
1726     return trotter_seq;
1727 }
1728
1729 static real energyNoseHoover(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ)
1730 {
1731     real energy = 0;
1732
1733     int nh = state->nhchainlength;
1734
1735     for (int i = 0; i < ir->opts.ngtc; i++)
1736     {
1737         const double* ixi   = &state->nosehoover_xi[i * nh];
1738         const double* ivxi  = &state->nosehoover_vxi[i * nh];
1739         const double* iQinv = &(MassQ->Qinv[i * nh]);
1740
1741         real nd   = ir->opts.nrdf[i];
1742         real reft = std::max<real>(ir->opts.ref_t[i], 0);
1743         real kT   = BOLTZ * reft;
1744
1745         if (nd > 0.0)
1746         {
1747             if (inputrecNvtTrotter(ir) || inputrecNptTrotter(ir))
1748             {
1749                 /* contribution from the thermal momenta of the NH chain */
1750                 for (int j = 0; j < nh; j++)
1751                 {
1752                     if (iQinv[j] > 0)
1753                     {
1754                         energy += 0.5 * gmx::square(ivxi[j]) / iQinv[j];
1755                         /* contribution from the thermal variable of the NH chain */
1756                         real ndj = 0;
1757                         if (j == 0)
1758                         {
1759                             ndj = nd;
1760                         }
1761                         else
1762                         {
1763                             ndj = 1;
1764                         }
1765                         energy += ndj * ixi[j] * kT;
1766                     }
1767                 }
1768             }
1769             else /* Other non Trotter temperature NH control  -- no chains yet. */
1770             {
1771                 energy += 0.5 * BOLTZ * nd * gmx::square(ivxi[0]) / iQinv[0];
1772                 energy += nd * ixi[0] * kT;
1773             }
1774         }
1775     }
1776
1777     return energy;
1778 }
1779
1780 /* Returns the energy from the barostat thermostat chain */
1781 static real energyPressureMTTK(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ)
1782 {
1783     real energy = 0;
1784
1785     int nh = state->nhchainlength;
1786
1787     for (int i = 0; i < state->nnhpres; i++)
1788     {
1789         /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1790         real reft = std::max<real>(ir->opts.ref_t[0], 0.0); /* using 'System' temperature */
1791         real kT   = BOLTZ * reft;
1792
1793         for (int j = 0; j < nh; j++)
1794         {
1795             double iQinv = MassQ->QPinv[i * nh + j];
1796             if (iQinv > 0)
1797             {
1798                 energy += 0.5 * gmx::square(state->nhpres_vxi[i * nh + j]) / iQinv;
1799                 /* contribution from the thermal variable of the NH chain */
1800                 energy += state->nhpres_xi[i * nh + j] * kT;
1801             }
1802             if (debug)
1803             {
1804                 fprintf(debug,
1805                         "P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f",
1806                         i,
1807                         j,
1808                         state->nhpres_vxi[i * nh + j],
1809                         state->nhpres_xi[i * nh + j]);
1810             }
1811         }
1812     }
1813
1814     return energy;
1815 }
1816
1817 /* Returns the energy accumulated by the V-rescale or Berendsen thermostat */
1818 static real energyVrescale(const t_inputrec* ir, const t_state* state)
1819 {
1820     real energy = 0;
1821     for (int i = 0; i < ir->opts.ngtc; i++)
1822     {
1823         energy += state->therm_integral[i];
1824     }
1825
1826     return energy;
1827 }
1828
1829 real NPT_energy(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ)
1830 {
1831     real energyNPT = 0;
1832
1833     if (ir->epc != PressureCoupling::No)
1834     {
1835         /* Compute the contribution of the pressure to the conserved quantity*/
1836
1837         real vol = det(state->box);
1838
1839         switch (ir->epc)
1840         {
1841             case PressureCoupling::ParrinelloRahman:
1842             {
1843                 /* contribution from the pressure momenta */
1844                 tensor invMass;
1845                 calcParrinelloRahmanInvMass(ir, state->box, invMass);
1846                 for (int d = 0; d < DIM; d++)
1847                 {
1848                     for (int n = 0; n <= d; n++)
1849                     {
1850                         if (invMass[d][n] > 0)
1851                         {
1852                             energyNPT += 0.5 * gmx::square(state->boxv[d][n]) / (invMass[d][n] * PRESFAC);
1853                         }
1854                     }
1855                 }
1856
1857                 /* Contribution from the PV term.
1858                  * Not that with non-zero off-diagonal reference pressures,
1859                  * i.e. applied shear stresses, there are additional terms.
1860                  * We don't support this here, since that requires keeping
1861                  * track of unwrapped box diagonal elements. This case is
1862                  * excluded in integratorHasConservedEnergyQuantity().
1863                  */
1864                 energyNPT += vol * trace(ir->ref_p) / (DIM * PRESFAC);
1865                 break;
1866             }
1867             case PressureCoupling::Mttk:
1868                 /* contribution from the pressure momenta */
1869                 energyNPT += 0.5 * gmx::square(state->veta) / MassQ->Winv;
1870
1871                 /* contribution from the PV term */
1872                 energyNPT += vol * trace(ir->ref_p) / (DIM * PRESFAC);
1873
1874                 if (ir->epc == PressureCoupling::Mttk)
1875                 {
1876                     /* contribution from the MTTK chain */
1877                     energyNPT += energyPressureMTTK(ir, state, MassQ);
1878                 }
1879                 break;
1880             case PressureCoupling::Berendsen:
1881             case PressureCoupling::CRescale: energyNPT += state->baros_integral; break;
1882             default:
1883                 GMX_RELEASE_ASSERT(
1884                         false,
1885                         "Conserved energy quantity for pressure coupling is not handled. A case "
1886                         "should be added with either the conserved quantity added or nothing added "
1887                         "and an exclusion added to integratorHasConservedEnergyQuantity().");
1888         }
1889     }
1890
1891     switch (ir->etc)
1892     {
1893         case TemperatureCoupling::No: break;
1894         case TemperatureCoupling::VRescale:
1895         case TemperatureCoupling::Berendsen: energyNPT += energyVrescale(ir, state); break;
1896         case TemperatureCoupling::NoseHoover:
1897             energyNPT += energyNoseHoover(ir, state, MassQ);
1898             break;
1899         case TemperatureCoupling::Andersen:
1900         case TemperatureCoupling::AndersenMassive:
1901             // Not supported, excluded in integratorHasConservedEnergyQuantity()
1902             break;
1903         default:
1904             GMX_RELEASE_ASSERT(
1905                     false,
1906                     "Conserved energy quantity for temperature coupling is not handled. A case "
1907                     "should be added with either the conserved quantity added or nothing added and "
1908                     "an exclusion added to integratorHasConservedEnergyQuantity().");
1909     }
1910
1911     return energyNPT;
1912 }
1913
1914
1915 static real vrescale_sumnoises(real nn, gmx::ThreeFry2x64<>* rng, gmx::NormalDistribution<real>* normalDist)
1916 {
1917     /*
1918      * Returns the sum of nn independent gaussian noises squared
1919      * (i.e. equivalent to summing the square of the return values
1920      * of nn calls to a normal distribution).
1921      */
1922     const real                   ndeg_tol = 0.0001;
1923     real                         r;
1924     gmx::GammaDistribution<real> gammaDist(0.5 * nn, 1.0);
1925
1926     if (nn < 2 + ndeg_tol)
1927     {
1928         int  nn_int, i;
1929         real gauss;
1930
1931         nn_int = gmx::roundToInt(nn);
1932
1933         if (nn - nn_int < -ndeg_tol || nn - nn_int > ndeg_tol)
1934         {
1935             gmx_fatal(FARGS,
1936                       "The v-rescale thermostat was called with a group with #DOF=%f, but for "
1937                       "#DOF<3 only integer #DOF are supported",
1938                       nn + 1);
1939         }
1940
1941         r = 0;
1942         for (i = 0; i < nn_int; i++)
1943         {
1944             gauss = (*normalDist)(*rng);
1945             r += gauss * gauss;
1946         }
1947     }
1948     else
1949     {
1950         /* Use a gamma distribution for any real nn > 2 */
1951         r = 2.0 * gammaDist(*rng);
1952     }
1953
1954     return r;
1955 }
1956
1957 real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut, int64_t step, int64_t seed)
1958 {
1959     /*
1960      * Generates a new value for the kinetic energy,
1961      * according to Bussi et al JCP (2007), Eq. (A7)
1962      * kk:    present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1963      * sigma: target average value of the kinetic energy (ndeg k_b T/2)  (in the same units as kk)
1964      * ndeg:  number of degrees of freedom of the atoms to be thermalized
1965      * taut:  relaxation time of the thermostat, in units of 'how often this routine is called'
1966      */
1967     real                          factor, rr, ekin_new;
1968     gmx::ThreeFry2x64<64>         rng(seed, gmx::RandomDomain::Thermostat);
1969     gmx::NormalDistribution<real> normalDist;
1970
1971     if (taut > 0.1)
1972     {
1973         factor = exp(-1.0 / taut);
1974     }
1975     else
1976     {
1977         factor = 0.0;
1978     }
1979
1980     rng.restart(step, 0);
1981
1982     rr = normalDist(rng);
1983
1984     ekin_new = kk
1985                + (1.0 - factor)
1986                          * (sigma * (vrescale_sumnoises(ndeg - 1, &rng, &normalDist) + rr * rr) / ndeg - kk)
1987                + 2.0 * rr * std::sqrt(kk * sigma / ndeg * (1.0 - factor) * factor);
1988
1989     return ekin_new;
1990 }
1991
1992 void vrescale_tcoupl(const t_inputrec* ir, int64_t step, gmx_ekindata_t* ekind, real dt, double therm_integral[])
1993 {
1994     const t_grpopts* opts;
1995     int              i;
1996     real             Ek, Ek_ref1, Ek_ref, Ek_new;
1997
1998     opts = &ir->opts;
1999
2000     for (i = 0; (i < opts->ngtc); i++)
2001     {
2002         if (ir->eI == IntegrationAlgorithm::VV)
2003         {
2004             Ek = trace(ekind->tcstat[i].ekinf);
2005         }
2006         else
2007         {
2008             Ek = trace(ekind->tcstat[i].ekinh);
2009         }
2010
2011         if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0)
2012         {
2013             Ek_ref1 = 0.5 * opts->ref_t[i] * BOLTZ;
2014             Ek_ref  = Ek_ref1 * opts->nrdf[i];
2015
2016             Ek_new = vrescale_resamplekin(Ek, Ek_ref, opts->nrdf[i], opts->tau_t[i] / dt, step, ir->ld_seed);
2017
2018             /* Analytically Ek_new>=0, but we check for rounding errors */
2019             if (Ek_new <= 0)
2020             {
2021                 ekind->tcstat[i].lambda = 0.0;
2022             }
2023             else
2024             {
2025                 ekind->tcstat[i].lambda = std::sqrt(Ek_new / Ek);
2026             }
2027
2028             therm_integral[i] -= Ek_new - Ek;
2029
2030             if (debug)
2031             {
2032                 fprintf(debug,
2033                         "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
2034                         i,
2035                         Ek_ref,
2036                         Ek,
2037                         Ek_new,
2038                         ekind->tcstat[i].lambda);
2039             }
2040         }
2041         else
2042         {
2043             ekind->tcstat[i].lambda = 1.0;
2044         }
2045     }
2046 }
2047
2048 void rescale_velocities(const gmx_ekindata_t* ekind, const t_mdatoms* mdatoms, int start, int end, rvec v[])
2049 {
2050     const unsigned short*             cTC    = mdatoms->cTC;
2051     gmx::ArrayRef<const t_grp_tcstat> tcstat = ekind->tcstat;
2052
2053     for (int n = start; n < end; n++)
2054     {
2055         int gt = 0;
2056         if (cTC)
2057         {
2058             gt = cTC[n];
2059         }
2060         const real lg = tcstat[gt].lambda;
2061         for (int d = 0; d < DIM; d++)
2062         {
2063             v[n][d] *= lg;
2064         }
2065     }
2066 }
2067
2068 //! Check whether we're doing simulated annealing
2069 bool doSimulatedAnnealing(const t_inputrec* ir)
2070 {
2071     for (int i = 0; i < ir->opts.ngtc; i++)
2072     {
2073         /* set bSimAnn if any group is being annealed */
2074         if (ir->opts.annealing[i] != SimulatedAnnealing::No)
2075         {
2076             return true;
2077         }
2078     }
2079     return false;
2080 }
2081
2082 // TODO If we keep simulated annealing, make a proper module that
2083 // does not rely on changing inputrec.
2084 bool initSimulatedAnnealing(t_inputrec* ir, gmx::Update* upd)
2085 {
2086     bool doSimAnnealing = doSimulatedAnnealing(ir);
2087     if (doSimAnnealing)
2088     {
2089         update_annealing_target_temp(ir, ir->init_t, upd);
2090     }
2091     return doSimAnnealing;
2092 }
2093
2094 /* set target temperatures if we are annealing */
2095 void update_annealing_target_temp(t_inputrec* ir, real t, gmx::Update* upd)
2096 {
2097     int  i, j, n, npoints;
2098     real pert, thist = 0, x;
2099
2100     for (i = 0; i < ir->opts.ngtc; i++)
2101     {
2102         npoints = ir->opts.anneal_npoints[i];
2103         switch (ir->opts.annealing[i])
2104         {
2105             case SimulatedAnnealing::No: continue;
2106             case SimulatedAnnealing::Periodic:
2107                 /* calculate time modulo the period */
2108                 pert  = ir->opts.anneal_time[i][npoints - 1];
2109                 n     = static_cast<int>(t / pert);
2110                 thist = t - n * pert; /* modulo time */
2111                 /* Make sure rounding didn't get us outside the interval */
2112                 if (std::fabs(thist - pert) < GMX_REAL_EPS * 100)
2113                 {
2114                     thist = 0;
2115                 }
2116                 break;
2117             case SimulatedAnnealing::Single: thist = t; break;
2118             default:
2119                 gmx_fatal(FARGS,
2120                           "Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)",
2121                           i,
2122                           ir->opts.ngtc,
2123                           npoints);
2124         }
2125         /* We are doing annealing for this group if we got here,
2126          * and we have the (relative) time as thist.
2127          * calculate target temp */
2128         j = 0;
2129         while ((j < npoints - 1) && (thist > (ir->opts.anneal_time[i][j + 1])))
2130         {
2131             j++;
2132         }
2133         if (j < npoints - 1)
2134         {
2135             /* Found our position between points j and j+1.
2136              * Interpolate: x is the amount from j+1, (1-x) from point j
2137              * First treat possible jumps in temperature as a special case.
2138              */
2139             if ((ir->opts.anneal_time[i][j + 1] - ir->opts.anneal_time[i][j]) < GMX_REAL_EPS * 100)
2140             {
2141                 ir->opts.ref_t[i] = ir->opts.anneal_temp[i][j + 1];
2142             }
2143             else
2144             {
2145                 x = ((thist - ir->opts.anneal_time[i][j])
2146                      / (ir->opts.anneal_time[i][j + 1] - ir->opts.anneal_time[i][j]));
2147                 ir->opts.ref_t[i] =
2148                         x * ir->opts.anneal_temp[i][j + 1] + (1 - x) * ir->opts.anneal_temp[i][j];
2149             }
2150         }
2151         else
2152         {
2153             ir->opts.ref_t[i] = ir->opts.anneal_temp[i][npoints - 1];
2154         }
2155     }
2156
2157     upd->update_temperature_constants(*ir);
2158 }
2159
2160 void pleaseCiteCouplingAlgorithms(FILE* fplog, const t_inputrec& ir)
2161 {
2162     if (EI_DYNAMICS(ir.eI))
2163     {
2164         if (ir.etc == TemperatureCoupling::Berendsen)
2165         {
2166             please_cite(fplog, "Berendsen84a");
2167         }
2168         if (ir.etc == TemperatureCoupling::VRescale)
2169         {
2170             please_cite(fplog, "Bussi2007a");
2171         }
2172         if (ir.epc == PressureCoupling::CRescale)
2173         {
2174             please_cite(fplog, "Bernetti2020");
2175         }
2176         // TODO this is actually an integrator, not a coupling algorithm
2177         if (ir.eI == IntegrationAlgorithm::SD1)
2178         {
2179             please_cite(fplog, "Goga2012");
2180         }
2181     }
2182 }