Skip to content

Commit

Permalink
tutorials: Prevent particles from crossing ELC gap
Browse files Browse the repository at this point in the history 
The chosen time step allowed particles to move a distance of 15%
sigma per time step, which made the simulation unstable. Particles
would sometimes get too close to a wall and acquire more kinetic
energy than the thermostat could dissipate, and eventually move
too far per time step and croos the walls. The new time step
prevents particles from moving more than 5% sigma per time step.

Co-authored-by: Rudolf Weeber <weeber@icp.uni-stuttgart.de>
Co-authored-by: keerthirk <keerthiradhakrishnan1995@gmail.com>
  • Loading branch information
@jngrad @RudolfWeeber @keerthirk1995
3 people committed Mar 28, 2024
1 parent b51fc95 commit 0a5af4c
Show file tree
Hide file tree
Showing 3 changed files with 65 additions and 111 deletions.
4 changes: 2 additions & 2 deletions doc/tutorials/electrodes/CMakeLists.txt
View file
Expand Up @@ -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)
167 changes: 61 additions & 106 deletions doc/tutorials/electrodes/electrodes_part2.ipynb
View file
Expand Up @@ -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",
Expand All @@ -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",
Expand All @@ -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",
Expand All @@ -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)"
]
},
{
Expand Down Expand Up @@ -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, :]"
Expand Down Expand Up @@ -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",
Expand Down
5 changes: 2 additions & 3 deletions testsuite/scripts/tutorials/test_electrodes_2.py
View file
Expand Up @@ -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(
Expand Down Expand Up @@ -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__":
Expand Down

0 comments on commit 0a5af4c

Please sign in to comment.