src/mdlib/nm.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  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_nm_c = "$Id$";
  30
  31 #include "typedefs.h"
  32 #include "vec.h"
  33 #include "smalloc.h"
  34 #include "statutil.h"
  35 #include "vcm.h"
  36 #include "mdebin.h"
  37 #include "nrnb.h"
  38 #include "calcmu.h"
  39 #include "dummies.h"
  40 #include "pppm.h"
  41 #include "update.h"
  42 #include "mdrun.h"
  43
  44 time_t do_nm(FILE *log,t_commrec *cr,int nfile,t_filenm fnm[],
  45              bool bVerbose,bool bCompact,int stepout,
  46              t_parm *parm,t_groups *grps,
  47              t_topology *top,real ener[],
  48              rvec x[],rvec vold[],rvec v[],rvec vt[],rvec f[],
  49              rvec buf[],t_mdatoms *mdatoms,
  50              t_nsborder *nsb,t_nrnb nrnb[],
  51              t_graph *graph,t_edsamyn *edyn,
  52              t_forcerec *fr,rvec box_size)
  53 {
  54   t_mdebin   *mdebin;
  55   int        fp_ene,step,nre;
  56   time_t     start_t;
  57   real       t,lambda,t0,lam0;
  58   bool       bNS,bStopCM,bTYZ,bLR,bLJLR,bBHAM,b14,bBox;
  59   tensor     force_vir,shake_vir;
  60   t_nrnb     mynrnb;
  61   char       *mtx,*enerfile,wfile[80];
  62   int        i,m;
  63   rvec       vcm;
  64   rvec       *xx,*vv,*ff;
  65
  66   /* added with respect to mdrun */
  67
  68   int        idum,jdum,kdum;
  69   real       der_range=1.0e-6,fmax;
  70   rvec       *dfdx;
  71   snew(dfdx,top->atoms.nr);
  72
  73   vv=NULL;
  74   ff=NULL;
  75
  76   /* end nmrun differences */
  77
  78   /* Initial values */
  79   t0           = parm->ir.init_t;
  80   lam0         = parm->ir.init_lambda;
  81   t            = t0;
  82   lambda       = lam0;
  83
  84   /* Check Environment variables */
  85   bTYZ=getenv("TYZ") != NULL;
  86
  87   init_nrnb(&mynrnb);
  88
  89   bBox = (fr->eBox != ebtNONE);
  90   if (bBox) {
  91     calc_shifts(parm->box,box_size,fr->shift_vec,FALSE);
  92     fprintf(log,"Removing pbc first time\n");
  93     mk_mshift(log,graph,parm->box,x);
  94     shift_self(graph,fr->shift_vec,x);
  95   }
  96
  97   mtx=ftp2fn(efMTX,nfile,fnm);
  98
  99   /* Open the enrgy file */
 100   if (MASTER(cr))
 101     fp_ene=open_enx(ftp2fn(efENX,nfile,fnm),"w");
 102   else
 103     fp_ene=-1;
 104   set_pot_bools(&(parm->ir),top,&bLR,&bLJLR,&bBHAM,&b14);
 105   mdebin=init_mdebin(fp_ene,grps,&(top->atoms),&(top->idef),bLR,bLJLR,
 106                      bBHAM,b14,parm->ir.bPert,parm->ir.epc,
 107                      parm->ir.bDispCorr,cr);
 108
 109   /* Compute initial EKin for all.. */
 110   calc_ke_part(TRUE,0,top->atoms.nr,
 111                vold,v,vt,&(parm->ir.opts),
 112                mdatoms,grps,&mynrnb,
 113                lambda,&ener[F_DVDL]);
 114
 115   /* Calculate Temperature coupling parameters lambda */
 116   ener[F_TEMP]=sum_ekin(&(parm->ir.opts),grps,parm->ekin,bTYZ);
 117   tcoupl(parm->ir.btc,&(parm->ir.opts),grps,parm->ir.delta_t,lam0,0,
 118          parm->ir.ntcmemory);
 119   where();
 120
 121   /* Write start time and temperature */
 122   start_t=print_date_and_time(log,cr->pid,"Started nmrun");
 123   if (MASTER(cr)) {
 124     fprintf(stderr,"starting nmrun '%s'\n%d steps.\n\n",*(top->name),
 125             top->atoms.nr);
 126   }
 127
 128   clear_rvec(vcm);
 129
 130   /* Call do_force once to make pairlist */
 131
 132   clear_mat(force_vir);
 133
 134   bNS=TRUE;
 135   do_force(log,cr,parm,nsb,force_vir,0,&mynrnb,
 136            top,grps,x,v,f,buf,mdatoms,ener,bVerbose && !PAR(cr),
 137            lambda,graph,bNS,FALSE,fr);
 138   bNS=FALSE;
 139   if (bBox)
 140     /* Shift back the coordinates, since we're not calling update */
 141     unshift_self(graph,fr->shift_vec,x);
 142
 143
 144   /* if forces not small, warn user */
 145   fmax=f_max(cr->left,cr->right,nsb->nprocs,0,top->atoms.nr,f);
 146   fprintf(stderr,"Maximum force:%12.5e\n",fmax);
 147   if (fmax > 1.0e-3) {
 148     fprintf(stderr,"Maximum force probably not small enough to");
 149     fprintf(stderr," ensure that you are in a \nenergy well. ");
 150     fprintf(stderr,"Be aware that negative eigenvalues may occur");
 151     fprintf(stderr," when the\nresulting matrix is diagonalized.\n");
 152   }
 153
 154   /***********************************************************
 155    *
 156    *      Loop over all pairs in matrix
 157    *
 158    *      do_force called twice. Once with positive and
 159    *      once with negative displacement
 160    *
 161    ************************************************************/
 162
 163   /* fudge nr of steps to nr of atoms */
 164
 165   parm->ir.nsteps=top->atoms.nr;
 166
 167   for (step=0; (step<parm->ir.nsteps); step++) {
 168
 169     for (idum=0; (idum<DIM); idum++) {
 170
 171       x[step][idum]=x[step][idum]-der_range;
 172
 173       clear_mat(force_vir);
 174
 175       do_force(log,cr,parm,nsb,force_vir,2*(step*DIM+idum),&mynrnb,
 176                top,grps,x,v,f,buf,mdatoms,ener,bVerbose && !PAR(cr),
 177                lambda,graph,bNS,FALSE,fr);
 178       if (bBox)
 179         /* Shift back the coordinates, since we're not calling update */
 180         unshift_self(graph,fr->shift_vec,x);
 181
 182       for (jdum=0; (jdum<top->atoms.nr); jdum++) {
 183         for (kdum=0; (kdum<DIM); kdum++) {
 184           dfdx[jdum][kdum]=f[jdum][kdum];
 185         }
 186       }
 187
 188       x[step][idum]=x[step][idum]+2.0*der_range;
 189
 190       clear_mat(force_vir);
 191
 192       do_force(log,cr,parm,nsb,force_vir,2*(step*DIM+idum)+1,&mynrnb,
 193                top,grps,x,v,f,buf,mdatoms,ener,bVerbose && !PAR(cr),
 194                lambda,graph,bNS,FALSE,fr);
 195       if (bBox)
 196         /* Shift back the coordinates, since we're not calling update */
 197         unshift_self(graph,fr->shift_vec,x);
 198
 199       for (jdum=0; (jdum<top->atoms.nr); jdum++) {
 200         for (kdum=0; (kdum<DIM); kdum++) {
 201           dfdx[jdum][kdum]=(f[jdum][kdum]-dfdx[jdum][kdum])/2.0e-6;
 202         }
 203       }
 204
 205       /* store derivatives now, diagonalization later */
 206
 207       xx=dfdx;
 208       write_traj(log,cr,mtx,
 209                  nsb,step,t,lambda,nrnb,nsb->natoms,xx,vv,ff,parm->box);
 210
 211       /* x is restored to original */
 212
 213       x[step][idum]=x[step][idum]-der_range;
 214
 215       if (bVerbose)
 216         fflush(log);
 217
 218       /*if (MASTER(cr) && bVerbose && ((step % stepout)==0))
 219         print_time(stderr,start_t,step,&parm->ir);*/
 220     }
 221     /* write progress */
 222     if (MASTER(cr) && bVerbose) {
 223       fprintf(stderr,"\rFinished step %d out of %d",step+1,top->atoms.nr);
 224       fflush(stderr);
 225     }
 226   }
 227   t=t0+step*parm->ir.delta_t;
 228   lambda=lam0+step*parm->ir.delta_lambda;
 229
 230   if (MASTER(cr)) {
 231     print_ebin(-1,FALSE,FALSE,log,step,t,lambda,0.0,eprAVER,
 232                FALSE,mdebin,grps,&(top->atoms));
 233     print_ebin(-1,FALSE,FALSE,log,step,t,lambda,0.0,eprRMS,
 234                FALSE,mdebin,grps,&(top->atoms));
 235   }
 236
 237   /* Construct dummy particles, for last output frame */
 238   construct_dummies(log,x,&mynrnb,parm->ir.delta_t,v,&top->idef);
 239
 240   /*free_nslist(log);*/
 241
 242   return start_t;
 243 }