Fix unstable electrode tutorial (#4824)

Fixes #4798 Description of changes: - change steepest descent parameters to more reasonable values in electrode_part2 - re-enable test of electrode_part2 - Limit Grahame test to 1 point / shorten CI test duration
espressomd · Nov 24, 2023 · 2c14952 · 2c14952
2 parents 38a5338 + 11a92de
commit 2c14952
Show file tree

Hide file tree

Showing 3 changed files with 55 additions and 12 deletions.
diff --git a/doc/tutorials/electrodes/electrodes_part2.ipynb b/doc/tutorials/electrodes/electrodes_part2.ipynb
@@ -625,18 +625,58 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "system.integrator.set_steepest_descent(f_max=10, gamma=50.0,\n",
-    "                                       max_displacement=0.02)\n",
-    "system.integrator.run(250)\n",
-    "system.integrator.set_vv()\n",
-    "system.thermostat.set_langevin(kT=1.0, gamma=0.1, seed=42)"
+    "# suitable minimization parameters for this system\n",
+    "F_TOL = 1e-2\n",
+    "DAMPING = 30\n",
+    "MAX_STEPS = 10000\n",
+    "MAX_DISPLACEMENT = 0.01 * LJ_SIGMA\n",
+    "EM_STEP = 10"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "# Set up steepest descent integration\n",
+    "system.integrator.set_steepest_descent(f_max=0,  # use a relative convergence criterion only\n",
+    "                                       gamma=DAMPING,\n",
+    "                                       max_displacement=MAX_DISPLACEMENT)\n",
+    "\n",
+    "# Initialize integrator to obtain initial forces\n",
+    "system.integrator.run(0)\n",
+    "old_force = np.max(np.linalg.norm(system.part.all().f, axis=1))\n",
+    "\n",
+    "\n",
+    "while system.time / system.time_step < MAX_STEPS:\n",
+    "    system.integrator.run(EM_STEP)\n",
+    "    force = np.max(np.linalg.norm(system.part.all().f, axis=1))\n",
+    "    rel_force = np.abs((force - old_force) / old_force)\n",
+    "    print(f'rel. force change: {rel_force:.2e}')\n",
+    "    if rel_force < F_TOL:\n",
+    "        break\n",
+    "    old_force = force"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Equilibrate the ion distribution"
+    "### 2.2 Equilibrate the ion distribution"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Switch to velocity verlet integration using a Langevin thermostat\n",
+    "system.integrator.set_vv()\n",
+    "system.thermostat.set_langevin(kT=1.0, gamma=0.1, seed=42)"
    ]
   },
   {
@@ -1015,7 +1055,9 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "scrolled": true
+   },
    "outputs": [],
    "source": [
     "fig, ax = plt.subplots(figsize=(10, 6))\n",

diff --git a/testsuite/scripts/tutorials/CMakeLists.txt b/testsuite/scripts/tutorials/CMakeLists.txt
@@ -62,7 +62,7 @@ tutorial_test(FILE test_ferrofluid_2.py)
 tutorial_test(FILE test_ferrofluid_3.py)
 tutorial_test(FILE test_constant_pH__ideal.py)
 tutorial_test(FILE test_electrodes_1.py)
-# tutorial_test(FILE test_electrodes_2.py) # TODO: unstable, see issue #4798
+tutorial_test(FILE test_electrodes_2.py)
 tutorial_test(FILE test_constant_pH__interactions.py)
 tutorial_test(FILE test_widom_insertion.py)
 tutorial_test(FILE test_grand_canonical_monte_carlo.py)

diff --git a/testsuite/scripts/tutorials/test_electrodes_2.py b/testsuite/scripts/tutorials/test_electrodes_2.py
@@ -23,8 +23,8 @@
 from scipy import constants
 
 params = {'N_SAMPLES_EQUIL': 25, 'N_SAMPLES_PROD': 5,
-          'N_SAMPLES_EQUIL_CAP': 5, 'N_SAMPLES_CAP': 1,
-          'MIN_PHI': 1, 'MAX_PHI': 2.5, 'N_PHI': 4}
+          'N_SAMPLES_EQUIL_CAP': 0, 'N_SAMPLES_CAP': 5,
+          'MIN_PHI': 5, 'MAX_PHI': 5, 'N_PHI': 1}
 
 tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import(
     "@TUTORIALS_DIR@/electrodes/electrodes_part2.py",
@@ -60,12 +60,13 @@ def test_charge_profile(self):
 
     def test_capacitance(self):
         # For low potentials the capacitance should be in line with Grahame/DH
-        # equilibration performance limiting
+        # equilibration performance limiting, just use the system equilibrated
+        # in the first part
         grahame = -tutorial.sigma_vs_phi[:, 0] / (
             constants.elementary_charge / (constants.Boltzmann * tutorial.TEMPERATURE))
         msg = 'The capacitance at low potentials should be in line with Grahame/DH.'
         np.testing.assert_allclose(
-            grahame, tutorial.sigma_vs_phi[:, 1], atol=.015, err_msg=msg)
+            grahame, tutorial.sigma_vs_phi[:, 1], atol=.05, err_msg=msg)
 
 
 if __name__ == "__main__":