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