Fix constraint virial with multiple time stepping

With multiple time stepping the additional nstcalclr-1 force
contribution was constrained to remove it from the virial.
This procedure neglected the non-linear contribution due to
the rotation of constraints. Now the contribution of this force
component to the coordinate update is constrained instead and
the corresponding virial contribution is subtracted from the
constraint virial.

Fixes #1400

Change-Id: If3217f52808bf7491998324f8dc3161bc003ec1b

diff --git a/include/types/forcerec.h b/include/types/forcerec.h
index 920a81db5432cb34f9a367938c9332df2a79e750..7e14b7d5d4e5fa8ec571594cf3c321b553f629d0 100644 (file)
--- a/include/types/forcerec.h
+++ b/include/types/forcerec.h
@@ -332,6 +332,8 @@ typedef struct {
     gmx_bool bTwinRange;
     int      nlr;
     rvec    *f_twin;
+    /* Constraint virial correction for multiple time stepping */
+    tensor   vir_twin_constr;
 
     /* Forces that should not enter into the virial summation:
      * PPPM/PME/Ewald/posres

diff --git a/include/update.h b/include/update.h
index 8fde948258b6b0c8aada62a1c3e41cabb19ee23f..1d53f4a23a8400a97b51b7bbb6204b24cc9c1519 100644 (file)
--- a/include/update.h
+++ b/include/update.h
@@ -106,6 +106,7 @@ void update_coords(FILE             *fplog,
                    rvec             *f, /* forces on home particles */
                    gmx_bool          bDoLR,
                    rvec             *f_lr,
+                   tensor           *vir_lr_constr,
                    t_fcdata         *fcd,
                    gmx_ekindata_t   *ekind,
                    matrix            M,

diff --git a/src/kernel/md.c b/src/kernel/md.c
index 6b864b7a9e7106f8a4b652e928570b033298c084..da68da10f2bc219b4bc67902be786f81529dc0a4 100644 (file)
--- a/src/kernel/md.c
+++ b/src/kernel/md.c
@@ -1237,7 +1237,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
 
             update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC,
-                          f, bUpdateDoLR, fr->f_twin, fcd,
+                          f, bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
                           ekind, M, wcycle, upd, bInitStep, etrtVELOCITY1,
                           cr, nrnb, constr, &top->idef);
 
@@ -1288,6 +1288,12 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                                        cr, nrnb, wcycle, upd, constr,
                                        bInitStep, TRUE, bCalcVir, vetanew);
 
+                    if (bCalcVir && bUpdateDoLR && ir->nstcalclr > 1)
+                    {
+                        /* Correct the virial for multiple time stepping */
+                        m_sub(shake_vir, fr->vir_twin_constr, shake_vir);
+                    }
+
                     if (!bOK && !bFFscan)
                     {
                         gmx_fatal(FARGS, "Constraint error: Shake, Lincs or Settle could not solve the constrains");
@@ -1739,7 +1745,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
 
                     /* velocity half-step update */
                     update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
-                                  bUpdateDoLR, fr->f_twin, fcd,
+                                  bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
                                   ekind, M, wcycle, upd, FALSE, etrtVELOCITY2,
                                   cr, nrnb, constr, &top->idef);
                 }
@@ -1756,7 +1762,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                 bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
 
                 update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
-                              bUpdateDoLR, fr->f_twin, fcd,
+                              bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
                               ekind, M, wcycle, upd, bInitStep, etrtPOSITION, cr, nrnb, constr, &top->idef);
                 wallcycle_stop(wcycle, ewcUPDATE);
 
@@ -1766,6 +1772,12 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                                    cr, nrnb, wcycle, upd, constr,
                                    bInitStep, FALSE, bCalcVir, state->veta);
 
+                if (bCalcVir && bUpdateDoLR && ir->nstcalclr > 1)
+                {
+                    /* Correct the virial for multiple time stepping */
+                    m_sub(shake_vir, fr->vir_twin_constr, shake_vir);
+                }
+
                 if (ir->eI == eiVVAK)
                 {
                     /* erase F_EKIN and F_TEMP here? */
@@ -1784,7 +1796,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                     bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
 
                     update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
-                                  bUpdateDoLR, fr->f_twin, fcd,
+                                  bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
                                   ekind, M, wcycle, upd, bInitStep, etrtPOSITION, cr, nrnb, constr, &top->idef);
                     wallcycle_stop(wcycle, ewcUPDATE);
 

diff --git a/src/mdlib/update.c b/src/mdlib/update.c
index d1f851fbfcb9a1549c99ae4c85f9fc34a18b150e..d19d6015bb05e0d22b484e37e40eb0ea6a680254 100644 (file)
--- a/src/mdlib/update.c
+++ b/src/mdlib/update.c
@@ -1333,49 +1333,6 @@ static void deform(gmx_update_t upd,
     }
 }
 
