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