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