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