src/gromacs/mdlib/coupling.cpp

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