src/kernel/runner.c

   1 /*
   2  *
   3  *       This source code is part of
   4  *
   5  *        G   R   O   M   A   C   S
   6  *
   7  * GROningen MAchine for Chemical Simulations
   8  *
   9  *            VERSION 2.0
  10  *
  11  * Copyright (c) 1991-1997
  12  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  13  * University of Groningen, The Netherlands
  14  *
  15  * Please refer to:
  16  * GROMACS: A message-passing parallel molecular dynamics implementation
  17  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  18  * Comp. Phys. Comm. 91, 43-56 (1995)
  19  *
  20  * Also check out our WWW page:
  21  * http://rugmd0.chem.rug.nl/~gmx
  22  * or e-mail to:
  23  * gromacs@chem.rug.nl
  24  *
  25  * And Hey:
  26  * GRoups of Organic Molecules in ACtion for Science
  27  */
  28
  29 #include <string.h>
  30 #include <time.h>
  31 #include "typedefs.h"
  32 #include "sysstuff.h"
  33 #include "string2.h"
  34 #include "network.h"
  35 #include "confio.h"
  36 #include "copyrite.h"
  37 #include "smalloc.h"
  38 #include "main.h"
  39 #include "pbc.h"
  40 #include "force.h"
  41 #include "macros.h"
  42 #include "names.h"
  43 #include "fatal.h"
  44 #include "txtdump.h"
  45 #include "update.h"
  46 #include "random.h"
  47 #include "vec.h"
  48 #include "statutil.h"
  49 #include "tgroup.h"
  50 #include "nrnb.h"
  51 #include "disre.h"
  52 #include "mdatoms.h"
  53 #include "mdrun.h"
  54 #include "callf77.h"
  55 #include "pppm.h"
  56
  57 bool optRerunMDset (int nfile, t_filenm fnm[])
  58 {
  59   return opt2bSet("-rerun",nfile,fnm);
  60 }
  61
  62 void get_cmparm(t_inputrec *ir,int step,bool *bStopCM,bool *bStopRot)
  63 {
  64   if (ir->nstcomm == 0) {
  65     *bStopCM  = FALSE;
  66     *bStopRot = FALSE;
  67   } else if (ir->nstcomm > 0) {
  68     *bStopCM  = do_per_step(step,ir->nstcomm);
  69     *bStopRot = FALSE;
  70   } else {
  71     *bStopCM  = FALSE;
  72     *bStopRot = do_per_step(step,-ir->nstcomm);
  73   }
  74 }
  75
  76 void set_pot_bools(t_inputrec *ir,t_topology *top,
  77                    bool *bLR,bool *bLJLR,bool *bBHAM,bool *b14)
  78 {
  79   *bLR   = (ir->rcoulomb > ir->rlist) || EEL_LR(ir->coulombtype);
  80   *bLJLR = (ir->rvdw > ir->rlist);
  81   *bBHAM = (top->idef.functype[0]==F_BHAM);
  82   *b14   = (top->idef.il[F_LJ14].nr > 0);
  83 }
  84
  85 void finish_run(FILE *log,t_commrec *cr,
  86                 const char *confout,
  87                 t_nsborder *nsb,
  88                 t_topology *top,
  89                 t_parm *parm,
  90                 t_nrnb nrnb[],
  91                 double cputime,double realtime,int step,
  92                 bool bWriteStat)
  93 {
  94   int    i,j;
  95   t_nrnb ntot;
  96
  97   for(i=0; (i<eNRNB); i++)
  98     ntot.n[i]=0;
  99   for(i=0; (i<nsb->nprocs); i++)
 100     for(j=0; (j<eNRNB); j++)
 101       ntot.n[j]+=nrnb[i].n[j];
 102
 103   if (bWriteStat) {
 104     if (MASTER(cr)) {
 105       fprintf(stderr,"\n\n");
 106       print_perf(stderr,cputime,realtime,&ntot,nsb->nprocs);
 107     }
 108     else
 109       print_nrnb(log,&(nrnb[nsb->pid]));
 110   }
 111
 112   if (MASTER(cr)) {
 113     print_perf(log,cputime,realtime,&ntot,nsb->nprocs);
 114     if (nsb->nprocs > 1)
 115       pr_load(log,nsb->nprocs,nrnb);
 116   }
 117 }
 118
 119 void mdrunner(t_commrec *cr,int nfile,t_filenm fnm[],bool bVerbose,
 120               bool bCompact,int nDlb,bool bNM,int nstepout,t_edsamyn *edyn)
 121 {
 122   double     cputime=0,realtime;
 123   t_parm     *parm;
 124   rvec       *buf,*f,*vold,*v,*vt,*x,box_size;
 125   real       *ener;
 126   t_nrnb     *nrnb;
 127   t_nsborder *nsb;
 128   t_topology *top;
 129   t_groups   *grps;
 130   t_graph    *graph;
 131   t_mdatoms  *mdatoms;
 132   t_forcerec *fr;
 133   time_t     start_t=0;
 134   bool       bDummies;
 135   int        i,m;
 136
 137   /* Initiate everything (snew sets to zero!) */
 138   snew(ener,F_NRE);
 139   snew(nsb,1);
 140   snew(top,1);
 141   snew(grps,1);
 142   snew(parm,1);
 143   snew(nrnb,cr->nprocs);
 144
 145   /* Initiate invsqrt routines */
 146   init_lookup_table(stdlog);
 147
 148   if (bVerbose && MASTER(cr))
 149     fprintf(stderr,"Getting Loaded...\n");
 150
 151   if (PAR(cr)) {
 152     /* The master processor reads from disk, then passes everything
 153      * around the ring, and finally frees the stuff
 154      */
 155     if (MASTER(cr))
 156       distribute_parts(cr->left,cr->right,cr->pid,cr->nprocs,parm,
 157                        ftp2fn(efTPX,nfile,fnm),nDlb);
 158
 159     /* Every processor (including the master) reads the data from the ring */
 160     init_parts(stdlog,cr,
 161                parm,top,&x,&v,&mdatoms,nsb,
 162                MASTER(cr) ? LIST_SCALARS | LIST_PARM : 0);
 163   } else {
 164     /* Read it up... */
 165     init_single(stdlog,parm,ftp2fn(efTPX,nfile,fnm),top,&x,&v,&mdatoms,nsb);
 166   }
 167   snew(buf,nsb->natoms);
 168   snew(f,nsb->natoms);
 169   snew(vt,nsb->natoms);
 170   snew(vold,nsb->natoms);
 171
 172   if (bVerbose && MASTER(cr))
 173     fprintf(stderr,"Loaded with Money\n\n");
 174
 175   /* Index numbers for parallellism... */
 176   nsb->pid      = cr->pid;
 177   top->idef.pid = cr->pid;
 178
 179   /* Group stuff (energies etc) */
 180   init_groups(stdlog,mdatoms,&(parm->ir.opts),grps);
 181
 182   /* Periodicity stuff */
 183   graph=mk_graph(&(top->idef),top->atoms.nr,FALSE);
 184   if (debug)
 185     p_graph(debug,"Initial graph",graph);
 186
 187   /* Distance Restraints */
 188   init_disres(stdlog,top->idef.il[F_DISRES].nr,&(parm->ir));
 189
 190   /* check if there are dummies */
 191   bDummies=FALSE;
 192   for(i=0; (i<F_NRE) && !bDummies; i++)
 193     bDummies = ((interaction_function[i].flags & IF_DUMMY) &&
 194                 (top->idef.il[i].nr > 0));
 195
 196   /* Initiate forcerecord */
 197   fr = mk_forcerec();
 198   init_forcerec(stdlog,fr,&(parm->ir),&(top->blocks[ebMOLS]),cr,
 199                 &(top->blocks[ebCGS]),&(top->idef),mdatoms,parm->box,FALSE);
 200   /* Initiate box */
 201   for(m=0; (m<DIM); m++)
 202     box_size[m]=parm->box[m][m];
 203
 204   /* Initiate PPPM if necessary */
 205   if (fr->eeltype == eelPPPM)
 206     init_pppm(stdlog,cr,nsb,FALSE,TRUE,box_size,ftp2fn(efHAT,nfile,fnm),&parm->ir);
 207
 208   /* Now do whatever the user wants us to do (how flexible...) */
 209   if (bNM) {
 210     start_t=do_nm(stdlog,cr,nfile,fnm,
 211                   bVerbose,bCompact,nstepout,parm,grps,
 212                   top,ener,x,vold,v,vt,f,buf,
 213                   mdatoms,nsb,nrnb,graph,edyn,fr,box_size);
 214   }
 215   else {
 216     switch (parm->ir.eI) {
 217     case eiMD:
 218     case eiLD:
 219       start_t=do_md(stdlog,cr,nfile,fnm,
 220                     bVerbose,bCompact,bDummies,nstepout,parm,grps,
 221                     top,ener,x,vold,v,vt,f,buf,
 222                     mdatoms,nsb,nrnb,graph,edyn,fr,box_size);
 223       break;
 224     case eiCG:
 225       start_t=do_cg(stdlog,nfile,fnm,parm,top,grps,nsb,
 226                     x,f,buf,mdatoms,parm->ekin,ener,
 227                     nrnb,bVerbose,bDummies,cr,graph,fr,box_size);
 228       break;
 229     case eiSteep:
 230       start_t=do_steep(stdlog,nfile,fnm,parm,top,grps,nsb,
 231                        x,f,buf,mdatoms,parm->ekin,ener,
 232                        nrnb,bVerbose,bDummies,cr,graph,fr,box_size);
 233       break;
 234     default:
 235       fatal_error(0,"Invalid integrator (%d)...\n",parm->ir.eI);
 236     }
 237
 238     /* Some timing stats */
 239     if (MASTER(cr)) {
 240       print_time(stderr,start_t,parm->ir.nsteps,&parm->ir);
 241       realtime=difftime(time(NULL),start_t);
 242       if ((cputime=cpu_time()) == 0)
 243         cputime=realtime;
 244     }
 245     else
 246       realtime=0;
 247
 248     md2atoms(mdatoms,&(top->atoms),TRUE);
 249
 250     /* Finish up, write some stuff */
 251     {
 252       const char *gro=ftp2fn(efSTO,nfile,fnm);
 253
 254       /* if rerunMD, don't write last frame again */
 255       finish_run(stdlog,cr,gro,nsb,top,parm,
 256                  nrnb,cputime,realtime,parm->ir.nsteps,
 257                  (parm->ir.eI==eiMD) || (parm->ir.eI==eiLD));
 258     }
 259   }
 260
 261   /* Does what it says */
 262   print_date_and_time(stdlog,cr->pid,"Finished mdrun");
 263
 264   if (MASTER(cr)) {
 265     thanx(stderr);
 266   }
 267 }
 268
 269 void init_md(t_commrec *cr,t_inputrec *ir,real *t,real *t0,
 270              real *lambda,real *lam0,real *SAfactor,
 271              t_nrnb *mynrnb,bool *bTYZ,t_topology *top,
 272              int nfile,const t_filenm fnm[],char **traj,char **xtc_traj,
 273              ener_file_t *fp_ene, t_mdebin **mdebin,t_groups *grps,rvec vcm,
 274              tensor force_vir,tensor shake_vir,t_mdatoms *mdatoms)
 275 {
 276   bool bBHAM,b14,bLR,bLJLR;
 277
 278   /* Initial values */
 279   *t = *t0       = ir->init_t;
 280   if (ir->bPert) {
 281     *lambda = *lam0 = ir->init_lambda;
 282   }
 283   else {
 284     *lambda = *lam0   = 0.0;
 285   }
 286   if (ir->bSimAnn) {
 287     *SAfactor = 1.0 - *t0/ir->zero_temp_time;
 288     if (*SAfactor < 0)
 289       *SAfactor = 0;
 290   } else
 291     *SAfactor     = 1.0;
 292
 293   init_nrnb(mynrnb);
 294
 295   /* Check Environment variables & other booleans */
 296   *bTYZ=getenv("TYZ") != NULL;
 297   set_pot_bools(ir,top,&bLR,&bLJLR,&bBHAM,&b14);
 298
 299   /* Filenames */
 300   *traj     = ftp2fn(efTRN,nfile,fnm);
 301   *xtc_traj = ftp2fn(efXTC,nfile,fnm);
 302
 303   if (MASTER(cr)) {
 304     *fp_ene = open_enx(ftp2fn(efENX,nfile,fnm),"w");
 305     *mdebin = init_mdebin(*fp_ene,grps,&(top->atoms),&(top->idef),
 306                           bLR,bLJLR,bBHAM,b14,ir->bPert,ir->epc,
 307                           ir->bDispCorr);
 308   }
 309   else {
 310     *fp_ene = NULL;
 311     *mdebin = NULL;
 312   }
 313
 314   /* Initiate variables */
 315   clear_rvec(vcm);
 316   clear_mat(force_vir);
 317   clear_mat(shake_vir);
 318
 319   /* Set initial values for invmass etc. */
 320   init_mdatoms(mdatoms,*lambda,TRUE);
 321
 322   where();
 323 }
 324
 325 void do_pbc_first(FILE *log,t_parm *parm,rvec box_size,t_forcerec *fr,
 326                   t_graph *graph,rvec x[])
 327 {
 328   fprintf(log,"Removing pbc first time\n");
 329   calc_shifts(parm->box,box_size,fr->shift_vec,FALSE);
 330   mk_mshift(log,graph,parm->box,x);
 331   shift_self(graph,fr->shift_vec,x);
 332   fprintf(log,"Done rmpbc\n");
 333 }
 334