Reintroduce LJ Ewald tables

Fixes #3470 by reintroducing the LJ Ewald tables.

diff --git a/docs/release-notes/2020/2020.2.rst b/docs/release-notes/2020/2020.2.rst
index 2b2097684c6ec931afc4f5c8e0e8f393f7930728..341dfd445ce1235588e4e486f4339b699e6dcc8b 100644 (file)
--- a/docs/release-notes/2020/2020.2.rst
+++ b/docs/release-notes/2020/2020.2.rst
@@ -35,6 +35,15 @@ their original lambda state.
 
 :issue:`3465`
 
+Fixed free energy calculations with LJ PME
+""""""""""""""""""""""""""""""""""""""""""
+
+Fixed an issue that calculated wrong long-range corrections when using
+free energy perturbation with ``vdwtype = pme``. This affected forces,
+energies, lambda derivatives and foreign lambdas.
+
+:issue:`3470`
+
 Fixes for ``gmx`` tools
 ^^^^^^^^^^^^^^^^^^^^^^^
 

diff --git a/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp b/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp
index 8395379dad3f16295d38fe2df94a22242d05126c..19e5c660dd26bf9e9eaa0a07b502f15d26b376d0 100644 (file)
--- a/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp
+++ b/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017,2018,2019,2020, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -324,21 +324,32 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist,
     real rcutoff_max2 = std::max(ic->rcoulomb, ic->rvdw);
     rcutoff_max2      = rcutoff_max2 * rcutoff_max2;
 
-    const real* tab_ewald_F_lj = nullptr;
-    const real* tab_ewald_V_lj = nullptr;
-    const real* ewtab          = nullptr;
-    real        ewtabscale     = 0;
-    real        ewtabhalfspace = 0;
-    real        sh_ewald       = 0;
+    const real* tab_ewald_F_lj           = nullptr;
+    const real* tab_ewald_V_lj           = nullptr;
+    const real* ewtab                    = nullptr;
+    real        coulombTableScale        = 0;
+    real        coulombTableScaleInvHalf = 0;
+    real        vdwTableScale            = 0;
+    real        vdwTableScaleInvHalf     = 0;
+    real        sh_ewald                 = 0;
     if (elecInteractionTypeIsEwald || vdwInteractionTypeIsEwald)
     {
-        const auto& tables = *ic->coulombEwaldTables;
-        sh_ewald           = ic->sh_ewald;
-        ewtab              = tables.tableFDV0.data();
-        ewtabscale         = tables.scale;
-        ewtabhalfspace     = half / ewtabscale;
-        tab_ewald_F_lj     = tables.tableF.data();
-        tab_ewald_V_lj     = tables.tableV.data();
+        sh_ewald = ic->sh_ewald;
+    }
+    if (elecInteractionTypeIsEwald)
+    {
+        const auto& coulombTables = *ic->coulombEwaldTables;
+        ewtab                     = coulombTables.tableFDV0.data();
+        coulombTableScale         = coulombTables.scale;
+        coulombTableScaleInvHalf  = half / coulombTableScale;
+    }
+    if (vdwInteractionTypeIsEwald)
+    {
+        const auto& vdwTables = *ic->vdwEwaldTables;
+        tab_ewald_F_lj        = vdwTables.tableF.data();
+        tab_ewald_V_lj        = vdwTables.tableV.data();
+        vdwTableScale         = vdwTables.scale;
+        vdwTableScaleInvHalf  = half / vdwTableScale;
     }
 
     /* For Ewald/PME interactions we cannot easily apply the soft-core component to
@@ -710,12 +721,12 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist,
                  */
                 real v_lr, f_lr;
 
-                const real ewrt   = r * ewtabscale;
+                const real ewrt   = r * coulombTableScale;
                 int        ewitab = static_cast<int>(ewrt);
                 const real eweps  = ewrt - ewitab;
                 ewitab            = 4 * ewitab;
                 f_lr              = ewtab[ewitab] + eweps * ewtab[ewitab + 1];
-                v_lr = (ewtab[ewitab + 2] - ewtabhalfspace * eweps * (ewtab[ewitab] + f_lr));
+                v_lr = (ewtab[ewitab + 2] - coulombTableScaleInvHalf * eweps * (ewtab[ewitab] + f_lr));
                 f_lr *= rinv;
 
                 /* Note that any possible Ewald shift has already been applied in
@@ -755,7 +766,7 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist,
                  * r close to 0 for non-interacting pairs.
                  */
 
-                const real rs   = rsq * rinv * ewtabscale;
+                const real rs   = rsq * rinv * vdwTableScale;
                 const int  ri   = static_cast<int>(rs);
                 const real frac = rs - ri;
                 const real f_lr = (1 - frac) * tab_ewald_F_lj[ri] + frac * tab_ewald_F_lj[ri + 1];
@@ -763,7 +774,8 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist,
                  * the factor 1/6, we should add this.
                  */
                 const real FF = f_lr * rinv / six;
-                real VV = (tab_ewald_V_lj[ri] - ewtabhalfspace * frac * (tab_ewald_F_lj[ri] + f_lr)) / six;
+                real VV = (tab_ewald_V_lj[ri] - vdwTableScaleInvHalf * frac * (tab_ewald_F_lj[ri] + f_lr))
+                          / six;
 
                 if (ii == jnr)
                 {

diff --git a/src/gromacs/mdlib/forcerec.cpp b/src/gromacs/mdlib/forcerec.cpp
index 336162381160f0dbefd8beb21421193061b182e9..c7aeca1a8b7b1379716bda636cc67bf04cfb697d 100644 (file)
--- a/src/gromacs/mdlib/forcerec.cpp
+++ b/src/gromacs/mdlib/forcerec.cpp
@@ -760,12 +760,12 @@ static void init_ewald_f_table(const interaction_const_t& ic,
 
 void init_interaction_const_tables(FILE* fp, interaction_const_t* ic)
 {
-    if (EEL_PME_EWALD(ic->eeltype))
+    if (EEL_PME_EWALD(ic->eeltype) || EVDW_PME(ic->vdwtype))
     {
-        init_ewald_f_table(*ic, ic->coulombEwaldTables.get(), nullptr);
+        init_ewald_f_table(*ic, ic->coulombEwaldTables.get(), ic->vdwEwaldTables.get());
         if (fp != nullptr)
         {
-            fprintf(fp, "Initialized non-bonded Coulomb Ewald tables, spacing: %.2e size: %zu\n\n",
+            fprintf(fp, "Initialized non-bonded Ewald tables, spacing: %.2e size: %zu\n\n",
                     1 / ic->coulombEwaldTables->scale, ic->coulombEwaldTables->tableF.size());
         }
     }
@@ -826,6 +826,7 @@ static void init_interaction_const(FILE*                 fp,
     ic->cutoff_scheme = ir->cutoff_scheme;
 
     ic->coulombEwaldTables = std::make_unique<EwaldCorrectionTables>();
+    ic->vdwEwaldTables     = std::make_unique<EwaldCorrectionTables>();
 
     /* Lennard-Jones */
     ic->vdwtype         = ir->vdwtype;

diff --git a/src/gromacs/mdtypes/interaction_const.h b/src/gromacs/mdtypes/interaction_const.h
index 90b77649aea15c402d27eeefdbf95b0c0423098e..a034e7dd781ecefb99f7b987720608e4a65b2d9e 100644 (file)
--- a/src/gromacs/mdtypes/interaction_const.h
+++ b/src/gromacs/mdtypes/interaction_const.h
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2012,2013,2014,2015,2017,2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2012,2013,2014,2015,2017,2018,2019,2020, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -143,9 +143,8 @@ struct interaction_const_t
 
     // Coulomb Ewald correction table
     std::unique_ptr<EwaldCorrectionTables> coulombEwaldTables;
-    /* Note that a Van der Waals Ewald correction table
-     * of type EwaldCorrectionTables can be added here if wanted.
-     */
+    // Van der Waals Ewald correction table
+    std::unique_ptr<EwaldCorrectionTables> vdwEwaldTables;
 };
 
 #endif