Refactor exception mechanism in P3M/DP3M tuning functions (espressomd…

…#3869)

Fixes espressomd#3868

Description of changes:
- core: correctly propagate an error code in the P3M tuning function to avoid an infinite loop
- python: fix the broken `'metallic'` case of the P3M `epsilon` parameter
- 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

@@ -650,7 +650,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

@@ -215,7 +215,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_velocity=True, have_quaternion=False)
     # or
     system.virtual_sites = VirtualSitesOff()
@@ -440,7 +440,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)
 
@@ -455,7 +455,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))
 
@@ -484,7 +484,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})
 
@@ -511,7 +511,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)
     system.integrator.set_isotropic_npt(ext_pressure=1.0, piston=1.0)
 
@@ -195,7 +195,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:
 
@@ -125,7 +125,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.
@@ -146,7 +146,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)
 
@@ -175,7 +175,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()
 
@@ -257,10 +257,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`.
 
@@ -419,7 +417,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)
     system.integrator.set_isotropic_npt(ext_pressure=1.0, piston=1.0)
 
@@ -441,7 +439,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
@@ -475,7 +473,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)
 
@@ -491,7 +489,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]

src/core/electrostatics_magnetostatics/p3m.cpp

@@ -1804,6 +1804,8 @@ int p3m_adaptive_tune(char **log) {
         p3m_m_time(log, tmp_mesh, cao_min, cao_max, &tmp_cao, r_cut_iL_min,
                    r_cut_iL_max, &tmp_r_cut_iL, &tmp_alpha_L, &tmp_accuracy);
     /* some error occurred during the tuning force evaluation */
+    if (tmp_time == -P3M_TUNE_FAIL)
+      return ES_ERROR;
     /* this mesh does not work at all */
     if (tmp_time < 0.0)
       continue;

src/python/espressomd/electrostatics.pyx

@@ -263,11 +263,11 @@ IF P3M == 1:
                 raise ValueError("P3M accuracy has to be positive")
 
             if self._params["epsilon"] == "metallic":
-                self._params = 0.0
+                self._params["epsilon"] = 0.0
 
-            if not (is_valid_type(self._params["epsilon"], float)
-                    or self._params["epsilon"] == "metallic"):
-                raise ValueError("epsilon should be a double or 'metallic'")
+            check_type_or_throw_except(
+                self._params["epsilon"], 1, float,
+                "epsilon should be a double or 'metallic'")
 
             if not (self._params["inter"] == default_params["inter"]
                     or self._params["inter"] >= 0):
@@ -279,8 +279,7 @@ IF P3M == 1:
 
             if not (self._params["alpha"] == default_params["alpha"]
                     or self._params["alpha"] > 0):
-                raise ValueError(
-                    "alpha should be positive")
+                raise ValueError("alpha should be positive")
 
         def valid_keys(self):
             return ["mesh", "cao", "accuracy", "epsilon", "alpha", "r_cut",
@@ -326,13 +325,15 @@ IF P3M == 1:
 
         def _tune(self):
             set_prefactor(self._params["prefactor"])
+            p3m_set_eps(self._params["epsilon"])
             python_p3m_set_tune_params(self._params["r_cut"],
                                        self._params["mesh"],
                                        self._params["cao"],
                                        -1.0,
                                        self._params["accuracy"],
                                        self._params["inter"])
             resp = python_p3m_adaptive_tune()
+            handle_errors("P3M tuning failed")
             if resp:
                 raise Exception(
                     "failed to tune P3M parameters to required accuracy")
@@ -412,13 +413,12 @@ IF P3M == 1:
                 if not (self._params["accuracy"] >= 0):
                     raise ValueError("P3M accuracy has to be positive")
 
-                # if self._params["epsilon"] == "metallic":
-                #  self._params = 0.0
+                if self._params["epsilon"] == "metallic":
+                    self._params["epsilon"] = 0.0
 
-                if not (is_valid_type(self._params["epsilon"], float)
-                        or self._params["epsilon"] == "metallic"):
-                    raise ValueError(
-                        "epsilon should be a double or 'metallic'")
+                check_type_or_throw_except(
+                    self._params["epsilon"], 1, float,
+                    "epsilon should be a double or 'metallic'")
 
                 if not (self._params["inter"] == default_params["inter"]
                         or self._params["inter"] > 0):
@@ -456,13 +456,15 @@ IF P3M == 1:
 
             def _tune(self):
                 set_prefactor(self._params["prefactor"])
+                p3m_set_eps(self._params["epsilon"])
                 python_p3m_set_tune_params(self._params["r_cut"],
                                            self._params["mesh"],
                                            self._params["cao"],
                                            -1.0,
                                            self._params["accuracy"],
                                            self._params["inter"])
                 resp = python_p3m_adaptive_tune()
+                handle_errors("P3MGPU tuning failed")
                 if resp:
                     raise Exception(
                         "failed to tune P3M parameters to required accuracy")