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