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