src/mdlib/update.c

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