src/mdlib/sim_util.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.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * GROup of MAchos and Cynical Suckers
  31  */
  32 static char *SRCID_sim_util_c = "$Id$";
  33 #include <stdio.h>
  34 #include <time.h>
  35 #include "typedefs.h"
  36 #include "string2.h"
  37 #include "smalloc.h"
  38 #include "names.h"
  39 #include "confio.h"
  40 #include "mvdata.h"
  41 #include "txtdump.h"
  42 #include "pbc.h"
  43 #include "vec.h"
  44 #include "time.h"
  45 #include "nrnb.h"
  46 #include "mshift.h"
  47 #include "mdrun.h"
  48 #include "update.h"
  49 #include "physics.h"
  50 #include "main.h"
  51 #include "mdatoms.h"
  52 #include "pme.h"
  53 #include "pppm.h"
  54 #include "disre.h"
  55 #include "orires.h"
  56 #include "network.h"
  57 #include "calcmu.h"
  58 #include "constr.h"
  59 #include "xvgr.h"
  60
  61 #define difftime(end,start) ((double)(end)-(double)(start))
  62
  63 void print_time(FILE *out,time_t start,int step,t_inputrec *ir)
  64 {
  65   static real time_per_step;
  66   static time_t end;
  67   time_t finish;
  68   double dt;
  69   char buf[48];
  70
  71   fprintf(out,"\rstep %d",step);
  72   if ((step >= ir->nstlist)) {
  73     if ((ir->nstlist == 0) || ((step % ir->nstlist) == 0)) {
  74       /* We have done a full cycle let's update time_per_step */
  75       end=time(NULL);
  76       dt=difftime(end,start);
  77       time_per_step=dt/(step+1);
  78     }
  79     dt=(ir->nsteps-step)*time_per_step;
  80
  81     if (dt >= 300) {
  82       finish = end+(time_t)dt;
  83       sprintf(buf,"%s",ctime(&finish));
  84       buf[strlen(buf)-1]='\0';
  85       fprintf(out,", will finish at %s",buf);
  86     }
  87     else
  88       fprintf(out,", remaining runtime: %5d s               ",(int)dt);
  89   }
  90 #ifdef USE_MPI
  91   if (gmx_parallel)
  92     fprintf(out,"\n");
  93 #endif
  94   fflush(out);
  95 }
  96
  97 time_t print_date_and_time(FILE *log,int nodeid,char *title)
  98 {
  99   int i;
 100   char *ts,time_string[STRLEN];
 101   time_t now;
 102
 103   now=time(NULL);
 104   ts=ctime(&now);
 105   for (i=0; ts[i]>=' '; i++) time_string[i]=ts[i];
 106   time_string[i]='\0';
 107   fprintf(log,"%s on node %d %s\n",title,nodeid,time_string);
 108   return now;
 109 }
 110
 111 static void pr_commrec(FILE *log,t_commrec *cr)
 112 {
 113   fprintf(log,"commrec: nodeid=%d, nnodes=%d, left=%d, right=%d, threadid=%d, nthreads=%d\n",
 114           cr->nodeid,cr->nnodes,cr->left,cr->right,cr->threadid,cr->nthreads);
 115 }
 116
 117 static void sum_forces(int start,int end,rvec f[],rvec flr[])
 118 {
 119   int i;
 120
 121   if (debug) {
 122     pr_rvecs(debug,0,"fsr",f+start,end-start);
 123     pr_rvecs(debug,0,"flr",flr+start,end-start);
 124   }
 125   for(i=start; (i<end); i++)
 126     rvec_inc(f[i],flr[i]);
 127 }
 128
 129 static void reset_energies(t_grpopts *opts,t_groups *grp,
 130                            t_forcerec *fr,bool bNS,real epot[])
 131 {
 132   int   i,j;
 133
 134   /* First reset all energy components but the Long Range, except in
 135    * some special cases.
 136    */
 137   for(i=0; (i<egNR); i++)
 138     if (((i != egLR) && (i != egLJLR)) ||
 139         (fr->bTwinRange && bNS) || (!fr->bTwinRange))
 140       for(j=0; (j<grp->estat.nn); j++)
 141         grp->estat.ee[i][j]=0.0;
 142
 143   /* Normal potential energy components */
 144   for(i=0; (i<=F_EPOT); i++)
 145     epot[i] = 0.0;
 146   epot[F_DVDL]    = 0.0;
 147   epot[F_DVDLKIN] = 0.0;
 148 }
 149
 150 /*
 151  * calc_f_el calculates forces due to an electric field.
 152  *
 153  * force is kJ mol^-1 nm^-1 = e * kJ mol^-1 nm^-1 / e
 154  *
 155  * Et[] contains the parameters for the time dependent
 156  * part of the field (not yet used).
 157  * Ex[] contains the parameters for
 158  * the spatial dependent part of the field. You can have cool periodic
 159  * fields in principle, but only a constant field is supported
 160  * now.
 161  * The function should return the energy due to the electric field
 162  * (if any) but for now returns 0.
 163  */
 164
 165 static real calc_f_el(int start,int homenr,real charge[],rvec f[],t_cosines Ex[])
 166 {
 167   real Emu,fmu,strength;
 168   int  i,m;
 169
 170   Emu = 0;
 171   for(m=0; (m<DIM); m++)
 172     if (Ex[m].n) {
 173       /* Convert the field strength from V/nm to MD-units */
 174       strength = Ex[m].a[0]*FIELDFAC;
 175       for(i=start; (i<start+homenr); i++) {
 176         fmu      = charge[i]*strength;
 177         f[i][m] += fmu;
 178       }
 179     }
 180
 181   return Emu;
 182 }
 183
 184 void do_force(FILE *log,t_commrec *cr,t_commrec *mcr,
 185               t_parm *parm,t_nsborder *nsb,tensor vir_part,tensor pme_vir,
 186               int step,t_nrnb *nrnb,t_topology *top,t_groups *grps,
 187               rvec x[],rvec v[],rvec f[],rvec buf[],
 188               t_mdatoms *mdatoms,real ener[],t_fcdata *fcd,bool bVerbose,
 189               real lambda,t_graph *graph,
 190               bool bNS,bool bNBFonly,t_forcerec *fr, rvec mu_tot,
 191               bool bGatherOnly)
 192 {
 193   static rvec box_size;
 194   static real dvdl_lr = 0;
 195   int    cg0,cg1,i,j;
 196   int    start,homenr;
 197   static real mu_and_q[DIM+1];
 198   real   qsum;
 199
 200   start  = START(nsb);
 201   homenr = HOMENR(nsb);
 202   cg0    = CG0(nsb);
 203   cg1    = CG1(nsb);
 204
 205   update_forcerec(log,fr,parm->box);
 206
 207   /* Calculate total (local) dipole moment in a temporary common array.
 208    * This makes it possible to sum them over nodes faster.
 209    */
 210   calc_mu_and_q(nsb,x,mdatoms->chargeT,mu_and_q,mu_and_q+DIM);
 211
 212   if (fr->ePBC != epbcNONE) {
 213     /* Compute shift vectors every step, because of pressure coupling! */
 214     if (parm->ir.epc != epcNO)
 215       calc_shifts(parm->box,box_size,fr->shift_vec);
 216
 217     if (bNS) {
 218       put_charge_groups_in_box(log,cg0,cg1,parm->box,box_size,
 219                                &(top->blocks[ebCGS]),x,fr->cg_cm);
 220       inc_nrnb(nrnb,eNR_RESETX,homenr);
 221     } else if (parm->ir.eI==eiSteep || parm->ir.eI==eiCG)
 222       unshift_self(graph,parm->box,x);
 223
 224   }
 225   else if (bNS)
 226     calc_cgcm(log,cg0,cg1,&(top->blocks[ebCGS]),x,fr->cg_cm);
 227
 228   if (bNS) {
 229     inc_nrnb(nrnb,eNR_CGCM,cg1-cg0);
 230     if (PAR(cr))
 231       move_cgcm(log,cr,fr->cg_cm,nsb->workload);
 232     if (debug)
 233       pr_rvecs(debug,0,"cgcm",fr->cg_cm,nsb->cgtotal);
 234   }
 235
 236   /* Communicate coordinates and sum dipole and net charge if necessary */
 237   if (PAR(cr)) {
 238     move_x(log,cr->left,cr->right,x,nsb,nrnb);
 239     gmx_sum(DIM+1,mu_and_q,cr);
 240   }
 241   for(i=0;i<DIM;i++)
 242     mu_tot[i]=mu_and_q[i];
 243   qsum=mu_and_q[DIM];
 244
 245   /* Reset energies */
 246   reset_energies(&(parm->ir.opts),grps,fr,bNS,ener);
 247   if (bNS) {
 248     if (fr->ePBC != epbcNONE)
 249       /* Calculate intramolecular shift vectors to make molecules whole */
 250       mk_mshift(log,graph,parm->box,x);
 251
 252     /* Reset long range forces if necessary */
 253     if (fr->bTwinRange) {
 254       clear_rvecs(nsb->natoms,fr->f_twin);
 255       clear_rvecs(SHIFTS,fr->fshift_twin);
 256     }
 257     /* Do the actual neighbour searching and if twin range electrostatics
 258      * also do the calculation of long range forces and energies.
 259      */
 260     dvdl_lr = 0;
 261
 262     ns(log,fr,x,f,parm->box,grps,&(parm->ir.opts),top,mdatoms,
 263        cr,nrnb,nsb,step,lambda,&dvdl_lr);
 264   }
 265   /* Reset PME/Ewald forces if necessary */
 266   if (EEL_LR(fr->eeltype))
 267     clear_rvecs(homenr,fr->f_pme+start);
 268
 269   /* Copy long range forces into normal buffers */
 270   if (fr->bTwinRange) {
 271     for(i=0; i<nsb->natoms; i++)
 272       copy_rvec(fr->f_twin[i],f[i]);
 273     for(i=0; i<SHIFTS; i++)
 274       copy_rvec(fr->fshift_twin[i],fr->fshift[i]);
 275   }
 276   else {
 277     clear_rvecs(nsb->natoms,f);
 278     clear_rvecs(SHIFTS,fr->fshift);
 279   }
 280
 281   /* Compute the forces */
 282   force(log,step,fr,&(parm->ir),&(top->idef),nsb,cr,mcr,nrnb,grps,mdatoms,
 283         top->atoms.grps[egcENER].nr,&(parm->ir.opts),
 284         x,f,ener,fcd,bVerbose,parm->box,lambda,graph,&(top->atoms.excl),
 285         bNBFonly,pme_vir,mu_tot,qsum,bGatherOnly);
 286
 287   /* Take long range contribution to free energy into account */
 288   ener[F_DVDL] += dvdl_lr;
 289
 290 #ifdef DEBUG
 291   if (bNS)
 292     print_nrnb(log,nrnb);
 293 #endif
 294
 295   /* The short-range virial from surrounding boxes */
 296   clear_mat(vir_part);
 297   calc_vir(log,SHIFTS,fr->shift_vec,fr->fshift,vir_part);
 298   inc_nrnb(nrnb,eNR_VIRIAL,SHIFTS);
 299
 300   if (debug)
 301     pr_rvecs(debug,0,"vir_shifts",vir_part,DIM);
 302
 303   /* Compute forces due to electric field */
 304   calc_f_el(start,homenr,mdatoms->chargeT,f,parm->ir.ex);
 305
 306   /* When using PME/Ewald we compute the long range virial (pme_vir) there.
 307    * otherwise we do it based on long range forces from twin range
 308    * cut-off based calculation (or not at all).
 309    */
 310
 311   /* Communicate the forces */
 312   if (PAR(cr))
 313     move_f(log,cr->left,cr->right,f,buf,nsb,nrnb);
 314 }
 315
 316 void sum_lrforces(rvec f[],t_forcerec *fr,int start,int homenr)
 317 {
 318   /* Now add the forces from the PME calculation. Since this only produces
 319    * forces on the local atoms, this can be safely done after the
 320    * communication step.
 321    */
 322   if (EEL_LR(fr->eeltype))
 323     sum_forces(start,start+homenr,f,fr->f_pme);
 324 }
 325
 326 void calc_virial(FILE *log,int start,int homenr,rvec x[],rvec f[],
 327                  tensor vir_part,tensor pme_vir,
 328                  t_graph *graph,matrix box,
 329                  t_nrnb *nrnb,t_forcerec *fr,bool bTweak)
 330 {
 331   int i,j;
 332   tensor virtest;
 333
 334   /* Now it is time for the short range virial. At this timepoint vir_part
 335    * already contains the virial from surrounding boxes.
 336    * Calculate partial virial, for local atoms only, based on short range.
 337    * Total virial is computed in global_stat, called from do_md
 338    */
 339   f_calc_vir(log,start,start+homenr,x,f,vir_part,graph,box);
 340   inc_nrnb(nrnb,eNR_VIRIAL,homenr);
 341
 342   /* Add up the long range forces if necessary */
 343   /* if (!bTweak) {
 344     sum_lrforces(f,fr,start,homenr);
 345     }*/
 346
 347   /* Add up virial if necessary */
 348   if (EEL_LR(fr->eeltype) && (fr->eeltype != eelPPPM)) {
 349     if (debug && bTweak) {
 350       clear_mat(virtest);
 351       f_calc_vir(log,start,start+homenr,x,fr->f_pme,virtest,graph,box);
 352       pr_rvecs(debug,0,"virtest",virtest,DIM);
 353       pr_rvecs(debug,0,"pme_vir",pme_vir,DIM);
 354     }
 355     /* PPPM virial sucks */
 356     if (!bTweak)
 357       for(i=0; (i<DIM); i++)
 358         for(j=0; (j<DIM); j++)
 359           vir_part[i][j]+=pme_vir[i][j];
 360   }
 361   if (debug)
 362     pr_rvecs(debug,0,"vir_part",vir_part,DIM);
 363 }
 364
 365 #ifdef NO_CLOCK
 366 #define clock() -1
 367 #endif
 368 static double runtime=0;
 369 static clock_t cprev;
 370
 371 void start_time(void)
 372 {
 373   cprev   = clock();
 374   runtime = 0.0;
 375 }
 376
 377 void update_time(void)
 378 {
 379   clock_t c;
 380
 381   c        = clock();
 382   runtime += (c-cprev)/(double)CLOCKS_PER_SEC;
 383   cprev    = c;
 384 }
 385
 386 double node_time(void)
 387 {
 388   return runtime;
 389 }
 390
 391 void do_shakefirst(FILE *log,bool bTYZ,real lambda,real ener[],
 392                    t_parm *parm,t_nsborder *nsb,t_mdatoms *md,
 393                    rvec x[],rvec vold[],rvec buf[],rvec f[],
 394                    rvec v[],t_graph *graph,t_commrec *cr,t_nrnb *nrnb,
 395                    t_groups *grps,t_forcerec *fr,t_topology *top,
 396                    t_edsamyn *edyn,t_pull *pulldata)
 397 {
 398   int    i,m,start,homenr,end,step;
 399   tensor shake_vir;
 400   double mass,tmass,vcm[4];
 401   real   dt=parm->ir.delta_t;
 402   real   dt_1;
 403
 404   if (count_constraints(top,cr)) {
 405     start  = START(nsb);
 406     homenr = HOMENR(nsb);
 407     end    = start+homenr;
 408     if (debug)
 409       fprintf(debug,"vcm: start=%d, homenr=%d, end=%d\n",start,homenr,end);
 410     /* Do a first SHAKE to reset particles... */
 411     step = -2;
 412     fprintf(log,"\nConstraining the starting coordinates (step %d)\n",step);
 413     clear_mat(shake_vir);
 414     update(nsb->natoms,start,homenr,step,lambda,&ener[F_DVDL],
 415            parm,1.0,md,x,graph,
 416            NULL,NULL,vold,NULL,x,top,grps,shake_vir,cr,nrnb,bTYZ,
 417            FALSE,edyn,pulldata,FALSE);
 418     /* Compute coordinates at t=-dt, store them in buf */
 419     /* for(i=0; (i<nsb->natoms); i++) {*/
 420     for(i=start; (i<end); i++) {
 421       for(m=0; (m<DIM); m++) {
 422         f[i][m]=x[i][m];
 423         buf[i][m]=x[i][m]-dt*v[i][m];
 424       }
 425     }
 426
 427     /* Shake the positions at t=-dt with the positions at t=0
 428      * as reference coordinates.
 429      */
 430     step = -1;
 431     fprintf(log,"\nConstraining the coordinates at t0-dt (step %d)\n",step);
 432     clear_mat(shake_vir);
 433     update(nsb->natoms,start,homenr,
 434            step,lambda,&ener[F_DVDL],parm,1.0,md,f,graph,
 435            NULL,NULL,vold,NULL,buf,top,grps,shake_vir,cr,nrnb,bTYZ,FALSE,
 436            edyn,pulldata,FALSE);
 437
 438     /* Compute the velocities at t=-dt/2 using the coordinates at
 439      * t=-dt and t=0
 440      * Compute velocity of center of mass and total mass
 441      */
 442     for(m=0; (m<4); m++)
 443       vcm[m] = 0;
 444     dt_1=1.0/dt;
 445     for(i=start; (i<end); i++) {
 446       /*for(i=0; (i<nsb->natoms); i++) {*/
 447       mass = md->massA[i];
 448       for(m=0; (m<DIM); m++) {
 449         v[i][m]=(x[i][m]-f[i][m])*dt_1;
 450         vcm[m] += v[i][m]*mass;
 451       }
 452       vcm[3] += mass;
 453     }
 454     /* Compute the global sum of vcm */
 455     if (debug)
 456       fprintf(debug,"vcm: %8.3f  %8.3f  %8.3f,"
 457               " total mass = %12.5e\n",vcm[XX],vcm[YY],vcm[ZZ],vcm[3]);
 458     if (PAR(cr))
 459       gmx_sumd(4,vcm,cr);
 460     tmass = vcm[3];
 461     for(m=0; (m<DIM); m++)
 462       vcm[m] /= tmass;
 463     if (debug)
 464       fprintf(debug,"vcm: %8.3f  %8.3f  %8.3f,"
 465               " total mass = %12.5e\n",vcm[XX],vcm[YY],vcm[ZZ],tmass);
 466     /* Now we have the velocity of center of mass, let's remove it */
 467     for(i=start; (i<end); i++) {
 468       for(m=0; (m<DIM); m++)
 469         v[i][m] -= vcm[m];
 470     }
 471   }
 472 }
 473
 474 void calc_dispcorr(FILE *log,int eDispCorr,t_forcerec *fr,int natoms,
 475                    matrix box,tensor pres,tensor virial,real ener[])
 476 {
 477   static bool bFirst=TRUE;
 478   real vol,rc3,spres,svir;
 479   int  m;
 480
 481   ener[F_DISPCORR] = 0.0;
 482   ener[F_PRES]     = trace(pres)/3.0;
 483
 484   if (eDispCorr != edispcNO) {
 485     vol           = det(box);
 486     /* Forget the (small) effect of the shift on the LJ energy *
 487      * when fr->bLJShift = TRUE                                */
 488     rc3              = fr->rvdw*fr->rvdw*fr->rvdw;
 489     ener[F_DISPCORR] = -2.0*natoms*natoms*M_PI*fr->avcsix/(3.0*vol*rc3);
 490     if (eDispCorr == edispcEnerPres) {
 491       spres            = 2.0*ener[F_DISPCORR]*PRESFAC/vol;
 492       svir             = -6.0*ener[F_DISPCORR];
 493       ener[F_PRES]    += spres;
 494       for(m=0; (m<DIM); m++) {
 495         pres[m][m]    += spres;
 496         virial[m][m]  += svir;
 497       }
 498     }
 499     else {
 500       spres = 0;
 501       svir  = 0;
 502     }
 503     if (bFirst) {
 504       if (eDispCorr == edispcEner)
 505         fprintf(log,"Long Range LJ corr. to Epot: %10g\n",ener[F_DISPCORR]);
 506       else if (eDispCorr == edispcEnerPres)
 507         fprintf(log,"Long Range LJ corr. to Epot: %10g, Pres: %10g, Vir: %10g\n",
 508                 ener[F_DISPCORR],spres,svir);
 509       bFirst = FALSE;
 510     }
 511   }
 512   ener[F_EPOT] += ener[F_DISPCORR];
 513   ener[F_ETOT] += ener[F_DISPCORR];
 514 }
 515
 516 void do_pbc_first(FILE *log,t_parm *parm,rvec box_size,t_forcerec *fr,
 517                   t_graph *graph,rvec x[])
 518 {
 519   fprintf(log,"Removing pbc first time\n");
 520   calc_shifts(parm->box,box_size,fr->shift_vec);
 521   mk_mshift(log,graph,parm->box,x);
 522   if (getenv ("NOPBC") == NULL)
 523     shift_self(graph,parm->box,x);
 524   else
 525     fprintf(log,"Not doing first shift_self\n");
 526   fprintf(log,"Done rmpbc\n");
 527 }
 528
 529 void set_pot_bools(t_inputrec *ir,t_topology *top,
 530                    bool *bLR,bool *bLJLR,bool *bBHAM,bool *b14)
 531 {
 532   *bLR   = (ir->rcoulomb > ir->rlist) || EEL_LR(ir->coulombtype);
 533   *bLJLR = (ir->rvdw > ir->rlist);
 534   *bBHAM = (top->idef.functype[0]==F_BHAM);
 535   *b14   = (top->idef.il[F_LJ14].nr > 0);
 536 }
 537
 538 void finish_run(FILE *log,t_commrec *cr,char *confout,
 539                 t_nsborder *nsb,t_topology *top,t_parm *parm,
 540                 t_nrnb nrnb[],double nodetime,double realtime,int step,
 541                 bool bWriteStat)
 542 {
 543   int    i,j;
 544   t_nrnb ntot;
 545   real   runtime;
 546   for(i=0; (i<eNRNB); i++)
 547     ntot.n[i]=0;
 548   for(i=0; (i<nsb->nnodes); i++)
 549     for(j=0; (j<eNRNB); j++)
 550       ntot.n[j]+=nrnb[i].n[j];
 551   runtime=0;
 552   if (bWriteStat) {
 553     runtime=parm->ir.nsteps*parm->ir.delta_t;
 554     if (MASTER(cr)) {
 555       fprintf(stderr,"\n\n");
 556       print_perf(stderr,nodetime,realtime,runtime,&ntot,nsb->nnodes);
 557     }
 558     else
 559       print_nrnb(log,&(nrnb[nsb->nodeid]));
 560   }
 561
 562   if (MASTER(cr)) {
 563     print_perf(log,nodetime,realtime,runtime,&ntot,nsb->nnodes);
 564     if (nsb->nnodes > 1)
 565       pr_load(log,nsb->nnodes,nrnb);
 566   }
 567 }
 568
 569 void init_md(t_commrec *cr,t_inputrec *ir,tensor box,real *t,real *t0,
 570              real *lambda,real *lam0,real *SAfactor,
 571              t_nrnb *mynrnb,bool *bTYZ,t_topology *top,
 572              int nfile,t_filenm fnm[],char **traj,
 573              char **xtc_traj,int *fp_ene,
 574              FILE **fp_dgdl,t_mdebin **mdebin,t_groups *grps,
 575              tensor force_vir,tensor pme_vir,
 576              tensor shake_vir,t_mdatoms *mdatoms,rvec mu_tot,
 577              bool *bNEMD,t_vcm **vcm,t_nsborder *nsb)
 578 {
 579   bool bBHAM,b14,bLR,bLJLR;
 580   int  i;
 581
 582   /* Initial values */
 583   *t = *t0       = ir->init_t;
 584   if (ir->efep != efepNO) {
 585     *lambda = *lam0 = ir->init_lambda;
 586   }
 587   else {
 588     *lambda = *lam0   = 0.0;
 589   }
 590   if (ir->bSimAnn) {
 591     *SAfactor = 1.0 - *t0/ir->zero_temp_time;
 592     if (*SAfactor < 0)
 593       *SAfactor = 0;
 594   } else
 595     *SAfactor     = 1.0;
 596
 597   init_nrnb(mynrnb);
 598
 599   /* Check Environment variables & other booleans */
 600   *bTYZ=getenv("TYZ") != NULL;
 601   set_pot_bools(ir,top,&bLR,&bLJLR,&bBHAM,&b14);
 602
 603   if (nfile != -1) {
 604     /* Filenames */
 605     *traj     = ftp2fn(efTRN,nfile,fnm);
 606     *xtc_traj = ftp2fn(efXTC,nfile,fnm);
 607
 608     if (MASTER(cr)) {
 609       *fp_ene = open_enx(ftp2fn(efENX,nfile,fnm),"w");
 610       if ((fp_dgdl != NULL) && ir->efep!=efepNO)
 611         *fp_dgdl =
 612           xvgropen(opt2fn("-dgdl",nfile,fnm),
 613                    "dG/d\\8l\\4","Time (ps)",
 614                    "dG/d\\8l\\4 (kJ mol\\S-1\\N nm\\S-2\\N \\8l\\4\\S-1\\N)");
 615     } else
 616       *fp_ene = -1;
 617
 618     *mdebin = init_mdebin(*fp_ene,grps,&(top->atoms),&(top->idef),
 619                           bLR,bLJLR,bBHAM,b14,ir->efep!=efepNO,ir->epc,
 620                           ir->eDispCorr,(TRICLINIC(ir->compress) || TRICLINIC(box)),
 621                           (ir->etc==etcNOSEHOOVER),cr);
 622   }
 623
 624   /* Initiate variables */
 625   clear_mat(force_vir);
 626   clear_mat(pme_vir);
 627   clear_mat(shake_vir);
 628   clear_rvec(mu_tot);
 629
 630   /* Set initial values for invmass etc. */
 631   init_mdatoms(mdatoms,*lambda,TRUE);
 632
 633   *vcm = init_vcm(stdlog,top,cr,mdatoms,START(nsb),HOMENR(nsb),ir->nstcomm);
 634
 635   debug_gmx();
 636
 637   *bNEMD = (ir->opts.ngacc > 1) || (norm(ir->opts.acc[0]) > 0);
 638
 639   if (ir->eI == eiSD)
 640     init_sd_consts(ir->opts.ngtc,ir->opts.tau_t,ir->delta_t);
 641
 642 }
 643