Merge release-4-6 into master

[alexxy/gromacs.git] / src / programs / mdrun / md.c

diff --git a/src/programs/mdrun/md.c b/src/programs/mdrun/md.c
index 386b897e322ff58a7c38adb53ed2bcb644c979b3..da949ade755c10c105406c0d5d0a9a81a2fb5a99 100644 (file)
--- a/src/programs/mdrun/md.c
+++ b/src/programs/mdrun/md.c
@@ -570,7 +570,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
     {
         nstfep = ir->expandedvals->nstexpanded;
     }
-    if (repl_ex_nst > 0 && repl_ex_nst > nstfep)
+    if (repl_ex_nst > 0 && nstfep > repl_ex_nst)
     {
         nstfep = repl_ex_nst;
     }
@@ -1264,7 +1264,6 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
                                 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
                                 cglo_flags 
                                 | CGLO_ENERGY 
-                                | (bStopCM ? CGLO_STOPCM : 0)
                                 | (bTemp ? CGLO_TEMPERATURE:0) 
                                 | (bPres ? CGLO_PRESSURE : 0) 
                                 | (bPres ? CGLO_CONSTRAINT : 0)
@@ -1285,6 +1284,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
                 {
                     if (bTrotter)
                     {
+                        m_add(force_vir,shake_vir,total_vir); /* we need the un-dispersion corrected total vir here */
                         trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ2);
                     } 
                     else 
@@ -1301,9 +1301,6 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
                                             top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
                                             CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
                         }
-
-
-                        update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
                     }
                 }
                 
@@ -1350,7 +1347,10 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
                 saved_conserved_quantity -= enerd->term[F_DISPCORR];
             }
             /* sum up the foreign energy and dhdl terms for vv.  currently done every step so that dhdl is correct in the .edr */
-            sum_dhdl(enerd,state->lambda,ir->fepvals);
+            if (!bRerunMD)
+            {
+                sum_dhdl(enerd,state->lambda,ir->fepvals);
+            }
         }
         
         /* ########  END FIRST UPDATE STEP  ############## */
@@ -1557,20 +1557,26 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
             gs.sig[eglsCHKPT] = 1;
         }
   
-
-        /* at the start of step, randomize the velocities */
-        if (ETC_ANDERSEN(ir->etc) && EI_VV(ir->eI))
+        /* at the start of step, randomize or scale the velocities (trotter done elsewhere) */
+        if (EI_VV(ir->eI))
         {
-            gmx_bool bDoAndersenConstr;
-            bDoAndersenConstr = (constr && update_randomize_velocities(ir,step,mdatoms,state,upd,&top->idef,constr));
-            /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
-            if (bDoAndersenConstr)
+            if (!bInitStep)
             {
-                update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,
-                                   state,fr->bMolPBC,graph,f,
-                                   &top->idef,tmp_vir,NULL,
-                                   cr,nrnb,wcycle,upd,constr,
-                                   bInitStep,TRUE,bCalcVir,vetanew);
+                update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
+            }
+            if (ETC_ANDERSEN(ir->etc)) /* keep this outside of update_tcouple because of the extra info required to pass */
+            {
+                gmx_bool bIfRandomize;
+                bIfRandomize = update_randomize_velocities(ir,step,mdatoms,state,upd,&top->idef,constr);
+                /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
+                if (constr && bIfRandomize)
+                {
+                    update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,
+                                       state,fr->bMolPBC,graph,f,
+                                       &top->idef,tmp_vir,NULL,
+                                       cr,nrnb,wcycle,upd,constr,
+                                       bInitStep,TRUE,bCalcVir,vetanew);
+                }
             }
         }
 
@@ -1764,7 +1770,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
                                 lastbox,
                                 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
                                 cglo_flags 
-                                | (!EI_VV(ir->eI) ? CGLO_ENERGY : 0)
+                                | (!EI_VV(ir->eI) || bRerunMD ? CGLO_ENERGY : 0)
                                 | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
                                 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0) 
                                 | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0) 
@@ -1797,7 +1803,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
 
         /* only add constraint dvdl after constraints */
         enerd->term[F_DVDL_BONDED] += dvdl;
-        if (!bVV)
+        if (!bVV || bRerunMD)
         {
             /* sum up the foreign energy and dhdl terms for md and sd. currently done every step so that dhdl is correct in the .edr */
             sum_dhdl(enerd,state->lambda,ir->fepvals);
@@ -1940,7 +1946,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
             state->fep_state = lamnew;
             for (i=0;i<efptNR;i++)
             {
-                state->lambda[i] = ir->fepvals->all_lambda[i][lamnew];
+                state_global->lambda[i] = ir->fepvals->all_lambda[i][lamnew];
             }
         }
         /* Remaining runtime */