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 == etcANDERSEN)
 761             {
 762                 bRandomize = TRUE;
 763             }
 764             else
 765             {
 766                 double uniform[2];
 767
 768                 gmx_rng_cycle_2uniform(step*2, ng, ir->andersen_seed, RND_SEED_ANDERSEN, uniform);
 769                 bRandomize = (uniform[0] < rate);
 770             }
 771             if (bRandomize)
 772             {
 773                 real scal, gauss[3];
 774                 int  d;
 775
 776                 scal = sqrt(boltzfac[gc]*md->invmass[i]);
 777                 gmx_rng_cycle_3gaussian_table(step*2+1, ng, ir->andersen_seed, RND_SEED_ANDERSEN, gauss);
 778                 for (d = 0; d < DIM; d++)
 779                 {
 780                     state->v[i][d] = scal*gauss[d];
 781                 }
 782             }
 783         }
 784     }
 785 }
 786
 787
 788 void nosehoover_tcoupl(t_grpopts *opts, gmx_ekindata_t *ekind, real dt,
 789                        double xi[], double vxi[], t_extmass *MassQ)
 790 {
 791     int   i;
 792     real  reft, oldvxi;
 793
 794     /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
 795
 796     for (i = 0; (i < opts->ngtc); i++)
 797     {
 798         reft     = std::max<real>(0, opts->ref_t[i]);
 799         oldvxi   = vxi[i];
 800         vxi[i]  += dt*MassQ->Qinv[i]*(ekind->tcstat[i].Th - reft);
 801         xi[i]   += dt*(oldvxi + vxi[i])*0.5;
 802     }
 803 }
 804
 805 t_state *init_bufstate(const t_state *template_state)
 806 {
 807     t_state *state;
 808     int      nc = template_state->nhchainlength;
 809     snew(state, 1);
 810     snew(state->nosehoover_xi, nc*template_state->ngtc);
 811     snew(state->nosehoover_vxi, nc*template_state->ngtc);
 812     snew(state->therm_integral, template_state->ngtc);
 813     snew(state->nhpres_xi, nc*template_state->nnhpres);
 814     snew(state->nhpres_vxi, nc*template_state->nnhpres);
 815
 816     return state;
 817 }
 818
 819 void destroy_bufstate(t_state *state)
 820 {
 821     sfree(state->x);
 822     sfree(state->v);
 823     sfree(state->nosehoover_xi);
 824     sfree(state->nosehoover_vxi);
 825     sfree(state->therm_integral);
 826     sfree(state->nhpres_xi);
 827     sfree(state->nhpres_vxi);
 828     sfree(state);
 829 }
 830
 831 void trotter_update(t_inputrec *ir, gmx_int64_t step, gmx_ekindata_t *ekind,
 832                     gmx_enerdata_t *enerd, t_state *state,
 833                     tensor vir, t_mdatoms *md,
 834                     t_extmass *MassQ, int **trotter_seqlist, int trotter_seqno)
 835 {
 836
 837     int             n, i, d, ngtc, gc = 0;
 838     t_grp_tcstat   *tcstat;
 839     t_grpopts      *opts;
 840     gmx_int64_t     step_eff;
 841     real            dt;
 842     double         *scalefac, dtc;
 843     int            *trotter_seq;
 844     rvec            sumv = {0, 0, 0};
 845     gmx_bool        bCouple;
 846
 847     if (trotter_seqno <= ettTSEQ2)
 848     {
 849         step_eff = step-1;  /* the velocity verlet calls are actually out of order -- the first half step
 850                                is actually the last half step from the previous step.  Thus the first half step
 851                                actually corresponds to the n-1 step*/
 852
 853     }
 854     else
 855     {
 856         step_eff = step;
 857     }
 858
 859     bCouple = (ir->nsttcouple == 1 ||
 860                do_per_step(step_eff+ir->nsttcouple, ir->nsttcouple));
 861
 862     trotter_seq = trotter_seqlist[trotter_seqno];
 863
 864     if ((trotter_seq[0] == etrtSKIPALL) || (!bCouple))
 865     {
 866         return;
 867     }
 868     dtc  = ir->nsttcouple*ir->delta_t; /* This is OK for NPT, because nsttcouple == nstpcouple is enforcesd */
 869     opts = &(ir->opts);                /* just for ease of referencing */
 870     ngtc = opts->ngtc;
 871     assert(ngtc > 0);
 872     snew(scalefac, opts->ngtc);
 873     for (i = 0; i < ngtc; i++)
 874     {
 875         scalefac[i] = 1;
 876     }
 877     /* execute the series of trotter updates specified in the trotterpart array */
 878
 879     for (i = 0; i < NTROTTERPARTS; i++)
 880     {
 881         /* allow for doubled intgrators by doubling dt instead of making 2 calls */
 882         if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2) || (trotter_seq[i] == etrtNHC2))
 883         {
 884             dt = 2 * dtc;
 885         }
 886         else
 887         {
 888             dt = dtc;
 889         }
 890
 891         switch (trotter_seq[i])
 892         {
 893             case etrtBAROV:
 894             case etrtBAROV2:
 895                 boxv_trotter(ir, &(state->veta), dt, state->box, ekind, vir,
 896                              enerd->term[F_PDISPCORR], MassQ);
 897                 break;
 898             case etrtBARONHC:
 899             case etrtBARONHC2:
 900                 NHC_trotter(opts, state->nnhpres, ekind, dt, state->nhpres_xi,
 901                             state->nhpres_vxi, NULL, &(state->veta), MassQ, FALSE);
 902                 break;
 903             case etrtNHC:
 904             case etrtNHC2:
 905                 NHC_trotter(opts, opts->ngtc, ekind, dt, state->nosehoover_xi,
 906                             state->nosehoover_vxi, scalefac, NULL, MassQ, (ir->eI == eiVV));
 907                 /* need to rescale the kinetic energies and velocities here.  Could
 908                    scale the velocities later, but we need them scaled in order to
 909                    produce the correct outputs, so we'll scale them here. */
 910
 911                 for (i = 0; i < ngtc; i++)
 912                 {
 913                     tcstat                  = &ekind->tcstat[i];
 914                     tcstat->vscale_nhc      = scalefac[i];
 915                     tcstat->ekinscaleh_nhc *= (scalefac[i]*scalefac[i]);
 916                     tcstat->ekinscalef_nhc *= (scalefac[i]*scalefac[i]);
 917                 }
 918                 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
 919                 /* but do we actually need the total? */
 920
 921                 /* modify the velocities as well */
 922                 for (n = 0; n < md->homenr; n++)
 923                 {
 924                     if (md->cTC) /* does this conditional need to be here? is this always true?*/
 925                     {
 926                         gc = md->cTC[n];
 927                     }
 928                     for (d = 0; d < DIM; d++)
 929                     {
 930                         state->v[n][d] *= scalefac[gc];
 931                     }
 932
 933                     if (debug)
 934                     {
 935                         for (d = 0; d < DIM; d++)
 936                         {
 937                             sumv[d] += (state->v[n][d])/md->invmass[n];
 938                         }
 939                     }
 940                 }
 941                 break;
 942             default:
 943                 break;
 944         }
 945     }
 946     /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
 947 #if 0
 948     if (debug)
 949     {
 950         if (bFirstHalf)
 951         {
 952             for (d = 0; d < DIM; d++)
 953             {
 954                 consk[d] = sumv[d]*exp((1 + 1.0/opts->nrdf[0])*((1.0/DIM)*log(det(state->box)/state->vol0)) + state->nosehoover_xi[0]);
 955             }
 956             fprintf(debug, "Conserved kappa: %15.8f %15.8f %15.8f\n", consk[0], consk[1], consk[2]);
 957         }
 958     }
 959 #endif
 960     sfree(scalefac);
 961 }
 962
 963
 964 extern void init_npt_masses(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bInit)
 965 {
 966     int           n, i, j, d, ngtc, nh;
 967     t_grpopts    *opts;
 968     real          reft, kT, ndj, nd;
 969
 970     opts    = &(ir->opts); /* just for ease of referencing */
 971     ngtc    = ir->opts.ngtc;
 972     nh      = state->nhchainlength;
 973
 974     if (ir->eI == eiMD)
 975     {
 976         if (bInit)
 977         {
 978             snew(MassQ->Qinv, ngtc);
 979         }
 980         for (i = 0; (i < ngtc); i++)
 981         {
 982             if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
 983             {
 984                 MassQ->Qinv[i] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*opts->ref_t[i]);
 985             }
 986             else
 987             {
 988                 MassQ->Qinv[i] = 0.0;
 989             }
 990         }
 991     }
 992     else if (EI_VV(ir->eI))
 993     {
 994         /* Set pressure variables */
 995
 996         if (bInit)
 997         {
 998             if (state->vol0 == 0)
 999             {
1000                 state->vol0 = det(state->box);
1001                 /* because we start by defining a fixed
1002                    compressibility, we need the volume at this
1003                    compressibility to solve the problem. */
1004             }
1005         }
1006
1007         /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol  */
1008         /* Consider evaluating eventually if this the right mass to use.  All are correct, some might be more stable  */
1009         MassQ->Winv = (PRESFAC*trace(ir->compress)*BOLTZ*opts->ref_t[0])/(DIM*state->vol0*sqr(ir->tau_p/M_2PI));
1010         /* An alternate mass definition, from Tuckerman et al. */
1011         /* MassQ->Winv = 1.0/(sqr(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
1012         for (d = 0; d < DIM; d++)
1013         {
1014             for (n = 0; n < DIM; n++)
1015             {
1016                 MassQ->Winvm[d][n] = PRESFAC*ir->compress[d][n]/(state->vol0*sqr(ir->tau_p/M_2PI));
1017                 /* not clear this is correct yet for the anisotropic case. Will need to reevaluate
1018                    before using MTTK for anisotropic states.*/
1019             }
1020         }
1021         /* Allocate space for thermostat variables */
1022         if (bInit)
1023         {
1024             snew(MassQ->Qinv, ngtc*nh);
1025         }
1026
1027         /* now, set temperature variables */
1028         for (i = 0; i < ngtc; i++)
1029         {
1030             if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
1031             {
1032                 reft = std::max<real>(0, opts->ref_t[i]);
1033                 nd   = opts->nrdf[i];
1034                 kT   = BOLTZ*reft;
1035                 for (j = 0; j < nh; j++)
1036                 {
1037                     if (j == 0)
1038                     {
1039                         ndj = nd;
1040                     }
1041                     else
1042                     {
1043                         ndj = 1;
1044                     }
1045                     MassQ->Qinv[i*nh+j]   = 1.0/(sqr(opts->tau_t[i]/M_2PI)*ndj*kT);
1046                 }
1047             }
1048             else
1049             {
1050                 for (j = 0; j < nh; j++)
1051                 {
1052                     MassQ->Qinv[i*nh+j] = 0.0;
1053                 }
1054             }
1055         }
1056     }
1057 }
1058
1059 int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bTrotter)
1060 {
1061     int           i, j, nnhpres, nh;
1062     t_grpopts    *opts;
1063     real          bmass, qmass, reft, kT;
1064     int         **trotter_seq;
1065
1066     opts    = &(ir->opts); /* just for ease of referencing */
1067     nnhpres = state->nnhpres;
1068     nh      = state->nhchainlength;
1069
1070     init_npt_masses(ir, state, MassQ, TRUE);
1071
1072     /* first, initialize clear all the trotter calls */
1073     snew(trotter_seq, ettTSEQMAX);
1074     for (i = 0; i < ettTSEQMAX; i++)
1075     {
1076         snew(trotter_seq[i], NTROTTERPARTS);
1077         for (j = 0; j < NTROTTERPARTS; j++)
1078         {
1079             trotter_seq[i][j] = etrtNONE;
1080         }
1081         trotter_seq[i][0] = etrtSKIPALL;
1082     }
1083
1084     if (!bTrotter)
1085     {
1086         /* no trotter calls, so we never use the values in the array.
1087          * We access them (so we need to define them, but ignore
1088          * then.*/
1089
1090         return trotter_seq;
1091     }
1092
1093     /* compute the kinetic energy by using the half step velocities or
1094      * the kinetic energies, depending on the order of the trotter calls */
1095
1096     if (ir->eI == eiVV)
1097     {
1098         if (IR_NPT_TROTTER(ir))
1099         {
1100             /* This is the complicated version - there are 4 possible calls, depending on ordering.
1101                We start with the initial one. */
1102             /* first, a round that estimates veta. */
1103             trotter_seq[0][0] = etrtBAROV;
1104
1105             /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1106
1107             /* The first half trotter update */
1108             trotter_seq[2][0] = etrtBAROV;
1109             trotter_seq[2][1] = etrtNHC;
1110             trotter_seq[2][2] = etrtBARONHC;
1111
1112             /* The second half trotter update */
1113             trotter_seq[3][0] = etrtBARONHC;
1114             trotter_seq[3][1] = etrtNHC;
1115             trotter_seq[3][2] = etrtBAROV;
1116
1117             /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1118
1119         }
1120         else if (IR_NVT_TROTTER(ir))
1121         {
1122             /* This is the easy version - there are only two calls, both the same.
1123                Otherwise, even easier -- no calls  */
1124             trotter_seq[2][0] = etrtNHC;
1125             trotter_seq[3][0] = etrtNHC;
1126         }
1127         else if (IR_NPH_TROTTER(ir))
1128         {
1129             /* This is the complicated version - there are 4 possible calls, depending on ordering.
1130                We start with the initial one. */
1131             /* first, a round that estimates veta. */
1132             trotter_seq[0][0] = etrtBAROV;
1133
1134             /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1135
1136             /* The first half trotter update */
1137             trotter_seq[2][0] = etrtBAROV;
1138             trotter_seq[2][1] = etrtBARONHC;
1139
1140             /* The second half trotter update */
1141             trotter_seq[3][0] = etrtBARONHC;
1142             trotter_seq[3][1] = etrtBAROV;
1143
1144             /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1145         }
1146     }
1147     else if (ir->eI == eiVVAK)
1148     {
1149         if (IR_NPT_TROTTER(ir))
1150         {
1151             /* This is the complicated version - there are 4 possible calls, depending on ordering.
1152                We start with the initial one. */
1153             /* first, a round that estimates veta. */
1154             trotter_seq[0][0] = etrtBAROV;
1155
1156             /* The first half trotter update, part 1 -- double update, because it commutes */
1157             trotter_seq[1][0] = etrtNHC;
1158
1159             /* The first half trotter update, part 2 */
1160             trotter_seq[2][0] = etrtBAROV;
1161             trotter_seq[2][1] = etrtBARONHC;
1162
1163             /* The second half trotter update, part 1 */
1164             trotter_seq[3][0] = etrtBARONHC;
1165             trotter_seq[3][1] = etrtBAROV;
1166
1167             /* The second half trotter update */
1168             trotter_seq[4][0] = etrtNHC;
1169         }
1170         else if (IR_NVT_TROTTER(ir))
1171         {
1172             /* This is the easy version - there is only one call, both the same.
1173                Otherwise, even easier -- no calls  */
1174             trotter_seq[1][0] = etrtNHC;
1175             trotter_seq[4][0] = etrtNHC;
1176         }
1177         else if (IR_NPH_TROTTER(ir))
1178         {
1179             /* This is the complicated version - there are 4 possible calls, depending on ordering.
1180                We start with the initial one. */
1181             /* first, a round that estimates veta. */
1182             trotter_seq[0][0] = etrtBAROV;
1183
1184             /* The first half trotter update, part 1 -- leave zero */
1185             trotter_seq[1][0] = etrtNHC;
1186
1187             /* The first half trotter update, part 2 */
1188             trotter_seq[2][0] = etrtBAROV;
1189             trotter_seq[2][1] = etrtBARONHC;
1190
1191             /* The second half trotter update, part 1 */
1192             trotter_seq[3][0] = etrtBARONHC;
1193             trotter_seq[3][1] = etrtBAROV;
1194
1195             /* The second half trotter update -- blank for now */
1196         }
1197     }
1198
1199     switch (ir->epct)
1200     {
1201         case epctISOTROPIC:
1202         default:
1203             bmass = DIM*DIM; /* recommended mass parameters for isotropic barostat */
1204     }
1205
1206     snew(MassQ->QPinv, nnhpres*opts->nhchainlength);
1207
1208     /* barostat temperature */
1209     if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
1210     {
1211         reft = std::max<real>(0, opts->ref_t[0]);
1212         kT   = BOLTZ*reft;
1213         for (i = 0; i < nnhpres; i++)
1214         {
1215             for (j = 0; j < nh; j++)
1216             {
1217                 if (j == 0)
1218                 {
1219                     qmass = bmass;
1220                 }
1221                 else
1222                 {
1223                     qmass = 1;
1224                 }
1225                 MassQ->QPinv[i*opts->nhchainlength+j]   = 1.0/(sqr(opts->tau_t[0]/M_2PI)*qmass*kT);
1226             }
1227         }
1228     }
1229     else
1230     {
1231         for (i = 0; i < nnhpres; i++)
1232         {
1233             for (j = 0; j < nh; j++)
1234             {
1235                 MassQ->QPinv[i*nh+j] = 0.0;
1236             }
1237         }
1238     }
1239     return trotter_seq;
1240 }
1241
1242 real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ)
1243 {
1244     int     i, j;
1245     real    nd, ndj;
1246     real    ener_npt, reft, kT;
1247     double *ivxi, *ixi;
1248     double *iQinv;
1249     real    vol;
1250     int     nh = state->nhchainlength;
1251
1252     ener_npt = 0;
1253
1254     /* now we compute the contribution of the pressure to the conserved quantity*/
1255
1256     if (ir->epc == epcMTTK)
1257     {
1258         /* find the volume, and the kinetic energy of the volume */
1259
1260         switch (ir->epct)
1261         {
1262
1263             case epctISOTROPIC:
1264                 /* contribution from the pressure momenenta */
1265                 ener_npt += 0.5*sqr(state->veta)/MassQ->Winv;
1266
1267                 /* contribution from the PV term */
1268                 vol       = det(state->box);
1269                 ener_npt += vol*trace(ir->ref_p)/(DIM*PRESFAC);
1270
1271                 break;
1272             case epctANISOTROPIC:
1273
1274                 break;
1275
1276             case epctSURFACETENSION:
1277
1278                 break;
1279             case epctSEMIISOTROPIC:
1280
1281                 break;
1282             default:
1283                 break;
1284         }
1285     }
1286
1287     if (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir))
1288     {
1289         /* add the energy from the barostat thermostat chain */
1290         for (i = 0; i < state->nnhpres; i++)
1291         {
1292
1293             /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1294             ivxi  = &state->nhpres_vxi[i*nh];
1295             ixi   = &state->nhpres_xi[i*nh];
1296             iQinv = &(MassQ->QPinv[i*nh]);
1297             reft  = std::max<real>(ir->opts.ref_t[0], 0.0); /* using 'System' temperature */
1298             kT    = BOLTZ * reft;
1299
1300             for (j = 0; j < nh; j++)
1301             {
1302                 if (iQinv[j] > 0)
1303                 {
1304                     ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1305                     /* contribution from the thermal variable of the NH chain */
1306                     ener_npt += ixi[j]*kT;
1307                 }
1308                 if (debug)
1309                 {
1310                     fprintf(debug, "P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f", i, j, ivxi[j], ixi[j]);
1311                 }
1312             }
1313         }
1314     }
1315
1316     if (ir->etc)
1317     {
1318         for (i = 0; i < ir->opts.ngtc; i++)
1319         {
1320             ixi   = &state->nosehoover_xi[i*nh];
1321             ivxi  = &state->nosehoover_vxi[i*nh];
1322             iQinv = &(MassQ->Qinv[i*nh]);
1323
1324             nd   = ir->opts.nrdf[i];
1325             reft = std::max<real>(ir->opts.ref_t[i], 0);
1326             kT   = BOLTZ * reft;
1327
1328             if (nd > 0.0)
1329             {
1330                 if (IR_NVT_TROTTER(ir))
1331                 {
1332                     /* contribution from the thermal momenta of the NH chain */
1333                     for (j = 0; j < nh; j++)
1334                     {
1335                         if (iQinv[j] > 0)
1336                         {
1337                             ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1338                             /* contribution from the thermal variable of the NH chain */
1339                             if (j == 0)
1340                             {
1341                                 ndj = nd;
1342                             }
1343                             else
1344                             {
1345                                 ndj = 1.0;
1346                             }
1347                             ener_npt += ndj*ixi[j]*kT;
1348                         }
1349                     }
1350                 }
1351                 else  /* Other non Trotter temperature NH control  -- no chains yet. */
1352                 {
1353                     ener_npt += 0.5*BOLTZ*nd*sqr(ivxi[0])/iQinv[0];
1354                     ener_npt += nd*ixi[0]*kT;
1355                 }
1356             }
1357         }
1358     }
1359     return ener_npt;
1360 }
1361
1362 static real vrescale_gamdev(real ia,
1363                             gmx_int64_t step, gmx_int64_t *count,
1364                             gmx_int64_t seed1, gmx_int64_t seed2)
1365 /* Gamma distribution, adapted from numerical recipes */
1366 {
1367     real   am, e, s, v1, v2, x, y;
1368     double rnd[2];
1369
1370     assert(ia > 1);
1371
1372     do
1373     {
1374         do
1375         {
1376             do
1377             {
1378                 gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
1379                 v1 = rnd[0];
1380                 v2 = 2.0*rnd[1] - 1.0;
1381             }
1382             while (v1*v1 + v2*v2 > 1.0 ||
1383                    v1*v1*GMX_REAL_MAX < 3.0*ia);
1384             /* The last check above ensures that both x (3.0 > 2.0 in s)
1385              * and the pre-factor for e do not go out of range.
1386              */
1387             y  = v2/v1;
1388             am = ia - 1;
1389             s  = sqrt(2.0*am + 1.0);
1390             x  = s*y + am;
1391         }
1392         while (x <= 0.0);
1393
1394         e = (1.0 + y*y)*exp(am*log(x/am) - s*y);
1395
1396         gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
1397     }
1398     while (rnd[0] > e);
1399
1400     return x;
1401 }
1402
1403 static real gaussian_count(gmx_int64_t step, gmx_int64_t *count,
1404                            gmx_int64_t seed1, gmx_int64_t seed2)
1405 {
1406     double rnd[2], x, y, r;
1407
1408     do
1409     {
1410         gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
1411         x = 2.0*rnd[0] - 1.0;
1412         y = 2.0*rnd[1] - 1.0;
1413         r = x*x + y*y;
1414     }
1415     while (r > 1.0 || r == 0.0);
1416
1417     r = sqrt(-2.0*log(r)/r);
1418
1419     return x*r;
1420 }
1421
1422 static real vrescale_sumnoises(real nn,
1423                                gmx_int64_t step, gmx_int64_t *count,
1424                                gmx_int64_t seed1, gmx_int64_t seed2)
1425 {
1426 /*
1427  * Returns the sum of nn independent gaussian noises squared
1428  * (i.e. equivalent to summing the square of the return values
1429  * of nn calls to gmx_rng_gaussian_real).
1430  */
1431     const real ndeg_tol = 0.0001;
1432     real       r;
1433
1434     if (nn < 2 + ndeg_tol)
1435     {
1436         int  nn_int, i;
1437         real gauss;
1438
1439         nn_int = (int)(nn + 0.5);
1440
1441         if (nn - nn_int < -ndeg_tol || nn - nn_int > ndeg_tol)
1442         {
1443             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);
1444         }
1445
1446         r = 0;
1447         for (i = 0; i < nn_int; i++)
1448         {
1449             gauss = gaussian_count(step, count, seed1, seed2);
1450
1451             r += gauss*gauss;
1452         }
1453     }
1454     else
1455     {
1456         /* Use a gamma distribution for any real nn > 2 */
1457         r = 2.0*vrescale_gamdev(0.5*nn, step, count, seed1, seed2);
1458     }
1459
1460     return r;
1461 }
1462
1463 static real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut,
1464                                  gmx_int64_t step, gmx_int64_t seed)
1465 {
1466 /*
1467  * Generates a new value for the kinetic energy,
1468  * according to Bussi et al JCP (2007), Eq. (A7)
1469  * kk:    present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1470  * sigma: target average value of the kinetic energy (ndeg k_b T/2)  (in the same units as kk)
1471  * ndeg:  number of degrees of freedom of the atoms to be thermalized
1472  * taut:  relaxation time of the thermostat, in units of 'how often this routine is called'
1473  */
1474     /* rnd_count tracks the step-local state for the cycle random
1475      * number generator.
1476      */
1477     gmx_int64_t rnd_count = 0;
1478     real        factor, rr, ekin_new;
1479
1480     if (taut > 0.1)
1481     {
1482         factor = exp(-1.0/taut);
1483     }
1484     else
1485     {
1486         factor = 0.0;
1487     }
1488
1489     rr = gaussian_count(step, &rnd_count, seed, RND_SEED_VRESCALE);
1490
1491     ekin_new =
1492         kk +
1493         (1.0 - factor)*(sigma*(vrescale_sumnoises(ndeg-1, step, &rnd_count, seed, RND_SEED_VRESCALE) + rr*rr)/ndeg - kk) +
1494         2.0*rr*sqrt(kk*sigma/ndeg*(1.0 - factor)*factor);
1495
1496     return ekin_new;
1497 }
1498
1499 void vrescale_tcoupl(t_inputrec *ir, gmx_int64_t step,
1500                      gmx_ekindata_t *ekind, real dt,
1501                      double therm_integral[])
1502 {
1503     t_grpopts *opts;
1504     int        i;
1505     real       Ek, Ek_ref1, Ek_ref, Ek_new;
1506
1507     opts = &ir->opts;
1508
1509     for (i = 0; (i < opts->ngtc); i++)
1510     {
1511         if (ir->eI == eiVV)
1512         {
1513             Ek = trace(ekind->tcstat[i].ekinf);
1514         }
1515         else
1516         {
1517             Ek = trace(ekind->tcstat[i].ekinh);
1518         }
1519
1520         if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0)
1521         {
1522             Ek_ref1 = 0.5*opts->ref_t[i]*BOLTZ;
1523             Ek_ref  = Ek_ref1*opts->nrdf[i];
1524
1525             Ek_new  = vrescale_resamplekin(Ek, Ek_ref, opts->nrdf[i],
1526                                            opts->tau_t[i]/dt,
1527                                            step, ir->ld_seed);
1528
1529             /* Analytically Ek_new>=0, but we check for rounding errors */
1530             if (Ek_new <= 0)
1531             {
1532                 ekind->tcstat[i].lambda = 0.0;
1533             }
1534             else
1535             {
1536                 ekind->tcstat[i].lambda = sqrt(Ek_new/Ek);
1537             }
1538
1539             therm_integral[i] -= Ek_new - Ek;
1540
1541             if (debug)
1542             {
1543                 fprintf(debug, "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
1544                         i, Ek_ref, Ek, Ek_new, ekind->tcstat[i].lambda);
1545             }
1546         }
1547         else
1548         {
1549             ekind->tcstat[i].lambda = 1.0;
1550         }
1551     }
1552 }
1553
1554 real vrescale_energy(t_grpopts *opts, double therm_integral[])
1555 {
1556     int  i;
1557     real ener;
1558
1559     ener = 0;
1560     for (i = 0; i < opts->ngtc; i++)
1561     {
1562         ener += therm_integral[i];
1563     }
1564
1565     return ener;
1566 }
1567
1568 void rescale_velocities(gmx_ekindata_t *ekind, t_mdatoms *mdatoms,
1569                         int start, int end, rvec v[])
1570 {
1571     t_grp_acc      *gstat;
1572     t_grp_tcstat   *tcstat;
1573     unsigned short *cACC, *cTC;
1574     int             ga, gt, n, d;
1575     real            lg;
1576     rvec            vrel;
1577
1578     tcstat = ekind->tcstat;
1579     cTC    = mdatoms->cTC;
1580
1581     if (ekind->bNEMD)
1582     {
1583         gstat  = ekind->grpstat;
1584         cACC   = mdatoms->cACC;
1585
1586         ga = 0;
1587         gt = 0;
1588         for (n = start; n < end; n++)
1589         {
1590             if (cACC)
1591             {
1592                 ga   = cACC[n];
1593             }
1594             if (cTC)
1595             {
1596                 gt   = cTC[n];
1597             }
1598             /* Only scale the velocity component relative to the COM velocity */
1599             rvec_sub(v[n], gstat[ga].u, vrel);
1600             lg = tcstat[gt].lambda;
1601             for (d = 0; d < DIM; d++)
1602             {
1603                 v[n][d] = gstat[ga].u[d] + lg*vrel[d];
1604             }
1605         }
1606     }
1607     else
1608     {
1609         gt = 0;
1610         for (n = start; n < end; n++)
1611         {
1612             if (cTC)
1613             {
1614                 gt   = cTC[n];
1615             }
1616             lg = tcstat[gt].lambda;
1617             for (d = 0; d < DIM; d++)
1618             {
1619                 v[n][d] *= lg;
1620             }
1621         }
1622     }
1623 }
1624
1625
1626 /* set target temperatures if we are annealing */
1627 void update_annealing_target_temp(t_grpopts *opts, real t)
1628 {
1629     int  i, j, n, npoints;
1630     real pert, thist = 0, x;
1631
1632     for (i = 0; i < opts->ngtc; i++)
1633     {
1634         npoints = opts->anneal_npoints[i];
1635         switch (opts->annealing[i])
1636         {
1637             case eannNO:
1638                 continue;
1639             case  eannPERIODIC:
1640                 /* calculate time modulo the period */
1641                 pert  = opts->anneal_time[i][npoints-1];
1642                 n     = static_cast<int>(t / pert);
1643                 thist = t - n*pert; /* modulo time */
1644                 /* Make sure rounding didn't get us outside the interval */
1645                 if (fabs(thist-pert) < GMX_REAL_EPS*100)
1646                 {
1647                     thist = 0;
1648                 }
1649                 break;
1650             case eannSINGLE:
1651                 thist = t;
1652                 break;
1653             default:
1654                 gmx_fatal(FARGS, "Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)", i, opts->ngtc, npoints);
1655         }
1656         /* We are doing annealing for this group if we got here,
1657          * and we have the (relative) time as thist.
1658          * calculate target temp */
1659         j = 0;
1660         while ((j < npoints-1) && (thist > (opts->anneal_time[i][j+1])))
1661         {
1662             j++;
1663         }
1664         if (j < npoints-1)
1665         {
1666             /* Found our position between points j and j+1.
1667              * Interpolate: x is the amount from j+1, (1-x) from point j
1668              * First treat possible jumps in temperature as a special case.
1669              */
1670             if ((opts->anneal_time[i][j+1]-opts->anneal_time[i][j]) < GMX_REAL_EPS*100)
1671             {
1672                 opts->ref_t[i] = opts->anneal_temp[i][j+1];
1673             }
1674             else
1675             {
1676                 x = ((thist-opts->anneal_time[i][j])/
1677                      (opts->anneal_time[i][j+1]-opts->anneal_time[i][j]));
1678                 opts->ref_t[i] = x*opts->anneal_temp[i][j+1]+(1-x)*opts->anneal_temp[i][j];
1679             }
1680         }
1681         else
1682         {
1683             opts->ref_t[i] = opts->anneal_temp[i][npoints-1];
1684         }
1685     }
1686 }