src/gromacs/mdlib/coupling.cpp

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