src/contrib/glasmd.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  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_glasmd_c = "$Id$";
  30
  31 #include <stdio.h>
  32 #include <string.h>
  33 #include <time.h>
  34 #include "sysstuff.h"
  35 #include "string2.h"
  36 #include "led.h"
  37 #include "nrnb.h"
  38 #include "network.h"
  39 #include "confio.h"
  40 #include "binio.h"
  41 #include "copyrite.h"
  42 #include "smalloc.h"
  43 #include "main.h"
  44 #include "pbc.h"
  45 #include "force.h"
  46 #include "macros.h"
  47 #include "names.h"
  48 #include "stat.h"
  49 #include "fatal.h"
  50 #include "txtdump.h"
  51 #include "futil.h"
  52 #include "typedefs.h"
  53 #include "update.h"
  54 #include "random.h"
  55 #include "vec.h"
  56 #include "filenm.h"
  57 #include "statutil.h"
  58 #include "tgroup.h"
  59 #include "mdrun.h"
  60 #include "vcm.h"
  61 #include "ebin.h"
  62 #include "mdebin.h"
  63 #include "dummies.h"
  64 #include "mdrun.h"
  65 #include "physics.h"
  66 #include "glaasje.h"
  67
  68
  69 time_t do_md(FILE *log,t_commrec *cr,int nfile,t_filenm fnm[],
  70              bool bVerbose,bool bCompact,int stepout,
  71              t_parm *parm,t_groups *grps,
  72              t_topology *top,real ener[],
  73              rvec x[],rvec vold[],rvec v[],rvec vt[],rvec f[],
  74              rvec buf[],t_mdatoms *mdatoms,
  75              t_nsborder *nsb,t_nrnb nrnb[],
  76              t_graph *graph)
  77 {
  78   t_forcerec *fr;
  79   t_mdebin   *mdebin;
  80   FILE       *ene;
  81   int        step;
  82   time_t     start_t;
  83   real       t,lambda,t0,lam0,SAfactor;
  84   bool       bNS,bStopCM,bTYZ,bLR,bBHAM,b14;
  85   tensor     force_vir,shake_vir;
  86   t_nrnb     mynrnb;
  87   char       *traj,*enerfile;
  88   char       *xtc_traj; /* compressed trajectory filename */
  89   int        i,m;
  90   rvec       vcm,box_size;
  91   rvec       *xx,*vv,*ff;
  92
  93   /* Initial values */
  94   t0           = parm->ir.init_t;
  95   if (parm->ir.bPert) {
  96     lam0         = parm->ir.init_lambda;
  97     lambda       = lam0;
  98   }
  99   else {
 100     lam0   = 0.0;
 101     lambda = 0.0;
 102   }
 103   SAfactor     = 1.0;
 104
 105   /* Check Environment variables */
 106   bTYZ=getenv("TYZ") != NULL;
 107
 108   init_nrnb(&mynrnb);
 109
 110   fr=mk_forcerec();
 111   init_forcerec(log,fr,&(parm->ir),&(top->blocks[ebMOLS]),cr,
 112                 &(top->blocks[ebCGS]),&(top->idef),mdatoms,parm->box,FALSE);
 113   for(m=0; (m<DIM); m++)
 114     box_size[m]=parm->box[m][m];
 115   calc_shifts(parm->box,box_size,fr->shift_vec,FALSE);
 116
 117   fprintf(log,"Removing pbc first time\n");
 118   mk_mshift(log,graph,parm->box,x);
 119   shift_self(graph,fr->shift_vec,x);
 120   fprintf(log,"Done rmpbc\n");
 121
 122   traj=ftp2fn(efTRJ,nfile,fnm);
 123   xtc_traj=ftp2fn(efXTC,nfile,fnm);
 124   where();
 125
 126   if (MASTER(cr))
 127     ene=ftp2FILE(efENE,nfile,fnm,"w");
 128   else
 129     ene=NULL;
 130   where();
 131
 132   bLR=(parm->ir.rlong > parm->ir.rshort);
 133   bBHAM=(top->idef.functype[0]==F_BHAM);
 134   b14=(top->idef.il[F_LJ14].nr > 0);
 135   mdebin=init_mdebin(ene,grps,&(top->atoms),bLR,bBHAM,b14);
 136   where();
 137
 138   clear_rvec(vcm);
 139
 140   /* Set initial values for invmass etc. */
 141   init_mdatoms(mdatoms,lambda,TRUE);
 142   where();
 143
 144   if (parm->ir.bShakeFirst)
 145     do_shakefirst(log,bTYZ,lambda,ener,parm,nsb,mdatoms,x,vold,buf,f,v,
 146                   graph,cr,nrnb,grps,fr,top);
 147   where();
 148
 149   /* Compute initial EKin for all.. */
 150   calc_ke_part(TRUE,0,top->atoms.nr,
 151                vold,v,vt,&(parm->ir.opts),
 152                mdatoms,grps,&mynrnb,
 153                lambda,&ener[F_DVDLKIN]);
 154
 155   /* Calculate Temperature coupling parameters lambda */
 156   ener[F_TEMP]=sum_ekin(&(parm->ir.opts),grps,parm->ekin,bTYZ);
 157   tcoupl(parm->ir.btc,&(parm->ir.opts),grps,parm->ir.delta_t,SAfactor);
 158   where();
 159
 160   /* Write start time and temperature */
 161   start_t=print_date_and_time(log,cr->pid,"Started mdrun");
 162   if (MASTER(cr)) {
 163     fprintf(log,"Initial temperature: %g K\n",ener[F_TEMP]);
 164     fprintf(stderr,"starting mdrun '%s'\n%d steps, %8.1f ps.\n\n",*(top->name),
 165             parm->ir.nsteps,parm->ir.nsteps*parm->ir.delta_t);
 166   }
 167   /***********************************************************
 168    *
 169    *             Loop over MD steps
 170    *
 171    ************************************************************/
 172   for (step=0; (step<parm->ir.nsteps); step++) {
 173     /* Stop Center of Mass motion */
 174     bStopCM=do_per_step(step,parm->ir.nstcomm);
 175
 176     /* Determine whether or not to do Neighbour Searching */
 177     bNS=((parm->ir.nstlist && ((step % parm->ir.nstlist)==0)) || (step==0));
 178
 179     /* Construct dummy particles */
 180     construct_dummies(log,x,&mynrnb,parm->ir.delta_t,v,&top->idef);
 181
 182     /* Set values for invmass etc. */
 183     init_mdatoms(mdatoms,lambda,FALSE);
 184
 185     /*init_forcerec(log,fr,&(parm->ir),&(top->blocks[ebMOLS]),cr,
 186       &(top->blocks[ebCGS]),&(top->idef),mdatoms,parm->box,FALSE);
 187       */
 188     clear_mat(force_vir);
 189     do_force(log,cr,parm,nsb,force_vir,step,&mynrnb,
 190              top,grps,x,v,f,buf,mdatoms,ener,bVerbose && !PAR(cr),
 191              lambda,graph,bNS,FALSE,fr);
 192
 193     do_glas(log,START(nsb),HOMENR(nsb),x,f,fr,mdatoms,top->idef.atnr,
 194             &parm->ir);
 195
 196     /* Now we have the energies and forces corresponding to the
 197      * coordinates at time t. We must output all of this before
 198      * the update.
 199      */
 200     t        = t0   + step*parm->ir.delta_t;
 201     if (parm->ir.bPert)
 202       lambda   = lam0 + step*parm->ir.delta_lambda;
 203     SAfactor = 1.0  - step*parm->ir.delta_t*parm->ir.cooling_rate;
 204
 205     /* Spread the force on dummy particle to the other particles... */
 206     spread_dummy_f(log,x,f,&mynrnb,&top->idef);
 207
 208     if (do_per_step(step,parm->ir.nstxout)) xx=x; else xx=NULL;
 209     if (do_per_step(step,parm->ir.nstvout)) vv=v; else vv=NULL;
 210     if (do_per_step(step,parm->ir.nstfout)) ff=f; else ff=NULL;
 211     write_traj(log,cr,traj,
 212                nsb,step,t,lambda,nrnb,nsb->natoms,xx,vv,ff,parm->box);
 213     where();
 214
 215     if (do_per_step(step,parm->ir.nstxtcout)) {
 216       write_xtc_traj(log,cr,xtc_traj,nsb,
 217                      step,t,nsb->natoms,x,parm->box,parm->ir.xtcprec);
 218       where();
 219     }
 220
 221     where();
 222     clear_mat(shake_vir);
 223     update(nsb->natoms,START(nsb),HOMENR(nsb),step,lambda,&ener[F_DVDL],
 224            &(parm->ir),FALSE,mdatoms,x,graph,
 225            fr->shift_vec,f,buf,vold,v,vt,parm->pres,parm->box,
 226            top,grps,shake_vir,cr,&mynrnb,bTYZ,TRUE);
 227     if (PAR(cr))
 228       accumulate_u(cr,&(parm->ir.opts),grps);
 229
 230     where();
 231     /* Calculate partial Kinetic Energy (for this processor)
 232      * per group!
 233      */
 234     calc_ke_part(FALSE,START(nsb),HOMENR(nsb),
 235                  vold,v,vt,&(parm->ir.opts),
 236                  mdatoms,grps,&mynrnb,
 237                  lambda,&ener[F_DVDL]);
 238     where();
 239     if (bStopCM)
 240       calc_vcm(log,HOMENR(nsb),START(nsb),mdatoms->massT,v,vcm);
 241
 242     if (PAR(cr)) {
 243       global_stat(log,cr,ener,force_vir,shake_vir,
 244                   &(parm->ir.opts),grps,&mynrnb,nrnb,vcm);
 245       if (!bNS)
 246         for(i=0; (i<grps->estat.nn); i++)
 247           grps->estat.ee[egLR][i] /= cr->nprocs;
 248     }
 249     else
 250       cp_nrnb(&(nrnb[0]),&mynrnb);
 251
 252     if (bStopCM) {
 253       check_cm(log,vcm,mdatoms->tmass);
 254       do_stopcm(log,HOMENR(nsb),START(nsb),v,vcm,mdatoms->tmass);
 255       inc_nrnb(&mynrnb,eNR_STOPCM,HOMENR(nsb));
 256     }
 257
 258     /* Add force and shake contribution to the virial */
 259     m_add(force_vir,shake_vir,parm->vir);
 260
 261     /* Sum the potential energy terms from group contributions */
 262     sum_epot(&(parm->ir.opts),grps,ener);
 263
 264     /* Sum the kinetic energies of the groups & calc temp */
 265     ener[F_TEMP]=sum_ekin(&(parm->ir.opts),grps,parm->ekin,bTYZ);
 266     ener[F_EKIN]=trace(parm->ekin);
 267     ener[F_ETOT]=ener[F_EPOT]+ener[F_EKIN];
 268
 269     /* Calculate Temperature coupling parameters lambda */
 270     tcoupl(parm->ir.btc,&(parm->ir.opts),grps,parm->ir.delta_t,SAfactor);
 271
 272     /* Calculate pressure ! */
 273     calc_pres(parm->box,parm->ekin,parm->vir,parm->pres);
 274
 275     /* Calculate long range corrections to pressure and energy */
 276     calc_ljcorr(log,parm->ir.userint1,
 277                 fr,mdatoms->nr,parm->box,parm->pres,ener);
 278
 279     upd_mdebin(mdebin,mdatoms->tmass,step,ener,parm->box,shake_vir,
 280                force_vir,parm->vir,parm->pres,grps);
 281
 282     where();
 283     if ((MASTER(cr) && do_per_step(step,parm->ir.nstprint))) {
 284       print_ebin(ene,log,step,t,lambda,SAfactor,
 285                  eprNORMAL,bCompact,mdebin,grps,&(top->atoms));
 286     }
 287     if (bVerbose)
 288       fflush(log);
 289
 290     if ((step % 10) == 0)
 291       update_time();
 292     if (MASTER(cr) && bVerbose && ((step % stepout)==0))
 293       print_time(stderr,start_t,step,&parm->ir);
 294   }
 295   if (MASTER(cr)) {
 296     if (parm->ir.nstprint > 1)
 297       print_ebin(ene,log,step-1,t,lambda,SAfactor,
 298                  eprNORMAL,bCompact,mdebin,grps,&(top->atoms));
 299
 300     print_ebin(NULL,log,step,t,lambda,SAfactor,
 301                eprAVER,FALSE,mdebin,grps,&(top->atoms));
 302     print_ebin(NULL,log,step,t,lambda,SAfactor,
 303                eprRMS,FALSE,mdebin,grps,&(top->atoms));
 304   }
 305   if (ene)
 306     ffclose(ene);
 307
 308   /* Construct dummy particles, for last output frame */
 309   construct_dummies(log,x,&mynrnb,parm->ir.delta_t,v,&top->idef);
 310
 311   /*free_nslist(log);*/
 312
 313   return start_t;
 314 }
 315
 316 int main(int argc,char *argv[])
 317 {
 318   static char *desc[] = {
 319     "The mdrun program performs Molecular Dynamics simulations.",
 320     "It reads the binary topology (.tpb) file and distributes the",
 321     "topology over processors if needed. The coordinates are passed",
 322     "around, so that computations can begin.",
 323     "First a neighbourlist is made, then the forces are computed.",
 324     "The forces are globally summed, and the velocities and",
 325     "positions are updated. If necessary shake is performed to constrain",
 326     "bond lengths and/or bond angles.",
 327     "Temperature and Pressure can be controlled using weak coupling to a",
 328     "bath.[PAR]",
 329     "A number of environment variables can be set to influence the behaviour",
 330     "of the mdrun program. Most of these are for debugging purposes, but they",
 331     "sometimes come in handy when porting the software to an",
 332     "unsupported platform as well. These environment variables",
 333     "are listed elsewhere in the manual.[PAR]",
 334     "The mdrun produces three output file, plus one log file per processor.",
 335     "The first file is the trajectory, containing coordinates, velocities",
 336     "etc. The second file contains the coordinates and velocities at the end",
 337     "of the run plus the computational box. The third file contains energies.",
 338     "In the log file from processor 0 the energies, temperature, etc. are printed.[PAR]",
 339     "When run on a parallel computer or with PVM on a cluster of workstations",
 340     "the [BB]-np[bb] option must be given to indicate the number of",
 341     "processors. Note that at current PVM does work, but it may launch",
 342     "multiple processes on a single processor, which is not very sensible."
 343   };
 344   char         *lognm=NULL;
 345   t_commrec    *cr;
 346   static t_filenm fnm[] = {
 347     { efTPB, NULL, NULL,      ffREAD },
 348     { efTRJ, "-o", NULL,      ffWRITE },
 349     { efXTC, "-x", NULL,      ffOPTWR },
 350     { efGRO, "-c", "confout", ffWRITE },
 351     { efENE, "-e", "ener",    ffWRITE },
 352     { efLOG, "-g", "md",      ffWRITE },
 353   };
 354 #define NFILE asize(fnm)
 355
 356   /* Command line options ! */
 357   static bool bVerbose=FALSE,bCompact=TRUE;
 358   static int  nprocs=1,nDLB=0,nstepout=10;
 359   static t_pargs pa[] = {
 360     { "-np",      FALSE, etINT, &nprocs,
 361       "Number of processors, must be the same as used for grompp. THIS SHOULD BE THE FIRST ARGUMENT ON THE COMMAND LINE FOR MPI" },
 362     { "-v",       FALSE, etBOOL,&bVerbose, "Verbose mode" },
 363     { "-compact", FALSE, etBOOL,&bCompact,
 364       "Write a compact log file, i.e. do not write full virial and energy group matrix (these are also in the energy file, so this is redundant) " },
 365     { "-dlb",     FALSE, etINT, &nDLB,
 366       "Use dynamic load balancing every ... step. BUGGY do not use" },
 367     { "-stepout", FALSE, etINT, &nstepout,
 368       "Frequency of writing the remaining runtime" }
 369   };
 370
 371   get_pargs(argc,argv,asize(pa),pa);
 372   cr = init_par(nprocs,argv);
 373   bVerbose = bVerbose && MASTER(cr);
 374
 375   if (MASTER(cr)) {
 376     CopyRight(stderr,argv[0]);
 377     parse_common_args(&argc,argv,PCA_KEEP_ARGS,TRUE,
 378                       NFILE,fnm,0,NULL,asize(desc),desc,0,NULL);
 379     print_pargs(stderr,asize(pa),pa);
 380   }
 381   open_log(ftp2fn(efLOG,NFILE,fnm),cr);
 382
 383   if (MASTER(cr)) {
 384     CopyRight(stdlog,argv[0]);
 385     please_cite(stdlog,eCITEGMX);
 386   }
 387
 388   mdrunner(cr,NFILE,fnm,bVerbose,bCompact,nDLB,FALSE,nstepout);
 389
 390   exit(0);
 391
 392   return 0;
 393 }
 394