Fix ELC sample (#3468)

Fixes #3453

Description of changes:
- remove ELC from P3M sample
- create a stand-alone ELC sample with a gap and walls
- update documentation of electrostatics/magnetostatics methods

doc/sphinx/analysis.rst

@@ -376,12 +376,12 @@ where the first summand is the short ranged part and the second summand is the l
 
 The short ranged part is given by:
 
-.. math :: p^\text{Coulomb, P3M, dir}_{(k,l)}= \frac{1}{4\pi \epsilon_0 \epsilon_r} \frac{1}{2V} \sum_{\vec{n}}^* \sum_{i,j=1}^N q_i q_j \left( \frac{ \mathrm{erfc}(\beta |\vec{r}_j-\vec{r}_i+\vec{n}|)}{|\vec{r}_j-\vec{r}_i+\vec{n}|^3} + \\ \frac{2\beta \pi^{-1/2} \exp(-(\beta |\vec{r}_j-\vec{r}_i+\vec{n}|)^2)}{|\vec{r}_j-\vec{r}_i+\vec{n}|^2} \right) (\vec{r}_j-\vec{r}_i+\vec{n})_k (\vec{r}_j-\vec{r}_i+\vec{n})_l,
+.. math :: p^\text{Coulomb, P3M, dir}_{(k,l)}= \frac{1}{4\pi \varepsilon_0 \varepsilon_r} \frac{1}{2V} \sum_{\vec{n}}^* \sum_{i,j=1}^N q_i q_j \left( \frac{ \mathrm{erfc}(\beta |\vec{r}_j-\vec{r}_i+\vec{n}|)}{|\vec{r}_j-\vec{r}_i+\vec{n}|^3} + \\ \frac{2\beta \pi^{-1/2} \exp(-(\beta |\vec{r}_j-\vec{r}_i+\vec{n}|)^2)}{|\vec{r}_j-\vec{r}_i+\vec{n}|^2} \right) (\vec{r}_j-\vec{r}_i+\vec{n})_k (\vec{r}_j-\vec{r}_i+\vec{n})_l,
 
 where :math:`\beta` is the P3M splitting parameter, :math:`\vec{n}` identifies the periodic images, the asterisk denotes that terms with :math:`\vec{n}=\vec{0}` and i=j are omitted.
 The long ranged (k-space) part is given by:
 
-.. math :: p^\text{Coulomb, P3M, rec}_{(k,l)}= \frac{1}{4\pi \epsilon_0 \epsilon_r} \frac{1}{2 \pi V^2} \sum_{\vec{k} \neq \vec{0}} \frac{\exp(-\pi^2 \vec{k}^2/\beta^2)}{\vec{k}^2} |S(\vec{k})|^2 \cdot (\delta_{k,l}-2\frac{1+\pi^2\vec{k}^2/\beta^2}{\vec{k}^2} \vec{k}_k \vec{k}_l),
+.. math :: p^\text{Coulomb, P3M, rec}_{(k,l)}= \frac{1}{4\pi \varepsilon_0 \varepsilon_r} \frac{1}{2 \pi V^2} \sum_{\vec{k} \neq \vec{0}} \frac{\exp(-\pi^2 \vec{k}^2/\beta^2)}{\vec{k}^2} |S(\vec{k})|^2 \cdot (\delta_{k,l}-2\frac{1+\pi^2\vec{k}^2/\beta^2}{\vec{k}^2} \vec{k}_k \vec{k}_l),
 
 where :math:`S(\vec{k})` is the Fourier transformed charge density. Compared to Essmann we do not have the contribution :math:`p^\text{corr}_{k,l}` since we want to calculate the pressure that arises from all particles in the system.
 

doc/sphinx/appendix.rst

@@ -452,7 +452,7 @@ where
      &2u_xu_y\sum_{q>0}\frac{\cosh(2\pi f_q z)}{f_q(e^{2\pi f_q\lambda_z}
        - 1)}\cos(\omega_q y). \end{array}
 
-The implementation is very similar to MMM2d, except that the separation
+The implementation is very similar to MMM2D, except that the separation
 between slices close by, and above and below is not necessary.
 
 .. _Errors:

doc/sphinx/magnetostatics.rst

@@ -22,7 +22,8 @@ The prefactor :math:`D` is can be set by the user and is given by
    D =\frac{\mu_0 \mu}{4\pi}
    :label: dipolar_prefactor
 
-where :math:`\mu_0` and :math:`\mu` are the vacuum permittivity and the relative permittivity of the background material, respectively.
+where :math:`\mu_0` and :math:`\mu` are the vacuum permittivity and the
+relative permittivity of the background material, respectively.
 
 Magnetostatic interactions are activated via the actor framework::
 
@@ -38,48 +39,63 @@ some knowledge to use them properly. Uneducated use can result in
 completely unphysical simulations.
 
 
-
 .. _Dipolar P3M:
 
 Dipolar P3M
 ~~~~~~~~~~~
+
+:class:`espressomd.magnetostatics.DipolarP3M`
+
 This is the dipolar version of the P3M algorithm, described in :cite:`cerda08d`.
-It is interfaced via :class:`~espressomd.magnetostatics.DipolarP3M`.
 
 Make sure that you know the relevance of the P3M parameters before using
 P3M! If you are not sure, read the following references
 
 Note that dipolar P3M does not work with non-cubic boxes.
 
 
-The parameters of the dipolar P3M method can be tuned automatically, by providing ``accuracy=<TARGET_ACCURACY>`` to the method.
-It is also possible to pass a subset of the method parameters such as ``mesh``. In that case, only the omitted parameters are tuned::
-
+The parameters of the dipolar P3M method can be tuned automatically, by
+providing ``accuracy=<TARGET_ACCURACY>`` to the method. It is also possible to
+pass a subset of the method parameters such as ``mesh``. In that case, only
+the omitted parameters are tuned::
 
     import espressomd.magnetostatics as magnetostatics
     p3m = magnetostatics.DipolarP3M(prefactor=1, mesh=32, accuracy=1E-4)
     system.actors.add(p3m)
 
-It is important to note that the error estimates given in :cite:`cerda08d` used in the tuning contain assumptions about the system. In particular, a homogeneous system is assumed. If this is no longer the case during the simulation, actual force and torque errors can be significantly larger.
+It is important to note that the error estimates given in :cite:`cerda08d`
+used in the tuning contain assumptions about the system. In particular, a
+homogeneous system is assumed. If this is no longer the case during the
+simulation, actual force and torque errors can be significantly larger.
 
-.. _Dipolar Layer Correction (DLC):
 
+.. _Dipolar Layer Correction (DLC):
 
 Dipolar Layer Correction (DLC)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
 :class:`espressomd.magnetostatic_extensions.DLC`
 
-The dipolar layer correction (DLC) is used in conjunction with the dipolar P3M method to calculate dipolar interactions in a 2D-periodic system.
+The dipolar layer correction (DLC) is used in conjunction with the dipolar P3M
+method to calculate dipolar interactions in a 2D-periodic system.
 It is based on :cite:`brodka04a` and the dipolar version of
 :ref:`Electrostatic Layer Correction (ELC)`.
 
 Usage notes:
 
   * The non-periodic direction is always the **z-direction**.
 
-  * The method relies on a slab of the simulation box perpendicular to the z-direction not to contain particles. The size in z-direction of this slab is controlled by the ``gap_size`` parameter. The user has to ensure that no particles enter this region by means of constraints or by fixing the particles' z-coordinate. When there is no empty slab of the specified size, the method will silently produce wrong results.
+  * The method relies on a slab of the simulation box perpendicular to the
+    z-direction not to contain particles. The size in z-direction of this slab
+    is controlled by the ``gap_size`` parameter. The user has to ensure that
+    no particles enter this region by means of constraints or by fixing the
+    particles' z-coordinate. When there is no empty slab of the specified size,
+    the method will silently produce wrong results.
 
-  * The method can be tuned using the ``accuracy`` parameter. In contrast to the electrostatic method, it refers to the energy. Furthermore, it is assumed that all dipole moment are as large as the largest of the dipoles in the system.
+  * The method can be tuned using the ``accuracy`` parameter. In contrast to
+    the electrostatic method, it refers to the energy. Furthermore, it is
+    assumed that all dipole moment are as large as the largest of the dipoles
+    in the system.
 
 The method is used as follows::
 
@@ -92,8 +108,6 @@ The method is used as follows::
     system.actors.add(dlc)
 
 
-
-
 .. _Dipolar direct sum:
 
 Dipolar direct sum
@@ -114,7 +128,6 @@ Two methods are available:
 * :class:`~espressomd.magnetostatics.DipolarDirectSumCpu`
   performs the calculation in double precision on the Cpu.
 
-
 * :class:`~espressomd.magnetostatics.DipolarDirectSumGpu`
   performs the calculations in single precision on a Cuda-capable graphics card.
   The implementation is optimized for large systems of several thousand
@@ -123,13 +136,16 @@ Two methods are available:
   the rest of the gpu remains idle. Hence, the method will perform poorly
   for small systems.
 
-To use the methods, create an instance of either :class:`~espressomd.magnetostatics.DipolarDirectSumCpu` or :class:`~espressomd.magnetostatics.DipolarDirectSumGpu` and add it to the system's list of active actors. The only required parameter is the Prefactor :eq:`dipolar_prefactor`::
+To use the methods, create an instance of either
+:class:`~espressomd.magnetostatics.DipolarDirectSumCpu` or
+:class:`~espressomd.magnetostatics.DipolarDirectSumGpu` and add it to the
+system's list of active actors. The only required parameter is the Prefactor
+:eq:`dipolar_prefactor`::
 
   from espressomd.magnetostatics import DipolarDirectSumGpu
   dds = DipolarDirectSumGpu(bjerrum_length=1)
   system.actors.add(dds)
 
-
 For testing purposes, a variant of the dipolar direct sum is available which
 adds periodic copies to the system in periodic directions:
 :class:`~espressomd.magnetostatics.DipolarDirectSumWithReplicaCpu`.
@@ -141,12 +157,13 @@ rather to check the results you get from more efficient methods like P3M.
 do not support MPI parallelization.
 
 
+.. _Barnes-Hut octree sum on GPU:
 
-.. _Barnes-Hut octree sum on gpu:
-
-Barnes-Hut octree sum on gpu
+Barnes-Hut octree sum on GPU
 ----------------------------
 
+:class:`espressomd.magnetostatics.DipolarBarnesHutGpu`
+
 This interaction calculates energies and forces between dipoles by
 summing over the spatial octree cells (aka ``leaves``).
 Far enough cells are considered as a single dipole with a cumulative
@@ -157,17 +174,21 @@ an additive distance respectively. For the detailed description of the
 Barnes-Hut method application to the dipole-dipole interactions, please
 refer to :cite:`Polyakov2013`.
 
-To use the method, create an instance of :class:`~espressomd.magnetostatics.DipolarBarnesHutGpu` and add it to the system's list of active actors::
+To use the method, create an instance of :class:`~espressomd.magnetostatics.DipolarBarnesHutGpu`
+and add it to the system's list of active actors::
 
   from espressomd.magnetostatics import DipolarBarnesHutGpu
   bh = DipolarBarnesHutGpu(prefactor=pf_dds_gpu, epssq=200.0, itolsq=8.0)
   system.actors.add(bh)
 
-.. _ScaFaCoS Magnetostatics:
 
-ScaFaCoS Magnetostatics
+.. _ScaFaCoS magnetostatics:
+
+ScaFaCoS magnetostatics
 -----------------------
 
+:class:`espressomd.magnetostatics.Scafacos`
+
 |es| can use the methods from the ScaFaCoS *Scalable fast Coulomb solvers*
 library for dipoles, if the methods support dipolar calculations. The feature
 ``SCAFACOS_DIPOLES`` has to be added to :file:`myconfig.hpp` to activate this
@@ -176,15 +197,13 @@ the ScaFaCoS code.
 
 To use ScaFaCoS, create an instance of :class:`~espressomd.magnetostatics.Scafacos`
 and add it to the list of active actors. Three parameters have to be specified:
-
-* ``method_name``: name of the ScaFaCoS method being used.
-* ``method_params``: dictionary containing the method-specific parameters
-* ``prefactor``
-
-The method-specific parameters are described in the ScaFaCoS manual.
-In addition, methods supporting tuning have a parameter ``tolerance_field`` which sets the desired root mean square accuracy for the electric field
+``prefactor``, ``method_name``, ``method_params``. The method-specific
+parameters are described in the ScaFaCoS manual. In addition, methods
+supporting tuning have a parameter ``tolerance_field`` which sets the desired
+root mean square accuracy for the magnetic field.
 
 For details of the various methods and their parameters please refer to
-the ScaFaCoS manual. To use this feature, ScaFaCoS has to be built as a shared library. ScaFaCoS can be used only once, either for Coulomb or for dipolar interactions.
-
+the ScaFaCoS manual. To use this feature, ScaFaCoS has to be built as a
+shared library. ScaFaCoS can be used only once, either for Coulomb or for
+dipolar interactions.
 

doc/sphinx/under_the_hood.rst

@@ -16,13 +16,13 @@ Since basically all major parts of the main MD integration have to
 access the particle data, efficient access to the particle data is
 crucial for a fast MD code. Therefore the particle data needs some more
 elaborate organization, which will be presented here. A particle itself
-is represented by a structure (Particle) consisting of several
-substructures (e.g. ParticlePosition, ParticleForce or
-ParticleProperties), which in turn represent basic physical properties
+is represented by a structure (``Particle``) consisting of several
+substructures (e.g. ``ParticlePosition``, ``ParticleForce`` or
+``ParticleProperties``), which in turn represent basic physical properties
 such as position, force or charge. The particles are organized in one or
-more particle lists on each node, called Cell cells. The cells are
-arranged by several possible systems, the cellsystems as described
-above. A cell system defines a way the particles are stored in , i.e.
+more particle lists on each node, called ``CellPList``. The cells are
+arranged by several possible systems, as described in :ref:`Cellsystems`.
+A cell system defines a way the particles are stored in |es|, i.e.
 how they are distributed onto the processor nodes and how they are
 organized on each of them. Moreover a cell system also defines
 procedures to efficiently calculate the force, energy and pressure for

samples/p3m.py

@@ -18,9 +18,7 @@
 #
 """
 Simulate a Lennard-Jones liquid with charges. The P3M method is used to
-calculate electrostatic interactions. The ELC method can be optionally
-added to subtract the electrostatic contribution from the *z*-direction.
-For more details, see :ref:`Electrostatic Layer Correction (ELC)`.
+calculate electrostatic interactions.
 """
 import numpy as np
 import espressomd
@@ -29,7 +27,6 @@
 espressomd.assert_features(required_features)
 
 from espressomd import electrostatics
-from espressomd import electrostatic_extensions
 import argparse
 
 parser = argparse.ArgumentParser()
@@ -38,8 +35,6 @@
                    const="cpu", help="P3M on CPU", default="cpu")
 group.add_argument("--gpu", action="store_const", dest="mode",
                    const="gpu", help="P3M on GPU")
