Refactor new non-bonded interface

Interactions can now be reset. Skip inactive interactions when
reloading from a checkpoint file to avoid triggering range checks.
Re-introduce the original cutoff range checks.

doc/sphinx/inter_non-bonded.rst

@@ -18,23 +18,41 @@ explicitly set. A bonded interaction between a set of particles has to
 be specified explicitly by the command, while the command is used to
 define the interaction parameters.
 
-.. _Isotropic non-bonded interactions:
+Non-bonded interaction are configured via the
+:class:`espressomd.interactions.NonBondedInteraction` class,
+which is a member of :class:`espressomd.system.System`::
 
-Isotropic non-bonded interactions
----------------------------------
+    system.non_bonded_inter[type1, type2]
 
-Non-bonded interaction are configured via the :class:`espressomd.interactions.NonBondedInteraction` class, which is a member of :class:`espressomd.system.System`::
+This command defines an interaction between all particles of type ``type1``
+and ``type2``. All available interaction potentials and their parameters are
+listed below. For example, the following adds a WCA potential between
+particles of type 0::
 
-    system.non_bonded_inter[type1, type2]
+    system.non_bonded_inter[0, 0].wca.set_params(epsilon=1., sigma=2.)
+
+Each type pair can have multiple potentials active at the same time.
+To deactivate a specific potential for a given type pair, do::
+
+    system.non_bonded_inter[0, 0].wca.deactivate()
+
+To deactivate all potentials for a given type pair, do::
 
-This command defines an interaction between all particles of type ``type1`` and
-``type2``. Possible interaction types and their parameters are
-listed below.
+    system.non_bonded_inter[0, 0].reset()
+
+To deactivate all potentials between all type pairs, do::
+
+    system.non_bonded_inter.reset()
 
 For many non-bonded interactions, it is possible to artificially cap the
 forces, which often allows to equilibrate the system much faster. See
 the subsection :ref:`Capping the force during warmup` for more details.
 
+.. _Isotropic non-bonded interactions:
+
+Isotropic non-bonded interactions
+---------------------------------
+
 .. _Tabulated interaction:
 
 Tabulated interaction

src/core/nonbonded_interactions/bmhtf-nacl.cpp

@@ -43,6 +43,9 @@ BMHTF_Parameters::BMHTF_Parameters(double a, double b, double c, double d,
   if (d < 0.) {
     throw std::domain_error("BMHTF parameter 'd' has to be >= 0");
   }
+  if (cutoff < 0.) {
+    throw std::domain_error("BMHTF parameter 'cutoff' has to be >= 0");
+  }
   computed_shift = C / Utils::int_pow<6>(cut) + D / Utils::int_pow<8>(cut) -
                    A * std::exp(B * (sig - cut));
 }

src/core/nonbonded_interactions/buckingham.cpp

@@ -45,6 +45,9 @@ Buckingham_Parameters::Buckingham_Parameters(double a, double b, double c,
   if (d < 0.) {
     throw std::domain_error("Buckingham parameter 'd' has to be >= 0");
   }
+  if (cutoff < 0.) {
+    throw std::domain_error("Buckingham parameter 'cutoff' has to be >= 0");
+  }
 
   /* Replace the Buckingham potential for interatomic distance less
      than or equal to discontinuity by a straight line (F1+F2*r) */

src/core/nonbonded_interactions/gaussian.cpp

@@ -29,13 +29,16 @@
 
 #include <stdexcept>
 
-Gaussian_Parameters::Gaussian_Parameters(double eps, double sig, double cut)
-    : eps{eps}, sig{sig}, cut{cut} {
+Gaussian_Parameters::Gaussian_Parameters(double eps, double sig, double cutoff)
+    : eps{eps}, sig{sig}, cut{cutoff} {
   if (eps < 0.) {
     throw std::domain_error("Gaussian parameter 'eps' has to be >= 0");
   }
   if (sig < 0.) {
     throw std::domain_error("Gaussian parameter 'sig' has to be >= 0");
   }
+  if (cutoff < 0.) {
+    throw std::domain_error("Gaussian parameter 'cutoff' has to be >= 0");
+  }
 }
 #endif // GAUSSIAN

src/core/nonbonded_interactions/hat.cpp

@@ -34,6 +34,9 @@ Hat_Parameters::Hat_Parameters(double F_max, double cutoff)
   if (F_max < 0.) {
     throw std::domain_error("Hat parameter 'F_max' has to be >= 0");
   }
+  if (cutoff < 0.) {
+    throw std::domain_error("Hat parameter 'cutoff' has to be >= 0");
+  }
 }
 
 #endif // HAT

src/core/nonbonded_interactions/hertzian.cpp

@@ -34,6 +34,9 @@ Hertzian_Parameters::Hertzian_Parameters(double eps, double sig)
   if (eps < 0.) {
     throw std::domain_error("Hertzian parameter 'eps' has to be >= 0");
   }
+  if (sig < 0.) {
+    throw std::domain_error("Hertzian parameter 'sig' has to be >= 0");
+  }
 }
 
 #endif // HERTZIAN