-static void combine_forces(int nstcalclr,
-                           gmx_constr_t constr,
-                           t_inputrec *ir, t_mdatoms *md, t_idef *idef,
-                           t_commrec *cr,
-                           gmx_large_int_t step,
-                           t_state *state, gmx_bool bMolPBC,
-                           int start, int nrend,
-                           rvec f[], rvec f_lr[],
-                           t_nrnb *nrnb)
-{
-    int  i, d, nm1;
-
-    /* f contains the short-range forces + the long range forces
-     * which are stored separately in f_lr.
-     */
-
-    if (constr != NULL && !(ir->eConstrAlg == econtSHAKE && ir->epc == epcNO))
-    {
-        /* We need to constrain the LR forces separately,
-         * because due to the different pre-factor for the SR and LR
-         * forces in the update algorithm, we can not determine
-         * the constraint force for the coordinate constraining.
-         * Constrain only the additional LR part of the force.
-         */
-        /* MRS -- need to make sure this works with trotter integration -- the constraint calls may not be right.*/
-        constrain(NULL, FALSE, FALSE, constr, idef, ir, NULL, cr, step, 0, md,
-                  state->x, f_lr, f_lr, bMolPBC, state->box, state->lambda[efptBONDED], NULL,
-                  NULL, NULL, nrnb, econqForce, ir->epc == epcMTTK, state->veta, state->veta);
-    }
-
-    /* Add nstcalclr-1 times the LR force to the sum of both forces
-     * and store the result in forces_lr.
-     */
-    nm1 = nstcalclr - 1;
-    for (i = start; i < nrend; i++)
-    {
-        for (d = 0; d < DIM; d++)
-        {
-            f_lr[i][d] = f[i][d] + nm1*f_lr[i][d];
-        }
-    }
-}
-
 void update_tcouple(FILE             *fplog,
                     gmx_large_int_t   step,
                     t_inputrec       *inputrec,
@@ -1519,6 +1476,75 @@ static rvec *get_xprime(const t_state *state, gmx_update_t upd)
     return upd->xp;
 }
 
+static void combine_forces(gmx_update_t upd,
+                           int nstcalclr,
+                           gmx_constr_t constr,
+                           t_inputrec *ir, t_mdatoms *md, t_idef *idef,
+                           t_commrec *cr,
+                           gmx_large_int_t step,
+                           t_state *state, gmx_bool bMolPBC,
+                           int start, int nrend,
+                           rvec f[], rvec f_lr[],
+                           tensor *vir_lr_constr,
+                           t_nrnb *nrnb)
+{
+    int  i, d;
+
+    /* f contains the short-range forces + the long range forces
+     * which are stored separately in f_lr.
+     */
+
+    if (constr != NULL && vir_lr_constr != NULL &&
+        !(ir->eConstrAlg == econtSHAKE && ir->epc == epcNO))
+    {
+        /* We need to constrain the LR forces separately,
+         * because due to the different pre-factor for the SR and LR
+         * forces in the update algorithm, we have to correct
+         * the constraint virial for the nstcalclr-1 extra f_lr.
+         * Constrain only the additional LR part of the force.
+         */
+        /* MRS -- need to make sure this works with trotter integration -- the constraint calls may not be right.*/
+        rvec *xp;
+        real  fac;
+        int   gf = 0;
+
+        xp  = get_xprime(state, upd);
+
+        fac = (nstcalclr - 1)*ir->delta_t*ir->delta_t;
+
+        for (i = md->start; i < md->start + md->homenr; i++)
+        {
+            if (md->cFREEZE != NULL)
+            {
+                gf = md->cFREEZE[i];
+            }
+            for (d = 0; d < DIM; d++)
+            {
+                if ((md->ptype[i] != eptVSite) &&
+                    (md->ptype[i] != eptShell) &&
+                    !ir->opts.nFreeze[gf][d])
+                {
+                    xp[i][d] = state->x[i][d] + fac*f_lr[i][d]*md->invmass[i];
+                }
+            }
+        }
+        constrain(NULL, FALSE, FALSE, constr, idef, ir, NULL, cr, step, 0, md,
+                  state->x, xp, xp, bMolPBC, state->box, state->lambda[efptBONDED], NULL,
+                  NULL, vir_lr_constr, nrnb, econqCoord, ir->epc == epcMTTK, state->veta, state->veta);
+    }
+
+    /* Add nstcalclr-1 times the LR force to the sum of both forces
+     * and store the result in forces_lr.
+     */
+    for (i = start; i < nrend; i++)
+    {
+        for (d = 0; d < DIM; d++)
+        {
+            f_lr[i][d] = f[i][d] + (nstcalclr - 1)*f_lr[i][d];
+        }
+    }
+}
+
 void update_constraints(FILE             *fplog,
                         gmx_large_int_t   step,
                         real             *dvdlambda, /* the contribution to be added to the bonded interactions */
@@ -1837,6 +1863,7 @@ void update_coords(FILE             *fplog,
                    rvec             *f,    /* forces on home particles */
                    gmx_bool          bDoLR,
                    rvec             *f_lr,
+                   tensor           *vir_lr_constr,
                    t_fcdata         *fcd,
                    gmx_ekindata_t   *ekind,
                    matrix            M,
@@ -1898,9 +1925,10 @@ void update_coords(FILE             *fplog,
          * to produce twin time stepping.
          */
         /* is this correct in the new construction? MRS */
-        combine_forces(inputrec->nstcalclr, constr, inputrec, md, idef, cr,
+        combine_forces(upd,
+                       inputrec->nstcalclr, constr, inputrec, md, idef, cr,
                        step, state, bMolPBC,
-                       start, nrend, f, f_lr, nrnb);
+                       start, nrend, f, f_lr, vir_lr_constr, nrnb);
         force = f_lr;
     }
     else