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

diff --git a/include/types/iteratedconstraints.h b/include/types/iteratedconstraints.h
index 2af24618c8531fa1f50c8646c626409447b24634..3f75aab4b90e962d19f5c91d9defddc4d755acb8 100644 (file)
--- 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 e18a8d70cdd47cd5c6e59a705718e5ebb3e6eb2f..e349e139f22a18e1cabe09cd24589b33f4b6b46f 100644 (file)
--- 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 e22db5cacfcabd901c8917e73a0afa9604f3dd9b..a03c66e4627ccb17e4efdb74c6ea1fd115afaa40 100644 (file)
--- 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 2060205e70cf82882745b5a871ca915afb79d3b2..43bc09e013b37bd667e81bc6f958ff9cf140cb14 100644 (file)
--- 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 4c17f1077e7dbc670e1eb5f8a0e480c4eea307c8..b925b31c6796796e01bfdd5b49f2801f8186a0c7 100644 (file)
--- 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;i<MAXITERCONST+2;i++) 
     {
@@ -166,7 +166,7 @@ gmx_bool done_iterating(const t_commrec *cr,FILE *fplog, int nsteps, gmx_iterate
         
         if ((relerr < CONVERGEITER) || (err < CONVERGEITER) || (fom==0) || ((iterate->x == 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) {