Additional fixes to the nsttcouple and nstpcouple working correctly for trotter decom...

diff --git a/include/types/enums.h b/include/types/enums.h
index b6d7cef12159329ffedf18eb389e7d08be01bb83..bc5629e33cb78f5661835a7b344c4cc5f32d6926 100644 (file)
--- a/include/types/enums.h
+++ b/include/types/enums.h
@@ -61,6 +61,11 @@ enum {
   etrtVELOCITY1, etrtVELOCITY2, etrtPOSITION, etrtSKIPALL, etrtNR
 };
 
+/* sequenced parts of the trotter decomposition */
+enum {
+  ettTSEQ0,  ettTSEQ1,  ettTSEQ2,  ettTSEQ3,  ettTSEQ4, ettTSEQMAX
+};
+
 enum {
   epctISOTROPIC, epctSEMIISOTROPIC, epctANISOTROPIC,
   epctSURFACETENSION, epctNR

diff --git a/include/update.h b/include/update.h
index 0432c06005d89449a3a8fc1fa7779a108f59cb3c..382ce7237f857187af5c720e381336505f7381c5 100644 (file)
--- a/include/update.h
+++ b/include/update.h
@@ -194,7 +194,7 @@ extern void destroy_bufstate(t_state *state);
 
 extern void trotter_update(t_inputrec *ir, gmx_large_int_t step, gmx_ekindata_t *ekind, 
                           gmx_enerdata_t *enerd, t_state *state, tensor vir, t_mdatoms *md, 
-                          t_extmass *MassQ, int *trotter_seq);
+                          t_extmass *MassQ, int **trotter_seqlist, int trotter_seqno);
 
 extern int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *Mass, bool bTrotter); 
 

diff --git a/src/kernel/md.c b/src/kernel/md.c
index b44ed9f1a7273fb2ee61afbd63ec6c36ef6262ec..d980486b18c995edd6e3476ac8ef4afcbe9f8933 100644 (file)
--- a/src/kernel/md.c
+++ b/src/kernel/md.c
@@ -1955,7 +1955,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
                 copy_rvecn(state->v,cbuf,0,state->natoms); /* should make this better for parallelizing? */
             } else {
                 /* this is for NHC in the Ekin(t+dt/2) version of vv */
-                trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq[1]);            
+                trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ1);            
             }
 
             update_coords(fplog,step,ir,mdatoms,state,
@@ -1994,7 +1994,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
                            of veta.  */
                         
                         veta_save = state->veta;
-                        trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq[0]);
+                        trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ0);
                         vetanew = state->veta;
                         state->veta = veta_save;
                     } 
@@ -2059,7 +2059,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
                 {
                     if (bTrotter)
                     {
-                        trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq[2]);
+                        trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ2);
                     } 
                     else 
                     {
@@ -2376,7 +2376,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
                         m_add(force_vir,shake_vir,total_vir);
                         clear_mat(shake_vir);
                     }
-                    trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq[3]);
+                    trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ3);
                 /* We can only do Berendsen coupling after we have summed
                  * the kinetic energy or virial. Since the happens
                  * in global_state after update, we should only do it at
@@ -2429,7 +2429,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
                                     cglo_flags | CGLO_TEMPERATURE    
                         );
                     wallcycle_start(wcycle,ewcUPDATE);
-                    trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq[4]);            
+                    trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ4);            
                     /* now we know the scaling, we can compute the positions again again */
                     copy_rvecn(cbuf,state->x,0,state->natoms);
 

diff --git a/src/mdlib/coupling.c b/src/mdlib/coupling.c
index c67ebf3b94740cbb70184db2fad44868689f572d..b9ec60a7513d2c01e19f4ffce99760872ffd30c3 100644 (file)
--- a/src/mdlib/coupling.c
+++ b/src/mdlib/coupling.c
@@ -50,7 +50,6 @@
 #include "update.h"
 #include "mdrun.h"
 
-#define NTROTTERCALLS 5
 #define NTROTTERPARTS 3
 
 /* these integration routines are only referenced inside this file */
