src/gromacs/mdlib/update.c

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