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