src/gromacs/mdlib/coupling.cpp

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