Refactor exception mechanism in P3M/DP3M tuning functions (#3869)

Fixes #3868

Description of changes:
- core: correctly propagate errors codes in the P3M/DP3M tuning functions to avoid infinite loops
- core: rely on `runtimeErrorMsg()` instead of printing error messages to stderr or char arrays
- python: catch errors from the P3M/DP3M tuning functions in the `tune()` method of the relevant python Actors
- python: fix the broken `'metallic'` case of the P3M `epsilon` parameter
- testsuite: add tests for the `'metallic'` case and for the new exception mechanism
- python: fix the non-metallic epsilon case (epsilon was always set to 0 for `P3M ` and `P3MGPU` in the core since 4.0.0)
- documentation: fix broken code examples in the user guide
- documentation: mention P3M doesn't support non-cubic boxes for non-metallic epsilon

doc/sphinx/advanced_methods.rst

@@ -28,7 +28,7 @@ Several modes are available for different types of binding.
     import espressomd
     from espressomd.interactions import HarmonicBond
 
-    system = espressomd.System()
+    system = espressomd.System(box_l=[1, 1, 1])
     bond_centers = HarmonicBond(k=1000, r_0=<CUTOFF>)
     system.bonded_inter.add(bond_centers)
     system.collision_detection.set_params(mode="bind_centers", distance=<CUTOFF>,

doc/sphinx/constraints.rst

@@ -167,7 +167,7 @@ Available shapes
 Python syntax::
 
     import espressomd from espressomd.shapes import <SHAPE>
-    system = espressomd.System()
+    system = espressomd.System(box_l=[1, 1, 1])
 
 ``<SHAPE>`` can be any of the available shapes.
 

doc/sphinx/electrostatics.rst

@@ -35,20 +35,22 @@ Note that using the electrostatic interaction also requires assigning charges to
 the particles via the particle property
 :py:attr:`espressomd.particle_data.ParticleHandle.q`.
 
-This example shows the general usage of an electrostatic method ``<SOLVER>``.
-All of them need the Bjerrum length and a set of other required parameters.
-First, an instance of the solver is created and only after adding it to the actors
-list, it is activated. Internally the method calls a tuning routine on
-activation to achieve the given accuracy::
+All solvers need a prefactor and a set of other required parameters.
+This example shows the general usage of the electrostatic method ``P3M``.
+An instance of the solver is created and added to the actors list, at which
+point it will be automatically activated. This activation will internally
+call a tuning function to achieve the requested accuracy::
 
     import espressomd
-    from espressomd import electrostatics
+    import espressomd.electrostatics
 
-    system = espressomd.System()
-    solver = electrostatics.<SOLVER>(prefactor=C, <ADDITIONAL REQUIRED PARAMETERS>)
+    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)
     system.actors.add(solver)
 
-where the prefactor :math:`C` is defined as in Eqn. :eq:`coulomb_prefactor`
+where the prefactor is defined as :math:`C` in Eqn. :eq:`coulomb_prefactor`.
 
 .. _Coulomb P3M:
 
@@ -60,9 +62,10 @@ Coulomb P3M
 For this feature to work, you need to have the ``fftw3`` library
 installed on your system. In |es|, you can check if it is compiled in by
 checking for the feature ``FFTW`` with ``espressomd.features``.
-P3M requires full periodicity (1 1 1). Make sure that you know the relevance of the
-P3M parameters before using P3M! If you are not sure, read the following
-references:
+P3M requires full periodicity (1 1 1). When using a non-metallic dielectric
+constant (``epsilon != 0.0``), the box must be cubic.
+Make sure that you know the relevance of the P3M parameters before using P3M!
+If you are not sure, read the following references:
 :cite:`ewald21,hockney88,kolafa92,deserno98a,deserno98b,deserno00,deserno00a,cerda08d`.
 
 .. _Tuning Coulomb P3M:

doc/sphinx/inter_non-bonded.rst

@@ -620,7 +620,7 @@ that the Thole correction acts between all dipoles, intra- and intermolecular.
 Again, the accuracy is related to the P3M accuracy and the split between
 short-range and long-range electrostatics interaction. It is configured by::
 
-    system = espressomd.System()
+    system = espressomd.System(box_l=[1, 1, 1])
     system.non_bonded_inter[type_1,type_2].thole.set_params(scaling_coeff=<float>, q1q2=<float>)
 
 with parameters:

doc/sphinx/io.rst

@@ -211,7 +211,7 @@ capabilities. The usage is quite simple:
 .. code:: python
 
     from espressomd.io.mppiio import mpiio
-    system = espressomd.System()
+    system = espressomd.System(box_l=[1, 1, 1])
     # ... add particles here
     mpiio.write("/tmp/mydata", positions=True, velocities=True, types=True, bonds=True)
 

doc/sphinx/magnetostatics.rst

@@ -27,9 +27,15 @@ relative permittivity of the background material, respectively.
 
 Magnetostatic interactions are activated via the actor framework::
 
-    from espressomd.magnetostatics import DipolarDirectSumCpu
+    import espressomd
+    import espressomd.magnetostatics
 
-    direct_sum = DipolarDirectSumCpu(prefactor=1)
+    system = espressomd.System(box_l=[10, 10, 10])
+    system.time_step = 0.01
+    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)
@@ -49,7 +55,8 @@ Dipolar P3M
 This is the dipolar version of the P3M algorithm, described in :cite:`cerda08d`.
 
 Make sure that you know the relevance of the P3M parameters before using
