+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];
+ }
+ }
+}
+