src/gromacs/mdlib/coupling.cpp

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