Fixed free energy with LJ-PME

First, Only the normal LJ part of the potential was
shifted for LJ-PME with potential-shift
modifier in the free energy kernel,
and not the Ewald component. This had an effect
of roughly 0.05% on the total energy. Second,
Berk found that the verlet nbnxn SIMD LJ-PME kernel
had an issue that caused the potential and force
to be slightly off. With this patch, both
potential and force match the non-free-energy
kernel perfectly with at lambda=0.

Change-Id: I4bc53ce1bdb96ba06656b86470dbf8dbd1e81972

diff --git a/src/gromacs/gmxlib/ewald_util.c b/src/gromacs/gmxlib/ewald_util.c
index b6af0d607b11b4d92bb2851d7ddc6183be666d0f..d3e5f23748e36f07df9c21c073ae01b182fbb41a 100644 (file)
--- a/src/gromacs/gmxlib/ewald_util.c
+++ b/src/gromacs/gmxlib/ewald_util.c
@@ -226,7 +226,7 @@ void ewald_LRcorrection(int start, int end,
     {
         clear_mat(dxdf_lj);
     }
-    if ((calc_excl_corr || dipole_coeff != 0 || EVDW_PME(fr->vdwtype)) && !bFreeEnergy)
+    if ((calc_excl_corr || dipole_coeff != 0) && !bFreeEnergy)
     {
         for (i = start; (i < end); i++)
         {
@@ -387,7 +387,7 @@ void ewald_LRcorrection(int start, int end,
             }
         }
     }
-    else if (calc_excl_corr || dipole_coeff != 0 || EVDW_PME(fr->vdwtype))
+    else if (calc_excl_corr || dipole_coeff != 0)
     {
         for (i = start; (i < end); i++)
         {

diff --git a/src/gromacs/gmxlib/nonbonded/nb_free_energy.c b/src/gromacs/gmxlib/nonbonded/nb_free_energy.c
index 058b2a87eb2df4fae11ecc0e39fbded877bb2803..2bf8f90e3f8bd8bc60baf5462c60e5cf98e73781 100644 (file)
--- a/src/gromacs/gmxlib/nonbonded/nb_free_energy.c
+++ b/src/gromacs/gmxlib/nonbonded/nb_free_energy.c
@@ -106,6 +106,8 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict    nlist,
     const real *  chargeB;
     real          sigma6_min, sigma6_def, lam_power, sc_power, sc_r_power;
     real          alpha_coul, alpha_vdw, lambda_coul, lambda_vdw, ewc_lj;
+    real          ewcljrsq, ewclj, ewclj2, exponent, poly, vvdw_disp, vvdw_rep, sh_lj_ewald;
+    real          ewclj6;
     const real *  nbfp, *nbfp_grid;
     real *        dvdl;
     real *        Vv;
@@ -178,6 +180,10 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict    nlist,
     sh_ewald            = fr->ic->sh_ewald;
     rvdw                = fr->rvdw;
     sh_invrc6           = fr->ic->sh_invrc6;
+    sh_lj_ewald         = fr->ic->sh_lj_ewald;
+    ewclj               = fr->ewaldcoeff_lj;
+    ewclj2              = ewclj*ewclj;
+    ewclj6              = ewclj2*ewclj2*ewclj2;
 
     if (fr->coulomb_modifier == eintmodPOTSWITCH)
     {
@@ -273,6 +279,14 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict    nlist,
      */
     bConvertLJEwaldToLJ6   = (bEwaldLJ && (fr->vdw_modifier   != eintmodPOTSWITCH));
 
+    /* We currently don't implement exclusion correction, needed with the Verlet cut-off scheme, without conversion */
+    if (fr->cutoff_scheme == ecutsVERLET &&
+        ((bEwald   && !bConvertEwaldToCoulomb) ||
+         (bEwaldLJ && !bConvertLJEwaldToLJ6)))
+    {
+        gmx_incons("Unimplemented non-bonded setup");
+    }
+
     /* fix compiler warnings */
     nj1   = 0;
     n1C   = n1V   = 0;
@@ -514,7 +528,7 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict    nlist,
                                     Vcoul[i]   = qq[i]*rinvC;
                                     FscalC[i]  = Vcoul[i];
                                     /* The shift for the Coulomb potential is stored in
-                                     * the RF parameter c_rf, which is 0 without shift
+                                     * the RF parameter c_rf, which is 0 without shift.
                                      */
                                     Vcoul[i]  -= qq[i]*fr->ic->c_rf;
                                     break;
@@ -602,7 +616,6 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict    nlist,
                             switch (ivdw)
                             {
                                 case GMX_NBKERNEL_VDW_LENNARDJONES:
-                                case GMX_NBKERNEL_VDW_LJEWALD:
                                     /* cutoff LJ */
                                     if (sc_r_power == 6.0)
                                     {
@@ -610,19 +623,14 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict    nlist,
                                     }
                                     else
                                     {
-                                        rinv6            = pow(rinvV, 6.0);
+                                        rinv6            = rinvV*rinvV;
+                                        rinv6            = rinv6*rinv6*rinv6;
                                     }
                                     Vvdw6            = c6[i]*rinv6;
                                     Vvdw12           = c12[i]*rinv6*rinv6;
-                                    if (fr->vdw_modifier == eintmodPOTSHIFT)
-                                    {
-                                        Vvdw[i]          = ( (Vvdw12-c12[i]*sh_invrc6*sh_invrc6)*(1.0/12.0)
-                                                             -(Vvdw6-c6[i]*sh_invrc6)*(1.0/6.0));
-                                    }
-                                    else
-                                    {
-                                        Vvdw[i]          = Vvdw12*(1.0/12.0) - Vvdw6*(1.0/6.0);
-                                    }
+
+                                    Vvdw[i]          = ( (Vvdw12 - c12[i]*sh_invrc6*sh_invrc6)*(1.0/12.0)
+                                                         - (Vvdw6 - c6[i]*sh_invrc6)*(1.0/6.0));
                                     FscalV[i]        = Vvdw12 - Vvdw6;
                                     break;
 
@@ -656,6 +664,41 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict    nlist,
                                     FscalV[i] -= c12[i]*tabscale*FF*rV;
                                     break;
 
+                                case GMX_NBKERNEL_VDW_LJEWALD:
+                                    if (sc_r_power == 6.0)
+                                    {
+                                        rinv6            = rpinvV;
+                                    }
+                                    else
+                                    {
+                                        rinv6            = rinvV*rinvV;
+                                        rinv6            = rinv6*rinv6*rinv6;
+                                    }
+                                    c6grid           = nbfp_grid[tj[i]];
+
+                                    if (bConvertLJEwaldToLJ6)
+                                    {
+                                        /* cutoff LJ */
+                                        Vvdw6            = c6[i]*rinv6;
+                                        Vvdw12           = c12[i]*rinv6*rinv6;
+
+                                        Vvdw[i]          = ( (Vvdw12 - c12[i]*sh_invrc6*sh_invrc6)*(1.0/12.0)
+                                                             - (Vvdw6 - c6[i]*sh_invrc6 - c6grid*sh_lj_ewald)*(1.0/6.0));
+                                        FscalV[i]        = Vvdw12 - Vvdw6;
+                                    }
+                                    else
+                                    {
+                                        /* Normal LJ-PME */
+                                        ewcljrsq         = ewclj2*rV*rV;
+                                        exponent         = exp(-ewcljrsq);
+                                        poly             = exponent*(1.0 + ewcljrsq + ewcljrsq*ewcljrsq*0.5);
+                                        vvdw_disp        = (c6[i]-c6grid*(1.0-poly))*rinv6;
+                                        vvdw_rep         = c12[i]*rinv6*rinv6;
+                                        FscalV[i]        = vvdw_rep - vvdw_disp - c6grid*(1.0/6.0)*exponent*ewclj6;
+                                        Vvdw[i]          = (vvdw_rep - c12[i]*sh_invrc6*sh_invrc6)/12.0 - (vvdw_disp - c6[i]*sh_invrc6 - c6grid*sh_lj_ewald)/6.0;
+                                    }
+                                    break;
+
                                 case GMX_NBKERNEL_VDW_NONE:
                                     Vvdw[i]    = 0.0;
                                     FscalV[i]  = 0.0;
@@ -748,6 +791,10 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict    nlist,
                 v_lr      = (ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+f_lr));
                 f_lr     *= rinv;
 
+                /* Note that any possible Ewald shift has already been applied in
+                 * the normal interaction part above.
+                 */
+
                 if (ii == jnr)
                 {
                     /* If we get here, the i particle (ii) has itself (jnr)
@@ -776,6 +823,10 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict    nlist,
                  * the softcore to the entire VdW interaction,
                  * including the reciprocal-space component.
                  */
+                /* We could also use the analytical form here
+                 * iso a table, but that can cause issues for
+                 * r close to 0 for non-interacting pairs.
+                 */
                 real rs, frac, f_lr;
                 int  ri;
 
@@ -783,8 +834,11 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict    nlist,
                 ri     = (int)rs;
                 frac   = rs - ri;
                 f_lr   = (1 - frac)*tab_ewald_F_lj[ri] + frac*tab_ewald_F_lj[ri+1];
-                FF     = f_lr*rinv;
-                VV     = tab_ewald_V_lj[ri] - ewtabhalfspace*frac*(tab_ewald_F_lj[ri] + f_lr);
+                /* TODO: Currently the Ewald LJ table does not contain
+                 * the factor 1/6, we should add this.
+                 */
+                FF     = f_lr*rinv/6.0;
+                VV     = (tab_ewald_V_lj[ri] - ewtabhalfspace*frac*(tab_ewald_F_lj[ri] + f_lr))/6.0;
 
                 if (ii == jnr)
                 {
@@ -799,11 +853,10 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict    nlist,
                 for (i = 0; i < NSTATES; i++)
                 {
                     c6grid      = nbfp_grid[tj[i]];
-                    vvtot      += LFV[i]*c6grid*VV*(1.0/6.0);
-                    Fscal      += LFV[i]*c6grid*FF*(1.0/6.0);
-                    dvdl_vdw   += (DLF[i]*c6grid)*VV*(1.0/6.0);
+                    vvtot      += LFV[i]*c6grid*VV;
+                    Fscal      += LFV[i]*c6grid*FF;
+                    dvdl_vdw   += (DLF[i]*c6grid)*VV;
                 }
-
             }
 
             if (bDoForces)

diff --git a/src/gromacs/mdlib/nbnxn_atomdata.c b/src/gromacs/mdlib/nbnxn_atomdata.c
index 68365412121494b5daeab85794965ca63608977c..d271a8ce7a165392094f12f730b97ce1d9b661e1 100644 (file)
--- a/src/gromacs/mdlib/nbnxn_atomdata.c
+++ b/src/gromacs/mdlib/nbnxn_atomdata.c
@@ -851,7 +851,7 @@ static void copy_lj_to_nbat_lj_comb_x8(const real *ljparam_type,
     }
 }
 
-/* Sets the atom type and LJ data in nbnxn_atomdata_t */
+/* Sets the atom type in nbnxn_atomdata_t */
 static void nbnxn_atomdata_set_atomtypes(nbnxn_atomdata_t    *nbat,
                                          int                  ngrid,
                                          const nbnxn_search_t nbs,
@@ -872,9 +872,30 @@ static void nbnxn_atomdata_set_atomtypes(nbnxn_atomdata_t    *nbat,
 
             copy_int_to_nbat_int(nbs->a+ash, grid->cxy_na[i], ncz*grid->na_sc,
                                  type, nbat->ntype-1, nbat->type+ash);
+        }
+    }
+}
+
+/* Sets the LJ combination rule parameters in nbnxn_atomdata_t */
+static void nbnxn_atomdata_set_ljcombparams(nbnxn_atomdata_t    *nbat,
+                                            int                  ngrid,
+                                            const nbnxn_search_t nbs)
+{
+    int                 g, i, ncz, ash;
+    const nbnxn_grid_t *grid;
 
-            if (nbat->comb_rule != ljcrNONE)
+    if (nbat->comb_rule != ljcrNONE)
+    {
+        for (g = 0; g < ngrid; g++)
+        {
+            grid = &nbs->grid[g];
+
+            /* Loop over all columns and copy and fill */
+            for (i = 0; i < grid->ncx*grid->ncy; i++)
             {
+                ncz = grid->cxy_ind[i+1] - grid->cxy_ind[i];
+                ash = (grid->cell0 + grid->cxy_ind[i])*grid->na_sc;
+
                 if (nbat->XFormat == nbatX4)
                 {
                     copy_lj_to_nbat_lj_comb_x4(nbat->nbfp_comb,
@@ -1096,6 +1117,9 @@ void nbnxn_atomdata_set(nbnxn_atomdata_t    *nbat,
         nbnxn_atomdata_mask_fep(nbat, ngrid, nbs);
     }
 
+    /* This must be done after masking types for FEP */
+    nbnxn_atomdata_set_ljcombparams(nbat, ngrid, nbs);
+
     nbnxn_atomdata_set_energygroups(nbat, ngrid, nbs, atinfo);
 }
 

diff --git a/src/gromacs/mdlib/nbnxn_search.c b/src/gromacs/mdlib/nbnxn_search.c
index c06a5d10f1701d3d9981617c23622f525a44a1d2..ae8b82caf72abdd63fae0e4288d58d1a173ed23c 100644 (file)
--- a/src/gromacs/mdlib/nbnxn_search.c
+++ b/src/gromacs/mdlib/nbnxn_search.c
@@ -3440,7 +3440,6 @@ static void make_fep_list(const nbnxn_search_t    nbs,
                              * Note that the charge has been set to zero,
                              * but we need to avoid 0/0, as perturbed atoms
                              * can be on top of each other.
-                             * (and the LJ parameters have not been zeroed)
                              */
                             nbl->cj[cj_ind].excl &= ~(1U << (i*nbl->na_cj + j));
                         }