src/gromacs/mdlib/update.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 #include "gmxpre.h"
  38
  39 #include "config.h"
  40
  41
  42 #include <stdio.h>
  43 #include <math.h>
  44
  45 #include "gromacs/legacyheaders/types/commrec.h"
  46 #include "gromacs/legacyheaders/typedefs.h"
  47 #include "gromacs/legacyheaders/nrnb.h"
  48 #include "gromacs/math/units.h"
  49 #include "gromacs/legacyheaders/macros.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/legacyheaders/update.h"
  52 #include "gromacs/random/random.h"
  53 #include "gromacs/legacyheaders/tgroup.h"
  54 #include "gromacs/legacyheaders/force.h"
  55 #include "gromacs/legacyheaders/names.h"
  56 #include "gromacs/legacyheaders/txtdump.h"
  57 #include "gromacs/legacyheaders/mdrun.h"
  58 #include "gromacs/legacyheaders/constr.h"
  59 #include "gromacs/legacyheaders/disre.h"
  60 #include "gromacs/legacyheaders/orires.h"
  61 #include "gromacs/legacyheaders/gmx_omp_nthreads.h"
  62
  63 #include "gromacs/fileio/confio.h"
  64 #include "gromacs/pbcutil/mshift.h"
  65 #include "gromacs/pbcutil/pbc.h"
  66 #include "gromacs/pulling/pull.h"
  67 #include "gromacs/timing/wallcycle.h"
  68 #include "gromacs/utility/futil.h"
  69 #include "gromacs/utility/gmxomp.h"
  70 #include "gromacs/utility/smalloc.h"
  71
  72 /*For debugging, start at v(-dt/2) for velolcity verlet -- uncomment next line */
  73 /*#define STARTFROMDT2*/
  74
  75 typedef struct {
  76     double gdt;
  77     double eph;
  78     double emh;
  79     double em;
  80     double b;
  81     double c;
  82     double d;
  83 } gmx_sd_const_t;
  84
  85 typedef struct {
  86     real V;
  87     real X;
  88     real Yv;
  89     real Yx;
  90 } gmx_sd_sigma_t;
  91
  92 typedef struct {
  93     /* BD stuff */
  94     real           *bd_rf;
  95     /* SD stuff */
  96     gmx_sd_const_t *sdc;
  97     gmx_sd_sigma_t *sdsig;
  98     rvec           *sd_V;
  99     int             sd_V_nalloc;
 100     /* andersen temperature control stuff */
 101     gmx_bool       *randomize_group;
 102     real           *boltzfac;
 103 } gmx_stochd_t;
 104
 105 typedef struct gmx_update
 106 {
 107     gmx_stochd_t *sd;
 108     /* xprime for constraint algorithms */
 109     rvec         *xp;
 110     int           xp_nalloc;
 111
 112     /* Variables for the deform algorithm */
 113     gmx_int64_t     deformref_step;
 114     matrix          deformref_box;
 115 } t_gmx_update;
 116
 117
 118 static void do_update_md(int start, int nrend, double dt,
 119                          t_grp_tcstat *tcstat,
 120                          double nh_vxi[],
 121                          gmx_bool bNEMD, t_grp_acc *gstat, rvec accel[],
 122                          ivec nFreeze[],
 123                          real invmass[],
 124                          unsigned short ptype[], unsigned short cFREEZE[],
 125                          unsigned short cACC[], unsigned short cTC[],
 126                          rvec x[], rvec xprime[], rvec v[],
 127                          rvec f[], matrix M,
 128                          gmx_bool bNH, gmx_bool bPR)
 129 {
 130     double imass, w_dt;
 131     int    gf = 0, ga = 0, gt = 0;
 132     rvec   vrel;
 133     real   vn, vv, va, vb, vnrel;
 134     real   lg, vxi = 0, u;
 135     int    n, d;
 136
 137     if (bNH || bPR)
 138     {
 139         /* Update with coupling to extended ensembles, used for
 140          * Nose-Hoover and Parrinello-Rahman coupling
 141          * Nose-Hoover uses the reversible leap-frog integrator from
 142          * Holian et al. Phys Rev E 52(3) : 2338, 1995
 143          */
 144         for (n = start; n < nrend; n++)
 145         {
 146             imass = invmass[n];
 147             if (cFREEZE)
 148             {
 149                 gf   = cFREEZE[n];
 150             }
 151             if (cACC)
 152             {
 153                 ga   = cACC[n];
 154             }
 155             if (cTC)
 156             {
 157                 gt   = cTC[n];
 158             }
 159             lg   = tcstat[gt].lambda;
 160             if (bNH)
 161             {
 162                 vxi   = nh_vxi[gt];
 163             }
 164             rvec_sub(v[n], gstat[ga].u, vrel);
 165
 166             for (d = 0; d < DIM; d++)
 167             {
 168                 if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
 169                 {
 170                     vnrel = (lg*vrel[d] + dt*(imass*f[n][d] - 0.5*vxi*vrel[d]
 171                                               - iprod(M[d], vrel)))/(1 + 0.5*vxi*dt);
 172                     /* do not scale the mean velocities u */
 173                     vn             = gstat[ga].u[d] + accel[ga][d]*dt + vnrel;
 174                     v[n][d]        = vn;
 175                     xprime[n][d]   = x[n][d]+vn*dt;
 176                 }
 177                 else
 178                 {
 179                     v[n][d]        = 0.0;
 180                     xprime[n][d]   = x[n][d];
 181                 }
 182             }
 183         }
 184     }
 185     else if (cFREEZE != NULL ||
 186              nFreeze[0][XX] || nFreeze[0][YY] || nFreeze[0][ZZ] ||
 187              bNEMD)
 188     {
 189         /* Update with Berendsen/v-rescale coupling and freeze or NEMD */
 190         for (n = start; n < nrend; n++)
 191         {
 192             w_dt = invmass[n]*dt;
 193             if (cFREEZE)
 194             {
 195                 gf   = cFREEZE[n];
 196             }
 197             if (cACC)
 198             {
 199                 ga   = cACC[n];
 200             }
 201             if (cTC)
 202             {
 203                 gt   = cTC[n];
 204             }
 205             lg   = tcstat[gt].lambda;
 206
 207             for (d = 0; d < DIM; d++)
 208             {
 209                 vn             = v[n][d];
 210                 if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
 211                 {
 212                     vv             = lg*vn + f[n][d]*w_dt;
 213
 214                     /* do not scale the mean velocities u */
 215                     u              = gstat[ga].u[d];
 216                     va             = vv + accel[ga][d]*dt;
 217                     vb             = va + (1.0-lg)*u;
 218                     v[n][d]        = vb;
 219                     xprime[n][d]   = x[n][d]+vb*dt;
 220                 }
 221                 else
 222                 {
 223                     v[n][d]        = 0.0;
 224                     xprime[n][d]   = x[n][d];
 225                 }
 226             }
 227         }
 228     }
 229     else
 230     {
 231         /* Plain update with Berendsen/v-rescale coupling */
 232         for (n = start; n < nrend; n++)
 233         {
 234             if ((ptype[n] != eptVSite) && (ptype[n] != eptShell))
 235             {
 236                 w_dt = invmass[n]*dt;
 237                 if (cTC)
 238                 {
 239                     gt = cTC[n];
 240                 }
 241                 lg = tcstat[gt].lambda;
 242
 243                 for (d = 0; d < DIM; d++)
 244                 {
 245                     vn           = lg*v[n][d] + f[n][d]*w_dt;
 246                     v[n][d]      = vn;
 247                     xprime[n][d] = x[n][d] + vn*dt;
 248                 }
 249             }
 250             else
 251             {
 252                 for (d = 0; d < DIM; d++)
 253                 {
 254                     v[n][d]        = 0.0;
 255                     xprime[n][d]   = x[n][d];
 256                 }
 257             }
 258         }
 259     }
 260 }
 261
 262 static void do_update_vv_vel(int start, int nrend, double dt,
 263                              rvec accel[], ivec nFreeze[], real invmass[],
 264                              unsigned short ptype[], unsigned short cFREEZE[],
 265                              unsigned short cACC[], rvec v[], rvec f[],
 266                              gmx_bool bExtended, real veta, real alpha)
 267 {
 268     double imass, w_dt;
 269     int    gf = 0, ga = 0;
 270     rvec   vrel;
 271     real   u, vn, vv, va, vb, vnrel;
 272     int    n, d;
 273     double g, mv1, mv2;
 274
 275     if (bExtended)
 276     {
 277         g        = 0.25*dt*veta*alpha;
 278         mv1      = exp(-g);
 279         mv2      = series_sinhx(g);
 280     }
 281     else
 282     {
 283         mv1      = 1.0;
 284         mv2      = 1.0;
 285     }
 286     for (n = start; n < nrend; n++)
 287     {
 288         w_dt = invmass[n]*dt;
 289         if (cFREEZE)
 290         {
 291             gf   = cFREEZE[n];
 292         }
 293         if (cACC)
 294         {
 295             ga   = cACC[n];
 296         }
 297
 298         for (d = 0; d < DIM; d++)
 299         {
 300             if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
 301             {
 302                 v[n][d]             = mv1*(mv1*v[n][d] + 0.5*(w_dt*mv2*f[n][d]))+0.5*accel[ga][d]*dt;
 303             }
 304             else
 305             {
 306                 v[n][d]        = 0.0;
 307             }
 308         }
 309     }
 310 } /* do_update_vv_vel */
 311
 312 static void do_update_vv_pos(int start, int nrend, double dt,
 313                              ivec nFreeze[],
 314                              unsigned short ptype[], unsigned short cFREEZE[],
 315                              rvec x[], rvec xprime[], rvec v[],
 316                              gmx_bool bExtended, real veta)
 317 {
 318     double imass, w_dt;
 319     int    gf = 0;
 320     int    n, d;
 321     double g, mr1, mr2;
 322
 323     /* Would it make more sense if Parrinello-Rahman was put here? */
 324     if (bExtended)
 325     {
 326         g        = 0.5*dt*veta;
 327         mr1      = exp(g);
 328         mr2      = series_sinhx(g);
 329     }
 330     else
 331     {
 332         mr1      = 1.0;
 333         mr2      = 1.0;
 334     }
 335
 336     for (n = start; n < nrend; n++)
 337     {
 338
 339         if (cFREEZE)
 340         {
 341             gf   = cFREEZE[n];
 342         }
 343
 344         for (d = 0; d < DIM; d++)
 345         {
 346             if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
 347             {
 348                 xprime[n][d]   = mr1*(mr1*x[n][d]+mr2*dt*v[n][d]);
 349             }
 350             else
 351             {
 352                 xprime[n][d]   = x[n][d];
 353             }
 354         }
 355     }
 356 } /* do_update_vv_pos */
 357
 358 static void do_update_visc(int start, int nrend, double dt,
 359                            t_grp_tcstat *tcstat,
 360                            double nh_vxi[],
 361                            real invmass[],
 362                            unsigned short ptype[], unsigned short cTC[],
 363                            rvec x[], rvec xprime[], rvec v[],
 364                            rvec f[], matrix M, matrix box, real
 365                            cos_accel, real vcos,
 366                            gmx_bool bNH, gmx_bool bPR)
 367 {
 368     double imass, w_dt;
 369     int    gt = 0;
 370     real   vn, vc;
 371     real   lg, vxi = 0, vv;
 372     real   fac, cosz;
 373     rvec   vrel;
 374     int    n, d;
 375
 376     fac = 2*M_PI/(box[ZZ][ZZ]);
 377
 378     if (bNH || bPR)
 379     {
 380         /* Update with coupling to extended ensembles, used for
 381          * Nose-Hoover and Parrinello-Rahman coupling
 382          */
 383         for (n = start; n < nrend; n++)
 384         {
 385             imass = invmass[n];
 386             if (cTC)
 387             {
 388                 gt   = cTC[n];
 389             }
 390             lg   = tcstat[gt].lambda;
 391             cosz = cos(fac*x[n][ZZ]);
 392
 393             copy_rvec(v[n], vrel);
 394
 395             vc            = cosz*vcos;
 396             vrel[XX]     -= vc;
 397             if (bNH)
 398             {
 399                 vxi        = nh_vxi[gt];
 400             }
 401             for (d = 0; d < DIM; d++)
 402             {
 403                 vn             = v[n][d];
 404
 405                 if ((ptype[n] != eptVSite) && (ptype[n] != eptShell))
 406                 {
 407                     vn  = (lg*vrel[d] + dt*(imass*f[n][d] - 0.5*vxi*vrel[d]
 408                                             - iprod(M[d], vrel)))/(1 + 0.5*vxi*dt);
 409                     if (d == XX)
 410                     {
 411                         vn += vc + dt*cosz*cos_accel;
 412                     }
 413                     v[n][d]        = vn;
 414                     xprime[n][d]   = x[n][d]+vn*dt;
 415                 }
 416                 else
 417                 {
 418                     xprime[n][d]   = x[n][d];
 419                 }
 420             }
 421         }
 422     }
 423     else
 424     {
 425         /* Classic version of update, used with berendsen coupling */
 426         for (n = start; n < nrend; n++)
 427         {
 428             w_dt = invmass[n]*dt;
 429             if (cTC)
 430             {
 431                 gt   = cTC[n];
 432             }
 433             lg   = tcstat[gt].lambda;
 434             cosz = cos(fac*x[n][ZZ]);
 435
 436             for (d = 0; d < DIM; d++)
 437             {
 438                 vn             = v[n][d];
 439
 440                 if ((ptype[n] != eptVSite) && (ptype[n] != eptShell))
 441                 {
 442                     if (d == XX)
 443                     {
 444                         vc           = cosz*vcos;
 445                         /* Do not scale the cosine velocity profile */
 446                         vv           = vc + lg*(vn - vc + f[n][d]*w_dt);
 447                         /* Add the cosine accelaration profile */
 448                         vv          += dt*cosz*cos_accel;
 449                     }
 450                     else
 451                     {
 452                         vv           = lg*(vn + f[n][d]*w_dt);
 453                     }
 454                     v[n][d]        = vv;
 455                     xprime[n][d]   = x[n][d]+vv*dt;
 456                 }
 457                 else
 458                 {
 459                     v[n][d]        = 0.0;
 460                     xprime[n][d]   = x[n][d];
 461                 }
 462             }
 463         }
 464     }
 465 }
 466
 467 static gmx_stochd_t *init_stochd(t_inputrec *ir)
 468 {
 469     gmx_stochd_t   *sd;
 470     gmx_sd_const_t *sdc;
 471     int             ngtc, n, th;
 472     real            y;
 473
 474     snew(sd, 1);
 475
 476     ngtc = ir->opts.ngtc;
 477
 478     if (ir->eI == eiBD)
 479     {
 480         snew(sd->bd_rf, ngtc);
 481     }
 482     else if (EI_SD(ir->eI))
 483     {
 484         snew(sd->sdc, ngtc);
 485         snew(sd->sdsig, ngtc);
 486
 487         sdc = sd->sdc;
 488         for (n = 0; n < ngtc; n++)
 489         {
 490             if (ir->opts.tau_t[n] > 0)
 491             {
 492                 sdc[n].gdt = ir->delta_t/ir->opts.tau_t[n];
 493                 sdc[n].eph = exp(sdc[n].gdt/2);
 494                 sdc[n].emh = exp(-sdc[n].gdt/2);
 495                 sdc[n].em  = exp(-sdc[n].gdt);
 496             }
 497             else
 498             {
 499                 /* No friction and noise on this group */
 500                 sdc[n].gdt = 0;
 501                 sdc[n].eph = 1;
 502                 sdc[n].emh = 1;
 503                 sdc[n].em  = 1;
 504             }
 505             if (sdc[n].gdt >= 0.05)
 506             {
 507                 sdc[n].b = sdc[n].gdt*(sdc[n].eph*sdc[n].eph - 1)
 508                     - 4*(sdc[n].eph - 1)*(sdc[n].eph - 1);
 509                 sdc[n].c = sdc[n].gdt - 3 + 4*sdc[n].emh - sdc[n].em;
 510                 sdc[n].d = 2 - sdc[n].eph - sdc[n].emh;
 511             }
 512             else
 513             {
 514                 y = sdc[n].gdt/2;
 515                 /* Seventh order expansions for small y */
 516                 sdc[n].b = y*y*y*y*(1/3.0+y*(1/3.0+y*(17/90.0+y*7/9.0)));
 517                 sdc[n].c = y*y*y*(2/3.0+y*(-1/2.0+y*(7/30.0+y*(-1/12.0+y*31/1260.0))));
 518                 sdc[n].d = y*y*(-1+y*y*(-1/12.0-y*y/360.0));
 519             }
 520             if (debug)
 521             {
 522                 fprintf(debug, "SD const tc-grp %d: b %g  c %g  d %g\n",
 523                         n, sdc[n].b, sdc[n].c, sdc[n].d);
 524             }
 525         }
 526     }
 527     else if (ETC_ANDERSEN(ir->etc))
 528     {
 529         int        ngtc;
 530         t_grpopts *opts;
 531         real       reft;
 532
 533         opts = &ir->opts;
 534         ngtc = opts->ngtc;
 535
 536         snew(sd->randomize_group, ngtc);
 537         snew(sd->boltzfac, ngtc);
 538
 539         /* for now, assume that all groups, if randomized, are randomized at the same rate, i.e. tau_t is the same. */
 540         /* since constraint groups don't necessarily match up with temperature groups! This is checked in readir.c */
 541
 542         for (n = 0; n < ngtc; n++)
 543         {
 544             reft = max(0.0, opts->ref_t[n]);
 545             if ((opts->tau_t[n] > 0) && (reft > 0))  /* tau_t or ref_t = 0 means that no randomization is done */
 546             {
 547                 sd->randomize_group[n] = TRUE;
 548                 sd->boltzfac[n]        = BOLTZ*opts->ref_t[n];
 549             }
 550             else
 551             {
 552                 sd->randomize_group[n] = FALSE;
 553             }
 554         }
 555     }
 556     return sd;
 557 }
 558
 559 gmx_update_t init_update(t_inputrec *ir)
 560 {
 561     t_gmx_update *upd;
 562
 563     snew(upd, 1);
 564
 565     if (ir->eI == eiBD || EI_SD(ir->eI) || ir->etc == etcVRESCALE || ETC_ANDERSEN(ir->etc))
 566     {
 567         upd->sd    = init_stochd(ir);
 568     }
 569
 570     upd->xp        = NULL;
 571     upd->xp_nalloc = 0;
 572
 573     return upd;
 574 }
 575
 576 static void do_update_sd1(gmx_stochd_t *sd,
 577                           int start, int nrend, double dt,
 578                           rvec accel[], ivec nFreeze[],
 579                           real invmass[], unsigned short ptype[],
 580                           unsigned short cFREEZE[], unsigned short cACC[],
 581                           unsigned short cTC[],
 582                           rvec x[], rvec xprime[], rvec v[], rvec f[],
 583                           int ngtc, real ref_t[],
 584                           gmx_bool bDoConstr,
 585                           gmx_bool bFirstHalfConstr,
 586                           gmx_int64_t step, int seed, int* gatindex)
 587 {
 588     gmx_sd_const_t *sdc;
 589     gmx_sd_sigma_t *sig;
 590     real            kT;
 591     int             gf = 0, ga = 0, gt = 0;
 592     real            ism;
 593     int             n, d;
 594
 595     sdc = sd->sdc;
 596     sig = sd->sdsig;
 597
 598     for (n = 0; n < ngtc; n++)
 599     {
 600         kT = BOLTZ*ref_t[n];
 601         /* The mass is encounted for later, since this differs per atom */
 602         sig[n].V  = sqrt(kT*(1 - sdc[n].em*sdc[n].em));
 603     }
 604
 605     if (!bDoConstr)
 606     {
 607         for (n = start; n < nrend; n++)
 608         {
 609             real rnd[3];
 610             int  ng = gatindex ? gatindex[n] : n;
 611
 612             ism = sqrt(invmass[n]);
 613             if (cFREEZE)
 614             {
 615                 gf  = cFREEZE[n];
 616             }
 617             if (cACC)
 618             {
 619                 ga  = cACC[n];
 620             }
 621             if (cTC)
 622             {
 623                 gt  = cTC[n];
 624             }
 625
 626             gmx_rng_cycle_3gaussian_table(step, ng, seed, RND_SEED_UPDATE, rnd);
 627
 628             for (d = 0; d < DIM; d++)
 629             {
 630                 if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
 631                 {
 632                     real sd_V, vn;
 633
 634                     sd_V         = ism*sig[gt].V*rnd[d];
 635                     vn           = v[n][d] + (invmass[n]*f[n][d] + accel[ga][d])*dt;
 636                     v[n][d]      = vn*sdc[gt].em + sd_V;
 637                     /* Here we include half of the friction+noise
 638                      * update of v into the integration of x.
 639                      */
 640                     xprime[n][d] = x[n][d] + 0.5*(vn + v[n][d])*dt;
 641                 }
 642                 else
 643                 {
 644                     v[n][d]      = 0.0;
 645                     xprime[n][d] = x[n][d];
 646                 }
 647             }
 648         }
 649     }
 650     else
 651     {
 652         /* We do have constraints */
 653         if (bFirstHalfConstr)
 654         {
 655             /* First update without friction and noise */
 656             real im;
 657
 658             for (n = start; n < nrend; n++)
 659             {
 660                 im = invmass[n];
 661
 662                 if (cFREEZE)
 663                 {
 664                     gf  = cFREEZE[n];
 665                 }
 666                 if (cACC)
 667                 {
 668                     ga  = cACC[n];
 669                 }
 670                 if (cTC)
 671                 {
 672                     gt  = cTC[n];
 673                 }
 674
 675                 for (d = 0; d < DIM; d++)
 676                 {
 677                     if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
 678                     {
 679                         v[n][d]      = v[n][d] + (im*f[n][d] + accel[ga][d])*dt;
 680                         xprime[n][d] = x[n][d] +  v[n][d]*dt;
 681                     }
 682                     else
 683                     {
 684                         v[n][d]      = 0.0;
 685                         xprime[n][d] = x[n][d];
 686                     }
 687                 }
 688             }
 689         }
 690         else
 691         {
 692             /* Update friction and noise only */
 693             for (n = start; n < nrend; n++)
 694             {
 695                 real rnd[3];
 696                 int  ng = gatindex ? gatindex[n] : n;
 697
 698                 ism = sqrt(invmass[n]);
 699                 if (cFREEZE)
 700                 {
 701                     gf  = cFREEZE[n];
 702                 }
 703                 if (cACC)
 704                 {
 705                     ga  = cACC[n];
 706                 }
 707                 if (cTC)
 708                 {
 709                     gt  = cTC[n];
 710                 }
 711
 712                 gmx_rng_cycle_3gaussian_table(step, ng, seed, RND_SEED_UPDATE, rnd);
 713
 714                 for (d = 0; d < DIM; d++)
 715                 {
 716                     if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
 717                     {
 718                         real sd_V, vn;
 719
 720                         sd_V         = ism*sig[gt].V*rnd[d];
 721                         vn           = v[n][d];
 722                         v[n][d]      = vn*sdc[gt].em + sd_V;
 723                         /* Add the friction and noise contribution only */
 724                         xprime[n][d] = xprime[n][d] + 0.5*(v[n][d] - vn)*dt;
 725                     }
 726                 }
 727             }
 728         }
 729     }
 730 }
 731
 732 static void check_sd2_work_data_allocation(gmx_stochd_t *sd, int nrend)
 733 {
 734     if (nrend > sd->sd_V_nalloc)
 735     {
 736         sd->sd_V_nalloc = over_alloc_dd(nrend);
 737         srenew(sd->sd_V, sd->sd_V_nalloc);
 738     }
 739 }
 740
 741 static void do_update_sd2_Tconsts(gmx_stochd_t *sd,
 742                                   int           ngtc,
 743                                   const real    tau_t[],
 744                                   const real    ref_t[])
 745 {
 746     /* This is separated from the update below, because it is single threaded */
 747     gmx_sd_const_t *sdc;
 748     gmx_sd_sigma_t *sig;
 749     int             gt;
 750     real            kT;
 751
 752     sdc = sd->sdc;
 753     sig = sd->sdsig;
 754
 755     for (gt = 0; gt < ngtc; gt++)
 756     {
 757         kT = BOLTZ*ref_t[gt];
 758         /* The mass is encounted for later, since this differs per atom */
 759         sig[gt].V  = sqrt(kT*(1-sdc[gt].em));
 760         sig[gt].X  = sqrt(kT*sqr(tau_t[gt])*sdc[gt].c);
 761         sig[gt].Yv = sqrt(kT*sdc[gt].b/sdc[gt].c);
 762         sig[gt].Yx = sqrt(kT*sqr(tau_t[gt])*sdc[gt].b/(1-sdc[gt].em));
 763     }
 764 }
 765
 766 static void do_update_sd2(gmx_stochd_t *sd,
 767                           gmx_bool bInitStep,
 768                           int start, int nrend,
 769                           rvec accel[], ivec nFreeze[],
 770                           real invmass[], unsigned short ptype[],
 771                           unsigned short cFREEZE[], unsigned short cACC[],
 772                           unsigned short cTC[],
 773                           rvec x[], rvec xprime[], rvec v[], rvec f[],
 774                           rvec sd_X[],
 775                           const real tau_t[],
 776                           gmx_bool bFirstHalf, gmx_int64_t step, int seed,
 777                           int* gatindex)
 778 {
 779     gmx_sd_const_t *sdc;
 780     gmx_sd_sigma_t *sig;
 781     /* The random part of the velocity update, generated in the first
 782      * half of the update, needs to be remembered for the second half.
 783      */
 784     rvec  *sd_V;
 785     real   kT;
 786     int    gf = 0, ga = 0, gt = 0;
 787     real   vn = 0, Vmh, Xmh;
 788     real   ism;
 789     int    n, d, ng;
 790
 791     sdc  = sd->sdc;
 792     sig  = sd->sdsig;
 793     sd_V = sd->sd_V;
 794
 795     for (n = start; n < nrend; n++)
 796     {
 797         real rnd[6], rndi[3];
 798         ng  = gatindex ? gatindex[n] : n;
 799         ism = sqrt(invmass[n]);
 800         if (cFREEZE)
 801         {
 802             gf  = cFREEZE[n];
 803         }
 804         if (cACC)
 805         {
 806             ga  = cACC[n];
 807         }
 808         if (cTC)
 809         {
 810             gt  = cTC[n];
 811         }
 812
 813         gmx_rng_cycle_6gaussian_table(step*2+(bFirstHalf ? 1 : 2), ng, seed, RND_SEED_UPDATE, rnd);
 814         if (bInitStep)
 815         {
 816             gmx_rng_cycle_3gaussian_table(step*2, ng, seed, RND_SEED_UPDATE, rndi);
 817         }
 818         for (d = 0; d < DIM; d++)
 819         {
 820             if (bFirstHalf)
 821             {
 822                 vn             = v[n][d];
 823             }
 824             if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
 825             {
 826                 if (bFirstHalf)
 827                 {
 828                     if (bInitStep)
 829                     {
 830                         sd_X[n][d] = ism*sig[gt].X*rndi[d];
 831                     }
 832                     Vmh = sd_X[n][d]*sdc[gt].d/(tau_t[gt]*sdc[gt].c)
 833                         + ism*sig[gt].Yv*rnd[d*2];
 834                     sd_V[n][d] = ism*sig[gt].V*rnd[d*2+1];
 835
 836                     v[n][d] = vn*sdc[gt].em
 837                         + (invmass[n]*f[n][d] + accel[ga][d])*tau_t[gt]*(1 - sdc[gt].em)
 838                         + sd_V[n][d] - sdc[gt].em*Vmh;
 839
 840                     xprime[n][d] = x[n][d] + v[n][d]*tau_t[gt]*(sdc[gt].eph - sdc[gt].emh);
 841                 }
 842                 else
 843                 {
 844                     /* Correct the velocities for the constraints.
 845                      * This operation introduces some inaccuracy,
 846                      * since the velocity is determined from differences in coordinates.
 847                      */
 848                     v[n][d] =
 849                         (xprime[n][d] - x[n][d])/(tau_t[gt]*(sdc[gt].eph - sdc[gt].emh));
 850
 851                     Xmh = sd_V[n][d]*tau_t[gt]*sdc[gt].d/(sdc[gt].em-1)
 852                         + ism*sig[gt].Yx*rnd[d*2];
 853                     sd_X[n][d] = ism*sig[gt].X*rnd[d*2+1];
 854
 855                     xprime[n][d] += sd_X[n][d] - Xmh;
 856
 857                 }
 858             }
 859             else
 860             {
 861                 if (bFirstHalf)
 862                 {
 863                     v[n][d]        = 0.0;
 864                     xprime[n][d]   = x[n][d];
 865                 }
 866             }
 867         }
 868     }
 869 }
 870
 871 static void do_update_bd_Tconsts(double dt, real friction_coefficient,
 872                                  int ngtc, const real ref_t[],
 873                                  real *rf)
 874 {
 875     /* This is separated from the update below, because it is single threaded */
 876     int gt;
 877
 878     if (friction_coefficient != 0)
 879     {
 880         for (gt = 0; gt < ngtc; gt++)
 881         {
 882             rf[gt] = sqrt(2.0*BOLTZ*ref_t[gt]/(friction_coefficient*dt));
 883         }
 884     }
 885     else
 886     {
 887         for (gt = 0; gt < ngtc; gt++)
 888         {
 889             rf[gt] = sqrt(2.0*BOLTZ*ref_t[gt]);
 890         }
 891     }
 892 }
 893
 894 static void do_update_bd(int start, int nrend, double dt,
 895                          ivec nFreeze[],
 896                          real invmass[], unsigned short ptype[],
 897                          unsigned short cFREEZE[], unsigned short cTC[],
 898                          rvec x[], rvec xprime[], rvec v[],
 899                          rvec f[], real friction_coefficient,
 900                          real *rf, gmx_int64_t step, int seed,
 901                          int* gatindex)
 902 {
 903     /* note -- these appear to be full step velocities . . .  */
 904     int    gf = 0, gt = 0;
 905     real   vn;
 906     real   invfr = 0;
 907     int    n, d;
 908
 909     if (friction_coefficient != 0)
 910     {
 911         invfr = 1.0/friction_coefficient;
 912     }
 913
 914     for (n = start; (n < nrend); n++)
 915     {
 916         real rnd[3];
 917         int  ng  = gatindex ? gatindex[n] : n;
 918
 919         if (cFREEZE)
 920         {
 921             gf = cFREEZE[n];
 922         }
 923         if (cTC)
 924         {
 925             gt = cTC[n];
 926         }
 927         gmx_rng_cycle_3gaussian_table(step, ng, seed, RND_SEED_UPDATE, rnd);
 928         for (d = 0; (d < DIM); d++)
 929         {
 930             if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
 931             {
 932                 if (friction_coefficient != 0)
 933                 {
 934                     vn = invfr*f[n][d] + rf[gt]*rnd[d];
 935                 }
 936                 else
 937                 {
 938                     /* NOTE: invmass = 2/(mass*friction_constant*dt) */
 939                     vn = 0.5*invmass[n]*f[n][d]*dt
 940                         + sqrt(0.5*invmass[n])*rf[gt]*rnd[d];
 941                 }
 942
 943                 v[n][d]      = vn;
 944                 xprime[n][d] = x[n][d]+vn*dt;
 945             }
 946             else
 947             {
 948                 v[n][d]      = 0.0;
 949                 xprime[n][d] = x[n][d];
 950             }
 951         }
 952     }
 953 }
 954
 955 static void dump_it_all(FILE gmx_unused *fp, const char gmx_unused *title,
 956                         int gmx_unused natoms, rvec gmx_unused x[], rvec gmx_unused xp[],
 957                         rvec gmx_unused v[], rvec gmx_unused f[])
 958 {
 959 #ifdef DEBUG
 960     if (fp)
 961     {
 962         fprintf(fp, "%s\n", title);
 963         pr_rvecs(fp, 0, "x", x, natoms);
 964         pr_rvecs(fp, 0, "xp", xp, natoms);
 965         pr_rvecs(fp, 0, "v", v, natoms);
 966         pr_rvecs(fp, 0, "f", f, natoms);
 967     }
 968 #endif
 969 }
 970
 971 static void calc_ke_part_normal(rvec v[], t_grpopts *opts, t_mdatoms *md,
 972                                 gmx_ekindata_t *ekind, t_nrnb *nrnb, gmx_bool bEkinAveVel,
 973                                 gmx_bool bSaveEkinOld)
 974 {
 975     int           g;
 976     t_grp_tcstat *tcstat  = ekind->tcstat;
 977     t_grp_acc    *grpstat = ekind->grpstat;
 978     int           nthread, thread;
 979
 980     /* three main: VV with AveVel, vv with AveEkin, leap with AveEkin.  Leap with AveVel is also
 981        an option, but not supported now.  Additionally, if we are doing iterations.
 982        bEkinAveVel: If TRUE, we sum into ekin, if FALSE, into ekinh.
 983        bSavEkinOld: If TRUE (in the case of iteration = bIterate is TRUE), we don't copy over the ekinh_old.
 984        If FALSE, we overrwrite it.
 985      */
 986
 987     /* group velocities are calculated in update_ekindata and
 988      * accumulated in acumulate_groups.
 989      * Now the partial global and groups ekin.
 990      */
 991     for (g = 0; (g < opts->ngtc); g++)
 992     {
 993
 994         if (!bSaveEkinOld)
 995         {
 996             copy_mat(tcstat[g].ekinh, tcstat[g].ekinh_old);
 997         }
 998         if (bEkinAveVel)
 999         {
1000             clear_mat(tcstat[g].ekinf);
1001         }
1002         else
1003         {
1004             clear_mat(tcstat[g].ekinh);
1005         }
1006         if (bEkinAveVel)
1007         {
1008             tcstat[g].ekinscalef_nhc = 1.0; /* need to clear this -- logic is complicated! */
1009         }
1010     }
1011     ekind->dekindl_old = ekind->dekindl;
1012
1013     nthread = gmx_omp_nthreads_get(emntUpdate);
1014
1015 #pragma omp parallel for num_threads(nthread) schedule(static)
1016     for (thread = 0; thread < nthread; thread++)
1017     {
1018         int     start_t, end_t, n;
1019         int     ga, gt;
1020         rvec    v_corrt;
1021         real    hm;
1022         int     d, m;
1023         matrix *ekin_sum;
1024         real   *dekindl_sum;
1025
1026         start_t = ((thread+0)*md->homenr)/nthread;
1027         end_t   = ((thread+1)*md->homenr)/nthread;
1028
1029         ekin_sum    = ekind->ekin_work[thread];
1030         dekindl_sum = ekind->dekindl_work[thread];
1031
1032         for (gt = 0; gt < opts->ngtc; gt++)
1033         {
1034             clear_mat(ekin_sum[gt]);
1035         }
1036         *dekindl_sum = 0.0;
1037
1038         ga = 0;
1039         gt = 0;
1040         for (n = start_t; n < end_t; n++)
1041         {
1042             if (md->cACC)
1043             {
1044                 ga = md->cACC[n];
1045             }
1046             if (md->cTC)
1047             {
1048                 gt = md->cTC[n];
1049             }
1050             hm   = 0.5*md->massT[n];
1051
1052             for (d = 0; (d < DIM); d++)
1053             {
1054                 v_corrt[d]  = v[n][d]  - grpstat[ga].u[d];
1055             }
1056             for (d = 0; (d < DIM); d++)
1057             {
1058                 for (m = 0; (m < DIM); m++)
1059                 {
1060                     /* if we're computing a full step velocity, v_corrt[d] has v(t).  Otherwise, v(t+dt/2) */
1061                     ekin_sum[gt][m][d] += hm*v_corrt[m]*v_corrt[d];
1062                 }
1063             }
1064             if (md->nMassPerturbed && md->bPerturbed[n])
1065             {
1066                 *dekindl_sum +=
1067                     0.5*(md->massB[n] - md->massA[n])*iprod(v_corrt, v_corrt);
1068             }
1069         }
1070     }
1071
1072     ekind->dekindl = 0;
1073     for (thread = 0; thread < nthread; thread++)
1074     {
1075         for (g = 0; g < opts->ngtc; g++)
1076         {
1077             if (bEkinAveVel)
1078             {
1079                 m_add(tcstat[g].ekinf, ekind->ekin_work[thread][g],
1080                       tcstat[g].ekinf);
1081             }
1082             else
1083             {
1084                 m_add(tcstat[g].ekinh, ekind->ekin_work[thread][g],
1085                       tcstat[g].ekinh);
1086             }
1087         }
1088
1089         ekind->dekindl += *ekind->dekindl_work[thread];
1090     }
1091
1092     inc_nrnb(nrnb, eNR_EKIN, md->homenr);
1093 }
1094
1095 static void calc_ke_part_visc(matrix box, rvec x[], rvec v[],
1096                               t_grpopts *opts, t_mdatoms *md,
1097                               gmx_ekindata_t *ekind,
1098                               t_nrnb *nrnb, gmx_bool bEkinAveVel)
1099 {
1100     int           start = 0, homenr = md->homenr;
1101     int           g, d, n, m, gt = 0;
1102     rvec          v_corrt;
1103     real          hm;
1104     t_grp_tcstat *tcstat = ekind->tcstat;
1105     t_cos_acc    *cosacc = &(ekind->cosacc);
1106     real          dekindl;
1107     real          fac, cosz;
1108     double        mvcos;
1109
1110     for (g = 0; g < opts->ngtc; g++)
1111     {
1112         copy_mat(ekind->tcstat[g].ekinh, ekind->tcstat[g].ekinh_old);
1113         clear_mat(ekind->tcstat[g].ekinh);
1114     }
1115     ekind->dekindl_old = ekind->dekindl;
1116
1117     fac     = 2*M_PI/box[ZZ][ZZ];
1118     mvcos   = 0;
1119     dekindl = 0;
1120     for (n = start; n < start+homenr; n++)
1121     {
1122         if (md->cTC)
1123         {
1124             gt = md->cTC[n];
1125         }
1126         hm   = 0.5*md->massT[n];
1127
1128         /* Note that the times of x and v differ by half a step */
1129         /* MRS -- would have to be changed for VV */
1130         cosz         = cos(fac*x[n][ZZ]);
1131         /* Calculate the amplitude of the new velocity profile */
1132         mvcos       += 2*cosz*md->massT[n]*v[n][XX];
1133
1134         copy_rvec(v[n], v_corrt);
1135         /* Subtract the profile for the kinetic energy */
1136         v_corrt[XX] -= cosz*cosacc->vcos;
1137         for (d = 0; (d < DIM); d++)
1138         {
1139             for (m = 0; (m < DIM); m++)
1140             {
1141                 /* if we're computing a full step velocity, v_corrt[d] has v(t).  Otherwise, v(t+dt/2) */
1142                 if (bEkinAveVel)
1143                 {
1144                     tcstat[gt].ekinf[m][d] += hm*v_corrt[m]*v_corrt[d];
1145                 }
1146                 else
1147                 {
1148                     tcstat[gt].ekinh[m][d] += hm*v_corrt[m]*v_corrt[d];
1149                 }
1150             }
1151         }
1152         if (md->nPerturbed && md->bPerturbed[n])
1153         {
1154             /* The minus sign here might be confusing.
1155              * The kinetic contribution from dH/dl doesn't come from
1156              * d m(l)/2 v^2 / dl, but rather from d p^2/2m(l) / dl,
1157              * where p are the momenta. The difference is only a minus sign.
1158              */
1159             dekindl -= 0.5*(md->massB[n] - md->massA[n])*iprod(v_corrt, v_corrt);
1160         }
1161     }
1162     ekind->dekindl = dekindl;
1163     cosacc->mvcos  = mvcos;
1164
1165     inc_nrnb(nrnb, eNR_EKIN, homenr);
1166 }
1167
1168 void calc_ke_part(t_state *state, t_grpopts *opts, t_mdatoms *md,
1169                   gmx_ekindata_t *ekind, t_nrnb *nrnb, gmx_bool bEkinAveVel, gmx_bool bSaveEkinOld)
1170 {
1171     if (ekind->cosacc.cos_accel == 0)
1172     {
1173         calc_ke_part_normal(state->v, opts, md, ekind, nrnb, bEkinAveVel, bSaveEkinOld);
1174     }
1175     else
1176     {
1177         calc_ke_part_visc(state->box, state->x, state->v, opts, md, ekind, nrnb, bEkinAveVel);
1178     }
1179 }
1180
1181 extern void init_ekinstate(ekinstate_t *ekinstate, const t_inputrec *ir)
1182 {
1183     ekinstate->ekin_n = ir->opts.ngtc;
1184     snew(ekinstate->ekinh, ekinstate->ekin_n);
1185     snew(ekinstate->ekinf, ekinstate->ekin_n);
1186     snew(ekinstate->ekinh_old, ekinstate->ekin_n);
1187     snew(ekinstate->ekinscalef_nhc, ekinstate->ekin_n);
1188     snew(ekinstate->ekinscaleh_nhc, ekinstate->ekin_n);
1189     snew(ekinstate->vscale_nhc, ekinstate->ekin_n);
1190     ekinstate->dekindl = 0;
1191     ekinstate->mvcos   = 0;
1192 }
1193
1194 void update_ekinstate(ekinstate_t *ekinstate, gmx_ekindata_t *ekind)
1195 {
1196     int i;
1197
1198     for (i = 0; i < ekinstate->ekin_n; i++)
1199     {
1200         copy_mat(ekind->tcstat[i].ekinh, ekinstate->ekinh[i]);
1201         copy_mat(ekind->tcstat[i].ekinf, ekinstate->ekinf[i]);
1202         copy_mat(ekind->tcstat[i].ekinh_old, ekinstate->ekinh_old[i]);
1203         ekinstate->ekinscalef_nhc[i] = ekind->tcstat[i].ekinscalef_nhc;
1204         ekinstate->ekinscaleh_nhc[i] = ekind->tcstat[i].ekinscaleh_nhc;
1205         ekinstate->vscale_nhc[i]     = ekind->tcstat[i].vscale_nhc;
1206     }
1207
1208     copy_mat(ekind->ekin, ekinstate->ekin_total);
1209     ekinstate->dekindl = ekind->dekindl;
1210     ekinstate->mvcos   = ekind->cosacc.mvcos;
1211
1212 }
1213
1214 void restore_ekinstate_from_state(t_commrec *cr,
1215                                   gmx_ekindata_t *ekind, ekinstate_t *ekinstate)
1216 {
1217     int i, n;
1218
1219     if (MASTER(cr))
1220     {
1221         for (i = 0; i < ekinstate->ekin_n; i++)
1222         {
1223             copy_mat(ekinstate->ekinh[i], ekind->tcstat[i].ekinh);
1224             copy_mat(ekinstate->ekinf[i], ekind->tcstat[i].ekinf);
1225             copy_mat(ekinstate->ekinh_old[i], ekind->tcstat[i].ekinh_old);
1226             ekind->tcstat[i].ekinscalef_nhc = ekinstate->ekinscalef_nhc[i];
1227             ekind->tcstat[i].ekinscaleh_nhc = ekinstate->ekinscaleh_nhc[i];
1228             ekind->tcstat[i].vscale_nhc     = ekinstate->vscale_nhc[i];
1229         }
1230
1231         copy_mat(ekinstate->ekin_total, ekind->ekin);
1232
1233         ekind->dekindl      = ekinstate->dekindl;
1234         ekind->cosacc.mvcos = ekinstate->mvcos;
1235         n                   = ekinstate->ekin_n;
1236     }
1237
1238     if (PAR(cr))
1239     {
1240         gmx_bcast(sizeof(n), &n, cr);
1241         for (i = 0; i < n; i++)
1242         {
1243             gmx_bcast(DIM*DIM*sizeof(ekind->tcstat[i].ekinh[0][0]),
1244                       ekind->tcstat[i].ekinh[0], cr);
1245             gmx_bcast(DIM*DIM*sizeof(ekind->tcstat[i].ekinf[0][0]),
1246                       ekind->tcstat[i].ekinf[0], cr);
1247             gmx_bcast(DIM*DIM*sizeof(ekind->tcstat[i].ekinh_old[0][0]),
1248                       ekind->tcstat[i].ekinh_old[0], cr);
1249
1250             gmx_bcast(sizeof(ekind->tcstat[i].ekinscalef_nhc),
1251                       &(ekind->tcstat[i].ekinscalef_nhc), cr);
1252             gmx_bcast(sizeof(ekind->tcstat[i].ekinscaleh_nhc),
1253                       &(ekind->tcstat[i].ekinscaleh_nhc), cr);
1254             gmx_bcast(sizeof(ekind->tcstat[i].vscale_nhc),
1255                       &(ekind->tcstat[i].vscale_nhc), cr);
1256         }
1257         gmx_bcast(DIM*DIM*sizeof(ekind->ekin[0][0]),
1258                   ekind->ekin[0], cr);
1259
1260         gmx_bcast(sizeof(ekind->dekindl), &ekind->dekindl, cr);
1261         gmx_bcast(sizeof(ekind->cosacc.mvcos), &ekind->cosacc.mvcos, cr);
1262     }
1263 }
1264
1265 void set_deform_reference_box(gmx_update_t upd, gmx_int64_t step, matrix box)
1266 {
1267     upd->deformref_step = step;
1268     copy_mat(box, upd->deformref_box);
1269 }
1270
1271 static void deform(gmx_update_t upd,
1272                    int start, int homenr, rvec x[], matrix box, matrix *scale_tot,
1273                    const t_inputrec *ir, gmx_int64_t step)
1274 {
1275     matrix bnew, invbox, mu;
1276     real   elapsed_time;
1277     int    i, j;
1278
1279     elapsed_time = (step + 1 - upd->deformref_step)*ir->delta_t;
1280     copy_mat(box, bnew);
1281     for (i = 0; i < DIM; i++)
1282     {
1283         for (j = 0; j < DIM; j++)
1284         {
1285             if (ir->deform[i][j] != 0)
1286             {
1287                 bnew[i][j] =
1288                     upd->deformref_box[i][j] + elapsed_time*ir->deform[i][j];
1289             }
1290         }
1291     }
1292     /* We correct the off-diagonal elements,
1293      * which can grow indefinitely during shearing,
1294      * so the shifts do not get messed up.
1295      */
1296     for (i = 1; i < DIM; i++)
1297     {
1298         for (j = i-1; j >= 0; j--)
1299         {
1300             while (bnew[i][j] - box[i][j] > 0.5*bnew[j][j])
1301             {
1302                 rvec_dec(bnew[i], bnew[j]);
1303             }
1304             while (bnew[i][j] - box[i][j] < -0.5*bnew[j][j])
1305             {
1306                 rvec_inc(bnew[i], bnew[j]);
1307             }
1308         }
1309     }
1310     m_inv_ur0(box, invbox);
1311     copy_mat(bnew, box);
1312     mmul_ur0(box, invbox, mu);
1313
1314     for (i = start; i < start+homenr; i++)
1315     {
1316         x[i][XX] = mu[XX][XX]*x[i][XX]+mu[YY][XX]*x[i][YY]+mu[ZZ][XX]*x[i][ZZ];
1317         x[i][YY] = mu[YY][YY]*x[i][YY]+mu[ZZ][YY]*x[i][ZZ];
1318         x[i][ZZ] = mu[ZZ][ZZ]*x[i][ZZ];
1319     }
1320     if (scale_tot != NULL)
1321     {
1322         /* The transposes of the scaling matrices are stored,
1323          * so we need to do matrix multiplication in the inverse order.
1324          */
1325         mmul_ur0(*scale_tot, mu, *scale_tot);
1326     }
1327 }
1328
1329 void update_tcouple(gmx_int64_t       step,
1330                     t_inputrec       *inputrec,
1331                     t_state          *state,
1332                     gmx_ekindata_t   *ekind,
1333                     t_extmass        *MassQ,
1334                     t_mdatoms        *md)
1335
1336 {
1337     gmx_bool   bTCouple = FALSE;
1338     real       dttc;
1339     int        i, start, end, homenr, offset;
1340
1341     /* if using vv with trotter decomposition methods, we do this elsewhere in the code */
1342     if (inputrec->etc != etcNO &&
1343         !(IR_NVT_TROTTER(inputrec) || IR_NPT_TROTTER(inputrec) || IR_NPH_TROTTER(inputrec)))
1344     {
1345         /* We should only couple after a step where energies were determined (for leapfrog versions)
1346            or the step energies are determined, for velocity verlet versions */
1347
1348         if (EI_VV(inputrec->eI))
1349         {
1350             offset = 0;
1351         }
1352         else
1353         {
1354             offset = 1;
1355         }
1356         bTCouple = (inputrec->nsttcouple == 1 ||
1357                     do_per_step(step+inputrec->nsttcouple-offset,
1358                                 inputrec->nsttcouple));
1359     }
1360
1361     if (bTCouple)
1362     {
1363         dttc = inputrec->nsttcouple*inputrec->delta_t;
1364
1365         switch (inputrec->etc)
1366         {
1367             case etcNO:
1368                 break;
1369             case etcBERENDSEN:
1370                 berendsen_tcoupl(inputrec, ekind, dttc);
1371                 break;
1372             case etcNOSEHOOVER:
1373                 nosehoover_tcoupl(&(inputrec->opts), ekind, dttc,
1374                                   state->nosehoover_xi, state->nosehoover_vxi, MassQ);
1375                 break;
1376             case etcVRESCALE:
1377                 vrescale_tcoupl(inputrec, step, ekind, dttc,
1378                                 state->therm_integral);
1379                 break;
1380         }
1381         /* rescale in place here */
1382         if (EI_VV(inputrec->eI))
1383         {
1384             rescale_velocities(ekind, md, 0, md->homenr, state->v);
1385         }
1386     }
1387     else
1388     {
1389         /* Set the T scaling lambda to 1 to have no scaling */
1390         for (i = 0; (i < inputrec->opts.ngtc); i++)
1391         {
1392             ekind->tcstat[i].lambda = 1.0;
1393         }
1394     }
1395 }
1396
1397 void update_pcouple(FILE             *fplog,
1398                     gmx_int64_t       step,
1399                     t_inputrec       *inputrec,
1400                     t_state          *state,
1401                     matrix            pcoupl_mu,
1402                     matrix            M,
1403                     gmx_bool          bInitStep)
1404 {
1405     gmx_bool   bPCouple = FALSE;
1406     real       dtpc     = 0;
1407     int        i;
1408
1409     /* if using Trotter pressure, we do this in coupling.c, so we leave it false. */
1410     if (inputrec->epc != epcNO && (!(IR_NPT_TROTTER(inputrec) || IR_NPH_TROTTER(inputrec))))
1411     {
1412         /* We should only couple after a step where energies were determined */
1413         bPCouple = (inputrec->nstpcouple == 1 ||
1414                     do_per_step(step+inputrec->nstpcouple-1,
1415                                 inputrec->nstpcouple));
1416     }
1417
1418     clear_mat(pcoupl_mu);
1419     for (i = 0; i < DIM; i++)
1420     {
1421         pcoupl_mu[i][i] = 1.0;
1422     }
1423
1424     clear_mat(M);
1425
1426     if (bPCouple)
1427     {
1428         dtpc = inputrec->nstpcouple*inputrec->delta_t;
1429
1430         switch (inputrec->epc)
1431         {
1432             /* We can always pcoupl, even if we did not sum the energies
1433              * the previous step, since state->pres_prev is only updated
1434              * when the energies have been summed.
1435              */
1436             case (epcNO):
1437                 break;
1438             case (epcBERENDSEN):
1439                 if (!bInitStep)
1440                 {
1441                     berendsen_pcoupl(fplog, step, inputrec, dtpc, state->pres_prev, state->box,
1442                                      pcoupl_mu);
1443                 }
1444                 break;
1445             case (epcPARRINELLORAHMAN):
1446                 parrinellorahman_pcoupl(fplog, step, inputrec, dtpc, state->pres_prev,
1447                                         state->box, state->box_rel, state->boxv,
1448                                         M, pcoupl_mu, bInitStep);
1449                 break;
1450             default:
1451                 break;
1452         }
1453     }
1454 }
1455
1456 static rvec *get_xprime(const t_state *state, gmx_update_t upd)
1457 {
1458     if (state->nalloc > upd->xp_nalloc)
1459     {
1460         upd->xp_nalloc = state->nalloc;
1461         srenew(upd->xp, upd->xp_nalloc);
1462     }
1463
1464     return upd->xp;
1465 }
1466
1467 static void combine_forces(gmx_update_t upd,
1468                            int nstcalclr,
1469                            gmx_constr_t constr,
1470                            t_inputrec *ir, t_mdatoms *md, t_idef *idef,
1471                            t_commrec *cr,
1472                            gmx_int64_t step,
1473                            t_state *state, gmx_bool bMolPBC,
1474                            int start, int nrend,
1475                            rvec f[], rvec f_lr[],
1476                            tensor *vir_lr_constr,
1477                            t_nrnb *nrnb)
1478 {
1479     int  i, d;
1480
1481     /* f contains the short-range forces + the long range forces
1482      * which are stored separately in f_lr.
1483      */
1484
1485     if (constr != NULL && vir_lr_constr != NULL &&
1486         !(ir->eConstrAlg == econtSHAKE && ir->epc == epcNO))
1487     {
1488         /* We need to constrain the LR forces separately,
1489          * because due to the different pre-factor for the SR and LR
1490          * forces in the update algorithm, we have to correct
1491          * the constraint virial for the nstcalclr-1 extra f_lr.
1492          * Constrain only the additional LR part of the force.
1493          */
1494         /* MRS -- need to make sure this works with trotter integration -- the constraint calls may not be right.*/
1495         rvec *xp;
1496         real  fac;
1497         int   gf = 0;
1498
1499         xp  = get_xprime(state, upd);
1500
1501         fac = (nstcalclr - 1)*ir->delta_t*ir->delta_t;
1502
1503         for (i = 0; i < md->homenr; i++)
1504         {
1505             if (md->cFREEZE != NULL)
1506             {
1507                 gf = md->cFREEZE[i];
1508             }
1509             for (d = 0; d < DIM; d++)
1510             {
1511                 if ((md->ptype[i] != eptVSite) &&
1512                     (md->ptype[i] != eptShell) &&
1513                     !ir->opts.nFreeze[gf][d])
1514                 {
1515                     xp[i][d] = state->x[i][d] + fac*f_lr[i][d]*md->invmass[i];
1516                 }
1517             }
1518         }
1519         constrain(NULL, FALSE, FALSE, constr, idef, ir, NULL, cr, step, 0, 1.0, md,
1520                   state->x, xp, xp, bMolPBC, state->box, state->lambda[efptBONDED], NULL,
1521                   NULL, vir_lr_constr, nrnb, econqCoord, ir->epc == epcMTTK, state->veta, state->veta);
1522     }
1523
1524     /* Add nstcalclr-1 times the LR force to the sum of both forces
1525      * and store the result in forces_lr.
1526      */
1527     for (i = start; i < nrend; i++)
1528     {
1529         for (d = 0; d < DIM; d++)
1530         {
1531             f_lr[i][d] = f[i][d] + (nstcalclr - 1)*f_lr[i][d];
1532         }
1533     }
1534 }
1535
1536 void update_constraints(FILE             *fplog,
1537                         gmx_int64_t       step,
1538                         real             *dvdlambda, /* the contribution to be added to the bonded interactions */
1539                         t_inputrec       *inputrec,  /* input record and box stuff      */
1540                         gmx_ekindata_t   *ekind,
1541                         t_mdatoms        *md,
1542                         t_state          *state,
1543                         gmx_bool          bMolPBC,
1544                         t_graph          *graph,
1545                         rvec              force[],   /* forces on home particles */
1546                         t_idef           *idef,
1547                         tensor            vir_part,
1548                         t_commrec        *cr,
1549                         t_nrnb           *nrnb,
1550                         gmx_wallcycle_t   wcycle,
1551                         gmx_update_t      upd,
1552                         gmx_constr_t      constr,
1553                         gmx_bool          bFirstHalf,
1554                         gmx_bool          bCalcVir,
1555                         real              vetanew)
1556 {
1557     gmx_bool             bExtended, bLastStep, bLog = FALSE, bEner = FALSE, bDoConstr = FALSE;
1558     double               dt;
1559     real                 dt_1;
1560     int                  start, homenr, nrend, i, n, m, g, d;
1561     tensor               vir_con;
1562     rvec                *vbuf, *xprime = NULL;
1563     int                  nth, th;
1564
1565     if (constr)
1566     {
1567         bDoConstr = TRUE;
1568     }
1569     if (bFirstHalf && !EI_VV(inputrec->eI))
1570     {
1571         bDoConstr = FALSE;
1572     }
1573
1574     /* for now, SD update is here -- though it really seems like it
1575        should be reformulated as a velocity verlet method, since it has two parts */
1576
1577     start  = 0;
1578     homenr = md->homenr;
1579     nrend  = start+homenr;
1580
1581     dt   = inputrec->delta_t;
1582     dt_1 = 1.0/dt;
1583
1584     /*
1585      *  Steps (7C, 8C)
1586      *  APPLY CONSTRAINTS:
1587      *  BLOCK SHAKE
1588
1589      * When doing PR pressure coupling we have to constrain the
1590      * bonds in each iteration. If we are only using Nose-Hoover tcoupling
1591      * it is enough to do this once though, since the relative velocities
1592      * after this will be normal to the bond vector
1593      */
1594
1595     if (bDoConstr)
1596     {
1597         /* clear out constraints before applying */
1598         clear_mat(vir_part);
1599
1600         xprime = get_xprime(state, upd);
1601
1602         bLastStep = (step == inputrec->init_step+inputrec->nsteps);
1603         bLog      = (do_per_step(step, inputrec->nstlog) || bLastStep || (step < 0));
1604         bEner     = (do_per_step(step, inputrec->nstenergy) || bLastStep);
1605         /* Constrain the coordinates xprime */
1606         wallcycle_start(wcycle, ewcCONSTR);
1607         if (EI_VV(inputrec->eI) && bFirstHalf)
1608         {
1609             constrain(NULL, bLog, bEner, constr, idef,
1610                       inputrec, ekind, cr, step, 1, 1.0, md,
1611                       state->x, state->v, state->v,
1612                       bMolPBC, state->box,
1613                       state->lambda[efptBONDED], dvdlambda,
1614                       NULL, bCalcVir ? &vir_con : NULL, nrnb, econqVeloc,
1615                       inputrec->epc == epcMTTK, state->veta, vetanew);
1616         }
1617         else
1618         {
1619             constrain(NULL, bLog, bEner, constr, idef,
1620                       inputrec, ekind, cr, step, 1, 1.0, md,
1621                       state->x, xprime, NULL,
1622                       bMolPBC, state->box,
1623                       state->lambda[efptBONDED], dvdlambda,
1624                       state->v, bCalcVir ? &vir_con : NULL, nrnb, econqCoord,
1625                       inputrec->epc == epcMTTK, state->veta, state->veta);
1626         }
1627         wallcycle_stop(wcycle, ewcCONSTR);
1628
1629         where();
1630
1631         dump_it_all(fplog, "After Shake",
1632                     state->natoms, state->x, xprime, state->v, force);
1633
1634         if (bCalcVir)
1635         {
1636             if (inputrec->eI == eiSD2)
1637             {
1638                 /* A correction factor eph is needed for the SD constraint force */
1639                 /* Here we can, unfortunately, not have proper corrections
1640                  * for different friction constants, so we use the first one.
1641                  */
1642                 for (i = 0; i < DIM; i++)
1643                 {
1644                     for (m = 0; m < DIM; m++)
1645                     {
1646                         vir_part[i][m] += upd->sd->sdc[0].eph*vir_con[i][m];
1647                     }
1648                 }
1649             }
1650             else
1651             {
1652                 m_add(vir_part, vir_con, vir_part);
1653             }
1654             if (debug)
1655             {
1656                 pr_rvecs(debug, 0, "constraint virial", vir_part, DIM);
1657             }
1658         }
1659     }
1660
1661     where();
1662
1663     if (inputrec->eI == eiSD1 && bDoConstr && !bFirstHalf)
1664     {
1665         wallcycle_start(wcycle, ewcUPDATE);
1666         xprime = get_xprime(state, upd);
1667
1668         nth = gmx_omp_nthreads_get(emntUpdate);
1669
1670 #pragma omp parallel for num_threads(nth) schedule(static)
1671
1672         for (th = 0; th < nth; th++)
1673         {
1674             int start_th, end_th;
1675
1676             start_th = start + ((nrend-start)* th   )/nth;
1677             end_th   = start + ((nrend-start)*(th+1))/nth;
1678
1679             /* The second part of the SD integration */
1680             do_update_sd1(upd->sd,
1681                           start_th, end_th, dt,
1682                           inputrec->opts.acc, inputrec->opts.nFreeze,
1683                           md->invmass, md->ptype,
1684                           md->cFREEZE, md->cACC, md->cTC,
1685                           state->x, xprime, state->v, force,
1686                           inputrec->opts.ngtc, inputrec->opts.ref_t,
1687                           bDoConstr, FALSE,
1688                           step, inputrec->ld_seed,
1689                           DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
1690         }
1691         inc_nrnb(nrnb, eNR_UPDATE, homenr);
1692         wallcycle_stop(wcycle, ewcUPDATE);
1693
1694         if (bDoConstr)
1695         {
1696             /* Constrain the coordinates xprime for half a time step */
1697             wallcycle_start(wcycle, ewcCONSTR);
1698
1699             constrain(NULL, bLog, bEner, constr, idef,
1700                       inputrec, NULL, cr, step, 1, 0.5, md,
1701                       state->x, xprime, NULL,
1702                       bMolPBC, state->box,
1703                       state->lambda[efptBONDED], dvdlambda,
1704                       state->v, NULL, nrnb, econqCoord, FALSE, 0, 0);
1705
1706             wallcycle_stop(wcycle, ewcCONSTR);
1707         }
1708     }
1709
1710     if ((inputrec->eI == eiSD2) && !(bFirstHalf))
1711     {
1712         wallcycle_start(wcycle, ewcUPDATE);
1713         xprime = get_xprime(state, upd);
1714
1715         nth = gmx_omp_nthreads_get(emntUpdate);
1716
1717 #pragma omp parallel for num_threads(nth) schedule(static)
1718         for (th = 0; th < nth; th++)
1719         {
1720             int start_th, end_th;
1721
1722             start_th = start + ((nrend-start)* th   )/nth;
1723             end_th   = start + ((nrend-start)*(th+1))/nth;
1724
1725             /* The second part of the SD integration */
1726             do_update_sd2(upd->sd,
1727                           FALSE, start_th, end_th,
1728                           inputrec->opts.acc, inputrec->opts.nFreeze,
1729                           md->invmass, md->ptype,
1730                           md->cFREEZE, md->cACC, md->cTC,
1731                           state->x, xprime, state->v, force, state->sd_X,
1732                           inputrec->opts.tau_t,
1733                           FALSE, step, inputrec->ld_seed,
1734                           DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
1735         }
1736         inc_nrnb(nrnb, eNR_UPDATE, homenr);
1737         wallcycle_stop(wcycle, ewcUPDATE);
1738
1739         if (bDoConstr)
1740         {
1741             /* Constrain the coordinates xprime */
1742             wallcycle_start(wcycle, ewcCONSTR);
1743             constrain(NULL, bLog, bEner, constr, idef,
1744                       inputrec, NULL, cr, step, 1, 1.0, md,
1745                       state->x, xprime, NULL,
1746                       bMolPBC, state->box,
1747                       state->lambda[efptBONDED], dvdlambda,
1748                       NULL, NULL, nrnb, econqCoord, FALSE, 0, 0);
1749             wallcycle_stop(wcycle, ewcCONSTR);
1750         }
1751     }
1752
1753
1754     /* We must always unshift after updating coordinates; if we did not shake
1755        x was shifted in do_force */
1756
1757     if (!(bFirstHalf)) /* in the first half of vv, no shift. */
1758     {
1759         if (graph && (graph->nnodes > 0))
1760         {
1761             unshift_x(graph, state->box, state->x, upd->xp);
1762             if (TRICLINIC(state->box))
1763             {
1764                 inc_nrnb(nrnb, eNR_SHIFTX, 2*graph->nnodes);
1765             }
1766             else
1767             {
1768                 inc_nrnb(nrnb, eNR_SHIFTX, graph->nnodes);
1769             }
1770         }
1771         else
1772         {
1773 #pragma omp parallel for num_threads(gmx_omp_nthreads_get(emntUpdate)) schedule(static)
1774             for (i = start; i < nrend; i++)
1775             {
1776                 copy_rvec(upd->xp[i], state->x[i]);
1777             }
1778         }
1779
1780         dump_it_all(fplog, "After unshift",
1781                     state->natoms, state->x, upd->xp, state->v, force);
1782     }
1783 /* ############# END the update of velocities and positions ######### */
1784 }
1785
1786 void update_box(FILE             *fplog,
1787                 gmx_int64_t       step,
1788                 t_inputrec       *inputrec,  /* input record and box stuff      */
1789                 t_mdatoms        *md,
1790                 t_state          *state,
1791                 rvec              force[],   /* forces on home particles */
1792                 matrix           *scale_tot,
1793                 matrix            pcoupl_mu,
1794                 t_nrnb           *nrnb,
1795                 gmx_update_t      upd)
1796 {
1797     gmx_bool             bExtended, bLastStep, bLog = FALSE, bEner = FALSE;
1798     double               dt;
1799     real                 dt_1;
1800     int                  start, homenr, nrend, i, n, m, g;
1801     tensor               vir_con;
1802
1803     start  = 0;
1804     homenr = md->homenr;
1805     nrend  = start+homenr;
1806
1807     bExtended =
1808         (inputrec->etc == etcNOSEHOOVER) ||
1809         (inputrec->epc == epcPARRINELLORAHMAN) ||
1810         (inputrec->epc == epcMTTK);
1811
1812     dt = inputrec->delta_t;
1813
1814     where();
1815
1816     /* now update boxes */
1817     switch (inputrec->epc)
1818     {
1819         case (epcNO):
1820             break;
1821         case (epcBERENDSEN):
1822             berendsen_pscale(inputrec, pcoupl_mu, state->box, state->box_rel,
1823                              start, homenr, state->x, md->cFREEZE, nrnb);
1824             break;
1825         case (epcPARRINELLORAHMAN):
1826             /* The box velocities were updated in do_pr_pcoupl in the update
1827              * iteration, but we dont change the box vectors until we get here
1828              * since we need to be able to shift/unshift above.
1829              */
1830             for (i = 0; i < DIM; i++)
1831             {
1832                 for (m = 0; m <= i; m++)
1833                 {
1834                     state->box[i][m] += dt*state->boxv[i][m];
1835                 }
1836             }
1837             preserve_box_shape(inputrec, state->box_rel, state->box);
1838
1839             /* Scale the coordinates */
1840             for (n = start; (n < start+homenr); n++)
1841             {
1842                 tmvmul_ur0(pcoupl_mu, state->x[n], state->x[n]);
1843             }
1844             break;
1845         case (epcMTTK):
1846             switch (inputrec->epct)
1847             {
1848                 case (epctISOTROPIC):
1849                     /* DIM * eta = ln V.  so DIM*eta_new = DIM*eta_old + DIM*dt*veta =>
1850                        ln V_new = ln V_old + 3*dt*veta => V_new = V_old*exp(3*dt*veta) =>
1851                        Side length scales as exp(veta*dt) */
1852
1853                     msmul(state->box, exp(state->veta*dt), state->box);
1854
1855                     /* Relate veta to boxv.  veta = d(eta)/dT = (1/DIM)*1/V dV/dT.
1856                        o               If we assume isotropic scaling, and box length scaling
1857                        factor L, then V = L^DIM (det(M)).  So dV/dt = DIM
1858                        L^(DIM-1) dL/dt det(M), and veta = (1/L) dL/dt.  The
1859                        determinant of B is L^DIM det(M), and the determinant
1860                        of dB/dt is (dL/dT)^DIM det (M).  veta will be
1861                        (det(dB/dT)/det(B))^(1/3).  Then since M =
1862                        B_new*(vol_new)^(1/3), dB/dT_new = (veta_new)*B(new). */
1863
1864                     msmul(state->box, state->veta, state->boxv);
1865                     break;
1866                 default:
1867                     break;
1868             }
1869             break;
1870         default:
1871             break;
1872     }
1873
1874     if ((!(IR_NPT_TROTTER(inputrec) || IR_NPH_TROTTER(inputrec))) && scale_tot)
1875     {
1876         /* The transposes of the scaling matrices are stored,
1877          * therefore we need to reverse the order in the multiplication.
1878          */
1879         mmul_ur0(*scale_tot, pcoupl_mu, *scale_tot);
1880     }
1881
1882     if (DEFORM(*inputrec))
1883     {
1884         deform(upd, start, homenr, state->x, state->box, scale_tot, inputrec, step);
1885     }
1886     where();
1887     dump_it_all(fplog, "After update",
1888                 state->natoms, state->x, upd->xp, state->v, force);
1889 }
1890
1891 void update_coords(FILE             *fplog,
1892                    gmx_int64_t       step,
1893                    t_inputrec       *inputrec,  /* input record and box stuff   */
1894                    t_mdatoms        *md,
1895                    t_state          *state,
1896                    gmx_bool          bMolPBC,
1897                    rvec             *f,    /* forces on home particles */
1898                    gmx_bool          bDoLR,
1899                    rvec             *f_lr,
1900                    tensor           *vir_lr_constr,
1901                    t_fcdata         *fcd,
1902                    gmx_ekindata_t   *ekind,
1903                    matrix            M,
1904                    gmx_update_t      upd,
1905                    gmx_bool          bInitStep,
1906                    int               UpdatePart,
1907                    t_commrec        *cr, /* these shouldn't be here -- need to think about it */
1908                    t_nrnb           *nrnb,
1909                    gmx_constr_t      constr,
1910                    t_idef           *idef)
1911 {
1912     gmx_bool          bNH, bPR, bLastStep, bLog = FALSE, bEner = FALSE, bDoConstr = FALSE;
1913     double            dt, alpha;
1914     real             *imass, *imassin;
1915     rvec             *force;
1916     real              dt_1;
1917     int               start, homenr, nrend, i, j, d, n, m, g;
1918     int               blen0, blen1, iatom, jatom, nshake, nsettle, nconstr, nexpand;
1919     int              *icom = NULL;
1920     tensor            vir_con;
1921     rvec             *vcom, *xcom, *vall, *xall, *xin, *vin, *forcein, *fall, *xpall, *xprimein, *xprime;
1922     int               nth, th;
1923
1924     bDoConstr = (NULL != constr);
1925
1926     /* Running the velocity half does nothing except for velocity verlet */
1927     if ((UpdatePart == etrtVELOCITY1 || UpdatePart == etrtVELOCITY2) &&
1928         !EI_VV(inputrec->eI))
1929     {
1930         gmx_incons("update_coords called for velocity without VV integrator");
1931     }
1932
1933     start  = 0;
1934     homenr = md->homenr;
1935     nrend  = start+homenr;
1936
1937     xprime = get_xprime(state, upd);
1938
1939     dt   = inputrec->delta_t;
1940     dt_1 = 1.0/dt;
1941
1942     /* We need to update the NMR restraint history when time averaging is used */
1943     if (state->flags & (1<<estDISRE_RM3TAV))
1944     {
1945         update_disres_history(fcd, &state->hist);
1946     }
1947     if (state->flags & (1<<estORIRE_DTAV))
1948     {
1949         update_orires_history(fcd, &state->hist);
1950     }
1951
1952
1953     bNH = inputrec->etc == etcNOSEHOOVER;
1954     bPR = ((inputrec->epc == epcPARRINELLORAHMAN) || (inputrec->epc == epcMTTK));
1955
1956     if (bDoLR && inputrec->nstcalclr > 1 && !EI_VV(inputrec->eI))  /* get this working with VV? */
1957     {
1958         /* Store the total force + nstcalclr-1 times the LR force
1959          * in forces_lr, so it can be used in a normal update algorithm
1960          * to produce twin time stepping.
1961          */
1962         /* is this correct in the new construction? MRS */
1963         combine_forces(upd,
1964                        inputrec->nstcalclr, constr, inputrec, md, idef, cr,
1965                        step, state, bMolPBC,
1966                        start, nrend, f, f_lr, vir_lr_constr, nrnb);
1967         force = f_lr;
1968     }
1969     else
1970     {
1971         force = f;
1972     }
1973
1974     /* ############# START The update of velocities and positions ######### */
1975     where();
1976     dump_it_all(fplog, "Before update",
1977                 state->natoms, state->x, xprime, state->v, force);
1978
1979     if (inputrec->eI == eiSD2)
1980     {
1981         check_sd2_work_data_allocation(upd->sd, nrend);
1982
1983         do_update_sd2_Tconsts(upd->sd,
1984                               inputrec->opts.ngtc,
1985                               inputrec->opts.tau_t,
1986                               inputrec->opts.ref_t);
1987     }
1988     if (inputrec->eI == eiBD)
1989     {
1990         do_update_bd_Tconsts(dt, inputrec->bd_fric,
1991                              inputrec->opts.ngtc, inputrec->opts.ref_t,
1992                              upd->sd->bd_rf);
1993     }
1994
1995     nth = gmx_omp_nthreads_get(emntUpdate);
1996
1997 #pragma omp parallel for num_threads(nth) schedule(static) private(alpha)
1998     for (th = 0; th < nth; th++)
1999     {
2000         int start_th, end_th;
2001
2002         start_th = start + ((nrend-start)* th   )/nth;
2003         end_th   = start + ((nrend-start)*(th+1))/nth;
2004
2005         switch (inputrec->eI)
2006         {
2007             case (eiMD):
2008                 if (ekind->cosacc.cos_accel == 0)
2009                 {
2010                     do_update_md(start_th, end_th, dt,
2011                                  ekind->tcstat, state->nosehoover_vxi,
2012                                  ekind->bNEMD, ekind->grpstat, inputrec->opts.acc,
2013                                  inputrec->opts.nFreeze,
2014                                  md->invmass, md->ptype,
2015                                  md->cFREEZE, md->cACC, md->cTC,
2016                                  state->x, xprime, state->v, force, M,
2017                                  bNH, bPR);
2018                 }
2019                 else
2020                 {
2021                     do_update_visc(start_th, end_th, dt,
2022                                    ekind->tcstat, state->nosehoover_vxi,
2023                                    md->invmass, md->ptype,
2024                                    md->cTC, state->x, xprime, state->v, force, M,
2025                                    state->box,
2026                                    ekind->cosacc.cos_accel,
2027                                    ekind->cosacc.vcos,
2028                                    bNH, bPR);
2029                 }
2030                 break;
2031             case (eiSD1):
2032                 /* With constraints, the SD1 update is done in 2 parts */
2033                 do_update_sd1(upd->sd,
2034                               start_th, end_th, dt,
2035                               inputrec->opts.acc, inputrec->opts.nFreeze,
2036                               md->invmass, md->ptype,
2037                               md->cFREEZE, md->cACC, md->cTC,
2038                               state->x, xprime, state->v, force,
2039                               inputrec->opts.ngtc, inputrec->opts.ref_t,
2040                               bDoConstr, TRUE,
2041                               step, inputrec->ld_seed, DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
2042                 break;
2043             case (eiSD2):
2044                 /* The SD2 update is always done in 2 parts,
2045                  * because an extra constraint step is needed
2046                  */
2047                 do_update_sd2(upd->sd,
2048                               bInitStep, start_th, end_th,
2049                               inputrec->opts.acc, inputrec->opts.nFreeze,
2050                               md->invmass, md->ptype,
2051                               md->cFREEZE, md->cACC, md->cTC,
2052                               state->x, xprime, state->v, force, state->sd_X,
2053                               inputrec->opts.tau_t,
2054                               TRUE, step, inputrec->ld_seed,
2055                               DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
2056                 break;
2057             case (eiBD):
2058                 do_update_bd(start_th, end_th, dt,
2059                              inputrec->opts.nFreeze, md->invmass, md->ptype,
2060                              md->cFREEZE, md->cTC,
2061                              state->x, xprime, state->v, force,
2062                              inputrec->bd_fric,
2063                              upd->sd->bd_rf,
2064                              step, inputrec->ld_seed, DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
2065                 break;
2066             case (eiVV):
2067             case (eiVVAK):
2068                 alpha = 1.0 + DIM/((double)inputrec->opts.nrdf[0]); /* assuming barostat coupled to group 0. */
2069                 switch (UpdatePart)
2070                 {
2071                     case etrtVELOCITY1:
2072                     case etrtVELOCITY2:
2073                         do_update_vv_vel(start_th, end_th, dt,
2074                                          inputrec->opts.acc, inputrec->opts.nFreeze,
2075                                          md->invmass, md->ptype,
2076                                          md->cFREEZE, md->cACC,
2077                                          state->v, force,
2078                                          (bNH || bPR), state->veta, alpha);
2079                         break;
2080                     case etrtPOSITION:
2081                         do_update_vv_pos(start_th, end_th, dt,
2082                                          inputrec->opts.nFreeze,
2083                                          md->ptype, md->cFREEZE,
2084                                          state->x, xprime, state->v,
2085                                          (bNH || bPR), state->veta);
2086                         break;
2087                 }
2088                 break;
2089             default:
2090                 gmx_fatal(FARGS, "Don't know how to update coordinates");
2091                 break;
2092         }
2093     }
2094
2095 }
2096
2097
2098 void correct_ekin(FILE *log, int start, int end, rvec v[], rvec vcm, real mass[],
2099                   real tmass, tensor ekin)
2100 {
2101     /*
2102      * This is a debugging routine. It should not be called for production code
2103      *
2104      * The kinetic energy should calculated according to:
2105      *   Ekin = 1/2 m (v-vcm)^2
2106      * However the correction is not always applied, since vcm may not be
2107      * known in time and we compute
2108      *   Ekin' = 1/2 m v^2 instead
2109      * This can be corrected afterwards by computing
2110      *   Ekin = Ekin' + 1/2 m ( -2 v vcm + vcm^2)
2111      * or in hsorthand:
2112      *   Ekin = Ekin' - m v vcm + 1/2 m vcm^2
2113      */
2114     int    i, j, k;
2115     real   m, tm;
2116     rvec   hvcm, mv;
2117     tensor dekin;
2118
2119     /* Local particles */
2120     clear_rvec(mv);
2121
2122     /* Processor dependent part. */
2123     tm = 0;
2124     for (i = start; (i < end); i++)
2125     {
2126         m      = mass[i];
2127         tm    += m;
2128         for (j = 0; (j < DIM); j++)
2129         {
2130             mv[j] += m*v[i][j];
2131         }
2132     }
2133     /* Shortcut */
2134     svmul(1/tmass, vcm, vcm);
2135     svmul(0.5, vcm, hvcm);
2136     clear_mat(dekin);
2137     for (j = 0; (j < DIM); j++)
2138     {
2139         for (k = 0; (k < DIM); k++)
2140         {
2141             dekin[j][k] += vcm[k]*(tm*hvcm[j]-mv[j]);
2142         }
2143     }
2144     pr_rvecs(log, 0, "dekin", dekin, DIM);
2145     pr_rvecs(log, 0, " ekin", ekin, DIM);
2146     fprintf(log, "dekin = %g, ekin = %g  vcm = (%8.4f %8.4f %8.4f)\n",
2147             trace(dekin), trace(ekin), vcm[XX], vcm[YY], vcm[ZZ]);
2148     fprintf(log, "mv = (%8.4f %8.4f %8.4f)\n",
2149             mv[XX], mv[YY], mv[ZZ]);
2150 }
2151
2152 extern gmx_bool update_randomize_velocities(t_inputrec *ir, gmx_int64_t step, const t_commrec *cr,
2153                                             t_mdatoms *md, t_state *state, gmx_update_t upd, gmx_constr_t constr)
2154 {
2155
2156     int  i;
2157     real rate = (ir->delta_t)/ir->opts.tau_t[0];
2158
2159     if (ir->etc == etcANDERSEN && constr != NULL)
2160     {
2161         gmx_fatal(FARGS, "Normal Andersen is currently not supported with constraints, use massive Andersen instead");
2162     }
2163
2164     /* proceed with andersen if 1) it's fixed probability per
2165        particle andersen or 2) it's massive andersen and it's tau_t/dt */
2166     if ((ir->etc == etcANDERSEN) || do_per_step(step, (int)(1.0/rate)))
2167     {
2168         andersen_tcoupl(ir, step, cr, md, state, rate,
2169                         upd->sd->randomize_group, upd->sd->boltzfac);
2170         return TRUE;
2171     }
2172     return FALSE;
2173 }