Make ELC an actor

doc/sphinx/electrostatics.rst

@@ -234,7 +234,7 @@ using it.
 Electrostatic Layer Correction (ELC)
 ------------------------------------
 
-:class:`espressomd.electrostatic_extensions.ELC`
+:class:`espressomd.electrostatics.ELC`
 
 *ELC* is an extension of the P3M electrostatics solver for explicit 2D periodic
 systems. It can account for different dielectric jumps on both sides of the
@@ -260,8 +260,9 @@ Usage notes:
 
 *ELC* is an |es| actor and is used with::
 
-    import espressomd.electrostatic_extensions
-    elc = electrostatic_extensions.ELC(gap_size=box_l * 0.2, maxPWerror=1e-3)
+    import espressomd.electrostatics
+    p3m = espressomd.electrostatics.P3M(prefactor=1, accuracy=1e-4)
+    elc = espressomd.electrostatics.ELC(p3m_actor=p3m, gap_size=box_l * 0.2, maxPWerror=1e-3)
     system.actors.add(elc)
 
 *ELC* can also be used to simulate 2D periodic systems with image charges,
@@ -273,8 +274,8 @@ simulation region (*middle*) to *bottom* (at :math:`z=0`) and from *middle* to
 are :math:`\Delta_t=\frac{\varepsilon_m-\varepsilon_t}{\varepsilon_m+\varepsilon_t}`
 and :math:`\Delta_b=\frac{\varepsilon_m-\varepsilon_b}{\varepsilon_m+\varepsilon_b}`::
 
-    elc = electrostatic_extensions.ELC(gap_size=box_l * 0.2, maxPWerror=1e-3,
-                                       delta_mid_top=0.9, delta_mid_bot=0.1)
+    elc = espressomd.electrostatics.ELC(p3m_actor=p3m, gap_size=box_l * 0.2, maxPWerror=1e-3,
+                                        delta_mid_top=0.9, delta_mid_bot=0.1)
 
 The fully metallic case :math:`\Delta_t=\Delta_b=-1` would lead to divergence
 of the forces/energies in *ELC* and is therefore only possible with the
@@ -283,8 +284,8 @@ of the forces/energies in *ELC* and is therefore only possible with the
 Toggle ``const_pot`` on to maintain a constant electric potential difference
 ``pot_diff`` between the xy-planes at :math:`z=0` and :math:`z = L_z - h`::
 
-    elc = electrostatic_extensions.ELC(gap_size=box_l * 0.2, maxPWerror=1e-3,
-                                       const_pot=True, delta_mid_bot=100.0)
+    elc = espressomd.electrostatics.ELC(p3m_actor=p3m, gap_size=box_l * 0.2, maxPWerror=1e-3,
+                                        const_pot=True, delta_mid_bot=100.0)
 
 This is done by countering the total dipole moment of the system with the
 electric field :math:`E_{\textrm{induced}}` and superposing a homogeneous

samples/visualization_elc.py

@@ -27,7 +27,6 @@
 import espressomd
 import espressomd.shapes
 from espressomd import electrostatics
-from espressomd import electrostatic_extensions
 from espressomd import visualization
 
 required_features = ["P3M", "WCA"]
@@ -75,11 +74,8 @@
 system.thermostat.set_langevin(kT=0.1, gamma=1.0, seed=42)
 
 p3m = electrostatics.P3M(prefactor=1.0, accuracy=1e-2)
-
-system.actors.add(p3m)
-
-elc = electrostatic_extensions.ELC(maxPWerror=1.0, gap_size=elc_gap,
-                                   const_pot=True, pot_diff=potential_diff)
+elc = electrostatics.ELC(p3m_actor=p3m, maxPWerror=1.0, gap_size=elc_gap,
+                         const_pot=True, pot_diff=potential_diff)
 system.actors.add(elc)
 
 visualizer.run(1)

src/core/electrostatics_magnetostatics/elc.cpp

