From 28b67232d41fd9632d699e3e4852315b42791750 Mon Sep 17 00:00:00 2001 From: "kodiakhq[bot]" <49736102+kodiakhq[bot]@users.noreply.github.com> Date: Mon, 24 Aug 2020 09:10:39 +0000 Subject: [PATCH] Refactor exception mechanism in P3M/DP3M tuning functions (#3869) Fixes #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 | 2 +- doc/sphinx/constraints.rst | 2 +- doc/sphinx/electrostatics.rst | 27 +++++++++-------- doc/sphinx/inter_non-bonded.rst | 2 +- doc/sphinx/io.rst | 2 +- doc/sphinx/magnetostatics.rst | 13 ++++++-- doc/sphinx/particles.rst | 12 ++++---- doc/sphinx/running.rst | 4 +-- doc/sphinx/system_manipulation.rst | 2 +- doc/sphinx/system_setup.rst | 30 +++++++++---------- doc/sphinx/visualization.rst | 2 +- .../electrostatics_magnetostatics/p3m.cpp | 2 ++ src/python/espressomd/electrostatics.pyx | 26 ++++++++-------- src/python/espressomd/magnetostatics.pyx | 13 ++++---- testsuite/python/electrostaticInteractions.py | 8 ++--- 15 files changed, 78 insertions(+), 69 deletions(-) diff --git a/doc/sphinx/advanced_methods.rst b/doc/sphinx/advanced_methods.rst index 2d8e42b364d..c6e6ea8fbcc 100644 --- a/doc/sphinx/advanced_methods.rst +++ b/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=) system.bonded_inter.add(bond_centers) system.collision_detection.set_params(mode="bind_centers", distance=, diff --git a/doc/sphinx/constraints.rst b/doc/sphinx/constraints.rst index e841d44552e..69eedff49c5 100644 --- a/doc/sphinx/constraints.rst +++ b/doc/sphinx/constraints.rst @@ -167,7 +167,7 @@ Available shapes Python syntax:: import espressomd from espressomd.shapes import - system = espressomd.System() + system = espressomd.System(box_l=[1, 1, 1]) ```` can be any of the available shapes. diff --git a/doc/sphinx/electrostatics.rst b/doc/sphinx/electrostatics.rst index 48abd4f314b..8d2ab952389 100644 --- a/doc/sphinx/electrostatics.rst +++ b/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 ````. -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.(prefactor=C, ) + 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: diff --git a/doc/sphinx/inter_non-bonded.rst b/doc/sphinx/inter_non-bonded.rst index ce83de98bf3..16e16393989 100644 --- a/doc/sphinx/inter_non-bonded.rst +++ b/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=, q1q2=) with parameters: diff --git a/doc/sphinx/io.rst b/doc/sphinx/io.rst index 7b99bb52dec..b1da29dd81a 100644 --- a/doc/sphinx/io.rst +++ b/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) diff --git a/doc/sphinx/magnetostatics.rst b/doc/sphinx/magnetostatics.rst index 417a182d030..52fc41df91d 100644 --- a/doc/sphinx/magnetostatics.rst +++ b/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. diff --git a/doc/sphinx/particles.rst b/doc/sphinx/particles.rst index 7f9654ad64c..024cf013778 100644 --- a/doc/sphinx/particles.rst +++ b/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}) diff --git a/doc/sphinx/running.rst b/doc/sphinx/running.rst index aa63fa8b66b..cf0d574c598 100644 --- a/doc/sphinx/running.rst +++ b/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: diff --git a/doc/sphinx/system_manipulation.rst b/doc/sphinx/system_manipulation.rst index 4188de7e1f1..d399a59b2be 100644 --- a/doc/sphinx/system_manipulation.rst +++ b/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 diff --git a/doc/sphinx/system_setup.rst b/doc/sphinx/system_setup.rst index 03a85bac116..ea8d5a9b2d6 100644 --- a/doc/sphinx/system_setup.rst +++ b/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() `: +Details about the cell system can be obtained by :meth:`espressomd.system.System.cell_system.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 diff --git a/doc/sphinx/visualization.rst b/doc/sphinx/visualization.rst index a9b6e752b5d..e3cbe7e7c2c 100644 --- a/doc/sphinx/visualization.rst +++ b/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] diff --git a/src/core/electrostatics_magnetostatics/p3m.cpp b/src/core/electrostatics_magnetostatics/p3m.cpp index 4bc883048bf..210dc596339 100644 --- a/src/core/electrostatics_magnetostatics/p3m.cpp +++ b/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; diff --git a/src/python/espressomd/electrostatics.pyx b/src/python/espressomd/electrostatics.pyx index f8692d50b0a..6d62a604509 100644 --- a/src/python/espressomd/electrostatics.pyx +++ b/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,6 +325,7 @@ 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"], @@ -333,6 +333,7 @@ IF P3M == 1: 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,6 +456,7 @@ 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"], @@ -463,6 +464,7 @@ IF P3M == 1: 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") diff --git a/src/python/espressomd/magnetostatics.pyx b/src/python/espressomd/magnetostatics.pyx index 3d425e3b3a0..198dd180b5d 100644 --- a/src/python/espressomd/magnetostatics.pyx +++ b/src/python/espressomd/magnetostatics.pyx @@ -126,9 +126,9 @@ IF DP3M == 1: 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): @@ -179,9 +179,7 @@ IF DP3M == 1: self._params["r_cut"], self._params["mesh"], self._params["cao"], -1., self._params["accuracy"], self._params["inter"]) resp, log = self.python_dp3m_adaptive_tune() - if resp: - raise Exception( - "failed to tune dipolar P3M parameters to required accuracy") + handle_errors("dipolar P3M tuning failed") print(to_str(log)) self._params.update(self._get_params_from_es_core()) @@ -208,8 +206,7 @@ IF DP3M == 1: cdef char * log = NULL cdef int response response = dp3m_adaptive_tune(& log) - handle_errors( - "dipolar P3M_init: k-space cutoff is larger than half of box dimension") + handle_errors("dipolar P3M tuning failed") return response, log def python_dp3m_set_params(self, p_r_cut, p_mesh, p_cao, p_alpha, diff --git a/testsuite/python/electrostaticInteractions.py b/testsuite/python/electrostaticInteractions.py index 2c74c721bec..fbdf5f1fd77 100644 --- a/testsuite/python/electrostaticInteractions.py +++ b/testsuite/python/electrostaticInteractions.py @@ -21,7 +21,7 @@ import numpy as np import espressomd -from espressomd import electrostatics +import espressomd.electrostatics @utx.skipIfMissingFeatures(["ELECTROSTATICS"]) @@ -39,6 +39,9 @@ def setUp(self): self.system.part.add( id=1, pos=(3.0, 2.0, 2.0), q=-1) + def tearDown(self): + self.system.actors.clear() + def calc_dh_potential(self, r, df_params): kT = 1.0 q1 = self.system.part[0].q @@ -100,7 +103,6 @@ def test_p3m(self): [p3m_force, 0, 0], atol=1E-5) np.testing.assert_allclose(np.copy(self.system.part[1].f), [-p3m_force, 0, 0], atol=1E-10) - self.system.actors.remove(p3m) def test_dh(self): dh_params = dict(prefactor=1.2, kappa=0.8, r_cut=2.0) @@ -131,7 +133,6 @@ def test_dh(self): np.testing.assert_allclose(u_dh_core, u_dh, atol=1e-7) np.testing.assert_allclose(f_dh_core, -f_dh, atol=1e-2) - self.system.actors.remove(dh) def test_rf(self): """Tests the ReactionField coulomb interaction by comparing the @@ -173,7 +174,6 @@ def test_rf(self): np.testing.assert_allclose(u_rf_core, u_rf, atol=1e-7) np.testing.assert_allclose(f_rf_core, -f_rf, atol=1e-2) - self.system.actors.remove(rf) if __name__ == "__main__":