src/gromacs/mdlib/coupling.cpp

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