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