src/gromacs/mdlib/coupling.c

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