Merge release-5-0 into master

[alexxy/gromacs.git] / src / gromacs / mdlib / update.c

diff --git a/src/gromacs/mdlib/update.c b/src/gromacs/mdlib/update.c
index e55d50fa9401bad5f05e457eeae7ddb1ec3c62b3..9bdde8a28cb98ef207c43318b9bd4897dadcdddc 100644 (file)
--- a/src/gromacs/mdlib/update.c
+++ b/src/gromacs/mdlib/update.c
@@ -1232,49 +1232,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_int64_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(gmx_int64_t       step,
                     t_inputrec       *inputrec,
                     t_state          *state,
@@ -1413,6 +1370,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_int64_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 = 0; i < 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_int64_t       step,
                         real             *dvdlambda, /* the contribution to be added to the bonded interactions */
@@ -1726,6 +1752,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,
@@ -1786,9 +1813,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