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