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