src/core/nonbonded_interactions/lj.cpp

@@ -27,31 +27,22 @@
 #ifdef LENNARD_JONES
 #include "nonbonded_interaction_data.hpp"
 
-#include <utils/math/int_pow.hpp>
-
 #include <algorithm>
 #include <stdexcept>
 
-LJ_Parameters::LJ_Parameters(double eps, double sig, double cut, double offset,
-                             double min)
-    : LJ_Parameters(eps, sig, cut, offset, min, 0.) {
-  if (cut != 0.) {
-    auto const sig_cut = sig / cut;
-    shift = Utils::int_pow<6>(sig_cut) - Utils::int_pow<12>(sig_cut);
-  }
-}
-
-LJ_Parameters::LJ_Parameters(double eps, double sig, double cut, double offset,
-                             double min, double shift)
-    : eps{eps}, sig{sig}, cut{cut}, shift{shift}, offset{offset}, min{std::max(
-                                                                      min,
-                                                                      0.)} {
-  if (eps < 0.) {
+LJ_Parameters::LJ_Parameters(double epsilon, double sigma, double cutoff,
+                             double offset, double min, double shift)
+    : eps{epsilon}, sig{sigma}, cut{cutoff}, shift{shift}, offset{offset},
+      min{std::max(min, 0.)} {
+  if (epsilon < 0.) {
     throw std::domain_error("LJ parameter 'epsilon' has to be >= 0");
   }
-  if (sig < 0.) {
+  if (sigma < 0.) {
     throw std::domain_error("LJ parameter 'sigma' has to be >= 0");
   }
+  if (cutoff < 0.) {
+    throw std::domain_error("LJ parameter 'cutoff' has to be >= 0");
+  }
 }
 
 #endif /* ifdef LENNARD_JONES */

src/core/nonbonded_interactions/lj.hpp

@@ -40,8 +40,7 @@
 /** Calculate Lennard-Jones force factor */
 inline double lj_pair_force_factor(IA_parameters const &ia_params,
                                    double dist) {
-  if ((dist < ia_params.lj.cut + ia_params.lj.offset) &&
-      (dist > ia_params.lj.min + ia_params.lj.offset)) {
+  if (dist < ia_params.lj.max_cutoff() and dist > ia_params.lj.min_cutoff()) {
     auto const r_off = dist - ia_params.lj.offset;
     auto const frac6 = Utils::int_pow<6>(ia_params.lj.sig / r_off);
     return 48.0 * ia_params.lj.eps * frac6 * (frac6 - 0.5) / (r_off * dist);
@@ -51,8 +50,7 @@ inline double lj_pair_force_factor(IA_parameters const &ia_params,
 
 /** Calculate Lennard-Jones energy */
 inline double lj_pair_energy(IA_parameters const &ia_params, double dist) {
-  if ((dist < ia_params.lj.cut + ia_params.lj.offset) &&
-      (dist > ia_params.lj.min + ia_params.lj.offset)) {
+  if (dist < ia_params.lj.max_cutoff() and dist > ia_params.lj.min_cutoff()) {
     auto const r_off = dist - ia_params.lj.offset;
     auto const frac6 = Utils::int_pow<6>(ia_params.lj.sig / r_off);
     return 4.0 * ia_params.lj.eps *

src/core/nonbonded_interactions/ljcos.cpp

@@ -32,15 +32,18 @@
 
 #include <stdexcept>
 
-LJcos_Parameters::LJcos_Parameters(double eps, double sig, double cut,
+LJcos_Parameters::LJcos_Parameters(double epsilon, double sigma, double cutoff,
                                    double offset)
-    : eps{eps}, sig{sig}, cut{cut}, offset{offset} {
-  if (eps < 0.) {
+    : eps{epsilon}, sig{sigma}, cut{cutoff}, offset{offset} {
+  if (epsilon < 0.) {
     throw std::domain_error("LJcos parameter 'epsilon' has to be >= 0");
   }
-  if (sig < 0.) {
+  if (sigma < 0.) {
     throw std::domain_error("LJcos parameter 'sigma' has to be >= 0");
   }
+  if (cutoff < 0.) {
+    throw std::domain_error("LJcos parameter 'cutoff' has to be >= 0");
+  }
   auto const facsq = Utils::cbrt_2() * Utils::sqr(sig);
 
   rmin = sqrt(Utils::cbrt_2()) * sig;

src/core/nonbonded_interactions/ljcos2.cpp

@@ -30,18 +30,21 @@
 #include <cmath>
 #include <stdexcept>
 
-LJcos2_Parameters::LJcos2_Parameters(double eps, double sig, double offset,
-                                     double w)
-    : eps{eps}, sig{sig}, offset{offset}, w{w} {
-  if (eps < 0.) {
+LJcos2_Parameters::LJcos2_Parameters(double epsilon, double sigma,
+                                     double offset, double width)
+    : eps{epsilon}, sig{sigma}, offset{offset}, w{width} {
+  if (epsilon < 0.) {
     throw std::domain_error("LJcos2 parameter 'epsilon' has to be >= 0");
   }
-  if (sig < 0.) {
+  if (sigma < 0.) {
     throw std::domain_error("LJcos2 parameter 'sigma' has to be >= 0");
   }
-  if (sig != 0.) {
-    rchange = std::pow(2., 1. / 6.) * sig;
-    cut = w + rchange;
+  if (width < 0.) {
+    throw std::domain_error("LJcos2 parameter 'width' has to be >= 0");
+  }
+  if (sigma != 0.) {
+    rchange = std::pow(2., 1. / 6.) * sigma;
+    cut = width + rchange;
   }
 }
 

src/core/nonbonded_interactions/ljgen.cpp

@@ -46,6 +46,9 @@ LJGen_Parameters::LJGen_Parameters(double epsilon, double sigma, double cutoff,
   if (sigma < 0.) {
     throw std::domain_error("Generic LJ parameter 'sigma' has to be >= 0");
   }
+  if (cutoff < 0.) {
+    throw std::domain_error("Generic LJ parameter 'cutoff' has to be >= 0");
+  }
 #ifdef LJGEN_SOFTCORE
   if (delta < 0.) {
     throw std::domain_error("Generic LJ parameter 'delta' has to be >= 0");

src/core/nonbonded_interactions/morse.cpp

@@ -36,6 +36,9 @@ Morse_Parameters::Morse_Parameters(double eps, double alpha, double rmin,
   if (eps < 0.) {
     throw std::domain_error("Morse parameter 'eps' has to be >= 0");
   }
+  if (cutoff < 0.) {
+    throw std::domain_error("Morse parameter 'cutoff' has to be >= 0");
+  }
   auto const add1 = std::exp(-2.0 * alpha * (cut - rmin));
   auto const add2 = 2.0 * std::exp(-alpha * (cut - rmin));
   rest = eps * (add1 - add2);

src/core/nonbonded_interactions/nonbonded_interaction_data.hpp

@@ -28,6 +28,7 @@
 #include "config.hpp"
 
 #include <utils/index.hpp>
+#include <utils/math/int_pow.hpp>
 
 #include <algorithm>
 #include <cassert>
@@ -49,10 +50,17 @@ struct LJ_Parameters {
   double offset = 0.0;
   double min = 0.0;
   LJ_Parameters() = default;
-  LJ_Parameters(double eps, double sig, double cut, double offset, double min);
-  LJ_Parameters(double eps, double sig, double cut, double offset, double min,
-                double shift);
+  LJ_Parameters(double epsilon, double sigma, double cutoff, double offset,
+                double min, double shift);
+  double get_auto_shift() const {
+    auto auto_shift = 0.;
+    if (cut != 0.) {
+      auto_shift = Utils::int_pow<6>(sig / cut) - Utils::int_pow<12>(sig / cut);
+    }
+    return auto_shift;
+  }
   double max_cutoff() const { return cut + offset; }
+  double min_cutoff() const { return min + offset; }
 };
 
 /** WCA potential */
@@ -61,7 +69,7 @@ struct WCA_Parameters {
   double sig = 0.0;
   double cut = INACTIVE_CUTOFF;
   WCA_Parameters() = default;
-  WCA_Parameters(double eps, double sig);
+  WCA_Parameters(double epsilon, double sigma);
   double max_cutoff() const { return cut; }
 };
 
@@ -124,7 +132,7 @@ struct Gaussian_Parameters {
   double sig = 1.0;
   double cut = INACTIVE_CUTOFF;
   Gaussian_Parameters() = default;
-  Gaussian_Parameters(double eps, double sig, double cut);
+  Gaussian_Parameters(double eps, double sig, double cutoff);
   double max_cutoff() const { return cut; }
 };
 
@@ -202,7 +210,7 @@ struct LJcos_Parameters {
   double beta = 0.0;
   double rmin = 0.0;
   LJcos_Parameters() = default;
-  LJcos_Parameters(double eps, double sig, double cut, double offset);
+  LJcos_Parameters(double epsilon, double sigma, double cutoff, double offset);
   double max_cutoff() const { return cut + offset; }
 };
 
@@ -215,7 +223,7 @@ struct LJcos2_Parameters {
   double w = 0.0;
   double rchange = 0.0;
   LJcos2_Parameters() = default;
-  LJcos2_Parameters(double eps, double sig, double offset, double w);
+  LJcos2_Parameters(double epsilon, double sigma, double offset, double width);
   double max_cutoff() const { return cut + offset; }
 };
 

src/core/nonbonded_interactions/smooth_step.cpp

@@ -36,6 +36,12 @@ SmoothStep_Parameters::SmoothStep_Parameters(double eps, double sig,
   if (eps < 0.) {
     throw std::domain_error("SmoothStep parameter 'eps' has to be >= 0");
   }
+  if (sig < 0.) {
+    throw std::domain_error("SmoothStep parameter 'sig' has to be >= 0");
+  }
+  if (cutoff < 0.) {
+    throw std::domain_error("SmoothStep parameter 'cutoff' has to be >= 0");
+  }
 }
 
 #endif // SMOOTH_STEP

src/core/nonbonded_interactions/soft_sphere.cpp

@@ -36,6 +36,9 @@ SoftSphere_Parameters::SoftSphere_Parameters(double a, double n, double cutoff,
   if (a < 0.) {
     throw std::domain_error("SoftSphere parameter 'a' has to be >= 0");
   }
+  if (cutoff < 0.) {
+    throw std::domain_error("SoftSphere parameter 'cutoff' has to be >= 0");
+  }
   if (offset < 0.) {
     throw std::domain_error("SoftSphere parameter 'offset' has to be >= 0");
   }

src/core/nonbonded_interactions/wca.cpp

@@ -28,15 +28,16 @@
 #include <cmath>
 #include <stdexcept>
 
-WCA_Parameters::WCA_Parameters(double eps, double sig) : eps{eps}, sig{sig} {
-  if (eps < 0.) {
+WCA_Parameters::WCA_Parameters(double epsilon, double sigma)
+    : eps{epsilon}, sig{sigma} {
+  if (epsilon < 0.) {
     throw std::domain_error("WCA parameter 'epsilon' has to be >= 0");
   }
-  if (sig < 0.) {
+  if (sigma < 0.) {
     throw std::domain_error("WCA parameter 'sigma' has to be >= 0");
   }
-  if (sig != 0.) {
-    cut = sig * std::pow(2., 1. / 6.);
+  if (sigma != 0.) {
+    cut = sigma * std::pow(2., 1. / 6.);
   }
 }