Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix unstable electrode tutorial #4824

Merged
merged 3 commits into from
Nov 24, 2023
Merged

Fix unstable electrode tutorial #4824

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
7202e96
Improve steepest descent and enable CI
schlaicha Nov 24, 2023
2068be5
Limit Grahame test to 1 point
schlaicha Nov 24, 2023
11a92de
revert RNG seed
schlaicha Nov 24, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
56 changes: 49 additions & 7 deletions doc/tutorials/electrodes/electrodes_part2.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)"
]
},
{
Expand Down Expand Up @@ -1015,7 +1055,9 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize=(10, 6))\n",
Expand Down
2 changes: 1 addition & 1 deletion testsuite/scripts/tutorials/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
9 changes: 5 additions & 4 deletions testsuite/scripts/tutorials/test_electrodes_2.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down Expand Up @@ -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__":
Expand Down