@@ -1081,15 +1081,6 @@ void ELC_init() {
       runtimeErrorMsg() << "ELC auto-retuning failed";
     }
   }
-  if (elc_params.dielectric_contrast_on) {
-    p3m.params.additional_mesh[0] = 0;
-    p3m.params.additional_mesh[1] = 0;
-    p3m.params.additional_mesh[2] = elc_params.space_layer;
-  } else {
-    p3m.params.additional_mesh[0] = 0;
-    p3m.params.additional_mesh[1] = 0;
-    p3m.params.additional_mesh[2] = 0;
-  }
 }
 
 int ELC_set_params(double maxPWerror, double gap_size, double far_cut,

src/core/electrostatics_magnetostatics/p3m-common.cpp

@@ -115,14 +115,17 @@ double p3m_analytic_cotangent_sum(int n, double mesh_i, int cao) {
 
 void p3m_calc_local_ca_mesh(p3m_local_mesh &local_mesh,
                             const P3MParameters &params,
-                            const LocalBox<double> &local_geo, double skin) {
+                            const LocalBox<double> &local_geo, double skin,
+                            double space_layer) {
   int i;
   int ind[3];
   /* total skin size */
   double full_skin[3];
 
   for (i = 0; i < 3; i++)
-    full_skin[i] = params.cao_cut[i] + skin + params.additional_mesh[i];
+    full_skin[i] = params.cao_cut[i] + skin;
+
+  full_skin[2] += space_layer;
 
   /* inner left down grid point (global index) */
   for (i = 0; i < 3; i++)

src/core/electrostatics_magnetostatics/p3m-common.hpp

@@ -147,14 +147,11 @@ typedef struct {
   /** number of points unto which a single charge is interpolated, i.e.
    *  p3m.cao^3 */
   int cao3 = 0;
-  /** additional points around the charge assignment mesh, for method like
-   *  dielectric ELC creating virtual charges. */
-  double additional_mesh[3] = {};
 
   template <typename Archive> void serialize(Archive &ar, long int) {
     ar &tuning &alpha_L &r_cut_iL &mesh;
     ar &mesh_off &cao &accuracy &epsilon &cao_cut;
-    ar &a &ai &alpha &r_cut &cao3 &additional_mesh;
+    ar &a &ai &alpha &r_cut &cao3;
   }
 
 } P3MParameters;
@@ -186,7 +183,8 @@ double p3m_analytic_cotangent_sum(int n, double mesh_i, int cao);
  */
 void p3m_calc_local_ca_mesh(p3m_local_mesh &local_mesh,
                             const P3MParameters &params,
-                            const LocalBox<double> &local_geo, double skin);
+                            const LocalBox<double> &local_geo, double skin,
+                            double space_layer);
 
 /** Calculate the spatial position of the left down mesh
  *  point of the local mesh, to be stored in

src/core/electrostatics_magnetostatics/p3m-dipolar.cpp

@@ -200,7 +200,7 @@ void dp3m_init() {
    * and the cutoff for charge assignment dp3m.params.cao_cut */
   dp3m_init_a_ai_cao_cut();
 
-  p3m_calc_local_ca_mesh(dp3m.local_mesh, dp3m.params, local_geo, skin);
+  p3m_calc_local_ca_mesh(dp3m.local_mesh, dp3m.params, local_geo, skin, 0.0);
 
   dp3m.sm.resize(comm_cart, dp3m.local_mesh);
 

src/core/electrostatics_magnetostatics/p3m.cpp

@@ -185,7 +185,13 @@ void p3m_init() {
     return;
   }
 
-  p3m_calc_local_ca_mesh(p3m.local_mesh, p3m.params, local_geo, skin);
+  double elc_layer = 0.0;
+  if (coulomb.method == COULOMB_ELC_P3M) {
+    elc_layer = elc_params.space_layer;
+  }
+
+  p3m_calc_local_ca_mesh(p3m.local_mesh, p3m.params, local_geo, skin,
+                         elc_layer);
 
   p3m.sm.resize(comm_cart, p3m.local_mesh);
 

src/python/espressomd/electrostatic_extensions.pxd

@@ -20,28 +20,10 @@
 include "myconfig.pxi"
 from .electrostatics cimport *
 from libcpp.vector cimport vector
-from libcpp cimport bool
 from .utils cimport Vector3d
 
 IF ELECTROSTATICS and P3M:
 
-    cdef extern from "electrostatics_magnetostatics/elc.hpp":
-        ctypedef struct ELC_struct:
-            double maxPWerror
-            double gap_size
-            double far_cut
-            bool neutralize
-            double delta_mid_top,
-            double delta_mid_bot,
-            bool const_pot,
-            double pot_diff
-
-        int ELC_set_params(double maxPWerror, double min_dist, double far_cut,
-                           bool neutralize, double delta_mid_top, double delta_mid_bot, bool const_pot, double pot_diff)
-
-        # links intern C-struct with python object
-        ELC_struct elc_params
-
     cdef extern from "electrostatics_magnetostatics/icc.hpp":
         ctypedef struct iccp3m_struct:
             int n_ic

src/python/espressomd/electrostatic_extensions.pyx

@@ -31,124 +31,6 @@ IF ELECTROSTATICS and P3M:
     cdef class ElectrostaticExtensions(actors.Actor):
         pass
 
-    cdef class ELC(ElectrostaticExtensions):
-        """
-        Electrostatics solver for systems with two periodic dimensions.
-        See :ref:`Electrostatic Layer Correction (ELC)` for more details.
-
-        Parameters
-        ----------
-        gap_size : :obj:`float`, required
-            The gap size gives the height :math:`h` of the empty region between
-            the system box and the neighboring artificial images. |es| does not
-            make sure that the gap is actually empty, this is the user's
-            responsibility. The method will run even if the condition is not
-            fulfilled, however, the error bound will not be reached. Therefore
-            you should really make sure that the gap region is empty (e.g.
-            with wall constraints).
-        maxPWerror : :obj:`float`, required
-            The maximal pairwise error sets the least upper bound (LUB) error
-            of the force between any two charges without prefactors (see the
-            papers). The algorithm tries to find parameters to meet this LUB
-            requirements or will throw an error if there are none.
-        delta_mid_top : :obj:`float`, optional
-            Dielectric contrast :math:`\\Delta_t` between the upper boundary
-            and the simulation box.
-        delta_mid_bottom : :obj:`float`, optional
-            Dielectric contrast :math:`\\Delta_b` between the lower boundary
-            and the simulation box.
-        const_pot : :obj:`bool`, optional
-            Activate a constant electric potential between the top and bottom
-            of the simulation box.
-        pot_diff : :obj:`float`, optional
-            If ``const_pot`` is enabled, this parameter controls the applied
-            voltage between the boundaries of the simulation box in the
-            *z*-direction (at :math:`z = 0` and :math:`z = L_z - h`).
-        neutralize : :obj:`bool`, optional
-            By default, *ELC* just as P3M adds a homogeneous neutralizing
-            background to the system in case of a net charge. However, unlike
-            in three dimensions, this background adds a parabolic potential
-            across the slab :cite:`ballenegger09a`. Therefore, under normal
-            circumstances, you will probably want to disable the neutralization
-            for non-neutral systems. This corresponds then to a formal
-            regularization of the forces and energies :cite:`ballenegger09a`.
-            Also, if you add neutralizing walls explicitly as constraints, you
-            have to disable the neutralization. When using a dielectric
-            contrast or full metallic walls (``delta_mid_top != 0`` or
-            ``delta_mid_bot != 0`` or ``const_pot=True``), ``neutralize`` is
-            overwritten and switched off internally. Note that the special
-            case of non-neutral systems with a *non-metallic* dielectric jump
-            (e.g. ``delta_mid_top`` or ``delta_mid_bot`` in ``]-1,1[``) is not
-            covered by the algorithm and will throw an error.
-        far_cut : :obj:`float`, optional
-            Cutoff radius, use with care, intended for testing purposes. When
-            setting the cutoff directly, the maximal pairwise error is ignored.
-        """
-
-        def validate_params(self):
-            default_params = self.default_params()
-            check_type_or_throw_except(
-                self._params["maxPWerror"], 1, float, "")
-            check_range_or_except(
-                self._params, "maxPWerror", 0, False, "inf", True)
-            check_type_or_throw_except(self._params["gap_size"], 1, float, "")
-            check_range_or_except(
-                self._params, "gap_size", 0, False, "inf", True)
-            check_type_or_throw_except(self._params["far_cut"], 1, float, "")
-            check_type_or_throw_except(
-                self._params["neutralize"], 1, type(True), "")
-
-        def valid_keys(self):
-            return ["maxPWerror", "gap_size", "far_cut", "neutralize",
-                    "delta_mid_top", "delta_mid_bot", "const_pot", "pot_diff",
-                    "check_neutrality"]
-
-        def required_keys(self):
-            return ["maxPWerror", "gap_size"]
-
-        def default_params(self):
-            return {"maxPWerror": -1,
-                    "gap_size": -1,
-                    "far_cut": -1,
-                    "delta_mid_top": 0,
-                    "delta_mid_bot": 0,
-                    "const_pot": False,
-                    "pot_diff": 0.0,
-                    "neutralize": True,
-                    "check_neutrality": True}
-
-        def _get_params_from_es_core(self):
-            params = {}
-            params.update(elc_params)
-            return params
-
-        def _set_params_in_es_core(self):
-            if coulomb.method == COULOMB_P3M_GPU:
-                raise Exception("ELC is not set up to work with the GPU P3M")
-
-            if self._params["const_pot"]:
-                self._params["delta_mid_top"] = -1
-                self._params["delta_mid_bot"] = -1
-
-            if ELC_set_params(
-                self._params["maxPWerror"],
-                self._params["gap_size"],
-                self._params["far_cut"],
-                self._params["neutralize"],
-                self._params["delta_mid_top"],
-                self._params["delta_mid_bot"],
-                self._params["const_pot"],
-                    self._params["pot_diff"]):
-                handle_errors("ELC tuning failed")
-
-        def _activate_method(self):
-            check_neutrality(self._params)
-            self._set_params_in_es_core()
-
-        def _deactivate_method(self):
-            raise Exception(
-                "Unable to remove ELC as the state of the underlying electrostatics method will remain unclear.")
-
     cdef class ICC(ElectrostaticExtensions):
         """
         Interface to the induced charge calculation scheme for dielectric

src/python/espressomd/electrostatics.pxd

@@ -83,6 +83,24 @@ IF ELECTROSTATICS:
             cdef extern from "electrostatics_magnetostatics/p3m_gpu.hpp":
                 void p3m_gpu_init(int cao, int * mesh, double alpha) except +
 
+        cdef extern from "electrostatics_magnetostatics/elc.hpp":
+            ctypedef struct ELC_struct:
+                double maxPWerror
+                double gap_size
+                double far_cut
+                bool neutralize
+                double delta_mid_top
+                double delta_mid_bot
+                bool const_pot
+                double pot_diff
+
+            int ELC_set_params(double maxPWerror, double min_dist, double far_cut,
+                               bool neutralize, double delta_mid_top,
+                               double delta_mid_bot, bool const_pot, double pot_diff)
+
+            # links intern C-struct with python object
+            ELC_struct elc_params
+
     cdef extern from "electrostatics_magnetostatics/debye_hueckel.hpp":
         ctypedef struct Debye_hueckel_params:
             double r_cut