src/mdlib/update.c

   1 /*
   2  * $Id$
   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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_update_c = "$Id$";
  30
  31 #include <stdio.h>
  32 #include <math.h>
  33
  34 #include "assert.h"
  35 #include "sysstuff.h"
  36 #include "smalloc.h"
  37 #include "typedefs.h"
  38 #include "nrnb.h"
  39 #include "physics.h"
  40 #include "macros.h"
  41 #include "vveclib.h"
  42 #include "vec.h"
  43 #include "main.h"
  44 #include "confio.h"
  45 #include "update.h"
  46 #include "random.h"
  47 #include "futil.h"
  48 #include "mshift.h"
  49 #include "tgroup.h"
  50 #include "force.h"
  51 #include "names.h"
  52 #include "txtdump.h"
  53 #include "mdrun.h"
  54 #include "copyrite.h"
  55 #include "constr.h"
  56 #include "edsam.h"
  57 #include "pull.h"
  58
  59 static void calc_g(rvec x_unc,rvec x_cons,rvec g,double mdt_2)
  60 {
  61   int d;
  62
  63   for(d=0; (d<DIM);d++)
  64     g[d]=(x_cons[d]-x_unc[d])*mdt_2;
  65 }
  66
  67 static void do_shake_corr(rvec xold,rvec x,rvec v,double dt_1)
  68 {
  69   int    d;
  70
  71   for(d=0; (d<DIM);d++)
  72     v[d]=((double) x[d]-(double) xold[d])*dt_1;
  73 }
  74
  75 static void do_both(rvec xold,rvec x_unc,rvec x,rvec g,
  76                     rvec v,real mdt_2,real dt_1)
  77 {
  78   real xx,yy,zz;
  79
  80   xx=x[XX];
  81   yy=x[YY];
  82   zz=x[ZZ];
  83   g[XX]=(xx-x_unc[XX])*mdt_2;
  84   g[YY]=(yy-x_unc[YY])*mdt_2;
  85   g[ZZ]=(zz-x_unc[ZZ])*mdt_2;
  86   v[XX]=(xx-xold [XX])*dt_1;
  87   v[YY]=(yy-xold [YY])*dt_1;
  88   v[ZZ]=(zz-xold [ZZ])*dt_1;
  89 }
  90
  91 static void do_update(int start,int homenr,double dt,
  92                       rvec lamb[],t_grp_acc gstat[],
  93                       rvec accel[],ivec nFreeze[],
  94                       real invmass[],unsigned short ptype[],
  95                       unsigned short cFREEZE[],unsigned short cACC[],unsigned short cTC[],
  96                       rvec x[],rvec xprime[],rvec v[],rvec vold[],rvec f[])
  97 {
  98   double w_dt;
  99   int    gf,ga,gt;
 100   real   vn,vv,va,vb;
 101   real   uold,lg;
 102   int    n,d;
 103
 104   for (n=start; n<start+homenr; n++) {
 105     w_dt = invmass[n]*dt;
 106     gf   = cFREEZE[n];
 107     ga   = cACC[n];
 108     gt   = cTC[n];
 109
 110     for (d=0; d<DIM; d++) {
 111       vn             = v[n][d];
 112       lg             = lamb[gt][d];
 113       vold[n][d]     = vn;
 114
 115       if ((ptype[n] != eptDummy) && (ptype[n] != eptShell) && !nFreeze[gf][d]) {
 116         vv             = lg*(vn + f[n][d]*w_dt);
 117
 118         /* do not scale the mean velocities u */
 119         uold           = gstat[ga].uold[d];
 120         va             = vv + accel[ga][d]*dt;
 121         vb             = va + (1.0-lg)*uold;
 122         v[n][d]        = vb;
 123         xprime[n][d]   = x[n][d]+vb*dt;
 124       } else
 125         xprime[n][d]   = x[n][d];
 126     }
 127   }
 128 }
 129
 130 static void do_update_lang(int start,int homenr,double dt,
 131                            ivec nFreeze[],unsigned short ptype[],unsigned short cFREEZE[],
 132                            rvec x[],rvec xprime[],rvec v[],rvec vold[],
 133                            rvec f[],real temp, real fr, int *seed)
 134 {
 135   const unsigned long im = 0xffff;
 136   const unsigned long ia = 1093;
 137   const unsigned long ic = 18257;
 138   int    gf;
 139   real   vn,vv;
 140   real   rfac,invfr,rhalf,jr;
 141   int    n,d;
 142   unsigned long  jran;
 143
 144   /* (r-0.5) n times:  var_n = n * var_1 = n/12
 145      n=4:  var_n = 1/3, so multiply with 3 */
 146
 147   rfac  = sqrt(3.0 * 2.0*BOLTZ*temp/(fr*dt));
 148   rhalf = 2.0*rfac;
 149   rfac  = rfac/(real)im;
 150   invfr = 1.0/fr;
 151
 152   jran = (unsigned long)((real)im*rando(seed));
 153
 154   for (n=start; (n<start+homenr); n++) {
 155     gf   = cFREEZE[n];
 156     for (d=0; (d<DIM); d++) {
 157       vn             = v[n][d];
 158       vold[n][d]     = vn;
 159       if ((ptype[n]!=eptDummy) && (ptype[n]!=eptShell) && !nFreeze[gf][d]) {
 160         jran = (jran*ia+ic) & im;
 161         jr = (real)jran;
 162         jran = (jran*ia+ic) & im;
 163         jr += (real)jran;
 164         jran = (jran*ia+ic) & im;
 165         jr += (real)jran;
 166         jran = (jran*ia+ic) & im;
 167         jr += (real)jran;
 168         vv             = invfr*f[n][d] + rfac * jr - rhalf;
 169         v[n][d]        = vv;
 170         xprime[n][d]   = x[n][d]+v[n][d]*dt;
 171       } else
 172         xprime[n][d]   = x[n][d];
 173     }
 174   }
 175 }
 176
 177 static void shake_calc_vir(FILE *log,int nxf,rvec x[],rvec f[],tensor vir,
 178                            t_commrec *cr)
 179 {
 180   int    i,m,n;
 181   matrix dvir;
 182
 183   clear_mat(dvir);
 184   for(i=0; (i<nxf); i++) {
 185     for(m=0; (m<DIM); m++)
 186       for(n=0; (n<DIM); n++)
 187         dvir[m][n]+=x[i][m]*f[i][n];
 188   }
 189
 190   for(m=0; (m<DIM); m++)
 191     for(n=0; (n<DIM); n++)
 192       vir[m][n]-=0.5*dvir[m][n];
 193 }
 194
 195 static void dump_it_all(FILE *fp,char *title,
 196                  int natoms,rvec x[],rvec xp[],rvec v[],
 197                  rvec vold[],rvec f[])
 198 {
 199 #ifdef DEBUG
 200   fprintf(fp,"%s\n",title);
 201   pr_rvecs(fp,0,"x",x,natoms);
 202   pr_rvecs(fp,0,"xp",xp,natoms);
 203   pr_rvecs(fp,0,"v",v,natoms);
 204   pr_rvecs(fp,0,"vold",vold,natoms);
 205   pr_rvecs(fp,0,"f",f,natoms);
 206 #endif
 207 }
 208
 209 void calc_ke_part(bool bFirstStep,int start,int homenr,
 210                   rvec vold[],rvec v[],rvec vt[],
 211                   t_grpopts *opts,t_mdatoms *md,t_groups *grps,
 212                   t_nrnb *nrnb,real lambda,real *dvdlambda)
 213 {
 214   int          g,d,n,ga,gt;
 215   rvec         v_corrt;
 216   real         hm,vvt,vct;
 217   t_grp_tcstat *tcstat=grps->tcstat;
 218   t_grp_acc    *grpstat=grps->grpstat;
 219   real         dvdl;
 220
 221   /* group velocities are calculated in update_grps and
 222    * accumulated in acumulate_groups.
 223    * Now the partial global and groups ekin.
 224    */
 225   for (g=0; (g<opts->ngtc); g++)
 226     clear_mat(grps->tcstat[g].ekin);
 227
 228   if (bFirstStep) {
 229     for(n=start; (n<start+homenr); n++) {
 230       copy_rvec(v[n],vold[n]);
 231     }
 232     for (g=0; (g<opts->ngacc); g++) {
 233       for(d=0; (d<DIM); d++)
 234         grps->grpstat[g].ut[d]=grps->grpstat[g].u[d];
 235     }
 236   }
 237   else {
 238     for (g=0; (g<opts->ngacc); g++) {
 239       for(d=0; (d<DIM); d++)
 240         grps->grpstat[g].ut[d]=0.5*(grps->grpstat[g].u[d]+
 241                                     grps->grpstat[g].uold[d]);
 242     }
 243   }
 244
 245   dvdl = 0;
 246   for(n=start; (n<start+homenr); n++) {
 247     ga   = md->cACC[n];
 248     gt   = md->cTC[n];
 249     hm   = 0.5*md->massT[n];
 250
 251     for(d=0; (d<DIM); d++) {
 252       vvt        = 0.5*(v[n][d]+vold[n][d]);
 253       vt[n][d]   = vvt;
 254       vct        = vvt - grpstat[ga].ut[d];
 255       v_corrt[d] = vct;
 256     }
 257     for(d=0; (d<DIM); d++) {
 258       tcstat[gt].ekin[XX][d]+=hm*v_corrt[XX]*v_corrt[d];
 259       tcstat[gt].ekin[YY][d]+=hm*v_corrt[YY]*v_corrt[d];
 260       tcstat[gt].ekin[ZZ][d]+=hm*v_corrt[ZZ]*v_corrt[d];
 261     }
 262     if (md->bPerturbed[n]) {
 263       dvdl-=0.5*(md->massB[n]-md->massA[n])*iprod(v_corrt,v_corrt);
 264     }
 265   }
 266   *dvdlambda += dvdl;
 267
 268 #ifdef DEBUG
 269   fprintf(stdlog,"ekin: U=(%12e,%12e,%12e)\n",
 270           grpstat[0].ut[XX],grpstat[0].ut[YY],grpstat[0].ut[ZZ]);
 271   fprintf(stdlog,"ekin: %12e\n",trace(tcstat[0].ekin));
 272 #endif
 273
 274   inc_nrnb(nrnb,eNR_EKIN,homenr);
 275 }
 276
 277 void update(int          natoms,        /* number of atoms in simulation */
 278             int          start,
 279             int          homenr,        /* number of home particles     */
 280             int          step,
 281             real         lambda,
 282             real         *dvdlambda, /* FEP stuff */
 283             t_inputrec   *ir,           /* input record with constants  */
 284             t_mdatoms    *md,
 285             rvec         x[],   /* coordinates of home particles */
 286             t_graph      *graph,
 287             rvec         shift_vec[],
 288             rvec         force[],       /* forces on home particles     */
 289             rvec         delta_f[],
 290             rvec         vold[],        /* Old velocities                  */
 291             rvec         v[],           /* velocities of home particles */
 292             rvec         vt[],          /* velocity at time t           */
 293             tensor       pressure,      /* instantaneous pressure tensor */
 294             tensor       box,           /* instantaneous box lengths    */
 295             t_topology   *top,
 296             t_groups     *grps,
 297             tensor       vir_part,
 298             t_commrec    *cr,
 299             t_nrnb       *nrnb,
 300             bool         bTYZ,
 301             bool         bDoUpdate,
 302             t_edsamyn    *edyn,
 303             t_pull       *pulldata)
 304 {
 305   static char      buf[256];
 306   static bool      bFirst=TRUE;
 307   static rvec      *xprime,*x_unc=NULL;
 308   static int       ngtc,ngacc;
 309   static rvec      *lamb;
 310   static t_edpar   edpar;
 311   static bool      bConstraints;
 312
 313   t_idef           *idef=&(top->idef);
 314   double           dt;
 315   real             dt_1,dt_2;
 316   int              i,n,m,g;
 317
 318   if (bFirst) {
 319     bConstraints = init_constraints(stdlog,top,ir,md,
 320                                     start,homenr);
 321     bConstraints = bConstraints || pulldata->bPull;
 322
 323     if (edyn->bEdsam)
 324       init_edsam(stdlog,top,md,start,homenr,x,box,
 325                  edyn,&edpar);
 326
 327     /* Allocate memory for xold, original atomic positions
 328      * and for xprime.
 329      */
 330     snew(xprime,natoms);
 331     snew(x_unc,homenr);
 332
 333     /* Copy the pointer to the external acceleration in the opts */
 334     ngacc=ir->opts.ngacc;
 335
 336     ngtc=ir->opts.ngtc;
 337     snew(lamb,ir->opts.ngtc);
 338
 339     /* done with initializing */
 340     bFirst=FALSE;
 341   }
 342
 343   dt   = ir->delta_t;
 344   dt_1 = 1.0/dt;
 345   dt_2 = 1.0/(dt*dt);
 346
 347   for(i=0; (i<ngtc); i++) {
 348     real l=grps->tcstat[i].lambda;
 349
 350     if (bTYZ)
 351       lamb[i][XX]=1;
 352     else
 353       lamb[i][XX]=l;
 354     lamb[i][YY]=l;
 355     lamb[i][ZZ]=l;
 356   }
 357
 358   if (bDoUpdate) {
 359     /* update mean velocities */
 360     for (g=0; (g<ngacc); g++) {
 361       copy_rvec(grps->grpstat[g].u,grps->grpstat[g].uold);
 362       clear_rvec(grps->grpstat[g].u);
 363     }
 364
 365     /* Now do the actual update of velocities and positions */
 366     where();
 367     dump_it_all(stdlog,"Before update",natoms,x,xprime,v,vold,force);
 368     if (ir->eI==eiMD)
 369       /* use normal version of update */
 370       do_update(start,homenr,dt,
 371                 lamb,grps->grpstat,
 372                 ir->opts.acc,ir->opts.nFreeze,
 373                 md->invmass,md->ptype,
 374                 md->cFREEZE,md->cACC,md->cTC,
 375                 x,xprime,v,vold,force);
 376     else if (ir->eI==eiLD)
 377       do_update_lang(start,homenr,dt,
 378                      ir->opts.nFreeze,md->ptype,md->cFREEZE,
 379                      x,xprime,v,vold,force,
 380                      ir->ld_temp,ir->ld_fric,&ir->ld_seed);
 381     else
 382       fatal_error(0,"Don't know how to update coordinates");
 383
 384     where();
 385     inc_nrnb(nrnb,eNR_UPDATE,homenr);
 386     dump_it_all(stdlog,"After update",natoms,x,xprime,v,vold,force);
 387   }
 388   else {
 389     /* If we're not updating we're doing shakefirst!
 390      * In this case the extra coordinates are passed in v array
 391      */
 392     for(n=start; (n<start+homenr); n++) {
 393       copy_rvec(v[n],xprime[n]);
 394     }
 395   }
 396
 397   /*
 398    *  Steps (7C, 8C)
 399    *  APPLY CONSTRAINTS:
 400    *  BLOCK SHAKE
 401    */
 402
 403   if (bConstraints) {
 404     /* Copy Unconstrained X to temp array */
 405     for(n=start; (n<start+homenr); n++)
 406       copy_rvec(xprime[n],x_unc[n-start]);
 407
 408     /* Constrain the coordinates xprime */
 409     constrain(stdlog,top,ir,step,md,start,homenr,x,xprime,box,
 410               lambda,dvdlambda,nrnb);
 411
 412     where();
 413
 414     dump_it_all(stdlog,"After Shake",natoms,x,xprime,v,vold,force);
 415
 416     /* apply Essential Dynamics constraints when required */
 417     if (edyn->bEdsam)
 418       do_edsam(stdlog,top,ir,step,md,start,homenr,xprime,x,
 419                x_unc,force,box,edyn,&edpar,bDoUpdate);
 420
 421     /* apply pull constraints when required. Act on xprime, the SHAKED
 422        coordinates.  Don't do anything to f */
 423     if (pulldata->bPull && pulldata->runtype != eAfm &&
 424         pulldata->runtype != eUmbrella &&
 425         pulldata->runtype != eTest)
 426       pull(pulldata,xprime,force,box,top,dt,step,homenr,md);
 427
 428     where();
 429
 430     if (bDoUpdate) {
 431       real mdt_2;
 432
 433       for(n=start; (n<start+homenr); n++) {
 434         mdt_2 = dt_2*md->massT[n];
 435         do_both(x[n],x_unc[n-start],xprime[n],delta_f[n],
 436                 v[n],mdt_2,dt_1);
 437       }
 438       where();
 439
 440       inc_nrnb(nrnb,eNR_SHAKE_V,homenr);
 441       dump_it_all(stdlog,"After Shake-V",natoms,x,xprime,v,vold,force);
 442       where();
 443
 444       /* Calculate virial due to shake (for this proc) */
 445       calc_vir(stdlog,homenr,&(x[start]),&(delta_f[start]),vir_part,cr);
 446       inc_nrnb(nrnb,eNR_SHAKE_VIR,homenr);
 447       where();
 448     }
 449   }
 450
 451
 452   /* We must always unshift here, also if we did not shake
 453    * x was shifted in do_force */
 454   where();
 455   if ((graph->nnodes > 0) && bDoUpdate && (ir->eBox != ebtNONE)) {
 456     unshift_x(graph,shift_vec,x,xprime);
 457     inc_nrnb(nrnb,eNR_SHIFTX,graph->nnodes);
 458     for(n=start; (n<graph->start); n++)
 459       copy_rvec(xprime[n],x[n]);
 460     for(n=graph->start+graph->nnodes; (n<start+homenr); n++)
 461       copy_rvec(xprime[n],x[n]);
 462   }
 463   else {
 464     for(n=start; (n<start+homenr); n++)
 465       copy_rvec(xprime[n],x[n]);
 466   }
 467   dump_it_all(stdlog,"After unshift",natoms,x,xprime,v,vold,force);
 468   where();
 469
 470   if (bDoUpdate) {
 471     update_grps(start,homenr,grps,&(ir->opts),v,md);
 472     do_pcoupl(ir,step,pressure,box,start,homenr,x,md->cFREEZE,nrnb,
 473               ir->opts.nFreeze);
 474     where();
 475   }
 476 }