@@ -646,22 +645,36 @@ void destroy_bufstate(t_state *state)
 void trotter_update(t_inputrec *ir,gmx_large_int_t step, gmx_ekindata_t *ekind, 
                     gmx_enerdata_t *enerd, t_state *state, 
                     tensor vir, t_mdatoms *md, 
-                    t_extmass *MassQ, int *trotter_seq) 
+                    t_extmass *MassQ, int **trotter_seqlist, int trotter_seqno) 
 {
     
     int n,i,j,d,ntgrp,ngtc,gc=0;
     t_grp_tcstat *tcstat;
     t_grpopts *opts;
+    gmx_large_int_t step_eff;
     real ecorr,pcorr,dvdlcorr;
     real bmass,qmass,reft,kT,dt,nd;
     tensor dumpres,dumvir;
     double *scalefac,dtc;
+    int *trotter_seq;
     rvec sumv,consk;
     bool bCouple;
 
+    if (trotter_seqno <= ettTSEQ2)
+    {
+        step_eff = step-1;  /* the velocity verlet calls are actually out of order -- the first half step
+                               is actually the last half step from the previous step.  Thus the first half step
+                               actually corresponds to the n-1 step*/
+                               
+    } else {
+        step_eff = step;
+    }
+
     bCouple = (ir->nsttcouple == 1 ||
-               do_per_step(step+ir->nsttcouple,ir->nsttcouple));
-    
+               do_per_step(step_eff+ir->nsttcouple,ir->nsttcouple));
+
+    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))
     {
@@ -849,8 +862,8 @@ int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *MassQ, bool bTrot
     }
     
     /* first, initialize clear all the trotter calls */
-    snew(trotter_seq,NTROTTERCALLS);
-    for (i=0;i<NTROTTERCALLS;i++) 
+    snew(trotter_seq,ettTSEQMAX);
+    for (i=0;i<ettTSEQMAX;i++) 
     {
         snew(trotter_seq[i],NTROTTERPARTS);
         for (j=0;j<NTROTTERPARTS;j++) {

diff --git a/src/tools/gmx_membed.c b/src/tools/gmx_membed.c
index b5162b97fc6dce64b74ef07fa1917298934c78d9..78f368df409e414147d3b8cae42c3e847f4167e5 100644 (file)
--- a/src/tools/gmx_membed.c
+++ b/src/tools/gmx_membed.c
@@ -2690,7 +2690,7 @@ double do_md_membed(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
                 /* this is for NHC in the Ekin(t+dt/2) version of vv */
                 if (!bInitStep)
                 {
-                    trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq[1]);
+                 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ2);
                 }
 
                if (ir->eI == eiVVAK)
@@ -2734,7 +2734,7 @@ double do_md_membed(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
                            of veta.  */
 
                         veta_save = state->veta;
-                        trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq[0]);
+                        trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ0);
                         vetanew = state->veta;
                         state->veta = veta_save;
                     }
@@ -2792,7 +2792,7 @@ double do_md_membed(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
                 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
                 if (bVV && !bInitStep)
                 {
-                    trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq[2]);
+                 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ, trotter_seq,ettTSEQ2);
                 }
 
                 if (bIterations &&
@@ -3106,7 +3106,7 @@ double do_md_membed(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
                         m_add(force_vir,shake_vir,total_vir);
                         clear_mat(shake_vir);
                     }
-                    trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq[3]);
+                    trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ, trotter_seq,ettTSEQ3);
                 }
                 /* We can only do Berendsen coupling after we have summed
                  * the kinetic energy or virial. Since the happens
@@ -3155,7 +3155,7 @@ double do_md_membed(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
                                     cglo_flags | CGLO_TEMPERATURE | CGLO_CONSTRAINT
                         );
                     wallcycle_start(wcycle,ewcUPDATE);
-                    trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq[4]);
+                    trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ, trotter_seq,ettTSEQ4);
                     /* now we know the scaling, we can compute the positions again again */
                     copy_rvecn(cbuf,state->x,0,state->natoms);