-group.add_argument("--elc", action="store_const", dest="mode",
-                   const="elc", help="P3M and ELC on CPU")
 args = parser.parse_args()
 
 
@@ -142,10 +137,6 @@
 for key in list(p3m_params.keys()):
     print("{} = {}".format(key, p3m_params[key]))
 
-if args.mode == "elc":
-    elc = electrostatic_extensions.ELC(maxPWerror=1.0, gap_size=1.0)
-    system.actors.add(elc)
-
 print(system.actors)
 
 #############################################################

samples/visualization_elc.py

@@ -0,0 +1,85 @@
+# Copyright (C) 2010-2019 The ESPResSo project
+#
+# This file is part of ESPResSo.
+#
+# ESPResSo is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# ESPResSo is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+Visualize charged particles confined between two plates of a capacitor with
+a potential difference. The system is periodic in the *xy*-plane but has a gap
+in the *z*-direction. The ELC method subtracts the electrostatic contribution
+from the periodic images in the *z*-direction. The system total charge is zero.
+For more details, see :ref:`Electrostatic Layer Correction (ELC)`.
+"""
+
+import numpy as np
+
+import espressomd
+import espressomd.shapes
+from espressomd.minimize_energy import steepest_descent
+from espressomd import electrostatics
+from espressomd import electrostatic_extensions
+from espressomd import visualization
+
+required_features = ["P3M", "WCA"]
+espressomd.assert_features(required_features)
+
+box_l = 20
+elc_gap = 10
+potential_diff = -3.
+system = espressomd.System(box_l=[box_l, box_l, box_l + elc_gap])
+system.set_random_state_PRNG()
+np.random.seed(seed=system.seed)
+visualizer = visualization.openGLLive(
+    system,
+    background_color=[1, 1, 1],
+    constraint_type_colors=[[0, 0, 0]],
+    camera_position=[70, 10, 15],
+    camera_right=[0, 0, -1])
+
+system.time_step = 0.02
+system.cell_system.skin = 0.4
+system.periodicity = [1, 1, 1]
+
+qion = 1
+for i in range(300):
+    rpos = np.random.random(3) * box_l
+    system.part.add(pos=rpos, type=0, q=qion)
+    qion *= -1
+
+system.constraints.add(shape=espressomd.shapes.Wall(
+    dist=0, normal=[0, 0, 1]), particle_type=1)
+system.constraints.add(shape=espressomd.shapes.Wall(
+    dist=-box_l, normal=[0, 0, -1]), particle_type=1)
+
+system.non_bonded_inter[0, 1].wca.set_params(epsilon=1.0, sigma=1.0)
+system.non_bonded_inter[0, 0].wca.set_params(epsilon=1.0, sigma=1.0)
+
+energy = system.analysis.energy()
+print("Before Minimization: E_total=", energy['total'])
+steepest_descent(system, f_max=10, gamma=50.0, max_steps=1000,
+                 max_displacement=0.2)
+energy = system.analysis.energy()
+print("After Minimization: E_total=", energy['total'])
+
+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)
+system.actors.add(elc)
+
+visualizer.run(1)

src/core/collision.cpp

@@ -544,7 +544,7 @@ void handle_collisions() {
       // that see both particles
 
       // If we cannot access both particles, both are ghosts,
-      // ore one is ghost and one is not accessible
+      // or one is ghost and one is not accessible
       // we only increase the counter for the ext id to use based on the
       // number of particles created by other nodes
       if (((!p1 or p1->l.ghost) and (!p2 or p2->l.ghost)) or !p1 or !p2) {

src/core/energy.hpp

@@ -55,7 +55,7 @@ void master_energy_calc();
    all nodes. */
 void energy_calc(double *result, double time);
 
-/** Calculate long range energies (P3M, MMM2d...). */
+/** Calculate long range energies (P3M, ...). */
 void calc_long_range_energies(const ParticleRange &particles);
 
 /** Calculate the total energy */

src/core/forces.hpp

@@ -63,7 +63,7 @@ void init_forces_ghosts(const ParticleRange &particles);
  */
 void force_calc(CellStructure &cell_structure);
 
-/** Calculate long range forces (P3M, MMM2d...). */
+/** Calculate long range forces (P3M, ...). */
 void calc_long_range_forces(const ParticleRange &particles);
 /*@}*/
 

src/core/pressure.cpp

@@ -66,7 +66,7 @@ nptiso_struct nptiso = {0.0,
 /* local prototypes                                         */
 /************************************************************/
 
-/** Calculate long range virials (P3M, MMM2d...). */
+/** Calculate long range virials (P3M, ...). */
 void calc_long_range_virials(const ParticleRange &particles);
 
 /** Initializes a virials Observable stat. */