Fix for the pressure in MTTK with constraints + dispersion

Fix for pressure in MTTK - broken at some point in the
merge process.  The MTTK+constraints code definitely needs to
be simplified for 5.0, it's too easy to break on code rearrangements
right now!

Testing with simple systems over long time shows that the pressure is
unbiased, while the PR pressure is off by a factor of 1/N (which is totally
fine for most realistic systems).

Fixes redmine bug #1061
Change-Id: Ie82cf8d6bdeb3e47db044ae34fc42c077ad73c75

diff --git a/src/kernel/md.c b/src/kernel/md.c
index c72aea114f6b9f8843dbe23a02252bf5e1ac6a5c..ceedc979e5b73aa7a45d9c9c9a46bb86f477613a 100644 (file)
--- a/src/kernel/md.c
+++ b/src/kernel/md.c
@@ -1294,6 +1294,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 

diff --git a/src/mdlib/coupling.c b/src/mdlib/coupling.c
index c9a83119a0c6a3c98bd67503b1cd276723a1dc30..2060205e70cf82882745b5a871ca915afb79d3b2 100644 (file)
--- a/src/mdlib/coupling.c
+++ b/src/mdlib/coupling.c
@@ -178,7 +178,7 @@ static void NHC_trotter(t_grpopts *opts,int nvar, gmx_ekindata_t *ekind,real dtf
 }
 
 static void boxv_trotter(t_inputrec *ir, real *veta, real dt, tensor box, 
-                         gmx_ekindata_t *ekind, tensor vir, real pcorr, real ecorr, t_extmass *MassQ)
+                         gmx_ekindata_t *ekind, tensor vir, real pcorr, t_extmass *MassQ)
 {
 
     real  pscal;
@@ -218,7 +218,7 @@ static void boxv_trotter(t_inputrec *ir, real *veta, real dt, tensor box,
     /* for now, we use Elr = 0, because if you want to get it right, you
        really should be using PME. Maybe print a warning? */
 
-    pscal   = calc_pres(ir->ePBC,nwall,box,ekinmod,vir,localpres);
+    pscal   = calc_pres(ir->ePBC,nwall,box,ekinmod,vir,localpres)+pcorr;
 
     vol = det(box);
     GW = (vol*(MassQ->Winv/PRESFAC))*(DIM*pscal - trace(ir->ref_p));   /* W is in ps^2 * bar * nm^3 */
@@ -905,7 +905,7 @@ void trotter_update(t_inputrec *ir,gmx_large_int_t step, gmx_ekindata_t *ekind,
         case etrtBAROV:
         case etrtBAROV2:
             boxv_trotter(ir,&(state->veta),dt,state->box,ekind,vir,
-                         enerd->term[F_PDISPCORR],enerd->term[F_DISPCORR],MassQ);
+                         enerd->term[F_PDISPCORR],MassQ);
             break;
         case etrtBARONHC:
         case etrtBARONHC2: