Non-P3M Coulomb methods maintenance

doc/doxygen/bibliography.bib

@@ -567,6 +567,17 @@ @Article{thole81a
   doi       = {10.1016/0301-0104(81)85176-2},
 }
 
+@Article{tironi95a,
+  author    = {Tironi, Ilario G. and Sperb, Ren\'{e} and Smith, Paul E. and van Gunsteren, Wilfred F.},
+  title     = {A generalized Reaction Field Method for Molecular-Dynamics Simulations},
+  journal   = {The Journal of Chemical Physics},
+  year      = {1995},
+  volume    = {102},
+  number    = {13},
+  pages     = {5451--5459},
+  doi       = {10.1063/1.469273},
+}
+
 @Article{troester05a,
   title = {{Wang-Landau sampling with self-adaptive range}},
   author = {Tr\"{o}ster, Andreas and Dellago, Christoph},

doc/sphinx/electrostatics.rst

@@ -138,7 +138,8 @@ The Debye-Hückel electrostatic potential is defined by
 
   .. math:: U^{C-DH} = C \cdot \frac{q_1 q_2 \exp(-\kappa r)}{r}\quad \mathrm{for}\quad r<r_{\mathrm{cut}}
 
-where :math:`C` is defined as in Eqn. :eq:`coulomb_prefactor`.
+where :math:`C` is defined as in Eqn. :eq:`coulomb_prefactor` and
+:math:`\kappa` is the ionic strength.
 The Debye-Hückel potential is an approximate method for calculating
 electrostatic interactions, but technically it is treated as other
 short-ranged non-bonding potentials. For :math:`r > r_{\textrm{cut}}` it is
@@ -147,6 +148,29 @@ fluctuations in energy.
 
 For :math:`\kappa = 0`, this corresponds to the plain Coulomb potential.
 
+.. _Reaction Field method:
+
+Reaction Field method
+---------------------
+
+:class:`espressomd.electrostatics.ReactionField`
+
+The Reaction Field electrostatic potential is defined by
+
+  .. math:: U^{C-RF} = C \cdot q_1 q_2 \left[\frac{1}{r} - \frac{B r^2}{2r_{\mathrm{cut}}^3} -\frac{1 - B/2}{r_{\mathrm{cut}}}\right] \quad \mathrm{for}\quad r<r_{\mathrm{cut}}
+
+where :math:`C` is defined as in Eqn. :eq:`coulomb_prefactor` and :math:`B`
+is defined as:
+
+  .. math:: B = \frac{2(\varepsilon_1 - \varepsilon_2)(1 + \kappa r_{\mathrm{cut}}) - \varepsilon_2 (\kappa r_{\mathrm{cut}})^2}{(\varepsilon_1 + 2\varepsilon_2)(1 + \kappa r_{\mathrm{cut}}) + \varepsilon_2 (\kappa r_{\mathrm{cut}})^2}
+
+with :math:`\kappa` the ionic strength, :math:`\varepsilon_1` the dielectric
+constant inside the cavity and :math:`\varepsilon_2` the dielectric constant
+outside the cavity :cite:`tironi95a`.
+
+The term in :math:`1 - B/2` is a correction to make the
+potential continuous at :math:`r = r_{\mathrm{cut}}`.
+
 
 .. _Dielectric interfaces with the ICC algorithm:
 

doc/sphinx/zrefs.bib

@@ -884,6 +884,17 @@ @ARTICLE{thompson09a
   pages = {154107},
 }
 
+@Article{tironi95a,
+  author    = {Tironi, Ilario G. and Sperb, Ren\'{e} and Smith, Paul E. and van Gunsteren, Wilfred F.},
+  title     = {A generalized Reaction Field Method for Molecular-Dynamics Simulations},
+  journal   = {The Journal of Chemical Physics},
+  year      = {1995},
+  volume    = {102},
+  number    = {13},
+  pages     = {5451--5459},
+  doi       = {10.1063/1.469273},
+}
+
 @article{turner2008simulation,
   title={Simulation of chemical reaction equilibria by the reaction ensemble {M}onte {C}arlo method: a review},
   author={Heath Turner, C and Brennan, John K and Lisal, Martin and Smith, William R and Karl Johnson, J and Gubbins, Keith E},

src/core/electrostatics_magnetostatics/coulomb.cpp

@@ -139,15 +139,10 @@ void deactivate() {
     break;
 #endif
   case COULOMB_DH:
-    dh_params.r_cut = 0.0;
-    dh_params.kappa = 0.0;
+    dh_params = {};
     break;
   case COULOMB_RF:
-    rf_params.kappa = 0.0;
-    rf_params.epsilon1 = 0.0;
-    rf_params.epsilon2 = 0.0;
-    rf_params.r_cut = 0.0;
-    rf_params.B = 0.0;
+    rf_params = {};
     break;
   case COULOMB_MMM1D:
     mmm1d_params.maxPWerror = 1e40;
@@ -447,7 +442,7 @@ void bcast_coulomb_params() {
 
 void set_prefactor(double prefactor) {
   if (prefactor < 0.0) {
-    throw std::invalid_argument("Coulomb prefactor has to be >= 0");
+    throw std::domain_error("Coulomb prefactor has to be >= 0");
   }
 
   coulomb.prefactor = prefactor;

src/core/electrostatics_magnetostatics/debye_hueckel.cpp

@@ -28,21 +28,19 @@
 #ifdef ELECTROSTATICS
 #include "electrostatics_magnetostatics/common.hpp"
 
-Debye_hueckel_params dh_params{};
+#include <stdexcept>
 
-int dh_set_params(double kappa, double r_cut) {
-  if (dh_params.kappa < 0.0)
-    return -1;
+Debye_hueckel_params dh_params{};
 
-  if (dh_params.r_cut < 0.0)
-    return -2;
+void dh_set_params(double kappa, double r_cut) {
+  if (kappa < 0.0)
+    throw std::domain_error("kappa should be a non-negative number");
+  if (r_cut < 0.0)
+    throw std::domain_error("r_cut should be a non-negative number");
 
-  dh_params.kappa = kappa;
-  dh_params.r_cut = r_cut;
+  dh_params = {r_cut, kappa};
 
   mpi_bcast_coulomb_params();
-
-  return 1;
 }
 
 #endif

src/core/electrostatics_magnetostatics/debye_hueckel.hpp

@@ -21,7 +21,7 @@
 #ifndef DEBYE_HUECKEL_H
 #define DEBYE_HUECKEL_H
 /** \file
- *  Routines to calculate the Debye-Hueckel energy and force
+ *  Routines to calculate the Debye-Hückel energy and force
  *  for a particle pair.
  */
 #include "config.hpp"
@@ -32,18 +32,18 @@
 
 #include <cmath>
 
-/** Structure to hold Debye-Hueckel parameters. */
-typedef struct {
-  /** Cutoff for Debye-Hueckel interaction. */
+/** Debye-Hückel parameters. */
+struct Debye_hueckel_params {
+  /** Interaction cutoff. */
   double r_cut;
-  /** Debye kappa (inverse Debye length) . */
+  /** Ionic strength. */
   double kappa;
-} Debye_hueckel_params;
+};
 
-/** Debye-Hueckel parameters. */
+/** Global state of the Debye-Hückel method. */
 extern Debye_hueckel_params dh_params;
 
-int dh_set_params(double kappa, double r_cut);
+void dh_set_params(double kappa, double r_cut);
 
 /** Compute the Debye-Hueckel pair force.
  *  @param[in]  q1q2      Product of the charges on p1 and p2.

src/core/electrostatics_magnetostatics/reaction_field.cpp

@@ -29,27 +29,23 @@
 
 Reaction_field_params rf_params{};
 
-int rf_set_params(double kappa, double epsilon1, double epsilon2,
-                  double r_cut) {
-  rf_params.kappa = kappa;
-  rf_params.epsilon1 = epsilon1;
-  rf_params.epsilon2 = epsilon2;
-  rf_params.r_cut = r_cut;
-  rf_params.B = (2 * (epsilon1 - epsilon2) * (1 + kappa * r_cut) -
-                 epsilon2 * kappa * kappa * r_cut * r_cut) /
-                ((epsilon1 + 2 * epsilon2) * (1 + kappa * r_cut) +
-                 epsilon2 * kappa * kappa * r_cut * r_cut);
-  if (rf_params.epsilon1 < 0.0)
-    return -1;
+void rf_set_params(double kappa, double epsilon1, double epsilon2,
+                   double r_cut) {
+  if (kappa < 0.0)
+    throw std::domain_error("kappa should be a non-negative number");
+  if (epsilon1 < 0.0)
+    throw std::domain_error("epsilon1 should be a non-negative number");
+  if (epsilon2 < 0.0)
+    throw std::domain_error("epsilon2 should be a non-negative number");
+  if (r_cut < 0.0)
+    throw std::domain_error("r_cut should be a non-negative number");
 
-  if (rf_params.epsilon2 < 0.0)
-    return -1;
-
-  if (rf_params.r_cut < 0.0)
-    return -2;
+  auto const B = (2 * (epsilon1 - epsilon2) * (1 + kappa * r_cut) -
+                  epsilon2 * kappa * kappa * r_cut * r_cut) /
+                 ((epsilon1 + 2 * epsilon2) * (1 + kappa * r_cut) +
+                  epsilon2 * kappa * kappa * r_cut * r_cut);
+  rf_params = {kappa, epsilon1, epsilon2, r_cut, B};
 
   mpi_bcast_coulomb_params();
-
-  return 1;
 }
 #endif

src/core/electrostatics_magnetostatics/reaction_field.hpp

@@ -22,7 +22,7 @@
 #define REACTION_FIELD_H
 /** \file
  *  Routines to calculate the Reaction Field energy or/and force
- *  for a particle pair @cite neumann85b.
+ *  for a particle pair @cite neumann85b, @cite tironi95a.
  *
  *  Implementation in \ref reaction_field.cpp
  */
@@ -32,71 +32,60 @@
 #ifdef ELECTROSTATICS
 
 #include <utils/Vector.hpp>
+#include <utils/math/int_pow.hpp>
 
-/** Structure to hold Reaction Field Parameters. */
-typedef struct {
-  /** ionic strength. */
+/** Reaction Field parameters. */
+struct Reaction_field_params {
+  /** Ionic strength. */
   double kappa;
-  /** epsilon1 (continuum dielectric constant inside). */
+  /** Continuum dielectric constant inside the cavity. */
   double epsilon1;
-  /** epsilon2 (continuum dielectric constant outside). */
+  /** Continuum dielectric constant outside the cavity. */
   double epsilon2;
-  /** Cutoff for Reaction Field interaction. */
+  /** Interaction cutoff. */
   double r_cut;
-  /** B important prefactor. */
+  /** Interaction prefactor. Corresponds to the quantity
+   *  @f$ 1 + B_1 @f$ from eq. 22 in @cite tironi95a.
+   */
   double B;
-} Reaction_field_params;
+};
 
-/** Structure containing the Reaction Field parameters. */
+/** Global state of the Reaction Field method. */
 extern Reaction_field_params rf_params;
 
-int rf_set_params(double kappa, double epsilon1, double epsilon2, double r_cut);
-
-inline void add_rf_coulomb_pair_force_no_cutoff(double const q1q2,
-                                                Utils::Vector3d const &d,
-                                                double const dist,
-                                                Utils::Vector3d &force) {
-  double fac;
-  fac = 1.0 / (dist * dist * dist) +
-        rf_params.B / (rf_params.r_cut * rf_params.r_cut * rf_params.r_cut);
-  fac *= q1q2;
-  force += fac * d;
-}
+void rf_set_params(double kappa, double epsilon1, double epsilon2,
+                   double r_cut);
 
 /** Compute the Reaction Field pair force.
- *  @param q1q2      Product of the charges on p1 and p2.
- *  @param d         Vector pointing from p1 to p2.
- *  @param dist      Distance between p1 and p2.
- *  @param force     returns the force on particle 1.
+ *  @param[in]  q1q2      Product of the charges on p1 and p2.
+ *  @param[in]  d         Vector pointing from p1 to p2.
+ *  @param[in]  dist      Distance between p1 and p2.
+ *  @param[out] force     Calculated force on p1.
  */
 inline void add_rf_coulomb_pair_force(double const q1q2,
                                       Utils::Vector3d const &d,
                                       double const dist,
                                       Utils::Vector3d &force) {
   if (dist < rf_params.r_cut) {
-    add_rf_coulomb_pair_force_no_cutoff(q1q2, d, dist, force);
+    auto fac = 1.0 / Utils::int_pow<3>(dist) +
+               rf_params.B / Utils::int_pow<3>(rf_params.r_cut);
+    fac *= q1q2;
+    force += fac * d;
   }
 }
 
 /** Compute the Reaction Field pair energy.
  *  @param q1q2      Product of the charges on p1 and p2.
  *  @param dist      Distance between p1 and p2.
  */
-inline double rf_coulomb_pair_energy_no_cutoff(double const q1q2,
-                                               double const dist) {
-  double fac;
-  fac = 1.0 / dist -
-        (rf_params.B * dist * dist) /
-            (2 * rf_params.r_cut * rf_params.r_cut * rf_params.r_cut);
-  // cut off part
-  fac -= (1 - rf_params.B / 2) / rf_params.r_cut;
-  fac *= q1q2;
-  return fac;
-}
-
 inline double rf_coulomb_pair_energy(double const q1q2, double const dist) {
   if (dist < rf_params.r_cut) {
-    return rf_coulomb_pair_energy_no_cutoff(q1q2, dist);
+    auto fac = 1.0 / dist - (rf_params.B * dist * dist) /
+                                (2 * Utils::int_pow<3>(rf_params.r_cut));
+    // remove discontinuity at dist = r_cut
+    fac -= (1 - rf_params.B / 2) / rf_params.r_cut;
+    fac *= q1q2;
+    return fac;
   }
   return 0.0;
 }

src/python/espressomd/electrostatics.pxd

@@ -60,7 +60,7 @@ IF ELECTROSTATICS:
 
     cdef extern from "electrostatics_magnetostatics/coulomb.hpp" namespace "Coulomb":
 
-        int set_prefactor(double prefactor) except+
+        int set_prefactor(double prefactor) except +
         void deactivate_method()
 
     IF P3M:
@@ -108,7 +108,7 @@ IF ELECTROSTATICS:
 
         cdef extern Debye_hueckel_params dh_params
 
-        int dh_set_params(double kappa, double r_cut)
+        void dh_set_params(double kappa, double r_cut) except +
 
     cdef extern from "electrostatics_magnetostatics/reaction_field.hpp":
         ctypedef struct Reaction_field_params:
@@ -119,8 +119,8 @@ IF ELECTROSTATICS:
 
         cdef extern Reaction_field_params rf_params
 
-        int rf_set_params(double kappa, double epsilon1, double epsilon2,
-                          double r_cut)
+        void rf_set_params(double kappa, double epsilon1, double epsilon2,
+                           double r_cut) except +
 
 IF ELECTROSTATICS:
     cdef extern from "electrostatics_magnetostatics/mmm1d.hpp":