Fix group-scheme bug with changing LJ parameters in FE

We don't optimize for the case when we have only changed one of charge
or type, so the other vector must always be valid even when it is not
changing. The logic of calling ewald_LRcorrection didn't do this
correctly, perhaps because the construction logic in md2atoms was
unclear.

Changed name, origin and logic for bFreeEnergy to
bHaveChargeOrTypePerturbed to better reflect the correct usage and
meaning. Avoided testing any pointers for NULL - we should use
explicit control-flow constructs.

Fixes #1596

Change-Id: I61172681048075d3022bd6c4b781c6c9153eeadd

diff --git a/src/gromacs/gmxlib/ewald_util.c b/src/gromacs/gmxlib/ewald_util.c
index d3e5f23748e36f07df9c21c073ae01b182fbb41a..a319f0073fbe5f12422d6c1d22fd0f90fd421d27 100644 (file)
--- a/src/gromacs/gmxlib/ewald_util.c
+++ b/src/gromacs/gmxlib/ewald_util.c
@@ -128,12 +128,23 @@ real calc_ewaldcoeff_lj(real rc, real dtol)
     return x;
 }
 
+/* There's nothing special to do here if just masses are perturbed,
+ * but if either charge or type is perturbed then the implementation
+ * requires that B states are defined for both charge and type, and
+ * does not optimize for the cases where only one changes.
+ *
+ * The parameter vectors for B states are left undefined in atoms2md()
+ * when either FEP is inactive, or when there are no mass/charge/type
+ * perturbations. The parameter vectors for LJ-PME are likewise
+ * undefined when LJ-PME is not active. This works because
+ * bHaveChargeOrTypePerturbed handles the control flow. */
 void ewald_LRcorrection(int start, int end,
                         t_commrec *cr, int thread, t_forcerec *fr,
                         real *chargeA, real *chargeB,
                         real *C6A, real *C6B,
                         real *sigmaA, real *sigmaB,
                         real *sigma3A, real *sigma3B,
+                        gmx_bool bHaveChargeOrTypePerturbed,
                         gmx_bool calc_excl_corr,
                         t_blocka *excl, rvec x[],
                         matrix box, rvec mu_tot[],
@@ -156,7 +167,6 @@ void ewald_LRcorrection(int start, int end,
     tensor      dxdf_q, dxdf_lj;
     real        vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
     real        L1_q, L1_lj, dipole_coeff, qqA, qqB, qqL, vr0_q, vr0_lj = 0;
-    gmx_bool    bFreeEnergy  = (chargeB != NULL);
     gmx_bool    bMolPBC      = fr->bMolPBC;
     gmx_bool    bDoingLBRule = (fr->ljpme_combination_rule == eljpmeLB);
 
@@ -226,7 +236,7 @@ void ewald_LRcorrection(int start, int end,
     {
         clear_mat(dxdf_lj);
     }
-    if ((calc_excl_corr || dipole_coeff != 0) && !bFreeEnergy)
+    if ((calc_excl_corr || dipole_coeff != 0) && !bHaveChargeOrTypePerturbed)
     {
         for (i = start; (i < end); i++)
         {
@@ -553,7 +563,7 @@ void ewald_LRcorrection(int start, int end,
     /* Global corrections only on master process */
     if (MASTER(cr) && thread == 0)
     {
-        for (q = 0; q < (bFreeEnergy ? 2 : 1); q++)
+        for (q = 0; q < (bHaveChargeOrTypePerturbed ? 2 : 1); q++)
         {
             if (calc_excl_corr)
             {
@@ -581,7 +591,7 @@ void ewald_LRcorrection(int start, int end,
             }
         }
     }
-    if (!bFreeEnergy)
+    if (!bHaveChargeOrTypePerturbed)
     {
         *Vcorr_q = Vdipole[0] - Vself_q[0] - Vexcl_q;
         if (EVDW_PME(fr->vdwtype))

diff --git a/src/gromacs/legacyheaders/coulomb.h b/src/gromacs/legacyheaders/coulomb.h
index 7d086ddb4244baef3b06225cbcd9bdc2a425ac9a..089b2d20fd2bcacf66011f778cd109a28d26dffb 100644 (file)
--- a/src/gromacs/legacyheaders/coulomb.h
+++ b/src/gromacs/legacyheaders/coulomb.h
@@ -77,6 +77,7 @@ ewald_LRcorrection(int start, int end,
                    real *C6A, real *C6B,
                    real *sigmaA, real *sigmaB,
                    real *sigma3A, real *sigma3B,
+                   gmx_bool bHaveChargeOrTypePerturbed,
                    gmx_bool calc_excl_corr,
                    t_blocka *excl, rvec x[],
                    matrix box, rvec mu_tot[],

diff --git a/src/gromacs/mdlib/force.c b/src/gromacs/mdlib/force.c
index d4e6445491e36fe2f10bff1cf360d7b6d8ae8050..5230983cbe17802a374739a500948b343d4af8ba 100644 (file)
--- a/src/gromacs/mdlib/force.c
+++ b/src/gromacs/mdlib/force.c
@@ -541,14 +541,11 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
 
                     ewald_LRcorrection(fr->excl_load[t], fr->excl_load[t+1],
                                        cr, t, fr,
-                                       md->chargeA,
-                                       md->nChargePerturbed ? md->chargeB : NULL,
-                                       md->sqrt_c6A,
-                                       md->nTypePerturbed ? md->sqrt_c6B : NULL,
-                                       md->sigmaA,
-                                       md->nTypePerturbed ? md->sigmaB : NULL,
-                                       md->sigma3A,
-                                       md->nTypePerturbed ? md->sigma3B : NULL,
+                                       md->chargeA, md->chargeB,
+                                       md->sqrt_c6A, md->sqrt_c6B,
+                                       md->sigmaA, md->sigmaB,
+                                       md->sigma3A, md->sigma3B,
+                                       md->nChargePerturbed || md->nTypePerturbed,
                                        ir->cutoff_scheme != ecutsVERLET,
                                        excl, x, bSB ? boxs : box, mu_tot,
                                        ir->ewald_geometry,

diff --git a/src/gromacs/mdlib/mdatom.c b/src/gromacs/mdlib/mdatom.c
index b8d51a76da80566855c2088a159bc843d3a78e4e..e4b698a018a9f5a000109c8b7fd25eddc816af4e 100644 (file)
--- a/src/gromacs/mdlib/mdatom.c
+++ b/src/gromacs/mdlib/mdatom.c
@@ -156,27 +156,24 @@ void atoms2md(gmx_mtop_t *mtop, t_inputrec *ir,
         srenew(md->massT, md->nalloc);
         srenew(md->invmass, md->nalloc);
         srenew(md->chargeA, md->nalloc);
+        srenew(md->typeA, md->nalloc);
+        if (md->nPerturbed)
+        {
+            srenew(md->chargeB, md->nalloc);
+            srenew(md->typeB, md->nalloc);
+        }
         if (bLJPME)
         {
             srenew(md->sqrt_c6A, md->nalloc);
             srenew(md->sigmaA, md->nalloc);
             srenew(md->sigma3A, md->nalloc);
-        }
-        if (md->nPerturbed)
-        {
-            srenew(md->chargeB, md->nalloc);
-            if (bLJPME)
+            if (md->nPerturbed)
             {
                 srenew(md->sqrt_c6B, md->nalloc);
                 srenew(md->sigmaB, md->nalloc);
                 srenew(md->sigma3B, md->nalloc);
             }
         }
-        srenew(md->typeA, md->nalloc);
-        if (md->nPerturbed)
-        {
-            srenew(md->typeB, md->nalloc);
-        }
         srenew(md->ptype, md->nalloc);
         if (opts->ngtc > 1)
         {