Fix nightly builds (#4894)

Fixes #4850 Closes #4849 Closes #4851 Closes #4857 Closes #4860 Closes #4862 Closes #4861 Closes #4863 Closes #4865 Closes #4868 Closes #4869 Closes #4871 Closes #4873 Closes #4875 Closes #4876 Closes #4877 Closes #4878 Closes #4879 Closes #4880 Closes #4881 Closes #4882 Closes #4883 Closes #4885 Closes #4886 Closes #4887 Closes #4888 Closes #4889 Closes #4890 Closes #4891 Closes #4892 Closes #4893 Closes #4895 Closes #4896 Description of changes: - fix particle crossing the ELC gap in the electrodes tutorial - disable ASAN checks due to broken dependencies on the host machines
espressomd · Mar 28, 2024 · f642a45 · f642a45
2 parents 7ad0534 + 0a5af4c
commit f642a45
Show file tree

Hide file tree

Showing 4 changed files with 66 additions and 112 deletions.
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
@@ -171,7 +171,7 @@ clang-sanitizer:
      with_coverage: 'false'
      with_static_analysis: 'true'
      check_skip_long: 'true'
-     with_asan: 'true'
+     with_asan: 'false'
      with_ubsan: 'true'
      with_scafacos: 'true'
      with_walberla: 'true'

diff --git a/doc/tutorials/electrodes/CMakeLists.txt b/doc/tutorials/electrodes/CMakeLists.txt
@@ -22,5 +22,5 @@ configure_tutorial_target(TARGET tutorial_electrodes DEPENDS
 
 nb_export(TARGET tutorial_electrodes SUFFIX "1" FILE "electrodes_part1.ipynb"
           HTML_RUN)
-# TODO: fix time step issues (#4798) before adding HTML_RUN back
-nb_export(TARGET tutorial_electrodes SUFFIX "2" FILE "electrodes_part2.ipynb")
+nb_export(TARGET tutorial_electrodes SUFFIX "2" FILE "electrodes_part2.ipynb"
+          HTML_RUN)
diff --git a/doc/tutorials/electrodes/electrodes_part2.ipynb b/doc/tutorials/electrodes/electrodes_part2.ipynb
@@ -570,13 +570,16 @@
    "outputs": [],
    "source": [
     "# SOLUTION CELL\n",
+    "\n",
     "def setup_electrostatic_solver(potential_diff):\n",
     "    delta_mid_top = -1.  # (Fully metallic case both -1)                 \n",
     "    delta_mid_bot = -1.\n",
-    "    accuracy = 1e-3\n",
+    "    p3m_accuracy = 1e-3\n",
     "    elc_accuracy = 1e-3\n",
     "    p3m = espressomd.electrostatics.P3M(prefactor=BJERRUM_LENGTH,\n",
-    "                                        accuracy=accuracy,\n",
+    "                                        accuracy=p3m_accuracy,\n",
+    "                                        mesh=[12, 12, 48], # pinned for tuning reproducibility\n",
+    "                                        cao=4, # pinned for tuning reproducibility\n",
     "                                        tune=True,\n",
     "                                        verbose=False)\n",
     "    elc = espressomd.electrostatics.ELC(actor=p3m,\n",
@@ -589,32 +592,6 @@
     "    return elc"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "efbf4cf9",
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "markdown",
-   "id": "03ab39a1",
-   "metadata": {},
-   "source": [
-    "Now add the solver to the system:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "25219528",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "system.electrostatics.solver = setup_electrostatic_solver(POTENTIAL_DIFF)"
-   ]
-  },
   {
    "cell_type": "markdown",
    "id": "5fed3232",
@@ -625,111 +602,68 @@
     "### 2.1 Steepest descent\n",
     "\n",
     "Before we can start the simulation, we need to remove the overlap between particles.\n",
-    "For this, we use the steepest descent integrator.\n",
-    "Afterwards, we switch to a Velocity Verlet integrator and set up a Langevin thermostat.\n",
-    "Note, that we only analyze static properties, thus the damping and temperature chosen\n",
-    "here only determine the simulation time towards the equilibrium distribution."
+    "For this, we use the steepest descent integrator."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "51a25228",
+   "id": "7da4b3f6",
    "metadata": {},
    "outputs": [],
-   "source": [
-    "# 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,
-   "id": "1a3cacd2",
-   "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"
+    "system.integrator.set_steepest_descent(f_max=0., gamma=30., max_displacement=0.01)\n",
+    "system.integrator.run(100)\n",
+    "particle_forces = np.linalg.norm(system.part.all().f, axis=1)\n",
+    "assert np.max(particle_forces) < 1."
    ]
   },
   {
    "cell_type": "markdown",
    "id": "abbfc272",
    "metadata": {},
    "source": [
-    "### 2.2 Equilibrate the ion distribution"
+    "### 2.2 Warmup\n",
+    "\n",
+    "We now switch to a Velocity Verlet integrator and set up a Langevin thermostat.\n",
+    "Note, that we only analyze static properties, thus the damping and temperature chosen\n",
+    "here only determine the simulation time towards the equilibrium distribution."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "45c444f5",
+   "id": "38137a83",
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Switch to velocity verlet integration using a Langevin thermostat\n",
+    "# Equilibration parameters\n",
+    "system.time_step = 0.003 # this time step limits particle movement to ~5% sigma per integration step\n",
     "system.integrator.set_vv()\n",
-    "system.thermostat.set_langevin(kT=1.0, gamma=0.1, seed=42)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "e9c7fe2f",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# Equlibration parameters\n",
-    "STEPS_PER_SAMPLE = 200\n",
-    "N_SAMPLES_EQUIL = 50\n",
-    "N_PART = 2 * N_IONPAIRS\n",
+    "system.thermostat.set_langevin(kT=1., gamma=2., seed=42)\n",
+    "system.electrostatics.solver = setup_electrostatic_solver(POTENTIAL_DIFF)\n",
     "\n",
-    "times = np.zeros(N_SAMPLES_EQUIL)\n",
-    "e_total = np.zeros_like(times)\n",
-    "e_kin = np.zeros_like(times)\n",
-    "\n",
-    "for i in tqdm.trange(N_SAMPLES_EQUIL):\n",
-    "    times[i] = system.time\n",
+    "times = []\n",
+    "e_kin = []\n",
+    "for i in tqdm.trange(100):\n",
     "    energy = system.analysis.energy()\n",
-    "    e_total[i] = energy['total']\n",
-    "    e_kin[i] = energy['kinetic']\n",
-    "    system.integrator.run(STEPS_PER_SAMPLE)"
+    "    e_kin.append(energy['kinetic'])\n",
+    "    times.append(system.time)\n",
+    "    system.integrator.run(10)"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "c45afc27",
+   "id": "970825d8-6087-4b38-a39e-f67e718cbd6e",
    "metadata": {},
    "outputs": [],
    "source": [
     "# Plot the convergence of the total energy\n",
     "plt.figure(figsize=(10, 6))\n",
-    "plt.plot(times, e_total, label='total')\n",
-    "plt.plot(times, e_kin, label='kinetic')\n",
+    "plt.plot(times, len(e_kin) * [N_IONPAIRS * 2 * 3. / 2.], label='heat bath')\n",
+    "plt.plot(times, e_kin, label='kinetic energy')\n",
     "plt.xlabel('Simulation time')\n",
     "plt.ylabel('Energy')\n",
     "plt.legend()\n",
@@ -741,7 +675,32 @@
    "id": "1f1f7892",
    "metadata": {},
    "source": [
-    "Convergence after $t\\sim50$ time units."
+    "Convergence after $t \\sim 5$ time units."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9a22f6b5-33f9-4723-8e1c-41717cf6bb03",
+   "metadata": {},
+   "source": [
+    "### 2.3 Equilibrate the ion distribution\n",
+    "\n",
+    "Now we let ions diffuse to the electrodes."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e35f031d-4188-4d70-8852-7ff3ce571904",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "N_SAMPLES_EQUIL = 20\n",
+    "STEPS_PER_EQUIL = 200\n",
+    "\n",
+    "system.thermostat.set_langevin(kT=1.0, gamma=0.1, seed=42)\n",
+    "for tm in tqdm.trange(N_SAMPLES_EQUIL):\n",
+    "    system.integrator.run(STEPS_PER_EQUIL)"
    ]
   },
   {
@@ -831,20 +790,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "N_SAMPLES_PROD = 10\n",
+    "N_SAMPLES_PROD = 20\n",
+    "STEPS_PER_SAMPLE = 200\n",
     "\n",
     "# Add the accumulators\n",
-    "system.auto_update_accumulators.clear()\n",
     "system.auto_update_accumulators.add(density_accumulator_cation)\n",
     "system.auto_update_accumulators.add(density_accumulator_anion)\n",
-    "    \n",
-    "times = []\n",
-    "e_total = []\n",
+    "\n",
     "for tm in tqdm.trange(N_SAMPLES_PROD):\n",
     "    system.integrator.run(STEPS_PER_SAMPLE)\n",
-    "    times.append(system.time)\n",
-    "    energy = system.analysis.energy()\n",
-    "    e_total.append(energy['total'])\n",
     "\n",
     "cation_profile_mean = density_accumulator_cation.mean()[0, 0, :]\n",
     "anion_profile_mean = density_accumulator_anion.mean()[0, 0, :]"
@@ -1061,11 +1015,12 @@
     "MIN_PHI = 0.5\n",
     "MAX_PHI = 10\n",
     "N_PHI = 7\n",
-    "N_SAMPLES_EQUIL_CAP = 5\n",
+    "N_SAMPLES_EQUIL_CAP = 15\n",
     "N_SAMPLES_CAP = 5\n",
     "\n",
     "# sample from high to low potential to improve sampling\n",
     "for potential_diff in tqdm.tqdm(np.linspace(MIN_PHI, MAX_PHI, N_PHI)[::-1]):\n",
+    "    system.auto_update_accumulators.clear()\n",
     "    system.electrostatics.solver = setup_electrostatic_solver(potential_diff)\n",
     "    system.integrator.run(N_SAMPLES_EQUIL_CAP * STEPS_PER_SAMPLE)\n",
     "    sigmas = []\n",

diff --git a/testsuite/scripts/tutorials/test_electrodes_2.py b/testsuite/scripts/tutorials/test_electrodes_2.py
@@ -22,8 +22,7 @@
 import numpy as np
 from scipy import constants
 
-params = {'N_SAMPLES_EQUIL': 25, 'N_SAMPLES_PROD': 5,
-          'N_SAMPLES_EQUIL_CAP': 0, 'N_SAMPLES_CAP': 5,
+params = {'N_SAMPLES_PROD': 15, 'N_SAMPLES_EQUIL_CAP': 5, 'N_SAMPLES_CAP': 5,
           'MIN_PHI': 5, 'MAX_PHI': 5, 'N_PHI': 1}
 
 tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import(
@@ -66,7 +65,7 @@ def test_capacitance(self):
             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=.05, err_msg=msg)
+            grahame, tutorial.sigma_vs_phi[:, 1], atol=.08, err_msg=msg)
 
 
 if __name__ == "__main__":