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