-P3M! If you are not sure, read the following references
+P3M! If you are not sure, read the following references:
+:cite:`ewald21,hockney88,kolafa92,deserno98a,deserno98b,deserno00,deserno00a,cerda08d`.
 
 Note that dipolar P3M does not work with non-cubic boxes.
 

doc/sphinx/particles.rst

@@ -25,7 +25,7 @@ In order to add particles to the system, call
 :meth:`espressomd.particle_data.ParticleList.add`::
 
     import espressomd
-    system = espressomd.System()
+    system = espressomd.System(box_l=[1, 1, 1])
     system.part.add(pos=[1.0, 1.0, 1.0], id=0, type=0)
 
 This command adds a single particle to the system with properties given
@@ -308,7 +308,7 @@ To switch the active scheme, the attribute :attr:`espressomd.system.System.virtu
     import espressomd
     from espressomd.virtual_sites import VirtualSitesOff, VirtualSitesRelative
 
-    system = espressomd.System()
+    system = espressomd.System(box_l=[1, 1, 1])
     system.virtual_sites = VirtualSitesRelative(have_quaternion=False)
     # or
     system.virtual_sites = VirtualSitesOff()
@@ -439,7 +439,7 @@ Particle ids can be stored in a map for each
 individual type::
 
     import espressomd
-    system = espressomd.System()
+    system = espressomd.System(box_l=[1, 1, 1])
     system.setup_type_map([_type])
     system.number_of_particles(_type)
 
@@ -454,7 +454,7 @@ The keyword
 particles which have the given type. For counting the number of particles of a given type you could also use :meth:`espressomd.particle_data.ParticleList.select` ::
 
     import espressomd
-    system = espressomd.System()
+    system = espressomd.System(box_l=[1, 1, 1])
     ...
     number_of_particles = len(system.part.select(type=type))
 
@@ -483,7 +483,7 @@ Langevin swimmers
 
     import espressomd
 
-    system = espressomd.System()
+    system = espressomd.System(box_l=[1, 1, 1])
 
     system.part.add(id=0, pos=[1, 0, 0], swimming={'f_swim': 0.03})
 
@@ -510,7 +510,7 @@ Lattice-Boltzmann (LB) swimmers
 
     import espressomd
 
-    system = espressomd.System()
+    system = espressomd.System(box_l=[1, 1, 1])
 
     system.part.add(id=1, pos=[2, 0, 0], rotation=[1, 1, 1], swimming={
         'f_swim': 0.01, 'mode': 'pusher', 'dipole_length': 2.0})

doc/sphinx/running.rst

@@ -105,7 +105,7 @@ A code snippet would look like::
 
     import espressomd
 
-    system = espressomd.System()
+    system = espressomd.System(box_l=[1, 1, 1])
     system.thermostat.set_npt(kT=1.0, gamma0=1.0, gammav=1.0, seed=42)
     system.integrator.set_isotropic_npt(ext_pressure=1.0, piston=1.0)
 
@@ -204,7 +204,7 @@ In this way, just compiling in the ``ROTATION`` feature no longer changes the ph
 The rotation of a particle is controlled via the :attr:`espressomd.particle_data.ParticleHandle.rotation` property. E.g., the following code adds a particle with rotation enabled on the x axis::
 
     import espressomd
-    system = espressomd.System()
+    system = espressomd.System(box_l=[1, 1, 1])
     system.part.add(pos=(0, 0, 0), rotation=(1, 0, 0))
 
 Notes:

doc/sphinx/system_manipulation.rst

@@ -96,7 +96,7 @@ Galilei Transform and Particle Velocity Manipulation
 The following class :class:`espressomd.galilei.GalileiTransform` may be useful
 in affecting the velocity of the system. ::
 
-    system = espressomd.System()
+    system = espressomd.System(box_l=[1, 1, 1])
     gt = system.galilei
 
 * Particle motion and rotation

doc/sphinx/system_setup.rst

@@ -71,17 +71,17 @@ Accessing module states
 
 Some variables like or are no longer directly available as attributes.
 In these cases they can be easily derived from the corresponding Python
-objects like
+objects like::
 
-``n_part = len(espressomd.System().part[:].pos)``
+    n_part = len(system.part[:].pos)
 
 or by calling the corresponding ``get_state()`` methods like::
 
-    temperature = espressomd.System().thermostat.get_state()[0]['kT']
+    temperature = system.thermostat.get_state()[0]['kT']
 
-    gamma = espressomd.System().thermostat.get_state()[0]['gamma']
+    gamma = system.thermostat.get_state()[0]['gamma']
 
-    gamma_rot = espressomd.System().thermostat.get_state()[0]['gamma_rotation']
+    gamma_rot = system.thermostat.get_state()[0]['gamma_rotation']
 
 .. _Cellsystems:
 
@@ -112,7 +112,7 @@ The properties of the cell system can be accessed by
 
     (float) Skin for the Verlet list. This value has to be set, otherwise the simulation will not start.
 
-Details about the cell system can be obtained by :meth:`espressomd.System().cell_system.get_state() <espressomd.cellsystem.CellSystem.get_state>`:
+Details about the cell system can be obtained by :meth:`espressomd.system.System.cell_system.get_state() <espressomd.cellsystem.CellSystem.get_state>`:
 
     * ``cell_grid``       Dimension of the inner cell grid.
     * ``cell_size``       Box-length of a cell.
@@ -132,7 +132,7 @@ selects the domain decomposition cell scheme, using Verlet lists
 for the calculation of the interactions. If you specify ``use_verlet_lists=False``, only the
 domain decomposition is used, but not the Verlet lists. ::
 
-    system = espressomd.System()
+    system = espressomd.System(box_l=[1, 1, 1])
 
     system.cell_system.set_domain_decomposition(use_verlet_lists=True)
 
@@ -161,7 +161,7 @@ particles, giving an unfavorable computation time scaling of
 interaction in the cell model require the calculation of all pair
 interactions. ::
 
-    system = espressomd.System()
+    system = espressomd.System(box_l=[1, 1, 1])
 
     system.cell_system.set_n_square()
 
@@ -221,10 +221,8 @@ class :class:`espressomd.thermostat.Thermostat` has to be invoked.
 Best explained in an example::
 
     import espressomd
-    system = espressomd.System()
-    therm = system.Thermostat()
-
-    therm.set_langevin(kT=1.0, gamma=1.0, seed=41)
+    system = espressomd.System(box_l=[1, 1, 1])
+    system.thermostat.set_langevin(kT=1.0, gamma=1.0, seed=41)
 
 As explained before the temperature is set as thermal energy :math:`k_\mathrm{B} T`.
 
@@ -383,7 +381,7 @@ For example::
 
     import espressomd
 
-    system = espressomd.System()
+    system = espressomd.System(box_l=[1, 1, 1])
     system.thermostat.set_npt(kT=1.0, gamma0=1.0, gammav=1.0, seed=41)
     system.integrator.set_isotropic_npt(ext_pressure=1.0, piston=1.0)
 
@@ -416,7 +414,7 @@ The system integrator should be also changed.
 Best explained in an example::
 
     import espressomd
-    system = espressomd.System()
+    system = espressomd.System(box_l=[1, 1, 1])
     system.thermostat.set_brownian(kT=1.0, gamma=1.0, seed=41)
     system.integrator.set_brownian_dynamics()
 
@@ -505,7 +503,7 @@ example ``lbfluid``. If you do not choose the GPU manually before that,
 CUDA internally chooses one, which is normally the most powerful GPU
 available, but load-independent. ::
 
-    system = espressomd.System()
+    system = espressomd.System(box_l=[1, 1, 1])
 
     dev = system.cuda_init_handle.device
     system.cuda_init_handle.device = dev
@@ -539,7 +537,7 @@ List available CUDA devices
 If you want to list available CUDA devices
 you should access :attr:`espressomd.cuda_init.CudaInitHandle.device_list`, e.g., ::
 
-    system = espressomd.System()
+    system = espressomd.System(box_l=[1, 1, 1])
 
     print(system.cuda_init_handle.device_list)
 
@@ -555,7 +553,7 @@ When you start ``pypresso`` your first GPU should be selected.
 If you wanted to use the second GPU, this can be done
 by setting :attr:`espressomd.cuda_init.CudaInitHandle.device` as follows::
 
-    system = espressomd.System()
+    system = espressomd.System(box_l=[1, 1, 1])
 
     system.cuda_init_handle.device = 1
 

doc/sphinx/visualization.rst

@@ -42,7 +42,7 @@ window with ``start()``. See the following minimal code example::
     from espressomd import visualization
     from threading import Thread
 
-    system = espressomd.System()
+    system = espressomd.System(box_l=[1, 1, 1])
     system.cell_system.skin = 0.4
     system.time_step = 0.01
     system.box_l = [10, 10, 10]