Refactor integrator framework to reduce coupling

doc/sphinx/electrostatics.rst

@@ -433,8 +433,8 @@ does not need to be specified as it is automatically determined from the
 particle distances and maximal pairwise error. The second tuning form
 just takes the maximal pairwise error and tries out a lot of switching
 radii to find out the fastest one. If this takes too long, you can
-change the value of the setmd variable ``timings``, which controls the number of
-test force calculations.
+change the value of the property :attr:`espressomd.system.System.timings`,
+which controls the number of test force calculations.
 
 ::
 

doc/sphinx/running.rst

@@ -9,43 +9,43 @@ To run the integrator call the method
     system.integrator.run(number_of_steps, recalc_forces=False, reuse_forces=False)
 
 where ``number_of_steps`` is the number of time steps the integrator
-should perform. The two main integration schemes of |es| are the Velocity Verlet algorithm 
-and an adaption of the Velocity Verlet algorithm to simulate an NPT ensemble.
+should perform. The two main integration schemes of |es| are the Velocity Verlet algorithm
+and an adaption of the Velocity Verlet algorithm to simulate an NpT ensemble.
 A steepest descent implementation is also available.
 
 .. _Velocity Verlet Algorithm:
 
 Velocity Verlet algorithm
 -------------------------
 
-:func:`espressomd.integrate.Integrator.set_vv`
+:meth:`espressomd.integrate.IntegratorHandle.set_vv`
 
 The equations of motion for the trajectory of point-like particles read
 
 .. math:: \dot v_i(t) = F_i(\{x_j\},v_i,t)/m_i \\ \dot x_i(t) = v_i(t),
 
 where :math:`x_i`, :math:`v_i`, :math:`m_i` are position, velocity and mass of
-particle :math:`i` and :math:`F_i(\{x_j\},v_i,t)` the forces acting on it. 
+particle :math:`i` and :math:`F_i(\{x_j\},v_i,t)` the forces acting on it.
 These forces comprise all interactions with other particles and external fields
-as well as non-deterministic contributions described in :ref:`Thermostats`. 
+as well as non-deterministic contributions described in :ref:`Thermostats`.
 
 For numerical integration, this equation is discretized to the following steps (:cite:`rapaport04` eqs. 3.5.8 - 3.5.10):
 
 1. Calculate the velocity at the half step
 
-.. math:: v(t+dt/2) = v(t) + \frac{F(x(t),v(t-dt/2),t)}{m} dt/2
+   .. math:: v(t+dt/2) = v(t) + \frac{F(x(t),v(t-dt/2),t)}{m} dt/2
 
 2. Calculate the new position
 
-.. math:: x(t+dt) = x(t) + v(t+dt/2) dt
+   .. math:: x(t+dt) = x(t) + v(t+dt/2) dt
 
 3. Calculate the force based on the new position
 
-.. math:: F = F(x(t+dt), v(t+dt/2), t+dt)
+   .. math:: F = F(x(t+dt), v(t+dt/2), t+dt)
 
 4. Calculate the new velocity
 
-.. math:: v(t+dt) = v(t+dt/2) + \frac{F(x(t+dt),t+dt)}{m} dt/2
+   .. math:: v(t+dt) = v(t+dt/2) + \frac{F(x(t+dt),t+dt)}{m} dt/2
 
 Note that this implementation of the Velocity Verlet algorithm reuses
 forces in step 1. That is, they are computed once in step 3,
@@ -82,24 +82,24 @@ enforce force recalculation.
 Isotropic NPT integrator
 ------------------------
 
-:py:func:`~espressomd.integrate.Integrator.set_isotropic_npt`
+:meth:`espressomd.integrate.IntegratorHandle.set_isotropic_npt`
 
-As the NPT thermostat alters the way the equations of motion are integrated, it is
+As the NpT thermostat alters the way the equations of motion are integrated, it is
 discussed here and only a brief summary is given in :ref:`Thermostats`.
 
