Fixed exit with FEP with PME nodes without LJ-PME

With FEP runs with separate PME nodes without LJ-PME, mdrun would
exit with a gmx_incons.
Reorganized the PME some flags and conditionals for clarity.

Change-Id: I04513a2d84682617b36c0c8ad4afe632ce3f7ee1

diff --git a/src/gromacs/legacyheaders/pme.h b/src/gromacs/legacyheaders/pme.h
index b30678c2faf84ed8b9273b27537877fd471ed0ac..4159bf94eb05eac7d5f0f8795839af845bc012eb 100644 (file)
--- a/src/gromacs/legacyheaders/pme.h
+++ b/src/gromacs/legacyheaders/pme.h
@@ -87,7 +87,6 @@ int gmx_pme_destroy(FILE *log, gmx_pme_t *pmedata);
  */
 #define GMX_PME_DO_COULOMB    (1<<13)
 #define GMX_PME_DO_LJ         (1<<14)
-#define GMX_PME_LJ_LB         (1<<15)
 
 #define GMX_PME_DO_ALL_F  (GMX_PME_SPREAD_Q | GMX_PME_SOLVE | GMX_PME_CALC_F)
 

diff --git a/src/gromacs/mdlib/force.c b/src/gromacs/mdlib/force.c
index 3001236712a22d7217ade74c72de5f3ff0274538..6aa988000370bd9ae58eb2aff31dd3d2e3d0d1b3 100644 (file)
--- a/src/gromacs/mdlib/force.c
+++ b/src/gromacs/mdlib/force.c
@@ -595,11 +595,6 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
                     if (EVDW_PME(fr->vdwtype))
                     {
                         pme_flags |= GMX_PME_DO_LJ;
-                        if (fr->ljpme_combination_rule == eljpmeLB)
-                        {
-                            /*Lorentz-Berthelot Comb. Rules in LJ-PME*/
-                            pme_flags |= GMX_PME_LJ_LB;
-                        }
                     }
                     if (flags & GMX_FORCE_FORCES)
                     {

diff --git a/src/gromacs/mdlib/pme.c b/src/gromacs/mdlib/pme.c
index 0a9faf23a27ba32cb1c25eed856cca42508c2a8f..56478f1bb7347c10f2705ae5dfa556ebdc04b13e 100644 (file)
--- a/src/gromacs/mdlib/pme.c
+++ b/src/gromacs/mdlib/pme.c
@@ -304,14 +304,19 @@ typedef struct gmx_pme {
     int        nthread;       /* The number of threads doing PME on our rank */
 
     gmx_bool   bPPnode;       /* Node also does particle-particle forces */
+
     gmx_bool   bFEP;          /* Compute Free energy contribution */
     gmx_bool   bFEP_q;
     gmx_bool   bFEP_lj;
+
     int        nkx, nky, nkz; /* Grid dimensions */
+
     gmx_bool   bP3M;          /* Do P3M: optimize the influence function */
     int        pme_order;
     real       epsilon_r;
 
+    int        ljpme_combination_rule;  /* Type of combination rule in LJ-PME */
+
     int        ngrids;                  /* number of grids we maintain for pmegrid, (c)fftgrid and pfft_setups*/
 
     pmegrids_t pmegrid[DO_Q_AND_LJ_LB]; /* Grids on which we do spreading/interpolation,
@@ -3502,8 +3507,13 @@ int gmx_pme_init(gmx_pme_t *         pmedata,
     pme->nkz         = ir->nkz;
     pme->bP3M        = (ir->coulombtype == eelP3M_AD || getenv("GMX_PME_P3M") != NULL);
     pme->pme_order   = ir->pme_order;
+
+    /* Always constant electrostatics parameters */
     pme->epsilon_r   = ir->epsilon_r;
 
+    /* Always constant LJ parameters */
+    pme->ljpme_combination_rule = ir->ljpme_combination_rule;
+
     /* If we violate restrictions, generate a fatal error here */
     gmx_pme_check_restrictions(pme->pme_order,
                                pme->nkx, pme->nky, pme->nkz,
@@ -3614,14 +3624,28 @@ int gmx_pme_init(gmx_pme_t *         pmedata,
     ndata[0]    = pme->nkx;
     ndata[1]    = pme->nky;
     ndata[2]    = pme->nkz;
-    pme->ngrids = ((ir->ljpme_combination_rule == eljpmeLB) ? DO_Q_AND_LJ_LB : DO_Q_AND_LJ);
+    /* It doesn't matter if we allocate too many grids here,
+     * we only allocate and use the ones we need.
+     */
+    if (EVDW_PME(ir->vdwtype))
+    {
+        pme->ngrids = ((ir->ljpme_combination_rule == eljpmeLB) ? DO_Q_AND_LJ_LB : DO_Q_AND_LJ);
+    }
+    else
+    {
+        pme->ngrids = DO_Q;
+    }
     snew(pme->fftgrid, pme->ngrids);
     snew(pme->cfftgrid, pme->ngrids);
     snew(pme->pfft_setup, pme->ngrids);
 
     for (i = 0; i < pme->ngrids; ++i)
     {
-        if (((ir->ljpme_combination_rule == eljpmeLB) && i >= 2) || i % 2 == 0 || bFreeEnergy_q || bFreeEnergy_lj)
+        if ((i <  DO_Q && EEL_PME(ir->coulombtype) && (i == 0 ||
+                                                       bFreeEnergy_q))|| 
+            (i >= DO_Q && EVDW_PME(ir->vdwtype) && (i == 2 ||
+                                                    bFreeEnergy_lj ||
+                                                    ir->ljpme_combination_rule == eljpmeLB)))
         {
             pmegrids_init(&pme->pmegrid[i],
                           pme->pmegrid_nx, pme->pmegrid_ny, pme->pmegrid_nz,
@@ -4388,7 +4412,7 @@ void gmx_pme_calc_energy(gmx_pme_t pme, int n, rvec *x, real *q, real *V)
     {
         gmx_incons("gmx_pme_calc_energy called in parallel");
     }
-    if (pme->bFEP > 1)
+    if (pme->bFEP_q > 1)
     {
         gmx_incons("gmx_pme_calc_energy with free energy");
     }
@@ -4748,7 +4772,7 @@ int gmx_pme_do(gmx_pme_t pme,
      */
 
     /* If we are doing LJ-PME with LB, we only do Q here */
-    qmax = (flags & GMX_PME_LJ_LB) ? DO_Q : DO_Q_AND_LJ;
+    qmax = (pme->ljpme_combination_rule == eljpmeLB) ? DO_Q : DO_Q_AND_LJ;
 
     for (q = 0; q < qmax; ++q)
     {
@@ -4757,10 +4781,10 @@ int gmx_pme_do(gmx_pme_t pme,
          * If q < 2 we should be doing electrostatic PME
          * If q >= 2 we should be doing LJ-PME
          */
-        if ((!pme->bFEP_q && q == 1)
-            || (!pme->bFEP_lj && q == 3)
-            || (!(flags & GMX_PME_DO_COULOMB) && q < 2)
-            || (!(flags & GMX_PME_DO_LJ) && q >= 2))
+        if ((q <  DO_Q && (!(flags & GMX_PME_DO_COULOMB) ||
+                           (q == 1 && !pme->bFEP_q))) ||
+            (q >= DO_Q && (!(flags & GMX_PME_DO_LJ) ||
+                           (q == 3 && !pme->bFEP_lj))))
         {
             continue;
         }
@@ -4990,7 +5014,7 @@ int gmx_pme_do(gmx_pme_t pme,
     /* For Lorentz-Berthelot combination rules in LJ-PME, we need to calculate
      * seven terms. */
 
-    if (flags & GMX_PME_LJ_LB)
+    if ((flags & GMX_PME_DO_LJ) && pme->ljpme_combination_rule == eljpmeLB)
     {
         real *local_c6 = NULL, *local_sigma = NULL;
         if (pme->nnodes == 1)
@@ -5212,7 +5236,7 @@ int gmx_pme_do(gmx_pme_t pme,
                 bFirst = FALSE;
             } /*for (q = 8; q >= 2; --q)*/
         }     /*if (bCalcF)*/
-    }         /*if (flags & GMX_PME_LJ_LB)*/
+    }         /*if ((flags & GMX_PME_DO_LJ) && pme->ljpme_combination_rule == eljpmeLB) */
 
     if (bCalcF && pme->nnodes > 1)
     {

diff --git a/src/gromacs/mdlib/pme_pp.c b/src/gromacs/mdlib/pme_pp.c
index 773ad76f8ea432a1b7ac04c5d874c8df12302654..b375abcef0608beb82b648c13ee94f31e009bb70 100644 (file)
--- a/src/gromacs/mdlib/pme_pp.c
+++ b/src/gromacs/mdlib/pme_pp.c
@@ -582,8 +582,10 @@ int gmx_pme_recv_params_coords(struct gmx_pme_pp          *pme_pp,
             /* The box, FE flag and lambda are sent along with the coordinates
              *  */
             copy_mat(cnb.box, box);
-            *bFreeEnergy_q  = (cnb.flags & PP_PME_FEP_Q);
-            *bFreeEnergy_lj = (cnb.flags & PP_PME_FEP_LJ);
+            *bFreeEnergy_q  = ((cnb.flags & GMX_PME_DO_COULOMB) &&
+                               (cnb.flags & PP_PME_FEP_Q));
+            *bFreeEnergy_lj = ((cnb.flags & GMX_PME_DO_LJ) &&
+                               (cnb.flags & PP_PME_FEP_LJ));
             *lambda_q       = cnb.lambda_q;
             *lambda_lj      = cnb.lambda_lj;
             *bEnerVir       = (cnb.flags & PP_PME_ENER_VIR);

diff --git a/src/gromacs/mdlib/sim_util.c b/src/gromacs/mdlib/sim_util.c
index 52ab9568fd4eb5002f0cf18404af51af047c231e..d8ad3515540f1e4770c62b38951ca83fdc574fbd 100644 (file)
--- a/src/gromacs/mdlib/sim_util.c
+++ b/src/gromacs/mdlib/sim_util.c
@@ -934,10 +934,6 @@ void do_force_cutsVERLET(FILE *fplog, t_commrec *cr,
         if (EVDW_PME(fr->vdwtype))
         {
             pme_flags |= GMX_PME_DO_LJ;
-            if (fr->ljpme_combination_rule == eljpmeLB)
-            {
-                pme_flags |= GMX_PME_LJ_LB;
-            }
         }
 
         gmx_pme_send_coordinates(cr, bBS ? boxs : box, x,
@@ -1704,10 +1700,6 @@ void do_force_cutsGROUP(FILE *fplog, t_commrec *cr,
         if (EVDW_PME(fr->vdwtype))
         {
             pme_flags |= GMX_PME_DO_LJ;
-            if (fr->ljpme_combination_rule == eljpmeLB)
-            {
-                pme_flags |= GMX_PME_LJ_LB;
-            }
         }
 
         gmx_pme_send_coordinates(cr, bBS ? boxs : box, x,