src/mdlib/sim_util.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_sim_util_c = "$Id$";
  30
  31 #include <stdio.h>
  32 #include <time.h>
  33 #include "typedefs.h"
  34 #include "string2.h"
  35 #include "smalloc.h"
  36 #include "names.h"
  37 #include "confio.h"
  38 #include "mvdata.h"
  39 #include "txtdump.h"
  40 #include "pbc.h"
  41 #include "vec.h"
  42 #include "time.h"
  43 #include "nrnb.h"
  44 #include "mshift.h"
  45 #include "mdrun.h"
  46 #include "update.h"
  47 #include "physics.h"
  48 #include "main.h"
  49 #include "network.h"
  50
  51 #define difftime(end,start) ((double)(end)-(double)(start))
  52
  53 void print_time(FILE *out,time_t start,int step,t_inputrec *ir)
  54 {
  55   static real time_per_step;
  56   static time_t end;
  57   time_t finish;
  58   double dt;
  59   char buf[48];
  60
  61   fprintf(out,"\rstep %d",step);
  62   if ((step >= ir->nstlist)) {
  63     if ((ir->nstlist == 0) || ((step % ir->nstlist) == 0)) {
  64       /* We have done a full cycle let's update time_per_step */
  65       end=time(NULL);
  66       dt=difftime(end,start);
  67       time_per_step=dt/(step+1);
  68     }
  69     dt=(ir->nsteps-step)*time_per_step;
  70
  71     if (dt >= 300) {
  72       finish = end+(time_t)dt;
  73       sprintf(buf,"%s",ctime(&finish));
  74       buf[strlen(buf)-1]='\0';
  75       fprintf(out,", will finish at %s",buf);
  76     }
  77     else
  78       fprintf(out,", remaining runtime: %5d s               ",(int)dt);
  79   }
  80 #ifdef USE_MPI
  81   if (gmx_parallel)
  82     fprintf(out,"\n");
  83 #endif
  84   fflush(out);
  85 }
  86
  87 time_t print_date_and_time(FILE *log,int pid,char *title)
  88 {
  89   int i;
  90   char *ts,time_string[STRLEN];
  91   time_t now;
  92
  93   now=time(NULL);
  94   ts=ctime(&now);
  95   for (i=0; ts[i]>=' '; i++) time_string[i]=ts[i];
  96   time_string[i]='\0';
  97   fprintf(log,"%s on processor %d %s\n",title,pid,time_string);
  98   return now;
  99 }
 100
 101 static void pr_commrec(FILE *log,t_commrec *cr)
 102 {
 103   fprintf(log,"commrec: pid=%d, nprocs=%d, left=%d, right=%d\n",
 104           cr->pid,cr->nprocs,cr->left,cr->right);
 105 }
 106
 107 static void sum_forces(int start,int end,rvec f[],rvec flr[])
 108 {
 109   int i;
 110
 111   for(i=start; (i<end); i++)
 112     rvec_inc(f[i],flr[i]);
 113 }
 114
 115 static void reset_energies(t_grpopts *opts,t_groups *grp,
 116                            t_forcerec *fr,bool bNS,real epot[])
 117 {
 118   int   i,j;
 119
 120   /* First reset all energy components but the Long Range, except in
 121    * some special cases.
 122    */
 123   for(i=0; (i<egNR); i++)
 124     if (((i != egLR) && (i != egLJLR)) ||
 125         (fr->bTwinRange && bNS) || (!fr->bTwinRange))
 126       for(j=0; (j<grp->estat.nn); j++)
 127         grp->estat.ee[i][j]=0.0;
 128
 129   /* Normal potential energy components */
 130   for(i=0; (i<=F_EPOT); i++)
 131     epot[i] = 0.0;
 132   epot[F_DVDL]    = 0.0;
 133   epot[F_DVDLKIN] = 0.0;
 134 }
 135
 136 /*
 137  * calc_f_el calculates forces due to an electric field.
 138  *
 139  * force is kJ mol^-1 nm^-1 = e * kJ mol^-1 nm^-1 / e
 140  *
 141  * Et[] contains the parameters for the time dependent
 142  * part of the field (not yet used).
 143  * Ex[] contains the parameters for
 144  * the spatial dependent part of the field. You can have cool periodic
 145  * fields in principle, but only a constant field is supported
 146  * now.
 147  * The function should return the energy due to the electric field
 148  * (if any) but for now returns 0.
 149  */
 150
 151 static real calc_f_el(int start,int homenr,real charge[],rvec f[],t_cosines Ex[])
 152 {
 153   real Emu,fmu,strength;
 154   int  i,m;
 155
 156   Emu = 0;
 157   for(m=0; (m<DIM); m++)
 158     if (Ex[m].n) {
 159       /* Convert the field strength from V/nm to MD-units */
 160       strength = Ex[m].a[0]*FIELDFAC;
 161       for(i=start; (i<start+homenr); i++) {
 162         fmu      = charge[i]*strength;
 163         f[i][m] += fmu;
 164       }
 165     }
 166
 167   return Emu;
 168 }
 169
 170 void do_force(FILE *log,t_commrec *cr,
 171               t_parm *parm,t_nsborder *nsb,tensor vir_part,
 172               int step,t_nrnb *nrnb,t_topology *top,t_groups *grps,
 173               rvec x[],rvec v[],rvec f[],rvec buf[],
 174               t_mdatoms *mdatoms,real ener[],bool bVerbose,
 175               real lambda,t_graph *graph,
 176               bool bNS,bool bNBFonly,t_forcerec *fr)
 177 {
 178   static rvec box_size;
 179   static real dvdl_lr = 0;
 180   int    pid,cg0,cg1,i,j;
 181   int    start,homenr;
 182   matrix lr_vir; /* used for PME long range virial */
 183
 184   pid    = cr->pid;
 185   start  = START(nsb);
 186   homenr = HOMENR(nsb);
 187   cg0    = CG0(nsb);
 188   cg1    = CG1(nsb);
 189
 190   update_forcerec(log,fr,parm->box);
 191
 192   if (fr->eBox != ebtNONE) {
 193     /* Compute shift vectors every step, because of pressure coupling! */
 194     if (parm->ir.epc != epcNO)
 195       calc_shifts(parm->box,box_size,fr->shift_vec,FALSE);
 196
 197     if (bNS) {
 198       put_charge_groups_in_box(log,cg0,cg1,FALSE,
 199                                parm->box,box_size,&(top->blocks[ebCGS]),x,
 200                                fr->shift_vec,fr->cg_cm);
 201       inc_nrnb(nrnb,eNR_RESETX,homenr);
 202     }
 203   }
 204   else if (bNS)
 205     calc_cgcm(log,cg0,cg1,&(top->blocks[ebCGS]),x,fr->cg_cm);
 206
 207   if (bNS) {
 208     inc_nrnb(nrnb,eNR_CGCM,cg1-cg0);
 209     if (PAR(cr))
 210       move_cgcm(log,cr,fr->cg_cm,nsb->workload);
 211     if (debug)
 212       pr_rvecs(debug,0,"cgcm",fr->cg_cm,nsb->cgtotal);
 213   }
 214
 215   /* Communicate coordinates if necessary */
 216   if (PAR(cr))
 217     move_x(log,cr->left,cr->right,x,nsb,nrnb);
 218
 219   /* Reset energies */
 220   reset_energies(&(parm->ir.opts),grps,fr,bNS,ener);
 221
 222   if (bNS) {
 223     if (fr->eBox != ebtNONE)
 224       /* Calculate intramolecular shift vectors to make molecules whole */
 225       mk_mshift(log,graph,parm->box,x);
 226
 227     /* Reset long range forces if necessary */
 228     if (fr->bTwinRange) {
 229       clear_rvecs(nsb->natoms,fr->flr);
 230       clear_rvecs(SHIFTS,fr->fshift_lr);
 231     }
 232     /* Do the actual neighbour searching and if twin range electrostatics
 233      * also do the calculation of long range forces and energies.
 234      */
 235     dvdl_lr = 0;
 236     ns(log,fr,x,f,parm->box,grps,&(parm->ir.opts),top,mdatoms,
 237        cr,nrnb,nsb,step,lambda,&dvdl_lr);
 238   }
 239
 240   /* Reset forces or copy them from long range forces */
 241   if (fr->bTwinRange) {
 242     for(i=0; i<nsb->natoms; i++)
 243       copy_rvec(fr->flr[i],f[i]);
 244     for(i=0; i<SHIFTS; i++)
 245       copy_rvec(fr->fshift_lr[i],fr->fshift[i]);
 246   } else {
 247     if (fr->eeltype == eelPPPM ||
 248         fr->eeltype == eelPME ||
 249         fr->eeltype == eelEWALD)
 250       clear_rvecs(homenr,fr->flr+start);
 251     clear_rvecs(nsb->natoms,f);
 252     clear_rvecs(SHIFTS,fr->fshift);
 253   }
 254
 255   /* Compute the forces */
 256   force(log,step,fr,&(parm->ir),&(top->idef),nsb,cr,nrnb,grps,mdatoms,
 257         top->atoms.grps[egcENER].nr,&(parm->ir.opts),
 258         x,f,ener,bVerbose,parm->box,lambda,graph,&(top->atoms.excl),
 259         bNBFonly,lr_vir);
 260
 261   /* Take long range contribution to free energy into account */
 262   ener[F_DVDL] += dvdl_lr;
 263
 264   /* Compute forces due to electric field */
 265   calc_f_el(start,homenr,mdatoms->chargeT,f,parm->ir.ex);
 266
 267 #ifdef DEBUG
 268   if (bNS)
 269     print_nrnb(log,nrnb);
 270 #endif
 271
 272   /* The short-range virial from surrounding boxes */
 273   clear_mat(vir_part);
 274   calc_vir(log,SHIFTS,fr->shift_vec,fr->fshift,vir_part,cr);
 275   inc_nrnb(nrnb,eNR_VIRIAL,SHIFTS);
 276
 277 #ifdef DEBUG
 278   if (debug) {
 279     pr_rvecs(debug,0,"fr->fshift",fr->fshift,SHIFTS);
 280     pr_rvecs(debug,0,"in force.c vir",vir_part,DIM);
 281   }
 282 #endif
 283
 284   /* When using PME/Ewald we compute the long range virial (lr_vir) there.
 285    * otherwise we do it based on long range forces from twin range
 286    * cut-off based calculation (or not at all).
 287    */
 288
 289   /* Communicate the forces */
 290   if (PAR(cr))
 291     move_f(log,cr->left,cr->right,f,buf,nsb,nrnb);
 292
 293   /* Now it is time for the short range virial. At this timepoint vir_part
 294    * already contains the virial from surrounding boxes.
 295    * Calculate partial virial, for local atoms only, based on short range.
 296    * Total virial is computed in global_stat, called from do_md
 297    */
 298   f_calc_vir(log,start,start+homenr,x,f,vir_part,cr,graph,fr->shift_vec);
 299   inc_nrnb(nrnb,eNR_VIRIAL,homenr);
 300
 301   if ((fr->eeltype == eelPPPM) ||
 302       (fr->eeltype == eelPME)  ||
 303       (fr->eeltype == eelEWALD)) {
 304     /* Now add the forces from the PME calculation. Since this only produces
 305      * forces on the local atoms, this can be safely done after the
 306      * communication step. The same goes for the virial.
 307      */
 308     sum_forces(start,start+homenr,f,fr->flr);
 309     if (fr->eeltype != eelPPPM) /* PPPM virial sucks */
 310       for(i=0; (i<DIM); i++)
 311         for(j=0; (j<DIM); j++)
 312           vir_part[i][j]+=lr_vir[i][j];
 313   }
 314
 315   if (debug)
 316     pr_rvecs(debug,0,"vir_part",vir_part,DIM);
 317 }
 318
 319 #ifdef NO_CLOCK
 320 #define clock() -1
 321 #endif
 322 static double runtime=0;
 323 static clock_t cprev;
 324
 325 void start_time(void)
 326 {
 327   cprev   = clock();
 328   runtime = 0.0;
 329 }
 330
 331 void update_time(void)
 332 {
 333   clock_t c;
 334
 335   c        = clock();
 336   runtime += (c-cprev)/(double)CLOCKS_PER_SEC;
 337   cprev    = c;
 338 }
 339
 340 double cpu_time(void)
 341 {
 342   return runtime;
 343 }
 344
 345 void do_shakefirst(FILE *log,bool bTYZ,real lambda,real ener[],
 346                    t_parm *parm,t_nsborder *nsb,t_mdatoms *md,
 347                    rvec x[],rvec vold[],rvec buf[],rvec f[],
 348                    rvec v[],t_graph *graph,t_commrec *cr,t_nrnb *nrnb,
 349                    t_groups *grps,t_forcerec *fr,t_topology *top,
 350                    t_edsamyn *edyn,t_pull *pulldata)
 351 {
 352   int    i,m;
 353   tensor shake_vir;
 354   real   dt=parm->ir.delta_t;
 355   real   dt_1;
 356
 357   if ((top->idef.il[F_SHAKE].nr > 0)   ||
 358       (top->idef.il[F_SETTLE].nr > 0)) {
 359     /* Do a first SHAKE to reset particles... */
 360     clear_mat(shake_vir);
 361     update(nsb->natoms,START(nsb),HOMENR(nsb),
 362            -1,lambda,&ener[F_DVDL],
 363            &(parm->ir),md,x,graph,
 364            fr->shift_vec,NULL,NULL,vold,x,NULL,parm->pres,parm->box,
 365            top,grps,shake_vir,cr,nrnb,bTYZ,FALSE,edyn,pulldata);
 366     /* Compute coordinates at t=-dt, store them in buf */
 367     for(i=0; (i<nsb->natoms); i++) {
 368       for(m=0; (m<DIM); m++) {
 369         f[i][m]=x[i][m];
 370         buf[i][m]=x[i][m]-dt*v[i][m];
 371       }
 372     }
 373
 374     /* Shake the positions at t=-dt with the positions at t=0
 375      * as reference coordinates.
 376      */
 377     clear_mat(shake_vir);
 378     update(nsb->natoms,START(nsb),HOMENR(nsb),
 379            0,lambda,&ener[F_DVDL],&(parm->ir),md,f,graph,
 380            fr->shift_vec,NULL,NULL,vold,buf,NULL,parm->pres,parm->box,
 381            top,grps,shake_vir,cr,nrnb,bTYZ,FALSE,edyn,pulldata);
 382
 383     /* Compute the velocities at t=-dt/2 using the coordinates at
 384      * t=-dt and t=0
 385      */
 386     dt_1=1.0/dt;
 387     for(i=0; (i<nsb->natoms); i++) {
 388       for(m=0; (m<DIM); m++)
 389         v[i][m]=(x[i][m]-f[i][m])*dt_1;
 390     }
 391
 392     /* Shake the positions at t=-dt with the positions at t=0
 393      * as reference coordinates.
 394      */
 395     clear_mat(shake_vir);
 396     update(nsb->natoms,START(nsb),HOMENR(nsb),
 397            0,lambda,&ener[F_DVDL],&(parm->ir),md,f,graph,
 398            fr->shift_vec,NULL,NULL,vold,buf,NULL,parm->pres,parm->box,
 399            top,grps,shake_vir,cr,nrnb,bTYZ,FALSE,edyn,pulldata);
 400
 401     /* Compute the velocities at t=-dt/2 using the coordinates at
 402      * t=-dt and t=0
 403      */
 404     dt_1=1.0/dt;
 405     for(i=0; (i<nsb->natoms); i++) {
 406       for(m=0; (m<DIM); m++)
 407         v[i][m]=(x[i][m]-f[i][m])*dt_1;
 408     }
 409   }
 410 }
 411
 412 void calc_dispcorr(FILE *log,bool bDispCorr,t_forcerec *fr,int natoms,
 413                    matrix box,tensor pres,tensor virial,real ener[])
 414 {
 415   static bool bFirst=TRUE;
 416   real vol,rc3,spres,svir;
 417   int  m;
 418
 419   if (bDispCorr) {
 420     vol           = det(box);
 421     /* Forget the (small) effect of the shift on the LJ energy *
 422      * when fr->bLJShift = TRUE                                */
 423     rc3              = fr->rvdw*fr->rvdw*fr->rvdw;
 424     ener[F_DISPCORR] = -2.0*natoms*natoms*M_PI*fr->avcsix/(3.0*vol*rc3);
 425     spres            = 2.0*ener[F_DISPCORR]*PRESFAC/vol;
 426     svir             = -6.0*ener[F_DISPCORR];
 427     ener[F_PRES]     = trace(pres)/3.0+spres;
 428     for(m=0; (m<DIM); m++) {
 429       pres[m][m]    += spres;
 430       virial[m][m]  += svir;
 431     }
 432     if (bFirst) {
 433       fprintf(log,
 434               "Long Range LJ corrections: Epot=%10g, Pres=%10g, Vir=%10g\n",
 435               ener[F_DISPCORR],spres,svir);
 436       bFirst = FALSE;
 437     }
 438   }
 439   else {
 440     ener[F_DISPCORR] = 0.0;
 441     ener[F_PRES]     = trace(pres)/3.0;
 442   }
 443   ener[F_EPOT] += ener[F_DISPCORR];
 444   ener[F_ETOT] += ener[F_DISPCORR];
 445 }
 446