-To activate the NPT integrator, use :py:func:`~espressomd.integrate.Integrator.set_isotropic_npt`
+To activate the NpT integrator, use :meth:`~espressomd.integrate.IntegratorHandle.set_isotropic_npt`
 with parameters:
 
-    * ``ext_pressure``: The external pressure
-    * ``piston``: The mass of the applied piston
-    * ``direction``: Flags to enable/disable box dimensions to be subject to fluctuations. By default, all directions are enabled.
+* ``ext_pressure``: The external pressure
+* ``piston``: The mass of the applied piston
+* ``direction``: Flags to enable/disable box dimensions to be subject to fluctuations. By default, all directions are enabled.
 
-Additionally, a NPT thermostat has to be set by :py:func:`~espressomd.thermostat.Thermostat.set_npt()`
+Additionally, a NpT thermostat has to be set by :meth:`~espressomd.thermostat.Thermostat.set_npt()`
 with parameters:
 
-    * ``kT``: Thermal energy of the heat bath
-    * ``gamma0``: Friction coefficient of the bath
-    * ``gammav``: Artificial friction coefficient for the volume fluctuations.
+* ``kT``: Thermal energy of the heat bath
+* ``gamma0``: Friction coefficient of the bath
+* ``gammav``: Artificial friction coefficient for the volume fluctuations.
 
 A code snippet would look like::
 
@@ -111,66 +111,73 @@ A code snippet would look like::
 
 The physical meaning of these parameters is described below:
 
-The relaxation towards a desired pressure :math:`P` (parameter ``ext_pressure``) 
+The relaxation towards a desired pressure :math:`P` (parameter ``ext_pressure``)
 is enabled by treating the box
-volume :math:`V` as a degree of freedom with corresponding momentum :math:`\Pi = Q\dot{V}`, 
+volume :math:`V` as a degree of freedom with corresponding momentum :math:`\Pi = Q\dot{V}`,
 where :math:`Q` (parameter ``piston``) is an artificial piston mass.
 Which box dimensions are affected to change the volume can be controlled by a list of
 boolean flags for parameter ``direction``.
-An additional energy :math:`H_V = 1/(2Q)\Pi + PV` 
+An additional energy :math:`H_V = 1/(2Q)\Pi + PV`
 associated with the volume is postulated. This results in a "force" on the box such that
 
 .. math:: \dot{\Pi} = \mathcal{P} - P
 
-where 
+where
 
 .. math:: \mathcal{P} = \frac{1}{Vd} \sum_{i,j} f_{ij}x_{ij} + \frac{1}{Vd} \sum_i m_i v_i^2
 
-Here :math:`\mathcal{P}` is the instantaneous pressure, :math:`d` the dimension of the system (number of flags set by ``direction``), :math:`f_{ij}` the short range interaction force between particles :math:`i` and :math:`j` and :math:`x_{ij}= x_j - x_i`.
+Here :math:`\mathcal{P}` is the instantaneous pressure, :math:`d` the dimension
+of the system (number of flags set by ``direction``), :math:`f_{ij}` the
+short range interaction force between particles :math:`i` and :math:`j` and
+:math:`x_{ij}= x_j - x_i`.
 
-In addition to this deterministic force, a friction :math:`-\frac{\gamma^V}{Q}\Pi(t)` and noise :math:`\sqrt{k_B T \gamma^V} \eta(t)` are added for the box
+In addition to this deterministic force, a friction :math:`-\frac{\gamma^V}{Q}\Pi(t)`
+and noise :math:`\sqrt{k_B T \gamma^V} \eta(t)` are added for the box
 volume dynamics and the particle dynamics. This introduces three new parameters:
