/* * $Id: md.c,v 1.112.2.12 2007/09/20 08:56:03 spoel Exp $ * * This source code is part of * * G R O M A C S * * GROningen MAchine for Chemical Simulations * * VERSION 3.3.2 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others. * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2007, The GROMACS development team, * check out http://www.gromacs.org for more information. * This program is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License * as published by the Free Software Foundation; either version 2 * of the License, or (at your option) any later version. * * If you want to redistribute modifications, please consider that * scientific software is very special. Version control is crucial - * bugs must be traceable. We will be happy to consider code for * inclusion in the official distribution, but derived work must not * be called official GROMACS. Details are found in the README & COPYING * files - if they are missing, get the official version at www.gromacs.org. * * To help us fund GROMACS development, we humbly ask that you cite * the papers on the package - you can find them in the top README file. * * For more info, check our website at http://www.gromacs.org * * And Hey: * Groningen Machine for Chemical Simulation */ #ifdef HAVE_CONFIG_H #include #endif #include #include #include "typedefs.h" #include "smalloc.h" #include "sysstuff.h" #include "vec.h" #include "statutil.h" #include "vcm.h" #include "mdebin.h" #include "nrnb.h" #include "calcmu.h" #include "index.h" #include "vsite.h" #include "update.h" #include "ns.h" #include "trnio.h" #include "mdrun.h" #include "confio.h" #include "network.h" #include "pull.h" #include "xvgr.h" #include "physics.h" #include "names.h" #include "xmdrun.h" #include "disre.h" #include "orires.h" #include "dihre.h" #include "pppm.h" #include "pme.h" #include "mdatoms.h" #include "rmpbc.h" #include "repl_ex.h" #include "qmmm.h" #include "compute_io.h" #ifdef USE_MPE #include "mpe.h" #include "mpelogging.h" #endif volatile bool bGotTermSignal = FALSE, bGotUsr1Signal = FALSE; static RETSIGTYPE signal_handler(int n) { switch (n) { case SIGTERM: bGotTermSignal = TRUE; break; #ifndef WIN32 case SIGUSR1: bGotUsr1Signal = TRUE; break; #endif } } void mdrunner(t_commrec *cr,t_commrec *mcr,int nfile,t_filenm fnm[], bool bVerbose,bool bCompact, int nDlb,int nstepout, t_edsamyn *edyn,int repl_ex_nst,int repl_ex_seed, unsigned long Flags) { double nodetime=0,realtime; t_inputrec *inputrec; t_state *state; rvec *buf,*f,*vold,*vt; real tmpr1,tmpr2; real *ener; t_nrnb *nrnb; t_nsborder *nsb; t_topology *top; t_groups *grps; t_graph *graph; t_mdatoms *mdatoms; t_forcerec *fr; t_fcdata *fcd; time_t start_t=0; bool bVsites,bParVsites; t_comm_vsites vsitecomm; int i,m; char *gro; /* Initiate everything (snew sets to zero!) */ snew(ener,F_NRE); snew(fcd,1); snew(nsb,1); snew(top,1); snew(grps,1); snew(inputrec,1); snew(state,1); snew(nrnb,cr->nnodes); if (bVerbose && MASTER(cr)) fprintf(stderr,"Getting Loaded...\n"); if (PAR(cr)) { /* The master thread on the master node reads from disk, then passes everything * around the ring, and finally frees the stuff */ if (MASTER(cr)) distribute_parts(cr->left,cr->right,cr->nodeid,cr->nnodes,inputrec, ftp2fn(efTPX,nfile,fnm),nDlb); /* Every node (including the master) reads the data from the ring */ init_parts(stdlog,cr, inputrec,top,state,&mdatoms,nsb, MASTER(cr) ? LIST_SCALARS | LIST_INPUTREC : 0, &bParVsites,&vsitecomm); } else { /* Read it up... */ init_single(stdlog,inputrec,ftp2fn(efTPX,nfile,fnm),top,state,&mdatoms,nsb); bParVsites=FALSE; } if (inputrec->eI == eiSD) { /* Is not read from TPR yet, so we allocate space here */ snew(state->sd_X,nsb->natoms); } snew(buf,nsb->natoms); snew(f,nsb->natoms); snew(vt,nsb->natoms); snew(vold,nsb->natoms); if (bVerbose && MASTER(cr)) fprintf(stderr,"Loaded with Money\n\n"); /* Index numbers for parallellism... */ nsb->nodeid = cr->nodeid; top->idef.nodeid = cr->nodeid; /* Group stuff (energies etc) */ init_groups(stdlog,mdatoms,&(inputrec->opts),grps); /* Copy the cos acceleration to the groups struct */ grps->cosacc.cos_accel = inputrec->cos_accel; /* Periodicity stuff */ if (inputrec->ePBC == epbcXYZ) { graph=mk_graph(&(top->idef),top->atoms.nr,FALSE,FALSE); if (debug) p_graph(debug,"Initial graph",graph); } else graph = NULL; /* Distance Restraints */ init_disres(stdlog,top->idef.il[F_DISRES].nr,top->idef.il[F_DISRES].iatoms, top->idef.iparams,inputrec,mcr,fcd); /* Orientation restraints */ init_orires(stdlog,top->idef.il[F_ORIRES].nr,top->idef.il[F_ORIRES].iatoms, top->idef.iparams,state->x,mdatoms,inputrec,mcr, &(fcd->orires)); /* Dihedral Restraints */ init_dihres(stdlog,top->idef.il[F_DIHRES].nr,top->idef.il[F_DIHRES].iatoms, top->idef.iparams,inputrec,fcd); /* check if there are vsites */ bVsites=FALSE; for(i=0; (iidef.il[i].nr > 0)); /* Initiate forcerecord */ fr = mk_forcerec(); init_forcerec(stdlog,fr,inputrec,top,cr,mdatoms,nsb,state->box,FALSE, opt2fn("-table",nfile,fnm),opt2fn("-tablep",nfile,fnm),FALSE); fr->bSepDVDL = ((Flags & MD_SEPDVDL) == MD_SEPDVDL); /* Initialize QM-MM */ if(fr->bQMMM){ init_QMMMrec(cr,mdatoms,state->box,top,inputrec,fr); } /* Initiate PPPM if necessary */ if (fr->eeltype == eelPPPM) init_pppm(stdlog,cr,nsb,FALSE,TRUE,state->box,getenv("GMXGHAT"),inputrec); if ((fr->eeltype == eelPME) || (fr->eeltype == eelPMEUSER)) (void) init_pme(stdlog,cr,inputrec->nkx,inputrec->nky,inputrec->nkz, inputrec->pme_order, /*HOMENR(nsb),*/nsb->natoms, mdatoms->bChargePerturbed, inputrec->bOptFFT,inputrec->ewald_geometry); /* Make molecules whole at start of run */ if (fr->ePBC != epbcNONE) { do_pbc_first(stdlog,state->box,fr,graph,state->x); } /* Now do whatever the user wants us to do (how flexible...) */ switch (inputrec->eI) { case eiMD: case eiSD: case eiBD: start_t=do_md(stdlog,cr,mcr,nfile,fnm, bVerbose,bCompact,bVsites, bParVsites ? &vsitecomm : NULL, nstepout,inputrec,grps,top,ener,fcd,state,vold,vt,f,buf, mdatoms,nsb,nrnb,graph,edyn,fr, repl_ex_nst,repl_ex_seed, Flags); break; case eiCG: start_t=do_cg(stdlog,nfile,fnm,inputrec,top,grps,nsb, state,f,buf,mdatoms,ener,fcd, nrnb,bVerbose,bVsites, bParVsites ? &vsitecomm : NULL, cr,mcr,graph,fr); break; case eiLBFGS: start_t=do_lbfgs(stdlog,nfile,fnm,inputrec,top,grps,nsb, state,f,buf,mdatoms,ener,fcd, nrnb,bVerbose,bVsites, bParVsites ? &vsitecomm : NULL, cr,mcr,graph,fr); break; case eiSteep: start_t=do_steep(stdlog,nfile,fnm,inputrec,top,grps,nsb, state,f,buf,mdatoms,ener,fcd, nrnb,bVerbose,bVsites, bParVsites ? &vsitecomm : NULL, cr,mcr,graph,fr); break; case eiNM: start_t=do_nm(stdlog,cr,nfile,fnm, bVerbose,bCompact,nstepout,inputrec,grps, top,ener,fcd,state,vold,vt,f,buf, mdatoms,nsb,nrnb,graph,edyn,fr); break; case eiTPI: start_t=do_tpi(stdlog,nfile,fnm,inputrec,top,grps,nsb, state,f,buf,mdatoms,ener,fcd, nrnb,bVerbose, cr,mcr,graph,fr); break; default: gmx_fatal(FARGS,"Invalid integrator (%d)...\n",inputrec->eI); } /* Some timing stats */ if (MASTER(cr)) { realtime=difftime(time(NULL),start_t); if ((nodetime=node_time()) == 0) nodetime=realtime; } else realtime=0; /* Convert back the atoms */ md2atoms(mdatoms,&(top->atoms),TRUE); /* Finish up, write some stuff * if rerunMD, don't write last frame again */ finish_run(stdlog,cr,ftp2fn(efSTO,nfile,fnm), nsb,top,inputrec,nrnb,nodetime,realtime,inputrec->nsteps, EI_DYNAMICS(inputrec->eI)); /* Does what it says */ print_date_and_time(stdlog,cr->nodeid,"Finished mdrun"); } time_t do_md(FILE *log,t_commrec *cr,t_commrec *mcr,int nfile,t_filenm fnm[], bool bVerbose,bool bCompact, bool bVsites, t_comm_vsites *vsitecomm, int stepout,t_inputrec *inputrec,t_groups *grps,t_topology *top, real ener[],t_fcdata *fcd, t_state *state,rvec vold[],rvec vt[],rvec f[], rvec buf[],t_mdatoms *mdatoms,t_nsborder *nsb,t_nrnb nrnb[], t_graph *graph,t_edsamyn *edyn,t_forcerec *fr, int repl_ex_nst,int repl_ex_seed, unsigned long Flags) { t_mdebin *mdebin; int fp_ene=0,fp_trn=0,step,step_rel; FILE *fp_dgdl=NULL,*fp_field=NULL; time_t start_t; real t,t0,lam0; bool bNS,bSimAnn,bStopCM,bRerunMD,bNotLastFrame=FALSE, bFirstStep,bLastStep,bNEMD,do_log,bRerunWarnNoV=TRUE, bFullPBC,bForceUpdate=FALSE; tensor force_vir,shake_vir,total_vir,pres,ekin; t_nrnb mynrnb; char *traj,*xtc_traj; /* normal and compressed trajectory filename */ int i,m,status; rvec mu_tot; rvec *xx,*vv,*ff; t_vcm *vcm; t_trxframe rerun_fr; t_pull pulldata; /* for pull code */ gmx_repl_ex_t *repl_ex=NULL; /* A boolean (disguised as a real) to terminate mdrun */ real terminate=0; /* XMDRUN stuff: shell, general coupling etc. */ bool bFFscan; int nshell,nflexcon,nshell_flexcon_tot,count,nconverged=0; t_shell *shells=NULL; real timestep=0; double tcount=0; bool bShell_FlexCon,bIonize=FALSE,bGlas=FALSE; bool bTCR=FALSE,bConverged=TRUE,bOK,bExchanged; real temp0,mu_aver=0,fmax; int gnx,ii; atom_id *grpindex; char *grpname; t_coupl_rec *tcr=NULL; rvec *xcopy=NULL,*vcopy=NULL; matrix boxcopy,lastbox; /* End of XMDRUN stuff */ /* Turn on signal handling */ if (getenv("GMX_NO_TERM") == NULL) signal(SIGTERM,signal_handler); #ifndef WIN32 if (getenv("GMX_NO_USR1") == NULL) signal(SIGUSR1,signal_handler); #endif /* Check for special mdrun options */ bRerunMD = (Flags & MD_RERUN) == MD_RERUN; bIonize = (Flags & MD_IONIZE) == MD_IONIZE; bGlas = (Flags & MD_GLAS) == MD_GLAS; bFFscan = (Flags & MD_FFSCAN) == MD_FFSCAN; /* Initial values */ init_md(cr,inputrec,&t,&t0,&state->lambda,&lam0,&mynrnb,top, nfile,fnm,&traj,&xtc_traj,&fp_ene,&fp_dgdl,&fp_field,&mdebin,grps, force_vir,shake_vir,mdatoms,mu_tot,&bNEMD,&bSimAnn,&vcm,nsb); debug_gmx(); { double io = compute_io(inputrec,&top->atoms,mdebin->ebin->nener,1); if (io > 2000) fprintf(stderr, "\nWARNING: This run will generate roughly %.0f Mb of data\n\n", io); } /* init edsam, no effect if edyn->bEdsam==FALSE */ init_edsam(stdlog,top,inputrec,mdatoms,START(nsb),HOMENR(nsb),cr,edyn); /* Check for full periodicity calculations */ bFullPBC = (inputrec->ePBC == epbcFULL); /* Check for polarizable models */ shells = init_shells(log,START(nsb),HOMENR(nsb),&top->idef, mdatoms,&nshell); /* Check for flexible constraints */ nflexcon = count_flexible_constraints(log,fr,&top->idef); nshell_flexcon_tot = nshell + nflexcon; if (PAR(cr)) gmx_sumi(1,&nshell_flexcon_tot,cr); bShell_FlexCon = nshell_flexcon_tot > 0; gnx = top->blocks[ebMOLS].nr; snew(grpindex,gnx); for(i=0; (ix,mdatoms,inputrec->opts.nFreeze, state->box,START(nsb),HOMENR(nsb),cr); if (repl_ex_nst > 0) repl_ex = init_replica_exchange(log,mcr,state,inputrec, repl_ex_nst,repl_ex_seed); if (!inputrec->bUncStart && !bRerunMD) do_shakefirst(log,ener,inputrec,nsb,mdatoms,state,vold,buf,f, graph,cr,&mynrnb,grps,fr,top,edyn,&pulldata); debug_gmx(); /* Compute initial EKin for all.. */ if (grps->cosacc.cos_accel == 0) calc_ke_part(START(nsb),HOMENR(nsb),state->v,&(inputrec->opts), mdatoms,grps,&mynrnb,state->lambda); else calc_ke_part_visc(START(nsb),HOMENR(nsb), state->box,state->x,state->v,&(inputrec->opts), mdatoms,grps,&mynrnb,state->lambda); debug_gmx(); if (PAR(cr)) global_stat(log,cr,ener,force_vir,shake_vir, &(inputrec->opts),grps,&mynrnb,nrnb,vcm,&terminate); debug_gmx(); /* Calculate the initial half step temperature */ temp0 = sum_ekin(TRUE,&(inputrec->opts),grps,ekin,NULL); debug_gmx(); /* Initiate data for the special cases */ if (bFFscan) { snew(xcopy,nsb->natoms); snew(vcopy,nsb->natoms); for(ii=0; (iinatoms); ii++) { copy_rvec(state->x[ii],xcopy[ii]); copy_rvec(state->v[ii],vcopy[ii]); } copy_mat(state->box,boxcopy); } /* Write start time and temperature */ start_t=print_date_and_time(log,cr->nodeid,"Started mdrun"); if (MASTER(cr)) { fprintf(log,"Initial temperature: %g K\n",temp0); if (bRerunMD) { fprintf(stderr,"starting md rerun '%s', reading coordinates from" " input trajectory '%s'\n\n", *(top->name),opt2fn("-rerun",nfile,fnm)); if (bVerbose) fprintf(stderr,"Calculated time to finish depends on nsteps from " "run input file,\nwhich may not correspond to the time " "needed to process input trajectory.\n\n"); } else fprintf(stderr,"starting mdrun '%s'\n%d steps, %8.1f ps.\n\n", *(top->name),inputrec->nsteps,inputrec->nsteps*inputrec->delta_t); } /* Initialize values for invmass, etc. */ update_mdatoms(mdatoms,state->lambda,TRUE); /* Set the node time counter to 0 after initialisation */ start_time(); debug_gmx(); /*********************************************************** * * Loop over MD steps * ************************************************************/ /* if rerunMD then read coordinates and velocities from input trajectory */ if (bRerunMD) { if (getenv("GMX_FORCE_UPDATE")) bForceUpdate = TRUE; bNotLastFrame = read_first_frame(&status,opt2fn("-rerun",nfile,fnm), &rerun_fr,TRX_NEED_X | TRX_READ_V); if (rerun_fr.natoms != mdatoms->nr) gmx_fatal(FARGS,"Number of atoms in trajectory (%d) does not match the " "run input file (%d)\n",rerun_fr.natoms,mdatoms->nr); if (inputrec->ePBC != epbcNONE) { if (!rerun_fr.bBox) gmx_fatal(FARGS,"Rerun trajectory frame step %d time %f does not contain a box, while pbc is used",rerun_fr.step,rerun_fr.time); if (max_cutoff2(rerun_fr.box) < sqr(fr->rlistlong)) gmx_fatal(FARGS,"Rerun trajectory frame step %d time %f has too small box dimensions",rerun_fr.step,rerun_fr.time); } } /* loop over MD steps or if rerunMD to end of input trajectory */ bFirstStep = TRUE; bLastStep = FALSE; bExchanged = FALSE; step = inputrec->init_step; step_rel = 0; #ifdef USE_MPE MPE_Start_log( ); #endif while ((!bRerunMD && (step_rel <= inputrec->nsteps)) || (bRerunMD && bNotLastFrame)) { #ifdef USE_MPE /* MPI_Barrier( MPI_COMM_WORLD ); */ MPE_Log_event( ev_timestep1, 0, "time step starts" ); #endif if (bRerunMD) { if (rerun_fr.bStep) { step = rerun_fr.step; step_rel = step - inputrec->init_step; } if (rerun_fr.bTime) t = rerun_fr.time; else t = step; } else { bLastStep = (step_rel == inputrec->nsteps); t = t0 + step*inputrec->delta_t; } do_log = do_per_step(step,inputrec->nstlog) || bLastStep; if (inputrec->efep != efepNO) { if (bRerunMD && rerun_fr.bLambda && (inputrec->delta_lambda!=0)) state->lambda = rerun_fr.lambda; else state->lambda = lam0 + step*inputrec->delta_lambda; } if (MASTER(cr) && do_log && !bFFscan) print_ebin_header(log,step,t,state->lambda); if (bSimAnn) update_annealing_target_temp(&(inputrec->opts),t); if (bRerunMD) { for(i=0; inr; i++) copy_rvec(rerun_fr.x[i],state->x[i]); if (rerun_fr.bV) for(i=0; inr; i++) copy_rvec(rerun_fr.v[i],state->v[i]); else { for(i=0; inr; i++) clear_rvec(state->v[i]); if (bRerunWarnNoV) { fprintf(stderr,"\nWARNING: Some frames do not contain velocities.\n" " Ekin, temperature and pressure are incorrect,\n" " the virial will be incorrect when constraints are present.\n" "\n"); bRerunWarnNoV = FALSE; } } copy_mat(rerun_fr.box,state->box); /* for rerun MD always do Neighbour Searching */ bNS = ((inputrec->nstlist!=0) || bFirstStep); } else { /* Determine whether or not to do Neighbour Searching */ bNS = ((inputrec->nstlist && (step % inputrec->nstlist==0 || bExchanged)) || bFirstStep); } /* Stop Center of Mass motion */ bStopCM = do_per_step(step,abs(inputrec->nstcomm)); /* Copy back starting coordinates in case we're doing a forcefield scan */ if (bFFscan) { for(ii=0; (iinatoms); ii++) { copy_rvec(xcopy[ii],state->x[ii]); copy_rvec(vcopy[ii],state->v[ii]); } copy_mat(boxcopy,state->box); } if (bVsites) { if (graph) { /* Following is necessary because the graph may get out of sync * with the coordinates if we only have every N'th coordinate set */ if (bRerunMD || bExchanged) mk_mshift(log,graph,state->box,state->x); shift_self(graph,state->box,state->x); } construct_vsites(log,state->x,&mynrnb,inputrec->delta_t,state->v, &top->idef,graph,cr,fr->ePBC,state->box,vsitecomm); if (graph) unshift_self(graph,state->box,state->x); } debug_gmx(); /* Set values for invmass etc. This routine not parallellized, but hardly * ever used, only when doing free energy calculations. */ if(inputrec->efep != efepNO) update_mdatoms(mdatoms,state->lambda,FALSE); clear_mat(force_vir); /* Ionize the atoms if necessary */ if (bIonize) ionize(log,mdatoms,top->atoms.atomname,t,inputrec,state->x,state->v, START(nsb),START(nsb)+HOMENR(nsb),state->box,cr); /* Update force field in ffscan program */ if (bFFscan) { if (update_forcefield(nfile,fnm,fr,mdatoms->nr,state->x,state->box)) { if (gmx_parallel_env) gmx_finalize(cr); exit(0); } } if (bShell_FlexCon) { /* Now is the time to relax the shells */ count=relax_shells(log,cr,mcr,bVerbose,bFFscan ? step+1 : step, inputrec,bNS,bStopCM,top,ener,fcd, state,vold,vt,f,buf,mdatoms,nsb,&mynrnb,graph, grps, nshell,shells,nflexcon,fr,traj,t,mu_tot, nsb->natoms,&bConverged,bVsites,vsitecomm, fp_field); tcount+=count; if (bConverged) nconverged++; } else { /* The coordinates (x) are shifted (to get whole molecules) in do_force * This is parallellized as well, and does communication too. * Check comments in sim_util.c */ do_force(log,cr,mcr,inputrec,nsb,step,&mynrnb,top,grps, state->box,state->x,f,buf,mdatoms,ener,fcd,bVerbose && !PAR(cr), state->lambda,graph, TRUE,bNS,FALSE,TRUE,fr,mu_tot,FALSE,t,fp_field,edyn); } if (bTCR) mu_aver=calc_mu_aver(cr,nsb,state->x,mdatoms->chargeA,mu_tot,top,mdatoms, gnx,grpindex); if (bGlas) do_glas(log,START(nsb),HOMENR(nsb),state->x,f, fr,mdatoms,top->idef.atnr,inputrec,ener); if (bTCR && bFirstStep) { tcr=init_coupling(log,nfile,fnm,cr,fr,mdatoms,&(top->idef)); fprintf(log,"Done init_coupling\n"); fflush(log); } /* Now we have the energies and forces corresponding to the * coordinates at time t. We must output all of this before * the update. * for RerunMD t is read from input trajectory */ if (bVsites) spread_vsite_f(log,state->x,f,&mynrnb,&top->idef, fr,graph,state->box,vsitecomm,cr); /* Calculation of the virial must be done after vsites! */ /* Question: Is it correct to do the PME forces after this? */ calc_virial(log,START(nsb),HOMENR(nsb),state->x,f, force_vir,fr->vir_el_recip,graph,state->box,&mynrnb,fr); /* Spread the LR force on virtual sites to the other particles... * This is parallellized. MPI communication is performed * if the constructing atoms aren't local. */ if (bVsites && fr->bEwald) spread_vsite_f(log,state->x,fr->f_el_recip,&mynrnb,&top->idef, fr,graph,state->box,vsitecomm,cr); sum_lrforces(f,fr,START(nsb),HOMENR(nsb)); xx = (do_per_step(step,inputrec->nstxout) || bLastStep) ? state->x : NULL; vv = (do_per_step(step,inputrec->nstvout) || bLastStep) ? state->v : NULL; ff = (do_per_step(step,inputrec->nstfout)) ? f : NULL; fp_trn = write_traj(log,cr,traj,nsb,step,t,state->lambda, nrnb,nsb->natoms,xx,vv,ff,state->box, top); debug_gmx(); /* don't write xtc and last structure for rerunMD */ if (!bRerunMD && !bFFscan) { if (do_per_step(step,inputrec->nstxtcout)) write_xtc_traj(log,cr,xtc_traj,nsb,mdatoms, step,t,state->x,state->box,inputrec->xtcprec); if (bLastStep && MASTER(cr)) { fprintf(stderr,"Writing final coordinates.\n"); write_sto_conf(ftp2fn(efSTO,nfile,fnm), *top->name, &(top->atoms),state->x,state->v,state->box); } debug_gmx(); } clear_mat(shake_vir); /* Afm and Umbrella type pulling happens before the update, * other types in update */ if (pulldata.bPull && (pulldata.runtype == eAfm || pulldata.runtype == eUmbrella)) pull(&pulldata,state->x,f,force_vir,state->box, top,inputrec->delta_t,step,t, mdatoms,START(nsb),HOMENR(nsb),cr); if (bFFscan) clear_rvecs(nsb->natoms,buf); /* Box is changed in update() when we do pressure coupling, * but we should still use the old box for energy corrections and when * writing it to the energy file, so it matches the trajectory files for * the same timestep above. Make a copy in a separate array. */ copy_mat(state->box,lastbox); /* This is also parallellized, but check code in update.c */ /* bOK = update(nsb->natoms,START(nsb),HOMENR(nsb),step,state->lambda,&ener[F_DVDL], */ bOK = TRUE; if (!bRerunMD || rerun_fr.bV || bForceUpdate) update(nsb->natoms,START(nsb),HOMENR(nsb),step,&ener[F_DVDL], inputrec,mdatoms,state,graph,f,vold, top,grps,shake_vir,cr,&mynrnb,edyn,&pulldata,bNEMD, TRUE,bFirstStep,NULL,pres); else { /* Need to unshift here */ if ((inputrec->ePBC == epbcXYZ) && (graph->nnodes > 0)) unshift_self(graph,state->box,state->x); } if (!bOK && !bFFscan) gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains"); /* Correct the new box if it is too skewed */ if (DYNAMIC_BOX(*inputrec) && !bRerunMD) correct_box(state->box,fr,graph); /* (un)shifting should NOT be done after this, * since the box vectors might have changed */ /* Non-equilibrium MD: this is parallellized, but only does communication * when there really is NEMD. */ if (PAR(cr) && bNEMD) accumulate_u(cr,&(inputrec->opts),grps); debug_gmx(); if (grps->cosacc.cos_accel == 0) calc_ke_part(START(nsb),HOMENR(nsb),state->v,&(inputrec->opts), mdatoms,grps,&mynrnb,state->lambda); else calc_ke_part_visc(START(nsb),HOMENR(nsb), state->box,state->x,state->v,&(inputrec->opts), mdatoms,grps,&mynrnb,state->lambda); /* since we use the new coordinates in calc_ke_part_visc, we should use * the new box too. Still, won't this be offset by one timestep in the * energy file? / EL 20040121 */ debug_gmx(); /* Calculate center of mass velocity if necessary, also parallellized */ if (bStopCM && !bFFscan && !bRerunMD) calc_vcm_grp(log,START(nsb),HOMENR(nsb),mdatoms->massT, state->x,state->v,vcm); /* Check whether everything is still allright */ if (bGotTermSignal || bGotUsr1Signal) { if (bGotTermSignal) terminate = 1; else terminate = -1; fprintf(log,"\n\nReceived the %s signal\n\n", bGotTermSignal ? "TERM" : "USR1"); fflush(log); if (MASTER(cr)) { fprintf(stderr,"\n\nReceived the %s signal\n\n", bGotTermSignal ? "TERM" : "USR1"); fflush(stderr); } bGotTermSignal = FALSE; bGotUsr1Signal = FALSE; } if (PAR(cr)) { /* Globally (over all NODEs) sum energy, virial etc. * This includes communication */ global_stat(log,cr,ener,force_vir,shake_vir, &(inputrec->opts),grps,&mynrnb,nrnb,vcm,&terminate); /* Correct for double counting energies, should be moved to * global_stat */ if (fr->bTwinRange && !bNS) for(i=0; (iestat.nn); i++) { grps->estat.ee[egCOULLR][i] /= cr->nnodes; grps->estat.ee[egLJLR][i] /= cr->nnodes; } } else cp_nrnb(&(nrnb[0]),&mynrnb); /* This is just for testing. Nothing is actually done to Ekin * since that would require extra communication. */ if (!bNEMD && debug && (vcm->nr > 0)) correct_ekin(debug,START(nsb),START(nsb)+HOMENR(nsb),state->v, vcm->group_p[0], mdatoms->massT,mdatoms->tmass,ekin); if ((terminate != 0) && (step - inputrec->init_step < inputrec->nsteps)) { if (terminate<0 && inputrec->nstxout) /* this is the USR1 signal and we are writing x to trr, stop at next x frame in trr */ inputrec->nsteps = (step/inputrec->nstxout + 1) * inputrec->nstxout - inputrec->init_step; else inputrec->nsteps = step + 1 - inputrec->init_step; fprintf(log,"\nSetting nsteps to %d\n\n",inputrec->nsteps); fflush(log); if (MASTER(cr)) { fprintf(stderr,"\nSetting nsteps to %d\n\n",inputrec->nsteps); fflush(stderr); } /* erase the terminate signal */ terminate = 0; } /* Do center of mass motion removal */ if (bStopCM && !bFFscan && !bRerunMD) { /* The last parameter is the temperature of the center of mass. * In this case we put 1 there, to assure that the temperature of * COM motion for each group should be below 1 K. Obviously * one could implement it as a tolerance times the actual T as well * e.g. 1e-3 T, but it remains arbitrary anyway. DvdS 2007-07-14. */ check_cm_grp(log,vcm,1); do_stopcm_grp(log,START(nsb),HOMENR(nsb),state->x,state->v,vcm); inc_nrnb(&mynrnb,eNR_STOPCM,HOMENR(nsb)); /* calc_vcm_grp(log,START(nsb),HOMENR(nsb),mdatoms->massT,x,v,vcm); check_cm_grp(log,vcm); do_stopcm_grp(log,START(nsb),HOMENR(nsb),x,v,vcm); check_cm_grp(log,vcm); */ } /* Add force and shake contribution to the virial */ m_add(force_vir,shake_vir,total_vir); /* Sum the potential energy terms from group contributions */ sum_epot(&(inputrec->opts),grps,ener); /* Calculate the amplitude of the cosine velocity profile */ grps->cosacc.vcos = grps->cosacc.mvcos/mdatoms->tmass; /* Sum the kinetic energies of the groups & calc temp */ ener[F_TEMP]=sum_ekin(bRerunMD,&(inputrec->opts),grps,ekin, &(ener[F_DVDLKIN])); ener[F_EKIN]=trace(ekin); /* Calculate Temperature coupling parameters lambda and adjust * target temp when doing simulated annealing */ /* if(inputrec->etc==etcBERENDSEN) berendsen_tcoupl(&(inputrec->opts),grps,inputrec->delta_t); else if(inputrec->etc==etcNOSEHOOVER) nosehoover_tcoupl(&(inputrec->opts),grps,inputrec->delta_t); */ /* Calculate pressure and apply LR correction if PPPM is used. * Use the box from last timestep since we already called update(). */ calc_pres(fr->ePBC,lastbox,ekin,total_vir,pres, (fr->eeltype==eelPPPM) ? ener[F_COUL_RECIP] : 0.0); /* Calculate long range corrections to pressure and energy */ if (bTCR || bFFscan) set_avcsixtwelve(log,fr,mdatoms,&top->atoms.excl); /* Calculate long range corrections to pressure and energy */ calc_dispcorr(log,inputrec,fr,step,mdatoms->nr,lastbox,state->lambda, pres,total_vir,ener); ener[F_ETOT]=ener[F_EPOT]+ener[F_EKIN]; /* Check for excessively large energies */ if (bIonize) { #ifdef GMX_DOUBLE real etot_max = 1e200; #else real etot_max = 1e30; #endif if (fabs(ener[F_ETOT]) > etot_max) { fprintf(stderr,"Energy too large (%g), giving up\n",ener[F_ETOT]); break; } } /* The coordinates (x) were unshifted in update */ if (bFFscan && (!bShell_FlexCon || bConverged)) if (print_forcefield(log,ener,HOMENR(nsb),f,buf,xcopy, &(top->blocks[ebMOLS]),mdatoms->massT,pres)) { if (gmx_parallel_env) gmx_finalize(cr); fprintf(stderr,"\n"); exit(0); } if (bTCR) { /* Only do GCT when the relaxation of shells (minimization) has converged, * otherwise we might be coupling to bogus energies. * In parallel we must always do this, because the other sims might * update the FF. */ /* Since this is called with the new coordinates state->x, I assume * we want the new box state->box too. / EL 20040121 */ do_coupling(log,nfile,fnm,tcr,t,step,ener,fr, inputrec,MASTER(cr), mdatoms,&(top->idef),mu_aver, top->blocks[ebMOLS].nr,(mcr!=NULL) ? mcr : cr, state->box,total_vir,pres, mu_tot,state->x,f,bConverged); debug_gmx(); } /* Time for performance */ if (((step % stepout) == 0) || bLastStep) update_time(); /* Output stuff */ if (MASTER(cr)) { bool do_ene,do_dr,do_or,do_dihr; upd_mdebin(mdebin,fp_dgdl,mdatoms->tmass,step_rel,t,ener,state,lastbox, shake_vir,force_vir,total_vir,pres,grps,mu_tot); do_ene = do_per_step(step,inputrec->nstenergy) || bLastStep; do_dr = do_per_step(step,inputrec->nstdisreout) || bLastStep; do_or = do_per_step(step,inputrec->nstorireout) || bLastStep; do_dihr= do_per_step(step,inputrec->nstdihreout) || bLastStep; print_ebin(fp_ene,do_ene,do_dr,do_or,do_dihr,do_log?log:NULL, step,step_rel,t, eprNORMAL,bCompact,mdebin,fcd,&(top->atoms),&(inputrec->opts)); if (bVerbose) fflush(log); } /* Remaining runtime */ if (MASTER(cr) && bVerbose && ( ((step % stepout)==0) || bLastStep)) { if (bShell_FlexCon) fprintf(stderr,"\n"); print_time(stderr,start_t,step,inputrec); } bExchanged = FALSE; if ((repl_ex_nst > 0) && (step > 0) && !bLastStep && do_per_step(step,repl_ex_nst)) bExchanged = replica_exchange(log,mcr,repl_ex,state,ener[F_EPOT],step,t); bFirstStep = FALSE; if (bRerunMD) /* read next frame from input trajectory */ bNotLastFrame = read_next_frame(status,&rerun_fr); if (!bRerunMD || !rerun_fr.bStep) { /* increase the MD step number */ step++; step_rel++; } } /* End of main MD loop */ #ifdef USE_MPE MPE_Stop_log( ); #endif debug_gmx(); /* Dump the NODE time to the log file on each node */ if (PAR(cr)) { double *ct,ctmax,ctsum; snew(ct,cr->nnodes); ct[cr->nodeid] = node_time(); gmx_sumd(cr->nnodes,ct,cr); ctmax = ct[0]; ctsum = ct[0]; for(i=1; (inodeid); i++) { ctmax = max(ctmax,ct[i]); ctsum += ct[i]; } ctsum /= cr->nnodes; fprintf(log,"\nTotal NODE time on node %d: %g\n",cr->nodeid,ct[cr->nodeid]); fprintf(log,"Average NODE time: %g\n",ctsum); fprintf(log,"Load imbalance reduced performance to %3d%% of max\n", (int) (100.0*ctmax/ctsum)); sfree(ct); } if (bRerunMD) close_trj(status); if (MASTER(cr)) { print_ebin(fp_ene,FALSE,FALSE,FALSE,FALSE,log,step,step_rel,t, eprAVER,FALSE,mdebin,fcd,&(top->atoms),&(inputrec->opts)); print_ebin(fp_ene,FALSE,FALSE,FALSE,FALSE,log,step,step_rel,t, eprRMS,FALSE,mdebin,fcd,&(top->atoms),&(inputrec->opts)); close_enx(fp_ene); if (!bRerunMD && inputrec->nstxtcout) close_xtc_traj(); close_trn(fp_trn); if (fp_dgdl) fclose(fp_dgdl); if (fp_field) fclose(fp_field); } debug_gmx(); /* clean up edsam stuff, no effect if edyn->bEdsam == FALSE */ finish_edsam(stdlog,top,inputrec,mdatoms,START(nsb),HOMENR(nsb),cr,edyn); if (bShell_FlexCon) { fprintf(log,"Fraction of iterations that converged: %.2f %%\n", (nconverged*100.0)/(step_rel+1)); fprintf(log,"Average number of force evaluations per MD step: %.2f\n", tcount/(step_rel+1)); } if (repl_ex_nst > 0) print_replica_exchange_statistics(log,repl_ex); return start_t; }