From f6e7b9f55e96666f809baf41039a6afbce7cb3dd Mon Sep 17 00:00:00 2001 From: Michael Shirts Date: Sun, 13 Jan 2013 12:20:57 -0500 Subject: [PATCH] fixes issues with pressure control and infrequent evaluation Fixes a bug when pressure control in md-vv when nstcalcenergy is not a multiple of nstpcouple or nsttcouple. This bug results in boxes slowly expanding to unphysical sizes because the virial is neglected in the second half of the md-vv calculation. Also discovered that as part of the bug, global energies were being communicated where they did not need to be when nstpcouple and nsttcouple are > 1 in the case of md-vv, so redid some of the iteration counting and global communication to fix this all together. In the process, this simplified some of the interation counting. Change-Id: I3f77ba361fc88e03ccdfdde1b0be388ea46a8683 --- include/types/iteratedconstraints.h | 2 +- src/kernel/md.c | 108 +++++++++++++++++++----------------- src/kernel/readir.c | 13 ++--- src/mdlib/coupling.c | 6 +- src/mdlib/iteratedconstraints.c | 8 ++- 5 files changed, 71 insertions(+), 66 deletions(-) diff --git a/include/types/iteratedconstraints.h b/include/types/iteratedconstraints.h index 2af24618c8..3f75aab4b9 100644 --- a/include/types/iteratedconstraints.h +++ b/include/types/iteratedconstraints.h @@ -59,7 +59,7 @@ typedef struct { real f,fprev,x,xprev; int iter_i; - gmx_bool bIterate; + gmx_bool bIterationActive; real allrelerr[MAXITERCONST+2]; int num_close; /* number of "close" violations, caused by limited precision. */ } gmx_iterate_t; diff --git a/src/kernel/md.c b/src/kernel/md.c index e18a8d70cd..e349e139f2 100644 --- a/src/kernel/md.c +++ b/src/kernel/md.c @@ -204,7 +204,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[], gmx_bool bTCR=FALSE,bConverged=TRUE,bOK,bSumEkinhOld,bExchanged; gmx_bool bAppend; gmx_bool bResetCountersHalfMaxH=FALSE; - gmx_bool bVV,bIterations,bFirstIterate,bTemp,bPres,bTrotter; + gmx_bool bVV,bIterativeCase,bFirstIterate,bTemp,bPres,bTrotter; gmx_bool bUpdateDoLR; real mu_aver=0,dvdl; int a0,a1,gnx=0,ii; @@ -268,7 +268,11 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[], snew(cbuf,top_global->natoms); } /* all the iteratative cases - only if there are constraints */ - bIterations = ((IR_NPH_TROTTER(ir) || IR_NPT_TROTTER(ir)) && (constr) && (!bRerunMD)); + bIterativeCase = ((IR_NPH_TROTTER(ir) || IR_NPT_TROTTER(ir)) && (constr) && (!bRerunMD)); + gmx_iterate_init(&iterate,FALSE); /* The default value of iterate->bIterationActive is set to + false in this step. The correct value, true or false, + is set at each step, as it depends on the frequency of temperature + and pressure control.*/ bTrotter = (bVV && (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir) || IR_NVT_TROTTER(ir))); if (bRerunMD) @@ -620,7 +624,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[], } /* if using an iterative algorithm, we need to create a working directory for the state. */ - if (bIterations) + if (bIterativeCase) { bufstate = init_bufstate(state); } @@ -1067,15 +1071,21 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[], */ if (EI_VV(ir->eI) && (!bInitStep)) { - /* for vv, the first half actually corresponds to the last step */ + /* for vv, the first half of the integration actually corresponds + to the previous step. bCalcEner is only required to be evaluated on the 'next' step, + but the virial needs to be calculated on both the current step and the 'next' step. Future + reorganization may be able to get rid of one of the bCalcVir=TRUE steps. */ + bCalcEner = do_per_step(step-1,ir->nstcalcenergy); + bCalcVir = bCalcEner || + (ir->epc != epcNO && (do_per_step(step,ir->nstpcouple) || do_per_step(step-1,ir->nstpcouple))); } else { bCalcEner = do_per_step(step,ir->nstcalcenergy); + bCalcVir = bCalcEner || + (ir->epc != epcNO && do_per_step(step,ir->nstpcouple)); } - bCalcVir = bCalcEner || - (ir->epc != epcNO && do_per_step(step,ir->nstpcouple)); /* Do we need global communication ? */ bGStat = (bCalcVir || bCalcEner || bStopCM || @@ -1192,9 +1202,9 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[], ekind,M,wcycle,upd,bInitStep,etrtVELOCITY1, cr,nrnb,constr,&top->idef); - if (bIterations) + if (bIterativeCase && do_per_step(step-1,ir->nstpcouple) && !bInitStep) { - gmx_iterate_init(&iterate,bIterations && !bInitStep); + gmx_iterate_init(&iterate,TRUE); } /* for iterations, we save these vectors, as we will be self-consistently iterating the calculations */ @@ -1202,14 +1212,14 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[], /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */ /* save the state */ - if (bIterations && iterate.bIterate) { + if (iterate.bIterationActive) { copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts)); } bFirstIterate = TRUE; - while (bFirstIterate || (bIterations && iterate.bIterate)) + while (bFirstIterate || iterate.bIterationActive) { - if (bIterations && iterate.bIterate) + if (iterate.bIterationActive) { copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts)); if (bFirstIterate && bTrotter) @@ -1251,7 +1261,6 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[], unshift_self(graph,state->box,state->x); } - /* if VV, compute the pressure and constraints */ /* For VV2, we strictly only need this if using pressure * control, but we really would like to have accurate pressures @@ -1260,34 +1269,37 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[], * For now, keep this choice in comments. */ /*bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir)); */ - /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));*/ + /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));*/ bPres = TRUE; bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK)); if (bCalcEner && ir->eI==eiVVAK) /*MRS: 7/9/2010 -- this still doesn't fix it?*/ { bSumEkinhOld = TRUE; } - compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm, - wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot, - constr,NULL,FALSE,state->box, - top_global,&pcurr,top_global->natoms,&bSumEkinhOld, - cglo_flags - | CGLO_ENERGY - | (bTemp ? CGLO_TEMPERATURE:0) - | (bPres ? CGLO_PRESSURE : 0) - | (bPres ? CGLO_CONSTRAINT : 0) - | ((bIterations && iterate.bIterate) ? CGLO_ITERATE : 0) - | (bFirstIterate ? CGLO_FIRSTITERATE : 0) - | CGLO_SCALEEKIN - ); - /* explanation of above: - a) We compute Ekin at the full time step - if 1) we are using the AveVel Ekin, and it's not the - initial step, or 2) if we are using AveEkin, but need the full - time step kinetic energy for the pressure (always true now, since we want accurate statistics). - b) If we are using EkinAveEkin for the kinetic energy for the temperture control, we still feed in - EkinAveVel because it's needed for the pressure */ - + /* for vv, the first half of the integration actually corresponds to the previous step. + So we need information from the last step in the first half of the integration */ + if (bGStat || do_per_step(step-1,nstglobalcomm)) { + compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm, + wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot, + constr,NULL,FALSE,state->box, + top_global,&pcurr,top_global->natoms,&bSumEkinhOld, + cglo_flags + | CGLO_ENERGY + | (bTemp ? CGLO_TEMPERATURE:0) + | (bPres ? CGLO_PRESSURE : 0) + | (bPres ? CGLO_CONSTRAINT : 0) + | ((iterate.bIterationActive) ? CGLO_ITERATE : 0) + | (bFirstIterate ? CGLO_FIRSTITERATE : 0) + | CGLO_SCALEEKIN + ); + /* explanation of above: + a) We compute Ekin at the full time step + if 1) we are using the AveVel Ekin, and it's not the + initial step, or 2) if we are using AveEkin, but need the full + time step kinetic energy for the pressure (always true now, since we want accurate statistics). + b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in + EkinAveVel because it's needed for the pressure */ + } /* temperature scaling and pressure scaling to produce the extended variables at t+dt */ if (!bInitStep) { @@ -1313,7 +1325,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[], } } - if (bIterations && + if (iterate.bIterationActive && done_iterating(cr,fplog,step,&iterate,bFirstIterate, state->veta,&vetanew)) { @@ -1594,21 +1606,18 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[], } } - if (bIterations) - { - gmx_iterate_init(&iterate,bIterations); - } - - /* for iterations, we save these vectors, as we will be redoing the calculations */ - if (bIterations && iterate.bIterate) + if (bIterativeCase && do_per_step(step,ir->nstpcouple)) { + gmx_iterate_init(&iterate,TRUE); + /* for iterations, we save these vectors, as we will be redoing the calculations */ copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts)); } + bFirstIterate = TRUE; - while (bFirstIterate || (bIterations && iterate.bIterate)) + while (bFirstIterate || iterate.bIterationActive) { /* We now restore these vectors to redo the calculation with improved extended variables */ - if (bIterations) + if (iterate.bIterationActive) { copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts)); } @@ -1633,7 +1642,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[], /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */ if (bTrotter) { - if (bIterations && iterate.bIterate) + if (iterate.bIterationActive) { if (bFirstIterate) { @@ -1767,13 +1776,12 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[], wallcycle_stop(wcycle,ewcVSITECONSTR); } - /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */ + /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */ /* With Leap-Frog we can skip compute_globals at * non-communication steps, but we need to calculate * the kinetic energy one step before communication. */ - if (bGStat || do_per_step(step+1,nstglobalcomm) || - EI_VV(ir->eI)) + if (bGStat || (!EI_VV(ir->eI)&&do_per_step(step+1,nstglobalcomm))) { if (ir->nstlist == -1 && bFirstIterate) { @@ -1792,7 +1800,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[], | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0) | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0) | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0) - | (bIterations && iterate.bIterate ? CGLO_ITERATE : 0) + | (iterate.bIterationActive ? CGLO_ITERATE : 0) | (bFirstIterate ? CGLO_FIRSTITERATE : 0) | CGLO_CONSTRAINT ); @@ -1810,7 +1818,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[], but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could generate the new shake_vir, but test the veta value for convergence. This will take some thought. */ - if (bIterations && + if (iterate.bIterationActive && done_iterating(cr,fplog,step,&iterate,bFirstIterate, trace(shake_vir),&tracevir)) { diff --git a/src/kernel/readir.c b/src/kernel/readir.c index e22db5cacf..a03c66e462 100644 --- a/src/kernel/readir.c +++ b/src/kernel/readir.c @@ -2683,16 +2683,13 @@ void do_index(const char* mdparin, const char *ndx, } if ((ir->epc==epcMTTK) && (ir->etc>etcNO)) { - int mincouple; - mincouple = ir->nsttcouple; - if (ir->nstpcouple < mincouple) + if (ir->nstpcouple != ir->nsttcouple) { - mincouple = ir->nstpcouple; + int mincouple = min(ir->nstpcouple,ir->nsttcouple); + ir->nstpcouple = ir->nsttcouple = mincouple; + sprintf(warn_buf,"for current Trotter decomposition methods with vv, nsttcouple and nstpcouple must be equal. Both have been reset to min(nsttcouple,nstpcouple) = %d",mincouple); + warning_note(wi,warn_buf); } - ir->nstpcouple = mincouple; - ir->nsttcouple = mincouple; - sprintf(warn_buf,"for current Trotter decomposition methods with vv, nsttcouple and nstpcouple must be equal. Both have been reset to min(nsttcouple,nstpcouple) = %d",mincouple); - warning_note(wi,warn_buf); } } /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented diff --git a/src/mdlib/coupling.c b/src/mdlib/coupling.c index 2060205e70..43bc09e013 100644 --- a/src/mdlib/coupling.c +++ b/src/mdlib/coupling.c @@ -872,13 +872,11 @@ void trotter_update(t_inputrec *ir,gmx_large_int_t step, gmx_ekindata_t *ekind, trotter_seq = trotter_seqlist[trotter_seqno]; - /* signal we are returning if nothing is going to be done in this routine */ - if ((trotter_seq[0] == etrtSKIPALL) || !(bCouple)) + if ((trotter_seq[0] == etrtSKIPALL) || (!bCouple)) { return; } - - dtc = ir->nsttcouple*ir->delta_t; + dtc = ir->nsttcouple*ir->delta_t; /* This is OK for NPT, because nsttcouple == nstpcouple is enforcesd */ opts = &(ir->opts); /* just for ease of referencing */ ngtc = opts->ngtc; assert(ngtc>0); diff --git a/src/mdlib/iteratedconstraints.c b/src/mdlib/iteratedconstraints.c index 4c17f1077e..b925b31c67 100644 --- a/src/mdlib/iteratedconstraints.c +++ b/src/mdlib/iteratedconstraints.c @@ -64,12 +64,12 @@ /* maximum length of cyclic traps to check, emerging from limited numerical precision */ #define CYCLEMAX 20 -void gmx_iterate_init(gmx_iterate_t *iterate,gmx_bool bIterate) +void gmx_iterate_init(gmx_iterate_t *iterate,gmx_bool bSetIterationActive) { int i; iterate->iter_i = 0; - iterate->bIterate = bIterate; + iterate->bIterationActive = bSetIterationActive; iterate->num_close = 0; for (i=0;ix == iterate->xprev) && iterate->iter_i > 1)) { - iterate->bIterate = FALSE; + iterate->bIterationActive = FALSE; if (debug) { fprintf(debug,"Iterating NPT constraints: CONVERGED\n"); @@ -196,6 +196,7 @@ gmx_bool done_iterating(const t_commrec *cr,FILE *fplog, int nsteps, gmx_iterate Better to give up convergence here than have the simulation die. */ iterate->num_close++; + iterate->bIterationActive = FALSE; return TRUE; } else @@ -207,6 +208,7 @@ gmx_bool done_iterating(const t_commrec *cr,FILE *fplog, int nsteps, gmx_iterate { md_print_warn(cr,fplog,"Slight numerical convergence deviation with NPT at step %d, relative error only %10.5g, likely not a problem, continuing\n",nsteps,relerr); iterate->num_close++; + iterate->bIterationActive = FALSE; return TRUE; /* if more than a few, check the total fraction. If too high, die. */ } else if (iterate->num_close/(double)nsteps > FRACTION_CLOSE) { -- 2.11.4.GIT