-The friction coefficient for the box :math:`\gamma^V` (parameter ``gammav``), the friction coefficient of the particles :math:`\gamma^0` (parameter ``gamma0``) and the thermal energy :math:`k_BT` (parameter ``kT``).
-For a discussion of these terms and their discretisation, see :ref:`Langevin thermostat`, which uses the same approach, but only for particles.
-As a result of box geometry changes, the particle positions and velocities have to be rescaled 
+The friction coefficient for the box :math:`\gamma^V` (parameter ``gammav``),
+the friction coefficient of the particles :math:`\gamma^0` (parameter ``gamma0``)
+and the thermal energy :math:`k_BT` (parameter ``kT``).
+For a discussion of these terms and their discretisation, see :ref:`Langevin thermostat`,
+which uses the same approach, but only for particles.
+As a result of box geometry changes, the particle positions and velocities have to be rescaled
 during integration.
 
 The discretisation consists of the following steps (see :cite:`kolb99a` for a full derivation of the algorithm):
 
 1. Calculate the particle velocities at the half step
 
-.. math:: v'(t+dt/2) = v(t) + \frac{F(x(t),v(t-dt/2),t)}{m} dt/2 
+   .. math:: v'(t+dt/2) = v(t) + \frac{F(x(t),v(t-dt/2),t)}{m} dt/2
 
 2. Calculate the instantaneous pressure and "volume momentum"
 
