Rewrite long-range actors (#4506)

Fixes #4454

Description of changes:
- apply the SOLID principles
   - split electrostatics from magnetostatics
   - factor out code duplication
   - extract P3M/DP3M tuning code
   - hide ELC details in private class methods (the energy/force correction logic used to be exposed in `coulomb.cpp`)
   - hide Coulomb implementation details by no longer indirectly including `coulomb.hpp` in 30 unrelated files in the core
- write clear pre-conditions as event hooks
   - these events are triggered upon box size change, periodicity change, node grid change and cell structure change
      - partial fix for #4428
      - setting a `node_grid` incompatible with a long-range actor is no longer possible (the old `node_grid` is restored)
   - these preconditions always run once before the integration loop starts
   - write test cases to check for these events
   - bugfix: charge neutrality checks now run once before the integration loop starts, instead of once at actor activation
- encapsulate the state of long-range methods in classes managed by the script interface
   - bugfix: the cached `DipolarP3M` energy correction is now invalidated when the box length changes
   - remove global variables for CPU methods (partial fix for #2628) and introduce one global variable for each type of long-range actor: dipolar solver, Coulomb solver, Coulomb extension
  - convert Cython files for long-range actors to Python files using the `ScriptInterface` framework
     - long-range actors and their python object now have the same lifetime, regardless of whether the actor is active or not (partial fix for #4220)
  - the list of active long-range actors can no longer end up in an invalid state (partial fix for #4219):
     - updating an active actor with an invalid parameter now automatically restores the original parameter
     - inserting an invalid long-range actor in the list of actors now automatically removes the actor from the list of active actors, even when the exception happens during tuning
- API changes:
   - `ELC` and `DLC` are no longer extensions, instead they take another actor as argument
   - `P3MGPU` is now able to execute the CPU energy and pressure kernels using parameters tuned for the GPU force kernel
      - closes #239
      - P3M-based reaction scripts run with 2x speed-up when reaction steps take place every 1000 integration steps
   - `P3MGPU` is now supported by `ELC`
   - long-range actor parameters are now immutable, since changing one parameter (e.g. the prefactor) would usually invalidate the tuned parameters
   - the charge neutrality check sensitivity can now be modified via `actor.charge_neutrality_tolerance`; this is mostly relevant to actors coupled to `ICC`, whose induced charges can have value spanning several orders of magnitude
- under-the-hood changes:
   - CPU-based long-range actors no longer rely on the `MpiCallbacks` framework
   - the subtracted Coulomb bonded interaction and pair criteria no longer indirectly include `coulomb.hpp`

doc/sphinx/electrostatics.rst

@@ -47,11 +47,16 @@ call a tuning function to achieve the requested accuracy::
     system = espressomd.System(box_l=[10, 10, 10])
     system.time_step = 0.01
     system.part.add(pos=[[0, 0, 0], [1, 1, 1]], q=[-1, 1])
-    solver = espressomd.electrostatics.P3M(prefactor=2, accuracy=1e-3)
+    solver = espressomd.electrostatics.P3M(prefactor=2., accuracy=1e-3)
     system.actors.add(solver)
 
 where the prefactor is defined as :math:`C` in Eqn. :eq:`coulomb_prefactor`.
 
+The list of actors can be cleared with
+:meth:`system.actors.clear() <espressomd.actors.Actors.clear>` and
+:meth:`system.actors.remove(actor) <espressomd.actors.Actors.remove>`.
+
+
 .. _Coulomb P3M:
 
 Coulomb P3M
@@ -73,27 +78,16 @@ If you are not sure, read the following references:
 Tuning Coulomb P3M
 ~~~~~~~~~~~~~~~~~~
 
-The tuning method is called when the handle of the Coulomb P3M is added to the
-actor list. At this point, the system should already contain the charged
-particles. Set parameters are fixed and not changed by the tuning algorithm.
-This can be useful to speed up the tuning during testing or if the parameters
-are already known.
-
-To prevent the automatic tuning, set the ``tune`` parameter to ``False``.
-To manually tune or retune P3M, call :meth:`espressomd.electrostatics.P3M.tune
-<espressomd.electrostatics.ElectrostaticInteraction.tune>`.
-Note, however, that this is a method the P3M object inherited from
-:attr:`espressomd.electrostatics.ElectrostaticInteraction`.
-All parameters passed to the method are fixed in the tuning routine. If not
-specified in the ``tune()`` method, the parameters ``prefactor`` and
-``accuracy`` are reused.
-
 It is not easy to calculate the various parameters of the P3M method
 such that the method provides the desired accuracy at maximum speed. To
 simplify this, it provides a function to automatically tune the algorithm.
 Note that for this function to work properly, your system should already
-contain an initial configuration of charges and the correct initial box
-size. Also note that the provided tuning algorithms works very well on
+contain an initial configuration of charges and the correct initial box size.
+The tuning method is called when the handle of the Coulomb P3M is added to
+the actor list. Some parameters can be fixed (``r_cut``, ``cao``, ``mesh``)
+to speed up the tuning if the parameters are already known.
+
+Please note that the provided tuning algorithms works very well on
 homogeneous charge distributions, but might not achieve the requested
 precision for highly inhomogeneous or symmetric systems. For example,
 because of the nature of the P3M algorithm, systems are problematic
@@ -106,7 +100,7 @@ obtain sets of parameters that yield the desired accuracy, then it measures how
 long it takes to compute the Coulomb interaction using these parameter sets and
 chooses the set with the shortest run time.
 
-After execution the tuning routines report the tested parameter sets,
+During tuning, the algorithm reports the tested parameter sets,
 the corresponding k-space and real-space errors and the timings needed
 for force calculations. In the output, the timings are given in units of
 milliseconds, length scales are in units of inverse box lengths.
@@ -118,14 +112,16 @@ Coulomb P3M on GPU
 
 :class:`espressomd.electrostatics.P3MGPU`
 
-The GPU implementation of P3M calculates the far field portion on the GPU.
-It uses the same parameters and interface functionality as the CPU version of
-the solver. It should be noted that this does not always provide significant
+The GPU implementation of P3M calculates the far field contribution to the
+forces on the GPU. The near-field contribution to the forces, as well as the
+near- and far-field contributions to the energies are calculated on the CPU.
+It uses the same parameters
+and interface functionality as the CPU version of the solver.
+It should be noted that this does not always provide significant
 increase in performance. Furthermore it computes the far field interactions
-with only single precision which limits the maximum precision. The algorithm
-does not work in combination with the electrostatic extensions
-:ref:`Dielectric interfaces with the ICC* algorithm <Dielectric interfaces with the ICC algorithm>`
-and :ref:`Electrostatic Layer Correction (ELC)`.
+with only single precision which limits the maximum precision.
+The algorithm does not work in combination with the electrostatic extension
+:ref:`Dielectric interfaces with the ICC* algorithm <Dielectric interfaces with the ICC algorithm>`.
 
 The algorithm doesn't have kernels to compute energies, and will therefore
 contribute 0 to the long-range potential energy of the system. This can be
@@ -269,8 +265,6 @@ Electrostatic Layer Correction (ELC)
 systems. It can account for different dielectric jumps on both sides of the
 non-periodic direction. In more detail, it is a special procedure that
 converts a 3D electrostatic method to a 2D method in computational order N.
-Currently, it only supports P3M without GPU. This means,
-that you will first have to set up the P3M algorithm before using ELC.
 The periodicity has to be set to (1 1 1). *ELC* cancels the electrostatic
 contribution of the periodic replica in **z-direction**. Make sure that you
 read the papers on ELC (:cite:`arnold02c,arnold02d,tyagi08a`) before using it.
@@ -291,9 +285,15 @@ Usage notes:
 
     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)
+    elc = espressomd.electrostatics.ELC(actor=p3m, gap_size=box_l * 0.2, maxPWerror=1e-3)
     system.actors.add(elc)
 
+Although it is technically feasible to remove ``elc`` from the list of actors
+and then to add the ``p3m`` object, it is not recommended because the P3M
+parameters are mutated by *ELC*, e.g. the ``epsilon`` is made metallic.
+It is safer to instantiate a new P3M object instead of recycling one that
+has been adapted by *ELC*.
+
 *ELC* can also be used to simulate 2D periodic systems with image charges,
 specified by dielectric contrasts on the non-periodic boundaries
 (:cite:`tyagi08a`). This is achieved by setting the dielectric jump from the
@@ -303,7 +303,7 @@ 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 = espressomd.electrostatics.ELC(p3m_actor=p3m, gap_size=box_l * 0.2, maxPWerror=1e-3,
+    elc = espressomd.electrostatics.ELC(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
@@ -313,7 +313,7 @@ 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 = espressomd.electrostatics.ELC(p3m_actor=p3m, gap_size=box_l * 0.2, maxPWerror=1e-3,
+    elc = espressomd.electrostatics.ELC(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
@@ -396,7 +396,8 @@ ScaFaCoS electrostatics
 
 |es| can use the methods from the ScaFaCoS *Scalable fast Coulomb solvers*
 library. The specific methods available depend on the compile-time options of
-the library, and can be queried using :meth:`espressomd.scafacos.available_methods`.
+the library, and can be queried with
+:meth:`espressomd.electrostatics.Scafacos.get_available_methods`.
 
 To use ScaFaCoS, create an instance of :class:`~espressomd.electrostatics.Scafacos`
 and add it to the list of active actors. Three parameters have to be specified:

doc/sphinx/integration.rst

@@ -303,6 +303,9 @@ The algorithm can also be used for energy minimization::
     system.integrator.set_vv()
 
 Please note that not all features support energy calculation.
+For example :ref:`IBM <Immersed Boundary Method for soft elastic objects>`
+and :ref:`OIF <Object-in-fluid>` do not implement energy calculation for
+mesh surface deformation.
 
 .. _Brownian Dynamics:
 

doc/sphinx/magnetostatics.rst

@@ -35,14 +35,12 @@ Magnetostatic interactions are activated via the actor framework::
     system.part.add(pos=[[0, 0, 0], [1, 1, 1]],
                     rotation=2 * [(1, 1, 1)], dip=2 * [(1, 0, 0)])
 
-    direct_sum = espressomd.magnetostatics.DipolarDirectSumCpu(prefactor=1)
-    system.actors.add(direct_sum)
-    # ...
-    system.actors.remove(direct_sum)
+    actor = espressomd.magnetostatics.DipolarDirectSumCpu(prefactor=1.)
+    system.actors.add(actor)
 
-The magnetostatics algorithms for periodic boundary conditions require
-some knowledge to use them properly. Uneducated use can result in
-completely unphysical simulations.
+The list of actors can be cleared with
+:meth:`system.actors.clear() <espressomd.actors.Actors.clear>` and
+:meth:`system.actors.remove(actor) <espressomd.actors.Actors.remove>`.
 
 
 .. _Dipolar P3M:
@@ -81,7 +79,7 @@ simulation, actual force and torque errors can be significantly larger.
 Dipolar Layer Correction (DLC)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-:class:`espressomd.magnetostatic_extensions.DLC`
+:class:`espressomd.magnetostatics.DLC`
 
 The dipolar layer correction (DLC) is used in conjunction with the dipolar P3M
 method to calculate dipolar interactions in a 2D-periodic system.
@@ -104,15 +102,14 @@ Usage notes:
   assumed that all dipole moment are as large as the largest of the dipoles
   in the system.
 
-The method is used as follows::
+* When the base solver is not a P3M method, metallic epsilon is assumed.
 
-    import espressomd.magnetostatics as magnetostatics
-    import espressomd.magnetostatic_extensions as magnetostatic_extensions
+The method is used as follows::
 
-    p3m = magnetostatics.DipolarP3M(prefactor=1, accuracy=1E-4)
-    dlc = magnetostatic_extensions.DLC(maxPWerror=1E-5, gap_size=2.)
-    system.actors.add(p3m)
-    system.actors.add(dlc)
+    import espressomd.magnetostatics
+    dp3m = espressomd.magnetostatics.DipolarP3M(prefactor=1, accuracy=1E-4)
+    mdlc = espressomd.magnetostatics.DLC(actor=dp3m, maxPWerror=1E-5, gap_size=2.)
+    system.actors.add(mdlc)
 
 
 .. _Dipolar direct sum:
@@ -149,9 +146,9 @@ To use the methods, create an instance of either
 system's list of active actors. The only required parameter is the Prefactor
 :eq:`dipolar_prefactor`::
 
-  import espressomd.magnetostatics
-  dds = espressomd.magnetostatics.DipolarDirectSumGpu(bjerrum_length=1)
-  system.actors.add(dds)
+    import espressomd.magnetostatics
+    dds = espressomd.magnetostatics.DipolarDirectSumGpu(prefactor=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:
@@ -184,9 +181,9 @@ refer to :cite:`polyakov13a`.
 To use the method, create an instance of :class:`~espressomd.magnetostatics.DipolarBarnesHutGpu`
 and add it to the system's list of active actors::
 
-  import espressomd.magnetostatics
-  bh = espressomd.magnetostatics.DipolarBarnesHutGpu(prefactor=pf_dds_gpu, epssq=200.0, itolsq=8.0)
-  system.actors.add(bh)
+    import espressomd.magnetostatics
+    bh = espressomd.magnetostatics.DipolarBarnesHutGpu(prefactor=1., epssq=200.0, itolsq=8.0)
+    system.actors.add(bh)
 
 
 .. _ScaFaCoS magnetostatics:
@@ -199,8 +196,9 @@ ScaFaCoS magnetostatics
 |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
-feature. Dipolar calculations are only included in the ``dipolar`` branch of
-the ScaFaCoS code.
+feature. Dipolar calculations are only included in the ``dipoles`` branch of
+the ScaFaCoS code. The specific methods available can be queried with
+:meth:`espressomd.electrostatics.Scafacos.get_available_methods`.
 
 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:

doc/sphinx/reaction_methods.rst

@@ -16,7 +16,8 @@ Please keep in mind the following remarks:
 * All reaction methods uses Monte Carlo moves which require potential energies.
   Therefore reaction methods require support for energy calculations for all
   active interactions in the simulation. Some algorithms do not support energy
-  calculation, e.g. :ref:`Coulomb P3M on GPU`.
+  calculation, e.g. :ref:`OIF <Object-in-fluid>` and
+  :ref:`IBM <Immersed Boundary Method for soft elastic objects>`.
 
 * When modeling reactions that do not conserve the number of particles, the
   method has to create or delete particles from the system. This process can

doc/tutorials/charged_system/charged_system.ipynb

@@ -39,7 +39,7 @@
     "import espressomd.accumulators\n",
     "import espressomd.math\n",
     "\n",
-    "espressomd.assert_features(['WCA', 'ELECTROSTATICS'])\n",
+    "espressomd.assert_features(['ELECTROSTATICS', 'P3M', 'WCA'])\n",
     "\n",
     "import numpy as np\n",
     "import scipy.optimize\n",
@@ -503,11 +503,15 @@
     "* Write a function ``clear_system(system)`` that\n",
     "    * turns off the thermostat\n",
     "    * removes all particles\n",
+    "    * removes all actors\n",
     "    * removes all accumulators added to the auto-update-list\n",
     "    * resets the system clock\n",
     "\n",
     "**Hints:**\n",
-    "* The relevant parts of the documentation can be found here: [Thermostats](https://espressomd.github.io/doc/integration.html#thermostats), [ParticleList](https://espressomd.github.io/doc/espressomd.html#espressomd.particle_data.ParticleList),\n",
+    "* The relevant parts of the documentation can be found here:\n",
+    "[Thermostats](https://espressomd.github.io/doc/integration.html#thermostats),\n",
+    "[ParticleList](https://espressomd.github.io/doc/espressomd.html#espressomd.particle_data.ParticleList),\n",
+    "[Electrostatics](https://espressomd.github.io/doc/electrostatics.html),\n",
     "[AutoUpdateAccumulators](https://espressomd.github.io/doc/espressomd.html#espressomd.accumulators.AutoUpdateAccumulators),\n",
     "[System properties](https://espressomd.github.io/doc/espressomd.html#module-espressomd.system)"
    ]
@@ -522,6 +526,7 @@
     "def clear_system(system):\n",
     "    system.thermostat.turn_off()\n",
     "    system.part.clear()\n",
+    "    system.actors.clear()\n",
     "    system.auto_update_accumulators.clear()\n",
     "    system.time = 0.\n",
     "```"
@@ -597,7 +602,7 @@
     "\n",
     "**Hints:**\n",
     "* Don't forget to clear the system before setting up the system with a new set of parameters\n",
-    "* Don't forget to ``tune()`` the ``p3m`` instance after each change of parameters. If we reuse the p3m that was tuned before, likely the desired accuracy will not be achieved. \n",
+    "* Don't forget to add a new ``p3m`` instance after each change of parameters. If we reuse the p3m that was tuned before, likely the desired accuracy will not be achieved. \n",
     "* Extract the radial density profile from the accumulator via ``.mean()``"
    ]
   },
@@ -613,7 +618,8 @@
     "    setup_rod_and_counterions(\n",
     "        system, run['params']['counterion_valency'], COUNTERION_TYPE,\n",
     "        run['params']['rod_charge_dens'], N_rod_beads, ROD_TYPE)\n",
-    "    p3m.tune()\n",
+    "    p3m = espressomd.electrostatics.P3M(**p3m_params)\n",
+    "    system.actors.add(p3m)\n",
     "    remove_overlap(system, STEEPEST_DESCENT_PARAMS)\n",
     "    system.thermostat.set_langevin(**LANGEVIN_PARAMS)\n",
     "    print('starting warmup')\n",
@@ -840,7 +846,8 @@
     "anions, cations = add_salt(system, ANION_PARAMS, CATION_PARAMS)\n",
     "assert abs(sum(anions.q) + sum(cations.q)) < 1e-10\n",
     "\n",
-    "p3m.tune()\n",
+    "p3m = espressomd.electrostatics.P3M(**p3m_params)\n",
+    "system.actors.add(p3m)\n",
     "remove_overlap(system, STEEPEST_DESCENT_PARAMS)\n",
     "system.thermostat.set_langevin(**LANGEVIN_PARAMS)\n",
     "print('starting warmup')\n",

doc/tutorials/ferrofluid/ferrofluid_part1.ipynb

@@ -235,11 +235,10 @@
    "source": [
     "import espressomd\n",
     "import espressomd.magnetostatics\n",
-    "import espressomd.magnetostatic_extensions\n",
     "import espressomd.cluster_analysis\n",
     "import espressomd.pair_criteria\n",
     "\n",
-    "espressomd.assert_features('DIPOLES', 'LENNARD_JONES')\n",
+    "espressomd.assert_features('DIPOLES', 'DP3M', 'LENNARD_JONES')\n",
     "\n",
     "import numpy as np"
    ]
@@ -553,9 +552,8 @@
     "CI_DP3M_PARAMS = {} # debug variable for continuous integration, can be left empty\n",
     "# Setup dipolar P3M and dipolar layer correction\n",
     "dp3m = espressomd.magnetostatics.DipolarP3M(accuracy=5E-4, prefactor=DIP_LAMBDA * LJ_SIGMA**3 * KT, **CI_DP3M_PARAMS)\n",
-    "dlc = espressomd.magnetostatic_extensions.DLC(maxPWerror=1E-4, gap_size=BOX_SIZE - LJ_SIGMA)\n",
-    "system.actors.add(dp3m)\n",
-    "system.actors.add(dlc)\n",
+    "mdlc = espressomd.magnetostatics.DLC(actor=dp3m, maxPWerror=1E-4, gap_size=BOX_SIZE - LJ_SIGMA)\n",
+    "system.actors.add(mdlc)\n",
     "\n",
     "# tune verlet list skin\n",
     "system.cell_system.tune_skin(min_skin=0.4, max_skin=2., tol=0.2, int_steps=100)\n",

doc/tutorials/ferrofluid/ferrofluid_part2.ipynb

@@ -52,9 +52,8 @@
    "source": [
     "import espressomd\n",
     "import espressomd.magnetostatics\n",
-    "import espressomd.magnetostatic_extensions\n",
     "\n",
-    "espressomd.assert_features('DIPOLES', 'LENNARD_JONES')\n",
+    "espressomd.assert_features('DIPOLES', 'DP3M', 'LENNARD_JONES')\n",
     "\n",
     "import numpy as np"
    ]
@@ -170,9 +169,8 @@
     "\n",
     "# Setup dipolar P3M and dipolar layer correction (DLC)\n",
     "dp3m = espressomd.magnetostatics.DipolarP3M(accuracy=5E-4, prefactor=DIP_LAMBDA * LJ_SIGMA**3 * KT)\n",
-    "dlc = espressomd.magnetostatic_extensions.DLC(maxPWerror=1E-4, gap_size=box_size - LJ_SIGMA)\n",
-    "system.actors.add(dp3m)\n",
-    "system.actors.add(dlc)\n",
+    "mdlc = espressomd.magnetostatics.DLC(actor=dp3m, maxPWerror=1E-4, gap_size=box_size - LJ_SIGMA)\n",
+    "system.actors.add(mdlc)\n",
     "\n",
     "# tune verlet list skin again\n",
     "system.cell_system.tune_skin(min_skin=0.4, max_skin=2., tol=0.2, int_steps=100)\n",
@@ -499,6 +497,7 @@
    "source": [
     "# remove all particles\n",
     "system.part.clear()\n",
+    "system.actors.clear()\n",
     "system.thermostat.turn_off()\n",
     "\n",
     "# Random dipole moments\n",
@@ -530,14 +529,9 @@
     "system.cell_system.tune_skin(min_skin=0.4, max_skin=2., tol=0.2, int_steps=100)\n",
     "\n",
     "# Setup dipolar P3M and dipolar layer correction\n",
-    "system.actors.remove(dp3m)\n",
-    "system.actors.remove(dlc)\n",
-    "\n",
     "dp3m = espressomd.magnetostatics.DipolarP3M(accuracy=5E-4, prefactor=DIP_LAMBDA * LJ_SIGMA**3 * KT)\n",
-    "dlc = espressomd.magnetostatic_extensions.DLC(maxPWerror=1E-4, gap_size=box_size - LJ_SIGMA)\n",
-    "\n",
-    "system.actors.add(dp3m)\n",
-    "system.actors.add(dlc)\n",
+    "mdlc = espressomd.magnetostatics.DLC(actor=dp3m, maxPWerror=1E-4, gap_size=box_size - LJ_SIGMA)\n",
+    "system.actors.add(mdlc)\n",
     "\n",
     "# tune verlet list skin again\n",
     "system.cell_system.tune_skin(min_skin=0.4, max_skin=2., tol=0.2, int_steps=100)"

doc/tutorials/ferrofluid/ferrofluid_part3.ipynb

@@ -124,7 +124,7 @@
     "import espressomd\n",
     "import espressomd.magnetostatics\n",
     "\n",
-    "espressomd.assert_features('DIPOLES', 'LENNARD_JONES')\n",
+    "espressomd.assert_features('DIPOLES', 'DP3M', 'LENNARD_JONES')\n",
     "\n",
     "import numpy as np"
    ]