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