-.. math:: \mathcal{P} = \mathcal{P}(x(t),V(t),f(x(t)), v'(t+dt/2))
-.. math:: \Pi(t+dt/2) = \Pi(t) + (\mathcal{P}-P) dt/2 -\frac{\gamma^V}{Q}\Pi(t) dt/2  +  \sqrt{k_B T \gamma^V dt} \overline{\eta} 
+   .. math:: \mathcal{P} = \mathcal{P}(x(t),V(t),f(x(t)), v'(t+dt/2))
+   .. math:: \Pi(t+dt/2) = \Pi(t) + (\mathcal{P}-P) dt/2 -\frac{\gamma^V}{Q}\Pi(t) dt/2  +  \sqrt{k_B T \gamma^V dt} \overline{\eta}
 
 3. Calculate box volume and scaling parameter :math:`L` at half step and full step, scale the simulation box accordingly
 
-.. math:: V(t+dt/2) = V(t) + \frac{\Pi(t+dt/2)}{Q} dt/2
-.. math:: L(t+dt/2) = V(t+dt/2)^{1/d}
-.. math:: V(t+dt) = V(t+dt/2) + \frac{\Pi(t+dt/2)}{Q} dt/2
-.. math:: L(t+dt) = V(t+dt)^{1/d} 
+   .. math:: V(t+dt/2) = V(t) + \frac{\Pi(t+dt/2)}{Q} dt/2
+   .. math:: L(t+dt/2) = V(t+dt/2)^{1/d}
+   .. math:: V(t+dt) = V(t+dt/2) + \frac{\Pi(t+dt/2)}{Q} dt/2
+   .. math:: L(t+dt) = V(t+dt)^{1/d}
 
 4. Update particle positions and scale velocities
-    
-.. math:: x(t+dt) = \frac{L(t+dt)}{L(t)} \left[ x(t) + \frac{L^2(t)}{L^2(t+dt/2)} v(t+dt/2) dt \right]
-.. math:: v(t+dt/2) = \frac{L(t)}{L(t+dt)} v'(t+dt/2)
+
+   .. math:: x(t+dt) = \frac{L(t+dt)}{L(t)} \left[ x(t) + \frac{L^2(t)}{L^2(t+dt/2)} v(t+dt/2) dt \right]
+   .. math:: v(t+dt/2) = \frac{L(t)}{L(t+dt)} v'(t+dt/2)
 
 5. Calculate forces, instantaneous pressure and "volume momentum"
 
-.. math:: F = F(x(t+dt),v(t+dt/2),t)
-.. math:: \mathcal{P} = \mathcal{P}(x(t+dt),V(t+dt),f(x(t+dt)), v(t+dt/2))
-.. math:: \Pi(t+dt) = \Pi(t+dt/2) + (\mathcal{P}-P) dt/2 -\frac{\gamma^V}{Q}\Pi(t+dt/2) dt/2  +  \sqrt{k_B T \gamma^V dt} \overline{\eta}
+   .. math:: F = F(x(t+dt),v(t+dt/2),t)
+   .. math:: \mathcal{P} = \mathcal{P}(x(t+dt),V(t+dt),f(x(t+dt)), v(t+dt/2))
+   .. math:: \Pi(t+dt) = \Pi(t+dt/2) + (\mathcal{P}-P) dt/2 -\frac{\gamma^V}{Q}\Pi(t+dt/2) dt/2  +  \sqrt{k_B T \gamma^V dt} \overline{\eta}
 
 6. Update the velocities
 
-.. math:: v(t+dt) = v(t+dt/2) + \frac{F(t+dt)}{m} dt/2
+   .. math:: v(t+dt) = v(t+dt/2) + \frac{F(t+dt)}{m} dt/2
 
 Notes:
 
-* The NPT algorithm is only tested for all 3 directions enabled for scaling. Usage of ``direction`` is considered an experimental feature.
+* The NpT algorithm is only tested for all 3 directions enabled for scaling. Usage of ``direction`` is considered an experimental feature.
 * In step 4, only those coordinates are scaled for which ``direction`` is set.
 * For the instantaneous pressure, the same limitations of applicability hold as described in :ref:`Pressure`.
 * The particle forces :math:`F` include interactions as well as a friction and noise term analogous to the terms in the :ref:`Langevin thermostat`.
@@ -185,7 +192,7 @@ When the feature ``ROTATION`` is compiled in, particles not only have a position
 
 The rotational degrees of freedom are also integrated using a velocity Verlet scheme.
 The implementation is based on a quaternion representation of the particle orientation and described in :cite:`omelyan98` with quaternion components indexing made according to the formalism :math:`q = a + b\mathbf{i} + c\mathbf{j} + d\mathbf{k}` :cite:`allen2017`.
-   
+
 When the Langevin thermostat is enabled, the rotational degrees of freedom are also thermalized.
 
 Whether or not rotational degrees of freedom are propagated, is controlled on a per-particle and per-axis level, where the axes are the Cartesian axes of the particle in its body-fixed frame.
@@ -201,8 +208,8 @@ The rotation of a particle is controlled via the :attr:`espressomd.particle_data
 Notes:
 
 * The orientation of a particle is stored as a quaternion in the :attr:`espressomd.particle_data.ParticleHandle.quat` property. For a value of (1,0,0,0), the body and space frames coincide.
-* The space-frame direction of the particle's z-axis in its body frame is accessible through the ``espressomd.particle_data.ParticleHandle.director`` property.
-* Any other vector can be converted from body to space fixed frame using the ``espressomd.particle_data.ParticleHandle.convert_vector_body_to_space`` method.
+* The space-frame direction of the particle's z-axis in its body frame is accessible through the :attr:`espressomd.particle_data.ParticleHandle.director` property.
+* Any other vector can be converted from body to space fixed frame using the :meth:`espressomd.particle_data.ParticleHandle.convert_vector_body_to_space` method.
 * When ``DIPOLES`` are compiled in, the particles dipole moment is always co-aligned with the z-axis in the body-fixed frame.
 * Changing the particles dipole moment and director will re-orient the particle such that its z-axis in space frame is aligned parallel to the given vector. No guarantees are made for the other two axes after setting the director or the dipole moment.
 
@@ -226,7 +233,7 @@ The following particle properties are related to rotation:
 Steepest descent
 ----------------
 
-:func:`espressomd.integrate.Integrator.set_steepest_descent`
+:meth:`espressomd.integrate.IntegratorHandle.set_steepest_descent`
 
 This feature is used to propagate each particle by a small distance parallel to the force acting on it.
 When only conservative forces for which a potential exists are in use, this is equivalent to a steepest descent energy minimization.
@@ -248,7 +255,19 @@ coordinates that are fixed using the ``fix`` and ``rotation`` attribute of parti
 
 Usage example::
 
-        system.integrator.set_steepest_descent(
-            f_max=0, gamma=0.1, max_displacement=0.1)
-        system.integrator.run(20)
-        system.integrator.set_vv()  # to switch back to velocity verlet
+    system.integrator.set_steepest_descent(
+        f_max=0, gamma=0.1, max_displacement=0.1)
+    system.integrator.run(20)
+    system.integrator.set_vv()  # to switch back to velocity Verlet
+
+A wrapper function :func:`~espressomd.minimize_energy.steepest_descent` is
+available to set up the steepest descent integrator while preserving the
+original integrator. The snippet above can be rewritten to::
+
+    from espressomd.minimize_energy import steepest_descent
+    steepest_descent(system, f_max=0, gamma=0.1, max_displacement=0.1, max_steps=20)
+
+This convenience function only exists for historical reasons and remains available
+for backward compatibility. New scripts should setup the steepest descent
+integrator with the :meth:`~espressomd.integrate.IntegratorHandle.set_steepest_descent`
+handle directly.

doc/sphinx/system_setup.rst

@@ -27,7 +27,7 @@ for further calculations, you should explicitly make a copy e.g. via
     particle coordinates will remain the same, i.e., the particle stay in
     the same image box, but at the same relative position in their image
     box. If you want to scale the positions, use the command
-    :py:func:`~espressomd.system.System.change_volume_and_rescale_particles`.
+    :py:meth:`~espressomd.system.System.change_volume_and_rescale_particles`.
 
 * :py:attr:`~espressomd.system.System.periodicity`
 
@@ -411,8 +411,8 @@ Isotropic NPT thermostat
 
 This feature allows to simulate an (on average) homogeneous and isotropic system in the NPT ensemble.
 In order to use this feature, ``NPT`` has to be defined in the :file:`myconfig.hpp`.
-Activate the NPT thermostat with the command :py:func:`~espressomd.thermostat.Thermostat.set_npt`
-and setup the integrator for the NPT ensemble with :py:func:`~espressomd.integrate.Integrator.set_isotropic_npt`.
+Activate the NPT thermostat with the command :py:meth:`~espressomd.thermostat.Thermostat.set_npt`
+and setup the integrator for the NPT ensemble with :py:meth:`~espressomd.integrate.IntegratorHandle.set_isotropic_npt`.
 
 For example::
 

doc/tutorials/02-charged_system/02-charged_system-1.ipynb

@@ -40,6 +40,7 @@
    "outputs": [],
    "source": [
     "from espressomd import System, electrostatics\n",
+    "from espressomd.minimize_energy import steepest_descent\n",
     "import espressomd\n",
     "import numpy\n",
     "import matplotlib.pyplot as plt\n",
@@ -194,13 +195,12 @@
     "# WCA Equilibration\n",
     "max_sigma = max(wca_sigmas.values())\n",
     "min_dist = 0.0\n",
-    "system.minimize_energy.init(f_max=0, gamma=10.0, max_steps=10,\n",
-    "                            max_displacement=max_sigma * 0.01)\n",
     "\n",
     "while min_dist < max_sigma:\n",
     "    min_dist = system.analysis.min_dist([types[\"Anion\"], types[\"Cation\"]],\n",
     "                                        [types[\"Anion\"], types[\"Cation\"]])\n",
-    "    system.minimize_energy.minimize()\n",
+    "    steepest_descent(system, f_max=0, gamma=10.0, max_steps=10,\n",
+    "                     max_displacement=max_sigma * 0.01)\n",
     "\n",
     "# Set thermostat\n",
     "system.thermostat.set_langevin(kT=temp, gamma=gamma, seed=42)"

doc/tutorials/02-charged_system/02-charged_system-2.ipynb

@@ -69,6 +69,7 @@
    "source": [
     "from espressomd import System, electrostatics, electrostatic_extensions\n",
     "from espressomd.shapes import Wall\n",
+    "from espressomd.minimize_energy import steepest_descent\n",
     "import espressomd\n",
     "import numpy"
    ]
@@ -268,9 +269,8 @@
    "source": [
     "energy = system.analysis.energy()\n",
     "print(\"Before Minimization: E_total = {:.3e}\".format(energy['total']))\n",
-    "system.minimize_energy.init(f_max=10, gamma=10, max_steps=2000,\n",
-    "                            max_displacement=0.01)\n",
-    "system.minimize_energy.minimize()\n",
+    "steepest_descent(system, f_max=10, gamma=10, max_steps=2000,\n",
+    "                 max_displacement=0.01)\n",
     "energy = system.analysis.energy()\n",
     "print(\"After Minimization: E_total = {:.3e}\".format(energy['total']))\n",
     "\n",

doc/tutorials/02-charged_system/scripts/nacl.py

@@ -20,6 +20,7 @@
 #
 import espressomd
 from espressomd import assert_features, electrostatics
+from espressomd.minimize_energy import steepest_descent
 import numpy
 
 assert_features(["ELECTROSTATICS", "WCA"])
@@ -91,11 +92,10 @@ def combination_rule_sigma(rule, sig1, sig2):
 print("\n--->WCA Equilibration")
 max_sigma = max(wca_sigmas.values())
 min_dist = 0.0
-system.minimize_energy.init(f_max=0, gamma=10, max_steps=10,
-                            max_displacement=max_sigma * 0.01)
 
 while min_dist < max_sigma:
-    system.minimize_energy.minimize()
+    steepest_descent(system, f_max=0, gamma=10, max_steps=10,
+                     max_displacement=max_sigma * 0.01)
     min_dist = system.analysis.min_dist()
 
 # Set thermostat

doc/tutorials/02-charged_system/scripts/nacl_units_confined.py

@@ -21,6 +21,7 @@
 import espressomd
 from espressomd import electrostatics, electrostatic_extensions, assert_features
 from espressomd.shapes import Wall
+from espressomd.minimize_energy import steepest_descent
 import numpy
 
 assert_features(["ELECTROSTATICS", "MASS", "LENNARD_JONES"])
@@ -123,9 +124,8 @@ def combination_rule_sigma(rule, sig1, sig2):
 
 energy = system.analysis.energy()
 print("Before Minimization: E_total=", energy['total'])
-system.minimize_energy.init(
-    f_max=10, gamma=50.0, max_steps=1000, max_displacement=0.2)
-system.minimize_energy.minimize()
+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'])
 

doc/tutorials/02-charged_system/scripts/nacl_units_confined_vis.py

@@ -21,6 +21,7 @@
 import espressomd
 from espressomd import assert_features, electrostatics, electrostatic_extensions
 from espressomd.shapes import Wall
+from espressomd.minimize_energy import steepest_descent
 from espressomd import visualization_opengl
 import numpy
 from threading import Thread
@@ -126,9 +127,8 @@ def combination_rule_sigma(rule, sig1, sig2):
     system.non_bonded_inter[types[s[0]], types[s[1]]].lennard_jones.set_params(
         epsilon=lj_eps, sigma=lj_sig, cutoff=lj_cut, shift="auto")
 
-system.minimize_energy.init(
-    f_max=10, gamma=10, max_steps=2000, max_displacement=0.1)
-system.minimize_energy.minimize()
+steepest_descent(system, f_max=10, gamma=10, max_steps=2000,
+                 max_displacement=0.1)
 
 print("\n--->Tuning Electrostatics")
 p3m = electrostatics.P3M(prefactor=l_bjerrum, accuracy=1e-2)