From 712801035a5aabaf52074062acd864fcb80bce2e Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Mon, 31 Aug 2020 17:03:57 +0200 Subject: [PATCH 01/46] change indexing in the sampling loop now use to zero-based indexing in the loop, see #3818 --- doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index e2c79277544..09277ba0069 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -424,14 +424,14 @@ "instantaneous_temperature = np.zeros(sampling_iterations)\n", "etotal = np.zeros(sampling_iterations)\n", "\n", - "for i in range(1, sampling_iterations + 1):\n", + "for i in range(0, sampling_iterations):\n", " system.integrator.run(sampling_interval)\n", " # Measure energies\n", " energies = system.analysis.energy()\n", " kinetic_temperature = energies['kinetic'] / (1.5 * N_PART)\n", - " etotal[i - 1] = energies['total']\n", - " time[i - 1] = system.time\n", - " instantaneous_temperature[i - 1] = kinetic_temperature\n", + " etotal[i] = energies['total']\n", + " time[i] = system.time\n", + " instantaneous_temperature[i] = kinetic_temperature\n", "\n", "# Finalize the correlator and obtain the results\n", "msd_corr.finalize()\n", From a57181875df79de82b925dd4ad4e0b9e550aa70d Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Thu, 3 Sep 2020 21:24:56 +0200 Subject: [PATCH 02/46] omit default value in for loop --- doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index 09277ba0069..1b3f102795c 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -424,7 +424,7 @@ "instantaneous_temperature = np.zeros(sampling_iterations)\n", "etotal = np.zeros(sampling_iterations)\n", "\n", - "for i in range(0, sampling_iterations):\n", + "for i in range(sampling_iterations):\n", " system.integrator.run(sampling_interval)\n", " # Measure energies\n", " energies = system.analysis.energy()\n", From 0999091e422887009c3e5f814d80cbe578c08377 Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Thu, 3 Sep 2020 21:33:27 +0200 Subject: [PATCH 03/46] Add a first exercise to tutorial01 This follows the steps provided in https://github.com/espressomd/espresso/wiki/Documentation#tutorials --- .../01-lennard_jones/01-lennard_jones.ipynb | 36 ++++++++++++++++++- 1 file changed, 35 insertions(+), 1 deletion(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index 1b3f102795c..bd74175c388 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -170,13 +170,47 @@ "The next step would be to create an instance of the System class. This instance is used as a handle to the simulation system. At any time, only one instance of the System class can exist." ] }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "**Exercise:**\n", + "Create an instance of an espresso system and store it in a variable called `system`,\n", + "See [module documentaion](http://espressomd.org/html/doc/espressomd.html#module-espressomd.system)\n", + "Use `BOX_L` as box length." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "system = espressomd.System(box_l=BOX_L)\n", + "```" + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "code_folding": [] + }, + "outputs": [], "source": [ - "system = espressomd.System(box_l=BOX_L)" + "# Test solution of Exercise 1\n", + "assert(isinstance(system, espressomd.System))" ] }, { From 61abad67b0eaabce7ec389ff6e4f183bf1e83a30 Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Thu, 3 Sep 2020 21:58:22 +0200 Subject: [PATCH 04/46] Improve links to doc --- doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index bd74175c388..d8e816a9318 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -178,9 +178,11 @@ }, "source": [ "**Exercise:**\n", - "Create an instance of an espresso system and store it in a variable called `system`,\n", - "See [module documentaion](http://espressomd.org/html/doc/espressomd.html#module-espressomd.system)\n", - "Use `BOX_L` as box length." + "\n", + "* Create an instance of an espresso system and store it in a variable called `system`;\n", + " use `BOX_L` as box length.\n", + "\n", + "See [ESPResSo documentation](http://espressomd.org/html/doc/system_setup.html) and [module documentaion](http://espressomd.org/html/doc/espressomd.html#module-espressomd.system)." ] }, { From bb79a90cb6ed22112a301d351319dce4fb7d165f Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Thu, 3 Sep 2020 22:05:45 +0200 Subject: [PATCH 05/46] Convert irrelevant code cell to hint in markdown --- .../01-lennard_jones/01-lennard_jones.ipynb | 46 +++++++++++-------- 1 file changed, 27 insertions(+), 19 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index d8e816a9318..5a15c6d7e02 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -246,22 +246,12 @@ "\n", "Simulations can be carried out in different thermodynamic ensembles such as NVE (particle __N__umber, __V__olume, __E__nergy), NVT (particle __N__umber, __V__olume, __T__emperature) or NPT-isotropic (particle __N__umber, __P__ressure, __T__emperature).\n", "\n", - "The NVE ensemble is simulated without a thermostat. A previously enabled thermostat can be switched off as follows:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "system.thermostat.turn_off()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ + "The NVE ensemble is simulated without a thermostat. A previously enabled thermostat can be switched off as follows:\n", + "\n", + "``` python\n", + "system.thermostat.turn_off()\n", + "```\n", + "\n", "The NVT and NPT ensembles require a thermostat. In this tutorial, we use the Langevin thermostat.\n", "\n", "In ESPResSo, the thermostat is set as follows:" @@ -723,9 +713,9 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3.8.5 64-bit", "language": "python", - "name": "python3" + "name": "python38564bit31b7b39ba35a4507a8c61be7a2214877" }, "language_info": { "codemirror_mode": { @@ -737,7 +727,25 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.8" + "version": "3.8.5" + }, + "latex_envs": { + "LaTeX_envs_menu_present": true, + "autoclose": false, + "autocomplete": true, + "bibliofile": "biblio.bib", + "cite_by": "apalike", + "current_citInitial": 1, + "eqLabelWithNumbers": true, + "eqNumInitial": 1, + "hotkeys": { + "equation": "Ctrl-E", + "itemize": "Ctrl-I" + }, + "labels_anchors": false, + "latex_user_defs": false, + "report_style_numbering": false, + "user_envs_cfg": false } }, "nbformat": 4, From 405014db01b6ad573408d527b94d686597cc8e83 Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Thu, 3 Sep 2020 22:24:12 +0200 Subject: [PATCH 06/46] Exercise thermostat langevin --- .../01-lennard_jones/01-lennard_jones.ipynb | 30 ++++++++++++++----- 1 file changed, 23 insertions(+), 7 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index 5a15c6d7e02..fb39a1c04c8 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -258,21 +258,37 @@ ] }, { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, "source": [ - "system.thermostat.set_langevin(kT=TEMPERATURE, gamma=GAMMA, seed=42)" + "**Exercise:**\n", + "\n", + "* Use system.thermostat.set_langevin() to turn on the Langevin thermostat.\n", + "\n", + "Set the temperature to `temperature` and damping coefficient to `GAMMA`.\n", + "\n", + "For details see the [online documentation](http://espressomd.org/html/doc/espressomd.html#module-espressomd.thermostat)." ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "solution2": "hidden" + }, "source": [ - "Use a Langevin thermostat (NVT or NPT ensemble) with temperature set to `temperature` and damping coefficient to `GAMMA`." + "```python\n", + "system.thermostat.set_langevin(kT=TEMPERATURE, gamma=GAMMA, seed=42)\n", + "```" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [] + }, { "cell_type": "markdown", "metadata": {}, From 364c7c034c9ec5dc00c85fd936f641ca8e56c03d Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Thu, 3 Sep 2020 23:04:37 +0200 Subject: [PATCH 07/46] WIP: Add exercise to create particles Breaks tests: According to the documentation of Exercise2 (https://github.com/ipython-contrib/jupyter_contrib_nbextensions/tree/master/src/jupyter_contrib_nbextensions/nbextensions/exercise2) alternative solutions are consecutive cells in a solution block. However, converting using our CMake infrastructure *all* resulting code blocks are exectuded. --- .../01-lennard_jones/01-lennard_jones.ipynb | 72 +++++++++++++++++-- 1 file changed, 67 insertions(+), 5 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index fb39a1c04c8..8c57727ab28 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -300,16 +300,78 @@ "Particles can be grouped into several types, so that, e.g., a binary fluid can be simulated. Particle types are identified by integer ids, which are set via the particles' type attribute. If it is not specified, zero is implied." ] }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "**Exercise:**\n", + "\n", + "* Create N_PART particles at random positions.\n", + "\n", + " Use [system.part.add()](link).\n", + " \n", + " Either write a loop or use an (N_PART x 3) array for positions.\n", + " Use numpy.random.random() to generate random numbers." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "# Add particles to the simulation box at random positions\n", + "for i in range(N_PART):\n", + " system.part.add(type=0, pos=np.random.random(3) * system.box_l)\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "# Alternative solution:\n", + "system.part.add(type=[0]*N_PART, pos=np.random.random((N_PART,3)) * system.box_l)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Test that now we have indeed N_PART particles in the system\n", + "assert(len(system.part) == N_PART)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The particle properties can be accessed using standard numpy slicing syntax:" + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "# Add particles to the simulation box at random positions\n", - "for i in range(N_PART):\n", - " system.part.add(type=0, pos=np.random.random(3) * system.box_l)\n", - "\n", "# Access position of a single particle\n", "print(\"position of particle with id 0:\", system.part[0].pos)\n", "\n", @@ -327,7 +389,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Many objects in ESPResSo have a string representation, and thus can be displayed via python's print function:" + "Many objects in ESPResSo have a string representation, and thus can be displayed via python's `print` function:" ] }, { From d3f8514da305a7cdcc61467b3b528ae69722ddcc Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Thu, 3 Sep 2020 23:07:07 +0200 Subject: [PATCH 08/46] Init notes for tutor file --- doc/tutorials/01-lennard_jones/NotesForTutor.md | 4 ++++ 1 file changed, 4 insertions(+) create mode 100644 doc/tutorials/01-lennard_jones/NotesForTutor.md diff --git a/doc/tutorials/01-lennard_jones/NotesForTutor.md b/doc/tutorials/01-lennard_jones/NotesForTutor.md new file mode 100644 index 00000000000..18cee43b2b0 --- /dev/null +++ b/doc/tutorials/01-lennard_jones/NotesForTutor.md @@ -0,0 +1,4 @@ +# Hints for Tutors: Tutorial 01 - Lennard Jones + +* Re-explain Langevin equation +* Mention efficiency (loops and lists) vs. numpy arrays at the cost of some readability. \ No newline at end of file From 254812d9e99c99816b47b2dcfbe1a9fa22d1591c Mon Sep 17 00:00:00 2001 From: Rudolf Weeber Date: Fri, 4 Sep 2020 12:23:14 +0200 Subject: [PATCH 09/46] LJ Tutorial: Tutor's notes --- .../01-lennard_jones/NotesForTutor.md | 29 ++++++++++++++++++- 1 file changed, 28 insertions(+), 1 deletion(-) diff --git a/doc/tutorials/01-lennard_jones/NotesForTutor.md b/doc/tutorials/01-lennard_jones/NotesForTutor.md index 18cee43b2b0..bf8174675a1 100644 --- a/doc/tutorials/01-lennard_jones/NotesForTutor.md +++ b/doc/tutorials/01-lennard_jones/NotesForTutor.md @@ -1,4 +1,31 @@ # Hints for Tutors: Tutorial 01 - Lennard Jones +## Learning objectives (physics): + +After the tutorial, students should be able to + +* explain + * role of Lennard-Jones potential + * role of radial distribution function + * connection between mean squared displacement and diffusion coefficient + * how to use an auto correlation function to estimate correlation times and how that affects error estimation + + +## Learning objectives (Espresso) + +After the tutorial, students should be able to + +* instance a simulation object and set basic parameters like time step +* add particles to a simulation, +* access and change particle properties +* set up non-bonded interactions +* explain the role of particle type +* remove overlap using steepest descent +* setting up a thermostat +* integrating and collecting data +* explain the basic concept of observables and accumulators + +## Points to mention throughout the tutorial + * Re-explain Langevin equation -* Mention efficiency (loops and lists) vs. numpy arrays at the cost of some readability. \ No newline at end of file +* Mention efficiency (loops and lists) vs. numpy arrays at the cost of some readability. From 02013c51ede5efbb13caff420231f175a4ee18c1 Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Fri, 4 Sep 2020 11:22:12 +0200 Subject: [PATCH 10/46] Unify notation --- .../01-lennard_jones/01-lennard_jones.ipynb | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index 8c57727ab28..64088befa02 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -179,8 +179,8 @@ "source": [ "**Exercise:**\n", "\n", - "* Create an instance of an espresso system and store it in a variable called `system`;\n", - " use `BOX_L` as box length.\n", + "* Create an instance of an espresso system and store it in a variable called system;\n", + " use BOX_L as box length.\n", "\n", "See [ESPResSo documentation](http://espressomd.org/html/doc/system_setup.html) and [module documentaion](http://espressomd.org/html/doc/espressomd.html#module-espressomd.system)." ] @@ -266,9 +266,9 @@ "source": [ "**Exercise:**\n", "\n", - "* Use system.thermostat.set_langevin() to turn on the Langevin thermostat.\n", + "* Use system.thermostat.set_langevin() to turn on the Langevin thermostat.\n", "\n", - "Set the temperature to `temperature` and damping coefficient to `GAMMA`.\n", + " Set the temperature to TEMPERATURE and damping coefficient to GAMMA.\n", "\n", "For details see the [online documentation](http://espressomd.org/html/doc/espressomd.html#module-espressomd.thermostat)." ] @@ -309,12 +309,12 @@ "source": [ "**Exercise:**\n", "\n", - "* Create N_PART particles at random positions.\n", + "* Create N_PART particles at random positions.\n", "\n", " Use [system.part.add()](link).\n", " \n", - " Either write a loop or use an (N_PART x 3) array for positions.\n", - " Use numpy.random.random() to generate random numbers." + " Either write a loop or use an (N_PART x 3) array for positions.\n", + " Use numpy.random.random() to generate random numbers." ] }, { @@ -389,7 +389,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Many objects in ESPResSo have a string representation, and thus can be displayed via python's `print` function:" + "Many objects in ESPResSo have a string representation, and thus can be displayed via python's print function:" ] }, { From d90032495b0368c388ce4c236d59cd3369ef3ce5 Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Fri, 4 Sep 2020 11:25:02 +0200 Subject: [PATCH 11/46] Reformat alternative solution to avoid convertion to python code --- doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index 64088befa02..3f35c11462d 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -336,8 +336,8 @@ "solution2": "hidden" }, "source": [ + "***Alternative solution:***\n", "```python\n", - "# Alternative solution:\n", "system.part.add(type=[0]*N_PART, pos=np.random.random((N_PART,3)) * system.box_l)\n", "```" ] From 1c5e908606605d73af26bdc0c78233d7984f352f Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Fri, 4 Sep 2020 12:29:04 +0200 Subject: [PATCH 12/46] Replace force_cap warmup by steepest descent and add exercises --- .../01-lennard_jones/01-lennard_jones.ipynb | 98 +++++++++++++++---- 1 file changed, 77 insertions(+), 21 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index 3f35c11462d..541043d0df1 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -418,20 +418,49 @@ "source": [ "LJ_EPS = 1.0\n", "LJ_SIG = 1.0\n", - "LJ_CUT = 2.5 * LJ_SIG\n", - "LJ_CAP = 0.5\n", + "LJ_CUT = 2.5 * LJ_SIG" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "**Exercise:**\n", + "\n", + "* Setup a Lennard-Jones interaction with $\\epsilon=$LJ_EPS and $\\sigma=$LJ_SIG which are cut and shifted to zero at $r_\\mathrm{c}=$LJ_CUT$\\times\\sigma$.\n", + "\n", + " *Hint:* Use the [espressomd.system](http://espressomd.org/html/doc/espressomd.html#module-espressomd.system).non_bonded_inter handle object to instantiate a [espressomd.interactions.LennardJonesInteraction](http://espressomd.org/html/doc/espressomd.html#module-espressomd.interactions)." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", "system.non_bonded_inter[0, 0].lennard_jones.set_params(\n", " epsilon=LJ_EPS, sigma=LJ_SIG, cutoff=LJ_CUT, shift='auto')\n", - "system.force_cap = LJ_CAP" + "```" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Warmup\n", + "### Energy minimization\n", "\n", - "In many cases, including this tutorial, particles are initially placed randomly in the simulation box. It is therefore possible that particles overlap, resulting in a huge repulsive force between them. In this case, integrating the equations of motion would not be numerically stable. Hence, it is necessary to remove this overlap. This is done by limiting the maximum force between two particles, integrating the equations of motion, and increasing the force limit step by step as follows:" + "In many cases, including this tutorial, particles are initially placed randomly in the simulation box. It is therefore possible that particles overlap, resulting in a huge repulsive force between them. In this case, integrating the equations of motion would not be numerically stable. Hence, it is necessary to remove this overlap. This is typically done by performing a steepest descent minimization of the potential energy until a maximal force criterion is reached." ] }, { @@ -440,42 +469,69 @@ "metadata": {}, "outputs": [], "source": [ - "WARM_STEPS = 100\n", - "WARM_N_TIME = 2000\n", - "MIN_DIST = 0.87\n", - "\n", - "i = 0\n", - "act_min_dist = system.analysis.min_dist()\n", - "while i < WARM_N_TIME and act_min_dist < MIN_DIST:\n", - " system.integrator.run(WARM_STEPS)\n", - " act_min_dist = system.analysis.min_dist()\n", - " i += 1\n", - " LJ_CAP += 1.0\n", - " system.force_cap = LJ_CAP" + "F_MAX = 50\n", + "DAMPING = 30\n", + "MAX_STEPS = 10000\n", + "MAX_DISPLACEMENT = 0.01 * LJ_SIG" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, "source": [ - "### Integrating equations of motion and taking measurements\n", + "**Exercise:**\n", + "\n", + "* Use [espressomd.minimize_energy](http://espressomd.org/html/doc/espressomd.html#module-espressomd.minimize_energy) to make sure that there are no forces > F_MAX.\n", + " Limit the number of steepest descent steps to MAX_STEPS \n", + " and the maximal displacement to MAX_DISPLACEMENT.\n", + " A damping constant gamma = DAMPING usually is a good choice.\n", + "* Monitor the system total energy before and after the steepest descent to make sure the system is\n", + " properly relaxed. See the [espressomd.analyze](http://espressomd.org/html/doc/espressomd.html#module-espressomd.analyze) documentation for details." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "solution2": "hidden" + }, + "outputs": [], + "source": [ + "from espressomd.minimize_energy import steepest_descent\n", "\n", - "Once warmup is done, the force capping is switched off by setting it to zero." + "energy = system.analysis.energy()\n", + "print(\"Before Minimization: E_total = {}\".format(energy['total']))\n", + "steepest_descent(system, f_max=F_MAX, gamma=DAMPING, max_steps=MAX_STEPS, max_displacement=MAX_DISPLACEMENT)\n", + "energy = system.analysis.energy()\n", + "print(\"After Minimization: E_total = {}\".format(energy['total']))" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "system.force_cap = 0" + "# check that after the exercise the total energy is negative\n", + "assert(system.analysis.energy()['total'] < 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ + "### Integrating equations of motion and taking measurements\n", + "\n", "At this point, we have set the necessary environment and warmed up our system. Now, we integrate the equations of motion and take measurements. We first plot the radial distribution function which describes how the density varies as a function of distance from a tagged particle. The radial distribution function is averaged over several measurements to reduce noise.\n", "\n", "The potential and kinetic energies can be monitored using the analysis method system.analysis.energy().\n", From e8cb3b57d8645fb47fdc1031a68000ab59959f57 Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Fri, 4 Sep 2020 12:35:03 +0200 Subject: [PATCH 13/46] More points for tutors --- doc/tutorials/01-lennard_jones/NotesForTutor.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/doc/tutorials/01-lennard_jones/NotesForTutor.md b/doc/tutorials/01-lennard_jones/NotesForTutor.md index bf8174675a1..d1c6ba0ae90 100644 --- a/doc/tutorials/01-lennard_jones/NotesForTutor.md +++ b/doc/tutorials/01-lennard_jones/NotesForTutor.md @@ -29,3 +29,6 @@ After the tutorial, students should be able to * Re-explain Langevin equation * Mention efficiency (loops and lists) vs. numpy arrays at the cost of some readability. +* Non-bonded interactions are defined between pairs of particle types. + Mention that type=0 is implied when creating particles without specifying one. +* Explain Lennard-Jones cutoff and shift. \ No newline at end of file From 756214951c9fbf57a492d75f00a500f143b4fc66 Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Fri, 4 Sep 2020 12:39:11 +0200 Subject: [PATCH 14/46] Remove alternative solution path --- .../01-lennard_jones/01-lennard_jones.ipynb | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index 541043d0df1..faea0f775ae 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -325,19 +325,6 @@ "source": [ "```python\n", "# Add particles to the simulation box at random positions\n", - "for i in range(N_PART):\n", - " system.part.add(type=0, pos=np.random.random(3) * system.box_l)\n", - "```" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "solution2": "hidden" - }, - "source": [ - "***Alternative solution:***\n", - "```python\n", "system.part.add(type=[0]*N_PART, pos=np.random.random((N_PART,3)) * system.box_l)\n", "```" ] From 1c2f127499726ea790ea9a7923fb15c2b82dd322 Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Fri, 4 Sep 2020 14:10:16 +0200 Subject: [PATCH 15/46] Fix typo MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Jean-Noël Grad --- doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index faea0f775ae..dfb2c9acf8a 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -182,7 +182,7 @@ "* Create an instance of an espresso system and store it in a variable called system;\n", " use BOX_L as box length.\n", "\n", - "See [ESPResSo documentation](http://espressomd.org/html/doc/system_setup.html) and [module documentaion](http://espressomd.org/html/doc/espressomd.html#module-espressomd.system)." + "See [ESPResSo documentation](http://espressomd.org/html/doc/system_setup.html) and [module documentation](http://espressomd.org/html/doc/espressomd.html#module-espressomd.system)." ] }, { From e7a961b3a173e28302f38d9ec1715c7cba751a98 Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Fri, 4 Sep 2020 14:12:23 +0200 Subject: [PATCH 16/46] Fix typo MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Jean-Noël Grad --- doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index dfb2c9acf8a..1d3e711fadb 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -180,7 +180,7 @@ "**Exercise:**\n", "\n", "* Create an instance of an espresso system and store it in a variable called system;\n", - " use BOX_L as box length.\n", + " use BOX_L as box length.\n", "\n", "See [ESPResSo documentation](http://espressomd.org/html/doc/system_setup.html) and [module documentation](http://espressomd.org/html/doc/espressomd.html#module-espressomd.system)." ] From e42007ce91419777e3a404877f1cbd040df74ca0 Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Fri, 4 Sep 2020 14:12:43 +0200 Subject: [PATCH 17/46] Fix typo MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Jean-Noël Grad --- doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index 1d3e711fadb..a7223961690 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -471,7 +471,7 @@ "source": [ "**Exercise:**\n", "\n", - "* Use [espressomd.minimize_energy](http://espressomd.org/html/doc/espressomd.html#module-espressomd.minimize_energy) to make sure that there are no forces > F_MAX.\n", + "* Use [espressomd.minimize_energy](http://espressomd.org/html/doc/espressomd.html#module-espressomd.minimize_energy) to make sure that there are no forces > F_MAX.\n", " Limit the number of steepest descent steps to MAX_STEPS \n", " and the maximal displacement to MAX_DISPLACEMENT.\n", " A damping constant gamma = DAMPING usually is a good choice.\n", From 426f8e6ce3a6e4454d46c5bdba5fbf7f60afec90 Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Fri, 4 Sep 2020 14:13:34 +0200 Subject: [PATCH 18/46] Use ESPResSo in favor of Espresso MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Jean-Noël Grad --- doc/tutorials/01-lennard_jones/NotesForTutor.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/NotesForTutor.md b/doc/tutorials/01-lennard_jones/NotesForTutor.md index d1c6ba0ae90..f369cf965c3 100644 --- a/doc/tutorials/01-lennard_jones/NotesForTutor.md +++ b/doc/tutorials/01-lennard_jones/NotesForTutor.md @@ -11,7 +11,7 @@ After the tutorial, students should be able to * how to use an auto correlation function to estimate correlation times and how that affects error estimation -## Learning objectives (Espresso) +## Learning objectives (ESPResSo) After the tutorial, students should be able to @@ -31,4 +31,4 @@ After the tutorial, students should be able to * Mention efficiency (loops and lists) vs. numpy arrays at the cost of some readability. * Non-bonded interactions are defined between pairs of particle types. Mention that type=0 is implied when creating particles without specifying one. -* Explain Lennard-Jones cutoff and shift. \ No newline at end of file +* Explain Lennard-Jones cutoff and shift. From 4ba0765026931614f82ac5d9f17327f8b3240cea Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Fri, 4 Sep 2020 17:27:21 +0200 Subject: [PATCH 19/46] Use steepest descent integrator This also requires to shift the initialization of the thermostat to after the steepest descent. --- .../01-lennard_jones/01-lennard_jones.ipynb | 152 ++++++++++-------- .../01-lennard_jones/NotesForTutor.md | 1 + 2 files changed, 88 insertions(+), 65 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index a7223961690..df3d4a35efe 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -238,57 +238,6 @@ "system.cell_system.skin = SKIN" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Choosing the thermodynamic ensemble, thermostat\n", - "\n", - "Simulations can be carried out in different thermodynamic ensembles such as NVE (particle __N__umber, __V__olume, __E__nergy), NVT (particle __N__umber, __V__olume, __T__emperature) or NPT-isotropic (particle __N__umber, __P__ressure, __T__emperature).\n", - "\n", - "The NVE ensemble is simulated without a thermostat. A previously enabled thermostat can be switched off as follows:\n", - "\n", - "``` python\n", - "system.thermostat.turn_off()\n", - "```\n", - "\n", - "The NVT and NPT ensembles require a thermostat. In this tutorial, we use the Langevin thermostat.\n", - "\n", - "In ESPResSo, the thermostat is set as follows:" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "solution2": "hidden", - "solution2_first": true - }, - "source": [ - "**Exercise:**\n", - "\n", - "* Use system.thermostat.set_langevin() to turn on the Langevin thermostat.\n", - "\n", - " Set the temperature to TEMPERATURE and damping coefficient to GAMMA.\n", - "\n", - "For details see the [online documentation](http://espressomd.org/html/doc/espressomd.html#module-espressomd.thermostat)." - ] - }, - { - "cell_type": "markdown", - "metadata": { - "solution2": "hidden" - }, - "source": [ - "```python\n", - "system.thermostat.set_langevin(kT=TEMPERATURE, gamma=GAMMA, seed=42)\n", - "```" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [] - }, { "cell_type": "markdown", "metadata": {}, @@ -456,10 +405,13 @@ "metadata": {}, "outputs": [], "source": [ - "F_MAX = 50\n", + "F_TOL = 1e-2\n", + "E_TOL = 1e-5\n", + "F_MAX = 0 # Use a relative convergence criterion only\n", "DAMPING = 30\n", "MAX_STEPS = 10000\n", - "MAX_DISPLACEMENT = 0.01 * LJ_SIG" + "MAX_DISPLACEMENT = 0.01 * LJ_SIG\n", + "EMSTEP = 10" ] }, { @@ -471,12 +423,16 @@ "source": [ "**Exercise:**\n", "\n", - "* Use [espressomd.minimize_energy](http://espressomd.org/html/doc/espressomd.html#module-espressomd.minimize_energy) to make sure that there are no forces > F_MAX.\n", - " Limit the number of steepest descent steps to MAX_STEPS \n", - " and the maximal displacement to MAX_DISPLACEMENT.\n", + "* Use [espressomd.integrate](http://espressomd.org/html/doc/espressomd.html#module-espressomd.minimize_energy).set_steepest_descent to relax the initial configuration.\n", + " Use a maximal displacement to MAX_DISPLACEMENT.\n", " A damping constant gamma = DAMPING usually is a good choice.\n", - "* Monitor the system total energy before and after the steepest descent to make sure the system is\n", - " properly relaxed. See the [espressomd.analyze](http://espressomd.org/html/doc/espressomd.html#module-espressomd.analyze) documentation for details." + "* Use the relative change of the system total energy and the maximal force as a convergence criterion.\n", + " See the [espressomd.analyze](http://espressomd.org/html/doc/espressomd.html#module-espressomd.analyze) documentation on how to obtain energies and the [espressomd.particle_data module](http://espressomd.org/html/doc/espressomd.html#module-espressomd.particle_data) for the forces.\n", + " The steepest descent has converged if the relative force change < F_TOL\n", + "* Break the minimization loop after a maximal number of MAX_STEPS steps or if convergence is achieved.\n", + " Check for convergence every EMSTEP steps.\n", + " \n", + "***Hint:*** To obtain the initial fores one has to initialize integrator using integ_steps=0." ] }, { @@ -487,13 +443,26 @@ }, "outputs": [], "source": [ - "from espressomd.minimize_energy import steepest_descent\n", - "\n", - "energy = system.analysis.energy()\n", - "print(\"Before Minimization: E_total = {}\".format(energy['total']))\n", - "steepest_descent(system, f_max=F_MAX, gamma=DAMPING, max_steps=MAX_STEPS, max_displacement=MAX_DISPLACEMENT)\n", - "energy = system.analysis.energy()\n", - "print(\"After Minimization: E_total = {}\".format(energy['total']))" + "# Set up steepest descent integration\n", + "system.integrator.set_steepest_descent(f_max=F_MAX,\n", + " gamma=DAMPING,\n", + " max_displacement=MAX_DISPLACEMENT)\n", + "# Initialize integrator to obtain initial forces\n", + "system.integrator.run(0)\n", + "\n", + "i = 0\n", + "while i < MAX_STEPS//EMSTEP:\n", + " prev_maxforce = maxforce\n", + " prev_energy = system.analysis.energy()['total']\n", + " system.integrator.run(EMSTEP)\n", + " maxforce = np.max(np.abs(system.part[:].f))\n", + " relforce = np.abs((maxforce-prev_maxforce)/prev_maxforce)\n", + " energy = system.analysis.energy()['total']\n", + " relener = np.abs((energy-prev_energy)/prev_energy)\n", + " print(\"minimization step: {:4.0f}\\tmax. rel. force change:{:+3.3e}\\trel. energy change:{:+3.3e}\".format((i+1)*EMSTEP,relforce, relener))\n", + " if relforce < F_TOL or relener < E_TOL:\n", + " break\n", + " i += 1" ] }, { @@ -513,6 +482,59 @@ "assert(system.analysis.energy()['total'] < 0)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Choosing the thermodynamic ensemble, thermostat\n", + "\n", + "Simulations can be carried out in different thermodynamic ensembles such as NVE (particle __N__umber, __V__olume, __E__nergy), NVT (particle __N__umber, __V__olume, __T__emperature) or NPT-isotropic (particle __N__umber, __P__ressure, __T__emperature).\n", + "\n", + "The NVE ensemble is simulated without a thermostat. A previously enabled thermostat can be switched off as follows:\n", + "\n", + "``` python\n", + "system.thermostat.turn_off()\n", + "```\n", + "\n", + "The NVT and NPT ensembles require a thermostat. In this tutorial, we use the Langevin thermostat.\n", + "\n", + "In ESPResSo, the thermostat is set as follows:" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "**Exercise:**\n", + "\n", + "* Use system.thermostat.set_langevin() to turn on the Langevin thermostat.\n", + "\n", + " Set the temperature to TEMPERATURE and damping coefficient to GAMMA.\n", + "\n", + "For details see the [online documentation](http://espressomd.org/html/doc/espressomd.html#module-espressomd.thermostat)." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "system.thermostat.set_langevin(kT=TEMPERATURE, gamma=GAMMA, seed=42)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "markdown", "metadata": {}, diff --git a/doc/tutorials/01-lennard_jones/NotesForTutor.md b/doc/tutorials/01-lennard_jones/NotesForTutor.md index f369cf965c3..c401e145e5b 100644 --- a/doc/tutorials/01-lennard_jones/NotesForTutor.md +++ b/doc/tutorials/01-lennard_jones/NotesForTutor.md @@ -32,3 +32,4 @@ After the tutorial, students should be able to * Non-bonded interactions are defined between pairs of particle types. Mention that type=0 is implied when creating particles without specifying one. * Explain Lennard-Jones cutoff and shift. +* Explain steepest descent. \ No newline at end of file From 5104cdfecd5873f69402a8bbc1b499deda09664a Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Wed, 9 Sep 2020 20:53:51 +0200 Subject: [PATCH 20/46] Improve overlap removal --- .../01-lennard_jones/01-lennard_jones.ipynb | 12 ++++++++---- doc/tutorials/01-lennard_jones/NotesForTutor.md | 5 ++++- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index df3d4a35efe..f301d40c503 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -396,7 +396,12 @@ "source": [ "### Energy minimization\n", "\n", - "In many cases, including this tutorial, particles are initially placed randomly in the simulation box. It is therefore possible that particles overlap, resulting in a huge repulsive force between them. In this case, integrating the equations of motion would not be numerically stable. Hence, it is necessary to remove this overlap. This is typically done by performing a steepest descent minimization of the potential energy until a maximal force criterion is reached." + "In many cases, including this tutorial, particles are initially placed randomly in the simulation box. It is therefore possible that particles overlap, resulting in a huge repulsive force between them. In this case, integrating the equations of motion would not be numerically stable. Hence, it is necessary to remove this overlap. This is typically done by performing a steepest descent minimization of the potential energy until a maximal force criterion is reached.\n", + "\n", + "**Note:**\n", + "Making sure a system is well equilibrated highly depends on the system's details.\n", + "In most cases a relative convergence criterion on the forces and/or energies works well but you might have to make sure that the total force is smaller than a threshold value f_max at the end of the minimization.\n", + "Depending on the simulated system other strategies like simulations with small time step or capped forces might be necessary." ] }, { @@ -407,7 +412,6 @@ "source": [ "F_TOL = 1e-2\n", "E_TOL = 1e-5\n", - "F_MAX = 0 # Use a relative convergence criterion only\n", "DAMPING = 30\n", "MAX_STEPS = 10000\n", "MAX_DISPLACEMENT = 0.01 * LJ_SIG\n", @@ -444,7 +448,7 @@ "outputs": [], "source": [ "# Set up steepest descent integration\n", - "system.integrator.set_steepest_descent(f_max=F_MAX,\n", + "system.integrator.set_steepest_descent(f_max=0, # use a relative convergence criterion only\n", " gamma=DAMPING,\n", " max_displacement=MAX_DISPLACEMENT)\n", "# Initialize integrator to obtain initial forces\n", @@ -455,7 +459,7 @@ " prev_maxforce = maxforce\n", " prev_energy = system.analysis.energy()['total']\n", " system.integrator.run(EMSTEP)\n", - " maxforce = np.max(np.abs(system.part[:].f))\n", + " maxforce = np.max(np.linalg.norm(system.part[:].f, axis=1))\n", " relforce = np.abs((maxforce-prev_maxforce)/prev_maxforce)\n", " energy = system.analysis.energy()['total']\n", " relener = np.abs((energy-prev_energy)/prev_energy)\n", diff --git a/doc/tutorials/01-lennard_jones/NotesForTutor.md b/doc/tutorials/01-lennard_jones/NotesForTutor.md index c401e145e5b..71bbe02b97f 100644 --- a/doc/tutorials/01-lennard_jones/NotesForTutor.md +++ b/doc/tutorials/01-lennard_jones/NotesForTutor.md @@ -32,4 +32,7 @@ After the tutorial, students should be able to * Non-bonded interactions are defined between pairs of particle types. Mention that type=0 is implied when creating particles without specifying one. * Explain Lennard-Jones cutoff and shift. -* Explain steepest descent. \ No newline at end of file +* Explain steepest descent: In general overlap removal can + be done in multiple ways (warmup with force capping, + small time step, etc.). + Additionally the steepest descent algorithm can get trapped in local minima and the convergence criterion is system-dependent. \ No newline at end of file From 0ed1854bd35b23264131ad7ddb69f3c22d08102e Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Mon, 14 Sep 2020 18:07:35 +0200 Subject: [PATCH 21/46] Fix energy minimization --- .../01-lennard_jones/01-lennard_jones.ipynb | 32 +++++++++++++------ .../01-lennard_jones/NotesForTutor.md | 3 +- 2 files changed, 25 insertions(+), 10 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index f301d40c503..39a4b98be27 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -231,9 +231,6 @@ "SKIN = 0.4\n", "TIME_STEP = 0.01\n", "\n", - "TEMPERATURE = 0.728\n", - "GAMMA=1.0\n", - "\n", "system.time_step = TIME_STEP\n", "system.cell_system.skin = SKIN" ] @@ -440,24 +437,27 @@ ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": { + "code_folding": [], "solution2": "hidden" }, - "outputs": [], "source": [ + "```python\n", "# 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", + "maxforce = np.max(np.linalg.norm(system.part[:].f, axis=1))\n", + "energy = system.analysis.energy()['total']\n", "\n", "i = 0\n", "while i < MAX_STEPS//EMSTEP:\n", " prev_maxforce = maxforce\n", - " prev_energy = system.analysis.energy()['total']\n", + " prev_energy = energy\n", " system.integrator.run(EMSTEP)\n", " maxforce = np.max(np.linalg.norm(system.part[:].f, axis=1))\n", " relforce = np.abs((maxforce-prev_maxforce)/prev_maxforce)\n", @@ -466,7 +466,8 @@ " print(\"minimization step: {:4.0f}\\tmax. rel. force change:{:+3.3e}\\trel. energy change:{:+3.3e}\".format((i+1)*EMSTEP,relforce, relener))\n", " if relforce < F_TOL or relener < E_TOL:\n", " break\n", - " i += 1" + " i += 1\n", + "```" ] }, { @@ -505,6 +506,17 @@ "In ESPResSo, the thermostat is set as follows:" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Parameters for the Langevin thermostat\n", + "TEMPERATURE = 0.728\n", + "GAMMA=1.0" + ] + }, { "cell_type": "markdown", "metadata": { @@ -514,7 +526,8 @@ "source": [ "**Exercise:**\n", "\n", - "* Use system.thermostat.set_langevin() to turn on the Langevin thermostat.\n", + "* Use system.integrator.set_vv() to use a Velocity Verlet integration scheme and\n", + " system.thermostat.set_langevin() to turn on the Langevin thermostat.\n", "\n", " Set the temperature to TEMPERATURE and damping coefficient to GAMMA.\n", "\n", @@ -528,6 +541,7 @@ }, "source": [ "```python\n", + "system.integrator.set_vv()\n", "system.thermostat.set_langevin(kT=TEMPERATURE, gamma=GAMMA, seed=42)\n", "```" ] diff --git a/doc/tutorials/01-lennard_jones/NotesForTutor.md b/doc/tutorials/01-lennard_jones/NotesForTutor.md index 71bbe02b97f..6dc21cf2cea 100644 --- a/doc/tutorials/01-lennard_jones/NotesForTutor.md +++ b/doc/tutorials/01-lennard_jones/NotesForTutor.md @@ -35,4 +35,5 @@ After the tutorial, students should be able to * Explain steepest descent: In general overlap removal can be done in multiple ways (warmup with force capping, small time step, etc.). - Additionally the steepest descent algorithm can get trapped in local minima and the convergence criterion is system-dependent. \ No newline at end of file + Additionally the steepest descent algorithm can get trapped in local minima and the convergence criterion is system-dependent. +* Mention that accumulators are updated auomatically. \ No newline at end of file From 356b7fa41633180b8052dd960959c4c5960783bf Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Mon, 14 Sep 2020 18:18:20 +0200 Subject: [PATCH 22/46] Rewrite main integration loop - Energy and temperature observables - Calculation of the ACF and error estimate discussion --- .../01-lennard_jones/01-lennard_jones.ipynb | 199 +++++++++++++++--- 1 file changed, 170 insertions(+), 29 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index 39a4b98be27..11ba1374be0 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -559,11 +559,56 @@ "source": [ "### Integrating equations of motion and taking measurements\n", "\n", - "At this point, we have set the necessary environment and warmed up our system. Now, we integrate the equations of motion and take measurements. We first plot the radial distribution function which describes how the density varies as a function of distance from a tagged particle. The radial distribution function is averaged over several measurements to reduce noise.\n", + "At this point, we have set the necessary environment and warmed up our system. \n", + "Now, we integrate the equations of motion and take measurements and simultaneously\n", + "take measurements of relevant quantities.\n", + "\n", + "TODO: correlators-> accumulators. Mentoin correlator, mean variance calcualtor and time series" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "code_folding": [] + }, + "outputs": [], + "source": [ + "# Integration parameters\n", + "SAMPLING_INTERVAL = 100\n", + "SAMPLING_ITERATIONS = 200\n", + "\n", + "# Parameters for the radial distribution function\n", + "R_BINS = 100\n", + "R_MIN= 0.0\n", + "R_MAX = system.box_l[0] / 2.0" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "**Exercise:**\n", + "\n", + "* Set up the main integration loop\n", + " Take SAMPLING_ITERATIONS measurements every SAMPLING_INTERVAL integration steps.\n", + "* Monitor the potential and kinetic energies using the analysis method [system.analysis](http://espressomd.org/html/doc/espressomd.html#module-espressomd.analyze).energy().\n", + " The total energy time series should be stored in a numpy.array etotal\n", + " and the kinetic energy in ekin. Plot the timeseries of the total energy and\n", + " of the instantaneous temperature $T = 3/2 \\times E_\\mathrm{kin}$/N_PART.\n", + " Here $T$ refers to the measured temperature obtained from kinetic energy and the number\n", + " of degrees of freedom in the system. It should fluctuate around the preset temperature of the thermostat.\n", + " Also store the corresponding simulation time in time.\n", + "* Plot the instantaneous temperature and total energy versus time.\n", + "* Estimate the correlation time of the system from the autocorrelation function of the total energy.\n", + "* Extend the integration loop to include the calculation of the radial distribution function which describes how the density varies as a function of distance from a tagged particle. \n", + " The radial distribution function is averaged over several measurements to reduce noise.\n", + " To this end add an observable rdf_obs for the radial distirbution function \n", + "\n", "\n", - "The potential and kinetic energies can be monitored using the analysis method system.analysis.energy().\n", - "Here kinetic_temperature refers to the measured temperature obtained from kinetic energy and the number\n", - "of degrees of freedom in the system. It should fluctuate around the preset temperature of the thermostat.\n", "\n", "The mean square displacement of particle $i$ is given by:\n", "\n", @@ -574,15 +619,131 @@ "and can be calculated using \"observables and correlators\". An observable is an object which takes a measurement on the system. It can depend on parameters specified when the observable is instanced, such as the ids of the particles to be considered." ] }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "etotal = np.zeros(SAMPLING_ITERATIONS)\n", + "ekin = np.zeros(SAMPLING_ITERATIONS)\n", + "time = np.zeros(SAMPLING_ITERATIONS)\n", + "\n", + "for i in range(SAMPLING_ITERATIONS):\n", + " system.integrator.run(SAMPLING_INTERVAL)\n", + " \n", + " # Measure energies\n", + " energies = system.analysis.energy()\n", + " etotal[i] = energies['total']\n", + " ekin[i] = energies['kinetic']\n", + " time[i] = system.time\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "code_folding": [], + "solution2": "hidden" + }, + "source": [ + "```python\n", + "# Plot total energy and instantaneous temperature\n", + "import matplotlib.pyplot as plt\n", + "\n", + "fig1,ax = plt.subplots(num=None, figsize=(10, 6),\n", + " dpi=80, facecolor='w', edgecolor='k')\n", + "fig1.set_tight_layout(False)\n", + "\n", + "t = np.arange(0,SAMPLING_ITERATIONS*SAMPLING_INTERVAL,SAMPLING_INTERVAL)\n", + "ax.plot(t,etotal, lw=2, label='etotal')\n", + "ax.set_xlabel('time step', fontsize=20)\n", + "ax.set_ylabel('total energy', fontsize=20)\n", + "ax.legend(fontsize=16)\n", + "ax.tick_params(axis='both', labelsize=16)\n", + "ax2 = ax.twinx()\n", + "ax2.plot(t,ekin/ (1.5 * N_PART), color='r', ls='--', lw=2,\n", + " label='instantaneous temperature')\n", + "ax2.set_ylabel('instantaneous temperature', fontsize=20)\n", + "ax2.legend(fontsize=16)\n", + "ax2.tick_params(axis='both', labelsize=16)\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "Since the ensemble average $\\langle E_\\text{kin}\\rangle=3/2 N k_B T$ is related to the temperature,\n", + "we may compute the actual temperature of the system via $k_B T= 2/(3N) \\langle E_\\text{kin}\\rangle$.\n", + "The temperature is fixed and does not fluctuate in the NVT ensemble! The instantaneous temperature is\n", + "calculated via $2/(3N) E_\\text{kin}$ (without ensemble averaging), but it is not the temperature of the system." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "code_folding": [], + "solution2": "hidden" + }, + "source": [ + "``` python\n", + "# autocorrelation using numpy.correlate\n", + "def autocor(x):\n", + " mean=x.mean()\n", + " var=np.var(x)\n", + " xp=x-mean\n", + " corr=np.correlate(xp,xp,'full')[len(x)-1:]/var/len(x)\n", + " return corr\n", + "\n", + "fig2,ax = plt.subplots(num=None, figsize=(10, 6),\n", + " dpi=80, facecolor='w', edgecolor='k')\n", + "fig2.set_tight_layout(False)\n", + "\n", + "acf = autocor(etotal)\n", + "ax.plot(time-1,acf, lw=2, color='b')\n", + "ax.plot(time-1,-acf, lw=2, color='b', ls='--')\n", + "\n", + "from scipy.optimize import curve_fit\n", + "fitfn = lambda x,a,t: -a*x/t\n", + "popt,pcov = curve_fit(fitfn, time-1,acf)\n", + "tcor = popt[1]\n", + "x = np.linspace(0,15)\n", + "ax.plot(x,np.exp(-x/tcor), lw=2, color='r', ls='-.')\n", + "ax.set_yscale('log')\n", + "ax.set_ylim(bottom=1e-3)\n", + "ax.set_xlim(0,15)\n", + "ax.tick_params(axis='both', labelsize=16)\n", + "ax.set_ylabel('normalized ACF', fontsize=18)\n", + "ax.set_xlabel('time', fontsize=18)\n", + "\n", + "print (\"Estimated correlation time:\", np.round(tcor,2))\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "From the roughly exponential decay of the correlation function it dawns that only every second sample should be used to calculate statistically independent averages." + ] + }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "code_folding": [ + 0 + ] + }, "outputs": [], "source": [ - "# Integration parameters\n", - "sampling_interval = 100\n", - "sampling_iterations = 200\n", + "# Original solution\n", "\n", "from espressomd.observables import ParticlePositions, RDF\n", "from espressomd.accumulators import Correlator, MeanVarianceCalculator\n", @@ -629,22 +790,12 @@ "r = rdf_obs.bin_centers()" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We now use the plotting library matplotlib available in Python to visualize the measurements." - ] - }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "import matplotlib.pyplot as plt\n", - "plt.ion()" - ] + "source": [] }, { "cell_type": "markdown", @@ -692,16 +843,6 @@ "plt.show()" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Since the ensemble average $\\langle E_\\text{kin}\\rangle=3/2 N k_B T$ is related to the temperature,\n", - "we may compute the actual temperature of the system via $k_B T= 2/(3N) \\langle E_\\text{kin}\\rangle$.\n", - "The temperature is fixed and does not fluctuate in the NVT ensemble! The instantaneous temperature is\n", - "calculated via $2/(3N) E_\\text{kin}$ (without ensemble averaging), but it is not the temperature of the system." - ] - }, { "cell_type": "markdown", "metadata": {}, From 41d1e2d94f4c18b90af9bda097a4f3fb14e41582 Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Mon, 14 Sep 2020 20:54:49 +0200 Subject: [PATCH 23/46] Remove standard deviation estimate --- .../01-lennard_jones/01-lennard_jones.ipynb | 31 ------------------- 1 file changed, 31 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index 11ba1374be0..fcb497a50f8 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -890,37 +890,6 @@ "plt.show()" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Simple Error Estimation on Time Series Data\n", - "\n", - "A simple way to estimate the error of an observable is to use the standard error of the mean (SE) for $N$\n", - "_uncorrelated_ samples:\n", - "\n", - "\\begin{equation}\n", - " SE = \\sqrt{\\frac{\\sigma^2}{N}},\n", - "\\end{equation}\n", - "\n", - "where $\\sigma^2$ is the variance\n", - "\n", - "\\begin{equation}\n", - " \\sigma^2 = \\left\\langle x^2 - \\langle x\\rangle^2 \\right\\rangle\n", - "\\end{equation}" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# calculate the standard error of the mean of the total energy, assuming uncorrelatedness\n", - "standard_error_total_energy = np.sqrt(etotal.var()) / np.sqrt(sampling_iterations)\n", - "print(standard_error_total_energy)" - ] - }, { "cell_type": "markdown", "metadata": {}, From 5c470b6c5d622dbfd408b1137eb9542876e657ad Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Mon, 14 Sep 2020 21:39:18 +0200 Subject: [PATCH 24/46] (Re-)add RDF calculation --- .../01-lennard_jones/01-lennard_jones.ipynb | 153 ++++++++++-------- 1 file changed, 90 insertions(+), 63 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index fcb497a50f8..d545419a1c3 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -563,7 +563,13 @@ "Now, we integrate the equations of motion and take measurements and simultaneously\n", "take measurements of relevant quantities.\n", "\n", - "TODO: correlators-> accumulators. Mentoin correlator, mean variance calcualtor and time series" + "In general, observables extract properties of the particles and the LB fluid and \n", + "return either the raw data or a statistic derived from them.\n", + "We exemplarily calculate the [radial distribution function](https://en.wikipedia.org/wiki/Radial_distribution_function),\n", + "for which we create an [observable](http://espressomd.org/html/doc/espressomd.html#espressomd.observables.RDF) (i.e. a timeseries of radial distribution functions) that are averaged using ESPresSo's [mean-variance calculator](http://espressomd.org/html/doc/analysis.html#mean-variance-calculator).\n", + "\n", + "In order to measure the mean squared particle displacement (related to diffusion), we employ a\n", + "[correlator](http://espressomd.org/html/doc/analysis.html#correlations). The latter takes one or two observables, obtains values from them during the simulation and finally uses a fast correlation algorithm which enables efficient computation of correlation functions spanning many orders of magnitude in the lag time." ] }, { @@ -605,8 +611,12 @@ "* Plot the instantaneous temperature and total energy versus time.\n", "* Estimate the correlation time of the system from the autocorrelation function of the total energy.\n", "* Extend the integration loop to include the calculation of the radial distribution function which describes how the density varies as a function of distance from a tagged particle. \n", - " The radial distribution function is averaged over several measurements to reduce noise.\n", - " To this end add an observable rdf_obs for the radial distirbution function \n", + " The radial distribution function (RDF) is averaged over several measurements to reduce noise.\n", + " To this end add an observable rdf_obs for the radial distirbution function.\n", + " Use min_r=R_MIN, max_r=R_MAX, n_r_bins=R_BINS and create a corresponding [espressomd.accumulators.MeanVarianceCalculator](http://espressomd.org/html/doc/espressomd.html#espressomd.accumulators.MeanVarianceCalculator)\n", + " that gets updated every SAMPLING_INTERVAL timesteps.\n", + " Plot the RDF using the mean of the accumulator.\n", + " *Hint*: The radial coordinates can be obtained using the bin_centers() method of the espressomd.observables.RDF observable.\n", "\n", "\n", "\n", @@ -630,6 +640,15 @@ "ekin = np.zeros(SAMPLING_ITERATIONS)\n", "time = np.zeros(SAMPLING_ITERATIONS)\n", "\n", + "from espressomd.observables import ParticlePositions, RDF\n", + "from espressomd.accumulators import Correlator, MeanVarianceCalculator\n", + "\n", + "# Initialize RDF accumulator\n", + "rdf_obs = RDF(ids1=system.part[:].id, min_r=R_MIN, max_r=R_MAX, n_r_bins=R_BINS)\n", + "rdf_acc = MeanVarianceCalculator(obs=rdf_obs, delta_N=SAMPLING_INTERVAL)\n", + "# Accumulate results automatically during the integration\n", + "system.auto_update_accumulators.add(rdf_acc)\n", + "\n", "for i in range(SAMPLING_ITERATIONS):\n", " system.integrator.run(SAMPLING_INTERVAL)\n", " \n", @@ -638,6 +657,11 @@ " etotal[i] = energies['total']\n", " ekin[i] = energies['kinetic']\n", " time[i] = system.time\n", + " \n", + "# Finalize the accumulator and obtain the results\n", + "rdf = rdf_acc.get_mean()\n", + "r = rdf_obs.bin_centers()\n", + "\n", "```" ] }, @@ -717,8 +741,8 @@ "ax.set_ylim(bottom=1e-3)\n", "ax.set_xlim(0,15)\n", "ax.tick_params(axis='both', labelsize=16)\n", - "ax.set_ylabel('normalized ACF', fontsize=18)\n", - "ax.set_xlabel('time', fontsize=18)\n", + "ax.set_ylabel('normalized ACF', fontsize=20)\n", + "ax.set_xlabel('time', fontsize=20)\n", "\n", "print (\"Estimated correlation time:\", np.round(tcor,2))\n", "```" @@ -733,13 +757,71 @@ "From the roughly exponential decay of the correlation function it dawns that only every second sample should be used to calculate statistically independent averages." ] }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "Plot the experimental radial distribution.\n", + "Empirical radial distribution functions have been determined for pure fluids [5], mixtures [6] and confined fluids [7]. We will compare our distribution $g(r)$ to the theoretical distribution $g(r^*, \\rho^*, T^*)$ of a pure fluid [5]." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "fig3 = plt.figure(num=None, figsize=(10, 6),\n", + " dpi=80, facecolor='w', edgecolor='k')\n", + "fig3.set_tight_layout(False)\n", + "plt.plot(r, rdf, '-', color=\"#A60628\", linewidth=2, alpha=1, label='simulated')\n", + "\n", + "# Comparison to literature\n", + "# reduced units\n", + "T_star = TEMPERATURE / LJ_EPS\n", + "rho_star = DENSITY * LJ_SIG**3\n", + "# expression of the factors Pi from Equations 2-8 with coefficients qi from Table 1\n", + "P = lambda q1, q2, q3, q4, q5, q6, q7, q8, q9: q1 + q2 * np.exp(-q3 * T_star) + q4 * np.exp(-q5 * T_star) + q6 / rho_star + q7 / rho_star**2 + q8 * np.exp(-q3 * T_star) / rho_star**3 + q9 * np.exp(-q5 * T_star) / rho_star**4\n", + "a = P(9.24792, -2.64281, 0.133386, -1.35932, 1.25338, 0.45602, -0.326422, 0.045708, -0.0287681)\n", + "g = P(0.663161, -0.243089, 1.24749, -2.059, 0.04261, 1.65041, -0.343652, -0.037698, 0.008899)\n", + "P = lambda q1, q2, q3, q4, q5, q6, q7, q8: q1 + q2 * np.exp(-q3 * T_star) + q4 * rho_star + q5 * rho_star**2 + q6 * rho_star**3 + q7 * rho_star**4 + q8 * rho_star**5\n", + "c = P(-0.0677912, -1.39505, 0.512625, 36.9323, -36.8061, 21.7353, -7.76671, 1.36342)\n", + "k = P(16.4821, -0.300612, 0.0937844, -61.744, 145.285, -168.087, 98.2181, -23.0583)\n", + "P = lambda q1, q2, q3: q1 + q2 * np.exp(-q3 * rho_star)\n", + "b = P(-8.33289, 2.1714, 1.00063)\n", + "h = P(0.0325039, -1.28792, 2.5487)\n", + "P = lambda q1, q2, q3, q4: q1 + q2 * np.exp(-q3 * rho_star) + q4 * rho_star\n", + "d = P(-26.1615, 27.4846, 1.68124, 6.74296)\n", + "l = P(-6.7293, -59.5002, 10.2466, -0.43596)\n", + "P = lambda q1, q2, q3, q4, q5, q6, q7, q8: (q1 + q2 * rho_star + q3 / T_star + q4 / T_star**2 + q5 / T_star**3) / (q6 + q7 * rho_star + q8 * rho_star**2)\n", + "s = P(1.25225, -1.0179, 0.358564, -0.18533, 0.0482119, 1.27592, -1.78785, 0.634741)\n", + "P = lambda q1, q2, q3, q4, q5, q6: q1 + q2 * np.exp(-q3 * T_star) + q4 / T_star + q5 * rho_star + q6 * rho_star**2\n", + "m = P(-5.668, -3.62671, 0.680654, 0.294481, 0.186395, -0.286954)\n", + "P = lambda q1, q2, q3: q1 + q2 * np.exp(-q3 * T_star)\n", + "n = P(6.01325, 3.84098, 0.60793)\n", + "# theoretical curve\n", + "r_star = np.linspace(0, 3, 301)[1:]\n", + "r_star_cutoff = 1.02 # slightly more than 1 to smooth out the discontinuity in the range [1.0, 1.02]\n", + "theo_rdf = 1 + 1 / r_star**2 * (np.exp(-(a * r_star + b)) * np.sin(c * r_star + d) + np.exp(-(g * r_star + h)) * np.cos(k * r_star + l))\n", + "theo_rdf[np.nonzero(r_star <= r_star_cutoff)] = s * np.exp(-(m * r_star + n)**4)[np.nonzero(r_star <= r_star_cutoff)]\n", + "plt.plot(r_star, theo_rdf, '-', color=\"blue\", linewidth=2, alpha=1, ls='--', label=r'theoretical')\n", + "\n", + "\n", + "plt.xlabel('$r [\\sigma]$', fontsize=20)\n", + "plt.ylabel('$g(r)$', fontsize=20)\n", + "plt.gca().tick_params(axis='both', labelsize=16)\n", + "plt.show()\n", + "```" + ] + }, { "cell_type": "code", "execution_count": null, "metadata": { - "code_folding": [ - 0 - ] + "code_folding": [] }, "outputs": [], "source": [ @@ -911,61 +993,6 @@ "4. Plot the radial distribution functions for all three combinations of particle types. The mixed case will differ significantly, depending on your choice of lj_cut_mixed. Explain these differences." ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Comparison with the literature\n", - "\n", - "Empirical radial distribution functions have been determined for pure fluids [5], mixtures [6] and confined fluids [7]. We will compare our distribution $g(r)$ to the theoretical distribution $g(r^*, \\rho^*, T^*)$ of a pure fluid [5]." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "# reduced units\n", - "T_star = TEMPERATURE / LJ_EPS\n", - "rho_star = DENSITY * LJ_SIG**3\n", - "# expression of the factors Pi from Equations 2-8 with coefficients qi from Table 1\n", - "P = lambda q1, q2, q3, q4, q5, q6, q7, q8, q9: q1 + q2 * np.exp(-q3 * T_star) + q4 * np.exp(-q5 * T_star) + q6 / rho_star + q7 / rho_star**2 + q8 * np.exp(-q3 * T_star) / rho_star**3 + q9 * np.exp(-q5 * T_star) / rho_star**4\n", - "a = P(9.24792, -2.64281, 0.133386, -1.35932, 1.25338, 0.45602, -0.326422, 0.045708, -0.0287681)\n", - "g = P(0.663161, -0.243089, 1.24749, -2.059, 0.04261, 1.65041, -0.343652, -0.037698, 0.008899)\n", - "P = lambda q1, q2, q3, q4, q5, q6, q7, q8: q1 + q2 * np.exp(-q3 * T_star) + q4 * rho_star + q5 * rho_star**2 + q6 * rho_star**3 + q7 * rho_star**4 + q8 * rho_star**5\n", - "c = P(-0.0677912, -1.39505, 0.512625, 36.9323, -36.8061, 21.7353, -7.76671, 1.36342)\n", - "k = P(16.4821, -0.300612, 0.0937844, -61.744, 145.285, -168.087, 98.2181, -23.0583)\n", - "P = lambda q1, q2, q3: q1 + q2 * np.exp(-q3 * rho_star)\n", - "b = P(-8.33289, 2.1714, 1.00063)\n", - "h = P(0.0325039, -1.28792, 2.5487)\n", - "P = lambda q1, q2, q3, q4: q1 + q2 * np.exp(-q3 * rho_star) + q4 * rho_star\n", - "d = P(-26.1615, 27.4846, 1.68124, 6.74296)\n", - "l = P(-6.7293, -59.5002, 10.2466, -0.43596)\n", - "P = lambda q1, q2, q3, q4, q5, q6, q7, q8: (q1 + q2 * rho_star + q3 / T_star + q4 / T_star**2 + q5 / T_star**3) / (q6 + q7 * rho_star + q8 * rho_star**2)\n", - "s = P(1.25225, -1.0179, 0.358564, -0.18533, 0.0482119, 1.27592, -1.78785, 0.634741)\n", - "P = lambda q1, q2, q3, q4, q5, q6: q1 + q2 * np.exp(-q3 * T_star) + q4 / T_star + q5 * rho_star + q6 * rho_star**2\n", - "m = P(-5.668, -3.62671, 0.680654, 0.294481, 0.186395, -0.286954)\n", - "P = lambda q1, q2, q3: q1 + q2 * np.exp(-q3 * T_star)\n", - "n = P(6.01325, 3.84098, 0.60793)\n", - "# theoretical curve\n", - "r_star = np.linspace(0, 3, 301)[1:]\n", - "r_star_cutoff = 1.02 # slightly more than 1 to smooth out the discontinuity in the range [1.0, 1.02]\n", - "theo_rdf = 1 + 1 / r_star**2 * (np.exp(-(a * r_star + b)) * np.sin(c * r_star + d) + np.exp(-(g * r_star + h)) * np.cos(k * r_star + l))\n", - "theo_rdf[np.nonzero(r_star <= r_star_cutoff)] = s * np.exp(-(m * r_star + n)**4)[np.nonzero(r_star <= r_star_cutoff)]\n", - "# plot\n", - "fig = plt.figure(num=None, figsize=(10, 6), dpi=80, facecolor='w', edgecolor='k')\n", - "fig.set_tight_layout(False)\n", - "plt.plot(r_star, theo_rdf, '-', color=\"blue\", linewidth=2, alpha=1, label=r'theoretical')\n", - "plt.plot(r, rdf, '-', color=\"#A60628\", linewidth=2, alpha=1, label='simulation')\n", - "plt.xlabel('r $[\\sigma]$', fontsize='xx-large')\n", - "plt.ylabel('$g(r)$', fontsize='xx-large')\n", - "plt.legend(loc='upper right', fontsize='x-large')\n", - "plt.show()" - ] - }, { "cell_type": "markdown", "metadata": {}, From 697bda0e74a0973c4769d55d210a7cd61d8ba150 Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Wed, 23 Sep 2020 10:46:02 +0200 Subject: [PATCH 25/46] Make RDF comparison work * Don't use LJ shift * Use density/temperture from Verlet 1967 -> Successful comparison with fitted expression and LAMMPS --- .../01-lennard_jones/01-lennard_jones.ipynb | 179 ++++++++++-------- .../01-lennard_jones/NotesForTutor.md | 12 +- 2 files changed, 109 insertions(+), 82 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index d545419a1c3..e6e22959362 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -70,7 +70,7 @@ "where $\\epsilon$ is the depth of the potential well and $\\sigma$ is the (finite) distance at which the inter-particle potential is zero and $r$ is the distance between the particles. The $\\left(\\frac{1}{r}\\right)^{12}$ term describes repulsion and the $(\\frac{1}{r})^{6}$ term describes attraction. The Lennard-Jones potential is an\n", "approximation. The form of the repulsion term has no theoretical justification; the repulsion force should depend exponentially on the distance, but the repulsion term of the L-J formula is more convenient due to the ease and efficiency of computing $r^{12}$ as the square of $r^6$.\n", "\n", - "In practice, the L-J potential is cutoff beyond a specified distance $r_{c}$ and the potential at the cutoff distance is zero.\n", + "In practice, the L-J potential is typically cutoff beyond a specified distance $r_{c}$ and the potential at the cutoff distance is zero.\n", "\n", "
\n", "missing\n", @@ -155,10 +155,18 @@ "outputs": [], "source": [ "# Importing other relevant python modules\n", - "import numpy as np\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ "# System parameters\n", - "N_PART = 100\n", - "DENSITY = 0.5\n", + "N_PART = 200\n", + "DENSITY = 0.75\n", "\n", "BOX_L = np.power(N_PART / DENSITY, 1.0 / 3.0) * np.ones(3)" ] @@ -253,7 +261,7 @@ "solution2_first": true }, "source": [ - "**Exercise:**\n", + "\n", "\n", "* Create N_PART particles at random positions.\n", "\n", @@ -349,11 +357,36 @@ "metadata": {}, "outputs": [], "source": [ + "# use LJ units: EPS=SIG=1\n", "LJ_EPS = 1.0\n", "LJ_SIG = 1.0\n", "LJ_CUT = 2.5 * LJ_SIG" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In a periodic system it is in general not straight forward to calculate all non-bonded interactions. Due to the periodicity and to speed up calculations usually a cut-off for short-ranged potentials like Lennard-Jones is applied. The potential can be shifted to zero at the cutoff value to ensure consistent thermodynamics using the shift='auto' option of [espressomd.interactions.LennardJonesInteraction](http://espressomd.org/html/doc/espressomd.html#module-espressomd.interactions).\n", + "Alternatively, for an isotropic system one can assume that the density is homogeneous behind the cutoff, which allows to calculate the so-called long-range corrections to the energy and pressure,\n", + "$$V_\\mathrm{lr} = 12 N \\rho \\int_{r_\\mathrm{c}}^\\infty 4 \\pi r^2 g(r) V(r) \\,\\mathrm{d}r,$$\n", + "where $N=$N_PART is the particle number and $r_\\mathrm{c}=$LJ_CUT is the cutoff distance.\n", + "Using that the radial distribution function $g(r)=1$ for $r>r_\\mathrm{cut}$ and considering only the dispersion (attractive) part of the LJ potential, one obtains\n", + "$$V_\\mathrm{lr} = -\\frac{8}{3}\\pi N \\rho \\varepsilon \\frac{\\sigma^6}{r_\\mathrm{c}^3}.$$\n", + "Similarly, a long-range contribution to the pressure can be derived [5].\n", + "\n", + "To avoid spurious self-interactions of particles with their periodic images one usually forces that the shortest box length is at least twice the cutoff distance:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "assert((BOX_L-2*SKIN > LJ_CUT).all())" + ] + }, { "cell_type": "markdown", "metadata": { @@ -361,9 +394,11 @@ "solution2_first": true }, "source": [ - "**Exercise:**\n", + "\n", "\n", "* Setup a Lennard-Jones interaction with $\\epsilon=$LJ_EPS and $\\sigma=$LJ_SIG which are cut and shifted to zero at $r_\\mathrm{c}=$LJ_CUT$\\times\\sigma$.\n", + " To allow for comparison with the fundamental work on MD simulations of LJ systems [6] we don't shift the potential to zero at the cutoff and add $V_\\mathrm{lr}$ later. However, keep in mind that this choice will influence the thermodynamics of your simulation system!\n", + "\n", "\n", " *Hint:* Use the [espressomd.system](http://espressomd.org/html/doc/espressomd.html#module-espressomd.system).non_bonded_inter handle object to instantiate a [espressomd.interactions.LennardJonesInteraction](http://espressomd.org/html/doc/espressomd.html#module-espressomd.interactions)." ] @@ -376,7 +411,7 @@ "source": [ "```python\n", "system.non_bonded_inter[0, 0].lennard_jones.set_params(\n", - " epsilon=LJ_EPS, sigma=LJ_SIG, cutoff=LJ_CUT, shift='auto')\n", + " epsilon=LJ_EPS, sigma=LJ_SIG, cutoff=LJ_CUT, shift=0)\n", "```" ] }, @@ -513,7 +548,8 @@ "outputs": [], "source": [ "# Parameters for the Langevin thermostat\n", - "TEMPERATURE = 0.728\n", + "# reduced temperature T* = k_B T / LJ_EPS\n", + "TEMPERATURE = 0.827 # value from Tab. 1 in [6]\n", "GAMMA=1.0" ] }, @@ -689,6 +725,7 @@ "ax2 = ax.twinx()\n", "ax2.plot(t,ekin/ (1.5 * N_PART), color='r', ls='--', lw=2,\n", " label='instantaneous temperature')\n", + "ax2.plot(TEMPERATURE)\n", "ax2.set_ylabel('instantaneous temperature', fontsize=20)\n", "ax2.legend(fontsize=16)\n", "ax2.tick_params(axis='both', labelsize=16)\n", @@ -732,14 +769,14 @@ "ax.plot(time-1,-acf, lw=2, color='b', ls='--')\n", "\n", "from scipy.optimize import curve_fit\n", - "fitfn = lambda x,a,t: -a*x/t\n", + "fitfn = lambda x,a,t: a*np.exp(-x/t)\n", "popt,pcov = curve_fit(fitfn, time-1,acf)\n", "tcor = popt[1]\n", "x = np.linspace(0,15)\n", "ax.plot(x,np.exp(-x/tcor), lw=2, color='r', ls='-.')\n", "ax.set_yscale('log')\n", - "ax.set_ylim(bottom=1e-3)\n", - "ax.set_xlim(0,15)\n", + "# ax.set_ylim(bottom=1e-3)\n", + "# ax.set_xlim(0,15)\n", "ax.tick_params(axis='both', labelsize=16)\n", "ax.set_ylabel('normalized ACF', fontsize=20)\n", "ax.set_xlabel('time', fontsize=20)\n", @@ -764,12 +801,13 @@ }, "source": [ "Plot the experimental radial distribution.\n", - "Empirical radial distribution functions have been determined for pure fluids [5], mixtures [6] and confined fluids [7]. We will compare our distribution $g(r)$ to the theoretical distribution $g(r^*, \\rho^*, T^*)$ of a pure fluid [5]." + "Empirical radial distribution functions have been determined for pure fluids [7], mixtures [8] and confined fluids [9]. We will compare our distribution $g(r)$ to the theoretical distribution $g(r^*, \\rho^*, T^*)$ of a pure fluid [7]." ] }, { "cell_type": "markdown", "metadata": { + "code_folding": [], "solution2": "hidden" }, "source": [ @@ -777,38 +815,65 @@ "fig3 = plt.figure(num=None, figsize=(10, 6),\n", " dpi=80, facecolor='w', edgecolor='k')\n", "fig3.set_tight_layout(False)\n", - "plt.plot(r, rdf, '-', color=\"#A60628\", linewidth=2, alpha=1, label='simulated')\n", + "plt.plot(r, rdf, '-', color=\"#A60628\", linewidth=3, alpha=1, label='simulated')\n", "\n", "# Comparison to literature\n", + "# ========================\n", + "\n", "# reduced units\n", "T_star = TEMPERATURE / LJ_EPS\n", "rho_star = DENSITY * LJ_SIG**3\n", + "\n", "# expression of the factors Pi from Equations 2-8 with coefficients qi from Table 1\n", - "P = lambda q1, q2, q3, q4, q5, q6, q7, q8, q9: q1 + q2 * np.exp(-q3 * T_star) + q4 * np.exp(-q5 * T_star) + q6 / rho_star + q7 / rho_star**2 + q8 * np.exp(-q3 * T_star) / rho_star**3 + q9 * np.exp(-q5 * T_star) / rho_star**4\n", + "# expression for a,g\n", + "P = lambda q1, q2, q3, q4, q5, q6, q7, q8, q9: \\\n", + " q1 + q2 * np.exp(-q3 * T_star) + q4 * np.exp(-q5 * T_star) + q6 / rho_star \\\n", + " + q7 / rho_star**2 + q8 * np.exp(-q3 * T_star) / rho_star**3 + q9 * np.exp(-q5 * T_star) / rho_star**4\n", "a = P(9.24792, -2.64281, 0.133386, -1.35932, 1.25338, 0.45602, -0.326422, 0.045708, -0.0287681)\n", "g = P(0.663161, -0.243089, 1.24749, -2.059, 0.04261, 1.65041, -0.343652, -0.037698, 0.008899)\n", - "P = lambda q1, q2, q3, q4, q5, q6, q7, q8: q1 + q2 * np.exp(-q3 * T_star) + q4 * rho_star + q5 * rho_star**2 + q6 * rho_star**3 + q7 * rho_star**4 + q8 * rho_star**5\n", + "\n", + "# expression for c,k\n", + "P = lambda q1, q2, q3, q4, q5, q6, q7, q8: \\\n", + " q1 + q2 * np.exp(-q3 * T_star) + q4 * rho_star + q5 * rho_star**2 + q6 * rho_star**3 \\\n", + " + q7 * rho_star**4 + q8 * rho_star**5\n", "c = P(-0.0677912, -1.39505, 0.512625, 36.9323, -36.8061, 21.7353, -7.76671, 1.36342)\n", "k = P(16.4821, -0.300612, 0.0937844, -61.744, 145.285, -168.087, 98.2181, -23.0583)\n", + "\n", + "# expression for b,h\n", "P = lambda q1, q2, q3: q1 + q2 * np.exp(-q3 * rho_star)\n", "b = P(-8.33289, 2.1714, 1.00063)\n", "h = P(0.0325039, -1.28792, 2.5487)\n", + "\n", + "# expression for d,l\n", "P = lambda q1, q2, q3, q4: q1 + q2 * np.exp(-q3 * rho_star) + q4 * rho_star\n", "d = P(-26.1615, 27.4846, 1.68124, 6.74296)\n", "l = P(-6.7293, -59.5002, 10.2466, -0.43596)\n", - "P = lambda q1, q2, q3, q4, q5, q6, q7, q8: (q1 + q2 * rho_star + q3 / T_star + q4 / T_star**2 + q5 / T_star**3) / (q6 + q7 * rho_star + q8 * rho_star**2)\n", + "\n", + "# expression for s\n", + "P = lambda q1, q2, q3, q4, q5, q6, q7, q8: \\\n", + " (q1 + q2 * rho_star + q3 / T_star + q4 / T_star**2 + q5 / T_star**3) \\\n", + " / (q6 + q7 * rho_star + q8 * rho_star**2)\n", "s = P(1.25225, -1.0179, 0.358564, -0.18533, 0.0482119, 1.27592, -1.78785, 0.634741)\n", - "P = lambda q1, q2, q3, q4, q5, q6: q1 + q2 * np.exp(-q3 * T_star) + q4 / T_star + q5 * rho_star + q6 * rho_star**2\n", + "\n", + "# expression for m\n", + "P = lambda q1, q2, q3, q4, q5, q6: \\\n", + " q1 + q2 * np.exp(-q3 * T_star) + q4 / T_star + q5 * rho_star + q6 * rho_star**2\n", "m = P(-5.668, -3.62671, 0.680654, 0.294481, 0.186395, -0.286954)\n", + "\n", + "# expression for n\n", "P = lambda q1, q2, q3: q1 + q2 * np.exp(-q3 * T_star)\n", "n = P(6.01325, 3.84098, 0.60793)\n", - "# theoretical curve\n", - "r_star = np.linspace(0, 3, 301)[1:]\n", + "\n", + "# plot fitted expression (=theoretical curve)\n", + "r_star = np.linspace(0, 4, 301)[1:]\n", "r_star_cutoff = 1.02 # slightly more than 1 to smooth out the discontinuity in the range [1.0, 1.02]\n", - "theo_rdf = 1 + 1 / r_star**2 * (np.exp(-(a * r_star + b)) * np.sin(c * r_star + d) + np.exp(-(g * r_star + h)) * np.cos(k * r_star + l))\n", - "theo_rdf[np.nonzero(r_star <= r_star_cutoff)] = s * np.exp(-(m * r_star + n)**4)[np.nonzero(r_star <= r_star_cutoff)]\n", - "plt.plot(r_star, theo_rdf, '-', color=\"blue\", linewidth=2, alpha=1, ls='--', label=r'theoretical')\n", "\n", + "theo_rdf = 1 + 1 / r_star**2 * (np.exp(-(a * r_star + b)) * np.sin(c * r_star + d) \\\n", + " + np.exp(-(g * r_star + h)) * np.cos(k * r_star + l))\n", + "theo_rdf[np.nonzero(r_star <= r_star_cutoff)] = \\\n", + " s * np.exp(-(m * r_star + n)**4)[np.nonzero(r_star <= r_star_cutoff)]\n", + "\n", + "plt.plot(r_star, theo_rdf, '-', color=\"blue\", linewidth=2, alpha=1, ls='--', label=r'theoretical')\n", "\n", "plt.xlabel('$r [\\sigma]$', fontsize=20)\n", "plt.ylabel('$g(r)$', fontsize=20)\n", @@ -845,7 +910,7 @@ "r_max = system.box_l[0] / 2.0\n", "# Initialize RDF accumulator\n", "rdf_obs = RDF(ids1=system.part[:].id, min_r=r_min, max_r=r_max, n_r_bins=r_bins)\n", - "rdf_acc = MeanVarianceCalculator(obs=rdf_obs, delta_N=sampling_interval)\n", + "rdf_acc = MeanVarianceCalculator(obs=rdf_obs, delta_N=SAMPLING_INTERVAL)\n", "# Accumulate results automatically during the integration\n", "system.auto_update_accumulators.add(rdf_acc)\n", "\n", @@ -872,59 +937,6 @@ "r = rdf_obs.bin_centers()" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Plot the experimental radial distribution." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig1 = plt.figure(num=None, figsize=(10, 6), dpi=80, facecolor='w', edgecolor='k')\n", - "fig1.set_tight_layout(False)\n", - "plt.plot(r, rdf, '-', color=\"#A60628\", linewidth=2, alpha=1)\n", - "plt.xlabel('$r [\\sigma]$', fontsize=20)\n", - "plt.ylabel('$g(r)$', fontsize=20)\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now plot the instantaneous temperature." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig2 = plt.figure(num=None, figsize=(10, 6), dpi=80, facecolor='w', edgecolor='k')\n", - "fig2.set_tight_layout(False)\n", - "plt.plot(time, instantaneous_temperature, '-', color=\"red\", linewidth=2,\n", - " alpha=0.5, label='Instantaneous Temperature')\n", - "plt.plot([min(time), max(time)], [TEMPERATURE] * 2, '-', color=\"#348ABD\",\n", - " linewidth=2, alpha=1, label='Set Temperature')\n", - "plt.xlabel(r'Time [$\\delta t$]', fontsize=20)\n", - "plt.ylabel(r'$k_B$ Temperature [$k_B T$]', fontsize=20)\n", - "plt.legend(fontsize=16, loc=0)\n", - "plt.show()" - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -1003,10 +1015,19 @@ "[2] HJ Limbach, A. Arnold, and B. Mann. ESPResSo: An extensible simulation package for research on soft matter systems. *Computer Physics Communications*, 174(9):704–727, 2006. \n", "[3] A. Arnold, O. Lenz, S. Kesselheim, R. Weeber, F. Fahrenberger, D. Rohm, P. Kosovan, and C. Holm. ESPResSo 3.1 — molecular dynamics software for coarse-grained models. In M. Griebel and M. A. Schweitzer, editors, Meshfree Methods for Partial Differential Equations VI, volume 89 of Lecture Notes in Computational Science and Engineering, pages 1–23. Springer Berlin Heidelberg, 2013. \n", "[4] A. Arnold, BA Mann, HJ Limbach, and C. Holm. ESPResSo–An Extensible Simulation Package for Research on Soft Matter Systems. *Forschung und wissenschaftliches Rechnen*, 63:43–59, 2003. \n", - "[5] Morsali, Goharshadi, Mansoori, Abbaspour. An accurate expression for radial distribution function of the Lennard-Jones fluid. *Chemical Physics*, 310(1–3):11–15, 2005. DOI:10.1016/j.chemphys.2004.09.027 \n", - "[6] Matteoli. A simple expression for radial distribution functions of pure fluids and mixtures. *The Journal of Chemical Physics*, 103(11):4672, 1995. DOI:10.1063/1.470654 \n", - "[7] Abbaspour, Akbarzadeha, Abroodia. A new and accurate expression for the radial distribution function of confined Lennard-Jones fluid in carbon nanotubes. *RSC Advances*, 5(116): 95781–95787, 2015. DOI:10.1039/C5RA16151G " + "[5] W. P. Allen & D. J. Tildesley. Computer Simulation of Liquids. Oxford University Press, 2017. \n", + "[6] L. Verlet, “Computer ‘Experiments’ on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules, *Phys. Rev.*, 159(1):98–103, 1967. DOI:10.1103/PhysRev.159.98 \n", + "[6] Morsali, Goharshadi, Mansoori, Abbaspour. An accurate expression for radial distribution function of the Lennard-Jones fluid. *Chemical Physics*, 310(1–3):11–15, 2005. DOI:10.1016/j.chemphys.2004.09.027 \n", + "[7] Matteoli. A simple expression for radial distribution functions of pure fluids and mixtures. *The Journal of Chemical Physics*, 103(11):4672, 1995. DOI:10.1063/1.470654 \n", + "[8] Abbaspour, Akbarzadeha, Abroodia. A new and accurate expression for the radial distribution function of confined Lennard-Jones fluid in carbon nanotubes. *RSC Advances*, 5(116): 95781–95787, 2015. DOI:10.1039/C5RA16151G " ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/doc/tutorials/01-lennard_jones/NotesForTutor.md b/doc/tutorials/01-lennard_jones/NotesForTutor.md index 6dc21cf2cea..7c0ec15693b 100644 --- a/doc/tutorials/01-lennard_jones/NotesForTutor.md +++ b/doc/tutorials/01-lennard_jones/NotesForTutor.md @@ -5,8 +5,12 @@ After the tutorial, students should be able to * explain - * role of Lennard-Jones potential + * role of Lennard-Jones potential: + * noble gases + * first MD simulations by Verlet et al. (liquid argon) * role of radial distribution function + * relation to structure factor + * LJ agrees well with experiments for noble gases! * connection between mean squared displacement and diffusion coefficient * how to use an auto correlation function to estimate correlation times and how that affects error estimation @@ -31,9 +35,11 @@ After the tutorial, students should be able to * Mention efficiency (loops and lists) vs. numpy arrays at the cost of some readability. * Non-bonded interactions are defined between pairs of particle types. Mention that type=0 is implied when creating particles without specifying one. -* Explain Lennard-Jones cutoff and shift. +* Explain Lennard-Jones cutoff and shift - we don't shift here to make energies comparable to historic simulations (Verlet at al.). +* Explain reduced units - here: LJ units! * Explain steepest descent: In general overlap removal can be done in multiple ways (warmup with force capping, small time step, etc.). Additionally the steepest descent algorithm can get trapped in local minima and the convergence criterion is system-dependent. -* Mention that accumulators are updated auomatically. \ No newline at end of file +* Mention that accumulators are updated auomatically. +* Statistical independent measurements for t > correlation time. From c9e2817c09ec447a387ec79141edb9abb2cbb121 Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Wed, 23 Sep 2020 18:11:45 +0200 Subject: [PATCH 26/46] Remove MSD and redundant parts --- .../01-lennard_jones/01-lennard_jones.ipynb | 180 +++++------------- .../01-lennard_jones/NotesForTutor.md | 1 - 2 files changed, 43 insertions(+), 138 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index e6e22959362..148cacaef7e 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -550,7 +550,7 @@ "# Parameters for the Langevin thermostat\n", "# reduced temperature T* = k_B T / LJ_EPS\n", "TEMPERATURE = 0.827 # value from Tab. 1 in [6]\n", - "GAMMA=1.0" + "GAMMA = 1.0" ] }, { @@ -617,8 +617,8 @@ "outputs": [], "source": [ "# Integration parameters\n", - "SAMPLING_INTERVAL = 100\n", - "SAMPLING_ITERATIONS = 200\n", + "SAMPLING_INTERVAL = 20\n", + "SAMPLING_ITERATIONS = 100\n", "\n", "# Parameters for the radial distribution function\n", "R_BINS = 100\n", @@ -635,7 +635,8 @@ "source": [ "**Exercise:**\n", "\n", - "* Set up the main integration loop\n", + "**Part I:**\n", + "* Set up the main integration loop.\n", " Take SAMPLING_ITERATIONS measurements every SAMPLING_INTERVAL integration steps.\n", "* Monitor the potential and kinetic energies using the analysis method [system.analysis](http://espressomd.org/html/doc/espressomd.html#module-espressomd.analyze).energy().\n", " The total energy time series should be stored in a numpy.array etotal\n", @@ -643,31 +644,26 @@ " of the instantaneous temperature $T = 3/2 \\times E_\\mathrm{kin}$/N_PART.\n", " Here $T$ refers to the measured temperature obtained from kinetic energy and the number\n", " of degrees of freedom in the system. It should fluctuate around the preset temperature of the thermostat.\n", - " Also store the corresponding simulation time in time.\n", + " Use the corresponding simulation time stored in time.\n", "* Plot the instantaneous temperature and total energy versus time.\n", "* Estimate the correlation time of the system from the autocorrelation function of the total energy.\n", + "* Make sure to use only uncorrelated samples to measure the total energy per particle including the long-range correction.\n", + " Provide a confidence interval and compare your results to the value given in Ref. [6].\n", + "\n", + "**Part II:**\n", "* Extend the integration loop to include the calculation of the radial distribution function which describes how the density varies as a function of distance from a tagged particle. \n", " The radial distribution function (RDF) is averaged over several measurements to reduce noise.\n", " To this end add an observable rdf_obs for the radial distirbution function.\n", " Use min_r=R_MIN, max_r=R_MAX, n_r_bins=R_BINS and create a corresponding [espressomd.accumulators.MeanVarianceCalculator](http://espressomd.org/html/doc/espressomd.html#espressomd.accumulators.MeanVarianceCalculator)\n", " that gets updated every SAMPLING_INTERVAL timesteps.\n", " Plot the RDF using the mean of the accumulator.\n", - " *Hint*: The radial coordinates can be obtained using the bin_centers() method of the espressomd.observables.RDF observable.\n", - "\n", - "\n", - "\n", - "The mean square displacement of particle $i$ is given by:\n", - "\n", - "\\begin{equation}\n", - "\\mathrm{msd}_i(t) =\\langle (\\vec{x}_i(t_0+t) -\\vec{x}_i(t_0))^2\\rangle,\n", - "\\end{equation}\n", - "\n", - "and can be calculated using \"observables and correlators\". An observable is an object which takes a measurement on the system. It can depend on parameters specified when the observable is instanced, such as the ids of the particles to be considered." + " *Hint*: The radial coordinates can be obtained using the bin_centers() method of the espressomd.observables.RDF observable." ] }, { "cell_type": "markdown", "metadata": { + "code_folding": [], "solution2": "hidden" }, "source": [ @@ -686,18 +682,20 @@ "system.auto_update_accumulators.add(rdf_acc)\n", "\n", "for i in range(SAMPLING_ITERATIONS):\n", + " \n", + " time[i] = system.time\n", + " print(i,'/',SAMPLING_ITERATIONS,end=\"\\r\")\n", + " \n", " system.integrator.run(SAMPLING_INTERVAL)\n", " \n", " # Measure energies\n", " energies = system.analysis.energy()\n", " etotal[i] = energies['total']\n", " ekin[i] = energies['kinetic']\n", - " time[i] = system.time\n", " \n", "# Finalize the accumulator and obtain the results\n", "rdf = rdf_acc.get_mean()\n", "r = rdf_obs.bin_centers()\n", - "\n", "```" ] }, @@ -725,10 +723,12 @@ "ax2 = ax.twinx()\n", "ax2.plot(t,ekin/ (1.5 * N_PART), color='r', ls='--', lw=2,\n", " label='instantaneous temperature')\n", - "ax2.plot(TEMPERATURE)\n", + "ax2.axhline(TEMPERATURE, ls='-.', lw=3, color='k')\n", "ax2.set_ylabel('instantaneous temperature', fontsize=20)\n", "ax2.legend(fontsize=16)\n", "ax2.tick_params(axis='both', labelsize=16)\n", + "print(\"average simulation temperature: {:1.3f}\".format(np.round(ekin.mean()/ (1.5 * N_PART),3)),\n", + " \"\\t(which differs by {:1.3f}% from TEMPERATURE)\".format(np.round(((ekin.mean()/ (1.5 * N_PART) - TEMPERATURE) / TEMPERATURE)*100,3)))\n", "```" ] }, @@ -765,21 +765,22 @@ "fig2.set_tight_layout(False)\n", "\n", "acf = autocor(etotal)\n", - "ax.plot(time-1,acf, lw=2, color='b')\n", - "ax.plot(time-1,-acf, lw=2, color='b', ls='--')\n", + "ax.plot(time/(system.time_step*SAMPLING_INTERVAL), acf, lw=2, color='b')\n", + "# ax.plot(time,-acf, lw=2, color='b', ls='--')\n", "\n", "from scipy.optimize import curve_fit\n", "fitfn = lambda x,a,t: a*np.exp(-x/t)\n", - "popt,pcov = curve_fit(fitfn, time-1,acf)\n", + "popt,pcov = curve_fit(fitfn, time/(system.time_step*SAMPLING_INTERVAL), acf)\n", "tcor = popt[1]\n", - "x = np.linspace(0,15)\n", + "x = np.linspace(0,5)\n", "ax.plot(x,np.exp(-x/tcor), lw=2, color='r', ls='-.')\n", - "ax.set_yscale('log')\n", + "# ax.set_yscale('log')\n", "# ax.set_ylim(bottom=1e-3)\n", - "# ax.set_xlim(0,15)\n", + "ax.set_xlim(0,5)\n", "ax.tick_params(axis='both', labelsize=16)\n", "ax.set_ylabel('normalized ACF', fontsize=20)\n", - "ax.set_xlabel('time', fontsize=20)\n", + "ax.set_xlabel(r'time / $\\tau_\\mathrm{S}$', fontsize=20)\n", + "ax.axhline(1./np.exp(1))\n", "\n", "print (\"Estimated correlation time:\", np.round(tcor,2))\n", "```" @@ -794,13 +795,27 @@ "From the roughly exponential decay of the correlation function it dawns that only every second sample should be used to calculate statistically independent averages." ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "solution2": "hidden" + }, + "outputs": [], + "source": [ + "corrected_energy = (etotal-ekin)[::2]/N_PART - 2/3 * np.pi * DENSITY * (4* LJ_EPS * LJ_SIG**6) / LJ_CUT**3\n", + "print(\"The obtained potential energy {:1.3f} +/- {:1.3f} agrees excellently with\\\n", + "the value -5.38 given in [6]\".format(*np.round(np.array([corrected_energy.mean(),\n", + " corrected_energy.std()/np.sqrt(len(corrected_energy)-1)]),3)))" + ] + }, { "cell_type": "markdown", "metadata": { "solution2": "hidden" }, "source": [ - "Plot the experimental radial distribution.\n", + "We now plot the experimental radial distribution.\n", "Empirical radial distribution functions have been determined for pure fluids [7], mixtures [8] and confined fluids [9]. We will compare our distribution $g(r)$ to the theoretical distribution $g(r^*, \\rho^*, T^*)$ of a pure fluid [7]." ] }, @@ -882,113 +897,11 @@ "```" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "code_folding": [] - }, - "outputs": [], - "source": [ - "# Original solution\n", - "\n", - "from espressomd.observables import ParticlePositions, RDF\n", - "from espressomd.accumulators import Correlator, MeanVarianceCalculator\n", - "# Pass the ids of the particles to be tracked to the observable.\n", - "part_pos = ParticlePositions(ids=range(N_PART))\n", - "# Initialize MSD correlator\n", - "msd_corr = Correlator(obs1=part_pos,\n", - " tau_lin=10, delta_N=10,\n", - " tau_max=1000 * TIME_STEP,\n", - " corr_operation=\"square_distance_componentwise\")\n", - "# Calculate results automatically during the integration\n", - "system.auto_update_accumulators.add(msd_corr)\n", - "\n", - "# Set parameters for the radial distribution function\n", - "r_bins = 100\n", - "r_min = 0.0\n", - "r_max = system.box_l[0] / 2.0\n", - "# Initialize RDF accumulator\n", - "rdf_obs = RDF(ids1=system.part[:].id, min_r=r_min, max_r=r_max, n_r_bins=r_bins)\n", - "rdf_acc = MeanVarianceCalculator(obs=rdf_obs, delta_N=SAMPLING_INTERVAL)\n", - "# Accumulate results automatically during the integration\n", - "system.auto_update_accumulators.add(rdf_acc)\n", - "\n", - "# Take measurements\n", - "time = np.zeros(sampling_iterations)\n", - "instantaneous_temperature = np.zeros(sampling_iterations)\n", - "etotal = np.zeros(sampling_iterations)\n", - "\n", - "for i in range(sampling_iterations):\n", - " system.integrator.run(sampling_interval)\n", - " # Measure energies\n", - " energies = system.analysis.energy()\n", - " kinetic_temperature = energies['kinetic'] / (1.5 * N_PART)\n", - " etotal[i] = energies['total']\n", - " time[i] = system.time\n", - " instantaneous_temperature[i] = kinetic_temperature\n", - "\n", - "# Finalize the correlator and obtain the results\n", - "msd_corr.finalize()\n", - "msd = msd_corr.result()\n", - "\n", - "# Finalize the accumulator and obtain the results\n", - "rdf = rdf_acc.get_mean()\n", - "r = rdf_obs.bin_centers()" - ] - }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The correlator output is stored in the array `msd` and has the following shape:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print(msd.shape)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The first column of this array contains the lag time in units of the time step.\n", - "The second column contains the number of values used to perform the averaging of the correlation.\n", - "The next three columns contain the x, y and z mean squared displacement of the msd of the first particle.\n", - "The next three columns then contain the x, y, z mean squared displacement of the next particle..." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig3 = plt.figure(num=None, figsize=(10, 6), dpi=80, facecolor='w', edgecolor='k')\n", - "fig3.set_tight_layout(False)\n", - "lag_time = msd_corr.lag_times()\n", - "for i in range(0, N_PART, 30):\n", - " msd_particle_i = msd[:, i, 0] + msd[:, i, 1] + msd[:, i, 2]\n", - " plt.plot(lag_time, msd_particle_i,\n", - " 'o-', linewidth=2, label=\"particle id =\" + str(i))\n", - "plt.xlabel(r'Lag time $\\tau$ [$\\delta t$]', fontsize=20)\n", - "plt.ylabel(r'Mean squared displacement [$\\sigma^2$]', fontsize=20)\n", - "plt.xscale('log')\n", - "plt.yscale('log')\n", - "plt.legend()\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Exercises" + "## Futher Exercises" ] }, { @@ -1021,13 +934,6 @@ "[7] Matteoli. A simple expression for radial distribution functions of pure fluids and mixtures. *The Journal of Chemical Physics*, 103(11):4672, 1995. DOI:10.1063/1.470654 \n", "[8] Abbaspour, Akbarzadeha, Abroodia. A new and accurate expression for the radial distribution function of confined Lennard-Jones fluid in carbon nanotubes. *RSC Advances*, 5(116): 95781–95787, 2015. DOI:10.1039/C5RA16151G " ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { diff --git a/doc/tutorials/01-lennard_jones/NotesForTutor.md b/doc/tutorials/01-lennard_jones/NotesForTutor.md index 7c0ec15693b..44e81566367 100644 --- a/doc/tutorials/01-lennard_jones/NotesForTutor.md +++ b/doc/tutorials/01-lennard_jones/NotesForTutor.md @@ -11,7 +11,6 @@ After the tutorial, students should be able to * role of radial distribution function * relation to structure factor * LJ agrees well with experiments for noble gases! - * connection between mean squared displacement and diffusion coefficient * how to use an auto correlation function to estimate correlation times and how that affects error estimation From 0d7726a1a3c1aa217ea154922b84b4ab05be20e2 Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Wed, 23 Sep 2020 20:12:20 +0200 Subject: [PATCH 27/46] Remove obsolete test for global standard error --- testsuite/scripts/tutorials/test_01-lennard_jones.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/testsuite/scripts/tutorials/test_01-lennard_jones.py b/testsuite/scripts/tutorials/test_01-lennard_jones.py index ad7e060efd9..ef70bbe70ee 100644 --- a/testsuite/scripts/tutorials/test_01-lennard_jones.py +++ b/testsuite/scripts/tutorials/test_01-lennard_jones.py @@ -28,9 +28,6 @@ class Tutorial(ut.TestCase): system = tutorial.system - def test_global_variables(self): - self.assertLess(tutorial.standard_error_total_energy, 2.5) - def test_rdf_curve(self): self.assertGreater( np.corrcoef(tutorial.rdf, tutorial.theo_rdf)[1, 0], 0.985) From e259d6e3f2a035a51effa16aab3c31de534375ff Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Wed, 23 Sep 2020 20:28:20 +0200 Subject: [PATCH 28/46] Remove obsolete statement about the core in C --- doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index 148cacaef7e..fc4393fe143 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -99,8 +99,7 @@ "\n", "What is ESPResSo? It is an extensible, efficient Molecular Dynamics package specially powerful on simulating charged systems. In depth information about the package can be found in the relevant sources [1,4,2,3].\n", "\n", - "ESPResSo consists of two components. The simulation engine is written in C and C++ for the sake of computational efficiency. The steering or control\n", - "level is interfaced to the kernel via an interpreter of the Python scripting languages.\n", + "ESPResSo consists of two components. The simulation engine is written in C++ for the sake of computational efficiency. The steering or control level is interfaced to the kernel via an interpreter of the Python scripting languages.\n", "\n", "The kernel performs all computationally demanding tasks. Before all, integration of Newton's equations of motion, including calculation of energies and forces. It also takes care of internal organization of data, storing the data about particles, communication between different processors or cells of the cell-system.\n", "\n", From 97a4602d51e242cc33e5fdf605fe21d33814c0d2 Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Thu, 24 Sep 2020 10:25:29 +0200 Subject: [PATCH 29/46] Improve testing --- .../01-lennard_jones/01-lennard_jones.ipynb | 27 +++++++++---------- .../tutorials/test_01-lennard_jones.py | 10 ++++--- 2 files changed, 20 insertions(+), 17 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index fc4393fe143..ae67b052883 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -795,17 +795,18 @@ ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": { "solution2": "hidden" }, - "outputs": [], "source": [ + "```python\n", "corrected_energy = (etotal-ekin)[::2]/N_PART - 2/3 * np.pi * DENSITY * (4* LJ_EPS * LJ_SIG**6) / LJ_CUT**3\n", - "print(\"The obtained potential energy {:1.3f} +/- {:1.3f} agrees excellently with\\\n", - "the value -5.38 given in [6]\".format(*np.round(np.array([corrected_energy.mean(),\n", - " corrected_energy.std()/np.sqrt(len(corrected_energy)-1)]),3)))" + "sim_energy = corrected_energy.mean()\n", + "dsim_energy = corrected_energy.std()/np.sqrt(len(corrected_energy)-1)\n", + "print(\"The obtained potential energy per particle is {:1.3f} +/- {:1.3f}, which agrees well with\\\n", + "the value -5.38 given in [6]\".format(*np.round(np.array([sim_energy,dsim_energy]),3)))\n", + "```" ] }, { @@ -879,20 +880,18 @@ "n = P(6.01325, 3.84098, 0.60793)\n", "\n", "# plot fitted expression (=theoretical curve)\n", - "r_star = np.linspace(0, 4, 301)[1:]\n", - "r_star_cutoff = 1.02 # slightly more than 1 to smooth out the discontinuity in the range [1.0, 1.02]\n", + "theo_rdf_cutoff = 1.02 # slightly more than 1 to smooth out the discontinuity in the range [1.0, 1.02]\n", "\n", - "theo_rdf = 1 + 1 / r_star**2 * (np.exp(-(a * r_star + b)) * np.sin(c * r_star + d) \\\n", - " + np.exp(-(g * r_star + h)) * np.cos(k * r_star + l))\n", - "theo_rdf[np.nonzero(r_star <= r_star_cutoff)] = \\\n", - " s * np.exp(-(m * r_star + n)**4)[np.nonzero(r_star <= r_star_cutoff)]\n", + "theo_rdf = 1 + 1 / r**2 * (np.exp(-(a * r + b)) * np.sin(c * r + d) \\\n", + " + np.exp(-(g * r + h)) * np.cos(k * r + l))\n", + "theo_rdf[np.nonzero(r <= theo_rdf_cutoff)] = \\\n", + " s * np.exp(-(m * r + n)**4)[np.nonzero(r <= theo_rdf_cutoff)]\n", "\n", - "plt.plot(r_star, theo_rdf, '-', color=\"blue\", linewidth=2, alpha=1, ls='--', label=r'theoretical')\n", + "plt.plot(r, theo_rdf, '-', color=\"blue\", linewidth=2, alpha=1, ls='--', label=r'theoretical')\n", "\n", "plt.xlabel('$r [\\sigma]$', fontsize=20)\n", "plt.ylabel('$g(r)$', fontsize=20)\n", "plt.gca().tick_params(axis='both', labelsize=16)\n", - "plt.show()\n", "```" ] }, diff --git a/testsuite/scripts/tutorials/test_01-lennard_jones.py b/testsuite/scripts/tutorials/test_01-lennard_jones.py index ef70bbe70ee..48d67036d73 100644 --- a/testsuite/scripts/tutorials/test_01-lennard_jones.py +++ b/testsuite/scripts/tutorials/test_01-lennard_jones.py @@ -28,9 +28,13 @@ class Tutorial(ut.TestCase): system = tutorial.system - def test_rdf_curve(self): - self.assertGreater( - np.corrcoef(tutorial.rdf, tutorial.theo_rdf)[1, 0], 0.985) + def test_rdf(self): + rms = np.sqrt(((tutorial.rdf - tutorial.theo_rdf)**2).mean()) / len(tutorial.r) + self.assertLess(rms, 5e-4) + + def test_potential_energy(self): + # Test that the potential energy/particle agrees with the value from Verlet, Phys. Rev. 1967 + self.assertLess(np.abs(tutorial.sim_energy + 5.38)/5.38,0.01) if __name__ == "__main__": From 32836a80e901337ac4ed0c1830def8a644a34d5b Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Thu, 24 Sep 2020 12:09:44 +0200 Subject: [PATCH 30/46] Refine tests --- .../01-lennard_jones/01-lennard_jones.ipynb | 2 +- .../scripts/tutorials/test_01-lennard_jones.py | 15 +++++++-------- 2 files changed, 8 insertions(+), 9 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index ae67b052883..e086c9c10f7 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -617,7 +617,7 @@ "source": [ "# Integration parameters\n", "SAMPLING_INTERVAL = 20\n", - "SAMPLING_ITERATIONS = 100\n", + "SAMPLING_ITERATIONS = 200\n", "\n", "# Parameters for the radial distribution function\n", "R_BINS = 100\n", diff --git a/testsuite/scripts/tutorials/test_01-lennard_jones.py b/testsuite/scripts/tutorials/test_01-lennard_jones.py index 48d67036d73..56ad69a1b50 100644 --- a/testsuite/scripts/tutorials/test_01-lennard_jones.py +++ b/testsuite/scripts/tutorials/test_01-lennard_jones.py @@ -20,22 +20,21 @@ import numpy as np tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import( - "@TUTORIALS_DIR@/01-lennard_jones/01-lennard_jones.py", - substitutions=lambda code: code.replace('r_star = ', 'r_star = r #', 1)) + "@TUTORIALS_DIR@/01-lennard_jones/01-lennard_jones.py") @skipIfMissingFeatures class Tutorial(ut.TestCase): - system = tutorial.system def test_rdf(self): - rms = np.sqrt(((tutorial.rdf - tutorial.theo_rdf)**2).mean()) / len(tutorial.r) - self.assertLess(rms, 5e-4) + np.testing.assert_allclose(tutorial.rdf, tutorial.theo_rdf, + rtol=0.0, atol=0.1) def test_potential_energy(self): - # Test that the potential energy/particle agrees with the value from Verlet, Phys. Rev. 1967 - self.assertLess(np.abs(tutorial.sim_energy + 5.38)/5.38,0.01) - + # Test that the potential energy/particle agrees with + # the value from Verlet, Phys. Rev. 1967 + ref_energy = -5.38 + self.assertAlmostEqual(tutorial.sim_energy, ref_energy, 1) if __name__ == "__main__": ut.main() From 6f172c09aedf207bc5f7966f62a4d4e47af9d29c Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Thu, 24 Sep 2020 12:27:28 +0200 Subject: [PATCH 31/46] formatting change --- testsuite/scripts/tutorials/test_01-lennard_jones.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/testsuite/scripts/tutorials/test_01-lennard_jones.py b/testsuite/scripts/tutorials/test_01-lennard_jones.py index 56ad69a1b50..62e99085dc7 100644 --- a/testsuite/scripts/tutorials/test_01-lennard_jones.py +++ b/testsuite/scripts/tutorials/test_01-lennard_jones.py @@ -28,7 +28,7 @@ class Tutorial(ut.TestCase): def test_rdf(self): np.testing.assert_allclose(tutorial.rdf, tutorial.theo_rdf, - rtol=0.0, atol=0.1) + rtol=0.0, atol=0.1) def test_potential_energy(self): # Test that the potential energy/particle agrees with @@ -36,5 +36,6 @@ def test_potential_energy(self): ref_energy = -5.38 self.assertAlmostEqual(tutorial.sim_energy, ref_energy, 1) + if __name__ == "__main__": ut.main() From e3504338957dbe4c5f446435662851af9cae2d11 Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Mon, 28 Sep 2020 11:12:22 +0200 Subject: [PATCH 32/46] Adapt naming convention --- .../01-lennard_jones/01-lennard_jones.ipynb | 46 +++++++++---------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index e086c9c10f7..150b1bc028b 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -485,20 +485,20 @@ "\n", "# Initialize integrator to obtain initial forces\n", "system.integrator.run(0)\n", - "maxforce = np.max(np.linalg.norm(system.part[:].f, axis=1))\n", + "max_force = np.max(np.linalg.norm(system.part[:].f, axis=1))\n", "energy = system.analysis.energy()['total']\n", "\n", "i = 0\n", "while i < MAX_STEPS//EMSTEP:\n", - " prev_maxforce = maxforce\n", + " prev_max_force = max_force\n", " prev_energy = energy\n", " system.integrator.run(EMSTEP)\n", - " maxforce = np.max(np.linalg.norm(system.part[:].f, axis=1))\n", - " relforce = np.abs((maxforce-prev_maxforce)/prev_maxforce)\n", + " max_force = np.max(np.linalg.norm(system.part[:].f, axis=1))\n", + " rel_force = np.abs((max_force-prev_max_force)/prev_max_force)\n", " energy = system.analysis.energy()['total']\n", - " relener = np.abs((energy-prev_energy)/prev_energy)\n", - " print(\"minimization step: {:4.0f}\\tmax. rel. force change:{:+3.3e}\\trel. energy change:{:+3.3e}\".format((i+1)*EMSTEP,relforce, relener))\n", - " if relforce < F_TOL or relener < E_TOL:\n", + " rel_ener = np.abs((energy-prev_energy)/prev_energy)\n", + " print(\"minimization step: {:4.0f}\\tmax. rel. force change:{:+3.3e}\\trel. energy change:{:+3.3e}\".format((i+1)*EMSTEP,rel_force, rel_ener))\n", + " if rel_force < F_TOL or rel_ener < E_TOL:\n", " break\n", " i += 1\n", "```" @@ -527,7 +527,7 @@ "source": [ "### Choosing the thermodynamic ensemble, thermostat\n", "\n", - "Simulations can be carried out in different thermodynamic ensembles such as NVE (particle __N__umber, __V__olume, __E__nergy), NVT (particle __N__umber, __V__olume, __T__emperature) or NPT-isotropic (particle __N__umber, __P__ressure, __T__emperature).\n", + "Simulations can be carried out in different thermodynamic ensembles such as NVE (particle __N__umber, __V__olume, __E__nergy), NVT (particle __N__umber, __V__olume, __T__emperature) or NpT-isotropic (particle __N__umber, __p__ressure, __T__emperature).\n", "\n", "The NVE ensemble is simulated without a thermostat. A previously enabled thermostat can be switched off as follows:\n", "\n", @@ -535,7 +535,7 @@ "system.thermostat.turn_off()\n", "```\n", "\n", - "The NVT and NPT ensembles require a thermostat. In this tutorial, we use the Langevin thermostat.\n", + "The NVT and NpT ensembles require a thermostat. In this tutorial, we use the Langevin thermostat.\n", "\n", "In ESPResSo, the thermostat is set as follows:" ] @@ -638,8 +638,8 @@ "* Set up the main integration loop.\n", " Take SAMPLING_ITERATIONS measurements every SAMPLING_INTERVAL integration steps.\n", "* Monitor the potential and kinetic energies using the analysis method [system.analysis](http://espressomd.org/html/doc/espressomd.html#module-espressomd.analyze).energy().\n", - " The total energy time series should be stored in a numpy.array etotal\n", - " and the kinetic energy in ekin. Plot the timeseries of the total energy and\n", + " The total energy time series should be stored in a numpy.array e_total\n", + " and the kinetic energy in e_kin. Plot the timeseries of the total energy and\n", " of the instantaneous temperature $T = 3/2 \\times E_\\mathrm{kin}$/N_PART.\n", " Here $T$ refers to the measured temperature obtained from kinetic energy and the number\n", " of degrees of freedom in the system. It should fluctuate around the preset temperature of the thermostat.\n", @@ -667,11 +667,11 @@ }, "source": [ "```python\n", - "etotal = np.zeros(SAMPLING_ITERATIONS)\n", - "ekin = np.zeros(SAMPLING_ITERATIONS)\n", + "e_total = np.zeros(SAMPLING_ITERATIONS)\n", + "e_kin = np.zeros(SAMPLING_ITERATIONS)\n", "time = np.zeros(SAMPLING_ITERATIONS)\n", "\n", - "from espressomd.observables import ParticlePositions, RDF\n", + "from espressomd.observables import RDF\n", "from espressomd.accumulators import Correlator, MeanVarianceCalculator\n", "\n", "# Initialize RDF accumulator\n", @@ -689,8 +689,8 @@ " \n", " # Measure energies\n", " energies = system.analysis.energy()\n", - " etotal[i] = energies['total']\n", - " ekin[i] = energies['kinetic']\n", + " e_total[i] = energies['total']\n", + " e_kin[i] = energies['kinetic']\n", " \n", "# Finalize the accumulator and obtain the results\n", "rdf = rdf_acc.get_mean()\n", @@ -714,20 +714,20 @@ "fig1.set_tight_layout(False)\n", "\n", "t = np.arange(0,SAMPLING_ITERATIONS*SAMPLING_INTERVAL,SAMPLING_INTERVAL)\n", - "ax.plot(t,etotal, lw=2, label='etotal')\n", + "ax.plot(t,e_total, lw=2, label='e_total')\n", "ax.set_xlabel('time step', fontsize=20)\n", "ax.set_ylabel('total energy', fontsize=20)\n", "ax.legend(fontsize=16)\n", "ax.tick_params(axis='both', labelsize=16)\n", "ax2 = ax.twinx()\n", - "ax2.plot(t,ekin/ (1.5 * N_PART), color='r', ls='--', lw=2,\n", + "ax2.plot(t,e_kin/ (1.5 * N_PART), color='r', ls='--', lw=2,\n", " label='instantaneous temperature')\n", "ax2.axhline(TEMPERATURE, ls='-.', lw=3, color='k')\n", "ax2.set_ylabel('instantaneous temperature', fontsize=20)\n", "ax2.legend(fontsize=16)\n", "ax2.tick_params(axis='both', labelsize=16)\n", - "print(\"average simulation temperature: {:1.3f}\".format(np.round(ekin.mean()/ (1.5 * N_PART),3)),\n", - " \"\\t(which differs by {:1.3f}% from TEMPERATURE)\".format(np.round(((ekin.mean()/ (1.5 * N_PART) - TEMPERATURE) / TEMPERATURE)*100,3)))\n", + "print(\"average simulation temperature: {:1.3f}\".format(np.round(e_kin.mean()/ (1.5 * N_PART),3)),\n", + " \"\\t(which differs by {:1.3f}% from TEMPERATURE)\".format(np.round(((e_kin.mean()/ (1.5 * N_PART) - TEMPERATURE) / TEMPERATURE)*100,3)))\n", "```" ] }, @@ -763,7 +763,7 @@ " dpi=80, facecolor='w', edgecolor='k')\n", "fig2.set_tight_layout(False)\n", "\n", - "acf = autocor(etotal)\n", + "acf = autocor(e_total)\n", "ax.plot(time/(system.time_step*SAMPLING_INTERVAL), acf, lw=2, color='b')\n", "# ax.plot(time,-acf, lw=2, color='b', ls='--')\n", "\n", @@ -801,7 +801,7 @@ }, "source": [ "```python\n", - "corrected_energy = (etotal-ekin)[::2]/N_PART - 2/3 * np.pi * DENSITY * (4* LJ_EPS * LJ_SIG**6) / LJ_CUT**3\n", + "corrected_energy = (e_total-e_kin)[::2]/N_PART - 2/3 * np.pi * DENSITY * (4* LJ_EPS * LJ_SIG**6) / LJ_CUT**3\n", "sim_energy = corrected_energy.mean()\n", "dsim_energy = corrected_energy.std()/np.sqrt(len(corrected_energy)-1)\n", "print(\"The obtained potential energy per particle is {:1.3f} +/- {:1.3f}, which agrees well with\\\n", @@ -973,4 +973,4 @@ }, "nbformat": 4, "nbformat_minor": 1 -} +} \ No newline at end of file From 0daa83995df7e8b0f3af2b4983c53dec0bc13c55 Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Mon, 28 Sep 2020 11:20:24 +0200 Subject: [PATCH 33/46] Adapt to PEP naming convention --- doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index 150b1bc028b..73d99cc4139 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -446,7 +446,7 @@ "DAMPING = 30\n", "MAX_STEPS = 10000\n", "MAX_DISPLACEMENT = 0.01 * LJ_SIG\n", - "EMSTEP = 10" + "EM_STEP = 10" ] }, { @@ -489,15 +489,15 @@ "energy = system.analysis.energy()['total']\n", "\n", "i = 0\n", - "while i < MAX_STEPS//EMSTEP:\n", + "while i < MAX_STEPS//EM_STEP:\n", " prev_max_force = max_force\n", " prev_energy = energy\n", - " system.integrator.run(EMSTEP)\n", + " system.integrator.run(EM_STEP)\n", " max_force = np.max(np.linalg.norm(system.part[:].f, axis=1))\n", " rel_force = np.abs((max_force-prev_max_force)/prev_max_force)\n", " energy = system.analysis.energy()['total']\n", " rel_ener = np.abs((energy-prev_energy)/prev_energy)\n", - " print(\"minimization step: {:4.0f}\\tmax. rel. force change:{:+3.3e}\\trel. energy change:{:+3.3e}\".format((i+1)*EMSTEP,rel_force, rel_ener))\n", + " print(\"minimization step: {:4.0f}\\tmax. rel. force change:{:+3.3e}\\trel. energy change:{:+3.3e}\".format((i+1)*EM_STEP,rel_force, rel_ener))\n", " if rel_force < F_TOL or rel_ener < E_TOL:\n", " break\n", " i += 1\n", From f270798669fafd4e3187d5841acd8a6e58a47970 Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Mon, 28 Sep 2020 11:24:04 +0200 Subject: [PATCH 34/46] Remove energy convergence criterion for steepest descent --- doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index 73d99cc4139..1464476a0f9 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -442,7 +442,6 @@ "outputs": [], "source": [ "F_TOL = 1e-2\n", - "E_TOL = 1e-5\n", "DAMPING = 30\n", "MAX_STEPS = 10000\n", "MAX_DISPLACEMENT = 0.01 * LJ_SIG\n", @@ -486,19 +485,15 @@ "# Initialize integrator to obtain initial forces\n", "system.integrator.run(0)\n", "max_force = np.max(np.linalg.norm(system.part[:].f, axis=1))\n", - "energy = system.analysis.energy()['total']\n", "\n", "i = 0\n", "while i < MAX_STEPS//EM_STEP:\n", " prev_max_force = max_force\n", - " prev_energy = energy\n", " system.integrator.run(EM_STEP)\n", " max_force = np.max(np.linalg.norm(system.part[:].f, axis=1))\n", " rel_force = np.abs((max_force-prev_max_force)/prev_max_force)\n", - " energy = system.analysis.energy()['total']\n", - " rel_ener = np.abs((energy-prev_energy)/prev_energy)\n", - " print(\"minimization step: {:4.0f}\\tmax. rel. force change:{:+3.3e}\\trel. energy change:{:+3.3e}\".format((i+1)*EM_STEP,rel_force, rel_ener))\n", - " if rel_force < F_TOL or rel_ener < E_TOL:\n", + " print(\"minimization step: {:4.0f}\\tmax. rel. force change:{:+3.3e}\".format((i+1)*EM_STEP,rel_force))\n", + " if rel_force < F_TOL:\n", " break\n", " i += 1\n", "```" From 14202ce4bbb2fb70fd46ae2c34a56c05ba756b4d Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Mon, 28 Sep 2020 11:24:32 +0200 Subject: [PATCH 35/46] Add more explanation for initial force calculation --- doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index 1464476a0f9..dcbe3ce3742 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -466,7 +466,7 @@ "* Break the minimization loop after a maximal number of MAX_STEPS steps or if convergence is achieved.\n", " Check for convergence every EMSTEP steps.\n", " \n", - "***Hint:*** To obtain the initial fores one has to initialize integrator using integ_steps=0." + "***Hint:*** To obtain the initial fores one has to initialize integrator using integ_steps=0, i.e. call system.integrator.run(0) before the force array can be accessed." ] }, { From 2c3f04159ea82f319dfff1a87dfe5ad2a295ab07 Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Mon, 28 Sep 2020 11:31:26 +0200 Subject: [PATCH 36/46] Add forgotten link --- doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index dcbe3ce3742..0f9a391912d 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -264,7 +264,7 @@ "\n", "* Create N_PART particles at random positions.\n", "\n", - " Use [system.part.add()](link).\n", + " Use [system.part.add()](http://espressomd.org/html/doc/espressomd.html#espressomd.particle_data.ParticleList).\n", " \n", " Either write a loop or use an (N_PART x 3) array for positions.\n", " Use numpy.random.random() to generate random numbers." From fecbdafa5d28ba67c5bf2e6b750cf93cc3f3fec3 Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Mon, 28 Sep 2020 14:44:17 +0200 Subject: [PATCH 37/46] Fix TOC to pass CI tests --- .../01-lennard_jones/01-lennard_jones.ipynb | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index 0f9a391912d..94f8720f9e5 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -19,16 +19,14 @@ "5. [First steps](#First-steps)\n", "6. [Overview of a simulation script](#Overview-of-a-simulation-script)\n", " 1. [System setup](#System-setup)\n", - " 2. [Choosing the thermodynamic ensemble, thermostat](#Choosing-the-thermodynamic-ensemble,-thermostat)\n", - " 3. [Placing and accessing particles](#Placing-and-accessing-particles)\n", - " 4. [Setting up non-bonded interactions](#Setting-up-non-bonded-interactions)\n", - " 5. [Warmup](#Warmup)\n", + " 2. [Placing and accessing particles](#Placing-and-accessing-particles)\n", + " 3. [Setting up non-bonded interactions](#Setting-up-non-bonded-interactions)\n", + " 4. [Energy minimization](#Energy-minimization)\n", + " 5. [Choosing the thermodynamic ensemble, thermostat](#Choosing-the-thermodynamic-ensemble,-thermostat)\n", " 6. [Integrating equations of motion and taking measurements](#Integrating-equations-of-motion-and-taking-measurements)\n", - " 7. [Simple Error Estimation on Time Series Data](#Simple-Error-Estimation-on-Time-Series-Data)\n", - "7. [Exercises](#Exercises)\n", + "7. [Futher Exercises](#Futher-Exercises)\n", " 1. [Binary Lennard-Jones Liquid](#Binary-Lennard-Jones-Liquid)\n", - "8. [Comparison with the literature](#Comparison-with-the-literature)\n", - "9. [References](#References)\n" + "8. [References](#References)\n" ] }, { From a537af4a83e96d1773b49998365710ec56d0283c Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Wed, 30 Sep 2020 09:24:32 +0200 Subject: [PATCH 38/46] Update doc/tutorials/01-lennard_jones/NotesForTutor.md MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Jean-Noël Grad --- doc/tutorials/01-lennard_jones/NotesForTutor.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/tutorials/01-lennard_jones/NotesForTutor.md b/doc/tutorials/01-lennard_jones/NotesForTutor.md index 44e81566367..0c469ad6f93 100644 --- a/doc/tutorials/01-lennard_jones/NotesForTutor.md +++ b/doc/tutorials/01-lennard_jones/NotesForTutor.md @@ -41,4 +41,4 @@ After the tutorial, students should be able to small time step, etc.). Additionally the steepest descent algorithm can get trapped in local minima and the convergence criterion is system-dependent. * Mention that accumulators are updated auomatically. -* Statistical independent measurements for t > correlation time. +* Statistically independent measurements for t > correlation time. From 01db3efe8af1c0a2278c8c03848f9a76367a886e Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Wed, 30 Sep 2020 09:25:46 +0200 Subject: [PATCH 39/46] Apply suggestions from code review MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Jean-Noël Grad --- doc/tutorials/01-lennard_jones/NotesForTutor.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/NotesForTutor.md b/doc/tutorials/01-lennard_jones/NotesForTutor.md index 0c469ad6f93..f84d39db9e6 100644 --- a/doc/tutorials/01-lennard_jones/NotesForTutor.md +++ b/doc/tutorials/01-lennard_jones/NotesForTutor.md @@ -1,6 +1,6 @@ # Hints for Tutors: Tutorial 01 - Lennard Jones -## Learning objectives (physics): +## Learning objectives (physics) After the tutorial, students should be able to @@ -18,8 +18,8 @@ After the tutorial, students should be able to After the tutorial, students should be able to -* instance a simulation object and set basic parameters like time step -* add particles to a simulation, +* instance a system object and set basic parameters like the time step +* add particles to a simulation * access and change particle properties * set up non-bonded interactions * explain the role of particle type From 2e0daad0d7ac33aabcb81d423c77d8cd17deee18 Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Wed, 30 Sep 2020 09:28:44 +0200 Subject: [PATCH 40/46] Apply suggestions from code review MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Jean-Noël Grad --- doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index 94f8720f9e5..3a4fd4e8715 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -393,7 +393,7 @@ "source": [ "\n", "\n", - "* Setup a Lennard-Jones interaction with $\\epsilon=$LJ_EPS and $\\sigma=$LJ_SIG which are cut and shifted to zero at $r_\\mathrm{c}=$LJ_CUT$\\times\\sigma$.\n", + "* Setup a Lennard-Jones interaction with $\\epsilon=$LJ_EPS and $\\sigma=$LJ_SIG which are cut and shifted to zero at $r_\\mathrm{c}=$LJ_CUT$\\times\\sigma$.\n", " To allow for comparison with the fundamental work on MD simulations of LJ systems [6] we don't shift the potential to zero at the cutoff and add $V_\\mathrm{lr}$ later. However, keep in mind that this choice will influence the thermodynamics of your simulation system!\n", "\n", "\n", @@ -455,7 +455,7 @@ "source": [ "**Exercise:**\n", "\n", - "* Use [espressomd.integrate](http://espressomd.org/html/doc/espressomd.html#module-espressomd.minimize_energy).set_steepest_descent to relax the initial configuration.\n", + "* Use [espressomd.integrate.set_steepest_descent](http://espressomd.org/html/doc/running.html#steepest-descent) to relax the initial configuration.\n", " Use a maximal displacement to MAX_DISPLACEMENT.\n", " A damping constant gamma = DAMPING usually is a good choice.\n", "* Use the relative change of the system total energy and the maximal force as a convergence criterion.\n", @@ -630,7 +630,7 @@ "**Part I:**\n", "* Set up the main integration loop.\n", " Take SAMPLING_ITERATIONS measurements every SAMPLING_INTERVAL integration steps.\n", - "* Monitor the potential and kinetic energies using the analysis method [system.analysis](http://espressomd.org/html/doc/espressomd.html#module-espressomd.analyze).energy().\n", + "* Monitor the potential and kinetic energies using the analysis method [system.analysis.energy()](http://espressomd.org/html/doc/espressomd.html#module-espressomd.analyze).\n", " The total energy time series should be stored in a numpy.array e_total\n", " and the kinetic energy in e_kin. Plot the timeseries of the total energy and\n", " of the instantaneous temperature $T = 3/2 \\times E_\\mathrm{kin}$/N_PART.\n", @@ -904,7 +904,7 @@ "A two-component Lennard-Jones liquid can be simulated by placing particles of two types (0 and 1) into the system. Depending on the Lennard-Jones parameters, the two components either mix or separate.\n", "\n", "1. Modify the code such that half of the particles are of type=1. Type 0 is implied for the remaining particles.\n", - "2. Specify Lennard-Jones interactions between type 0 particles with other type 0 particles, type 1 particles with other type 1 particles, and type 0 particles with type 1 particles (set parameters for system.non_bonded_inter[i,j].lennard_jones where {i,j} can be {0,0}, {1,1}, and {0,1}. Use the same Lennard-Jones parameters for interactions within a component, but use a different lj_cut_mixed parameter for the cutoff of the Lennard-Jones interaction between particles of type 0 and particles of type 1. Set this parameter to $2^{\\frac16}\\sigma$ to get de-mixing or to $2.5\\sigma$ to get mixing between the two components.\n", + "2. Specify Lennard-Jones interactions between type 0 particles with other type 0 particles, type 1 particles with other type 1 particles, and type 0 particles with type 1 particles (set parameters for system.non_bonded_inter[i,j].lennard_jones where {i,j} can be {0,0}, {1,1}, and {0,1}. Use the same Lennard-Jones parameters for interactions within a component, but use a different lj_cut_mixed parameter for the cutoff of the Lennard-Jones interaction between particles of type 0 and particles of type 1. Set this parameter to $2^{\\frac{1}{6}}\\sigma$ to get de-mixing or to $2.5\\sigma$ to get mixing between the two components.\n", "3. Record the radial distribution functions separately for particles of type 0 around particles of type 0, type 1 around particles of type 1, and type 0 around particles of type 1. This can be done by changing the ids1/ids2 arguments of the espressomd.observables.RDF command. You can record all three radial distribution functions in a single simulation. It is also possible to write them as several columns into a single file.\n", "4. Plot the radial distribution functions for all three combinations of particle types. The mixed case will differ significantly, depending on your choice of lj_cut_mixed. Explain these differences." ] @@ -966,4 +966,4 @@ }, "nbformat": 4, "nbformat_minor": 1 -} \ No newline at end of file +} From 675353c5ee75682d8a169692db0f5b4300547b9e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 30 Sep 2020 12:27:06 +0200 Subject: [PATCH 41/46] Formatting --- doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index 3a4fd4e8715..5c3d65854a5 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -674,12 +674,9 @@ "system.auto_update_accumulators.add(rdf_acc)\n", "\n", "for i in range(SAMPLING_ITERATIONS):\n", - " \n", " time[i] = system.time\n", - " print(i,'/',SAMPLING_ITERATIONS,end=\"\\r\")\n", - " \n", + " print(i + 1, '/', SAMPLING_ITERATIONS, end='\\r')\n", " system.integrator.run(SAMPLING_INTERVAL)\n", - " \n", " # Measure energies\n", " energies = system.analysis.energy()\n", " e_total[i] = energies['total']\n", @@ -835,8 +832,8 @@ "# expression of the factors Pi from Equations 2-8 with coefficients qi from Table 1\n", "# expression for a,g\n", "P = lambda q1, q2, q3, q4, q5, q6, q7, q8, q9: \\\n", - " q1 + q2 * np.exp(-q3 * T_star) + q4 * np.exp(-q5 * T_star) + q6 / rho_star \\\n", - " + q7 / rho_star**2 + q8 * np.exp(-q3 * T_star) / rho_star**3 + q9 * np.exp(-q5 * T_star) / rho_star**4\n", + " q1 + q2 * np.exp(-q3 * T_star) + q4 * np.exp(-q5 * T_star) + q6 / rho_star + q7 / rho_star**2 \\\n", + " + q8 * np.exp(-q3 * T_star) / rho_star**3 + q9 * np.exp(-q5 * T_star) / rho_star**4\n", "a = P(9.24792, -2.64281, 0.133386, -1.35932, 1.25338, 0.45602, -0.326422, 0.045708, -0.0287681)\n", "g = P(0.663161, -0.243089, 1.24749, -2.059, 0.04261, 1.65041, -0.343652, -0.037698, 0.008899)\n", "\n", From d35864918a9984be37ff196e8e7d3c9fa5994ff5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 30 Sep 2020 12:28:04 +0200 Subject: [PATCH 42/46] Remove custom Python kernel metadata Triggers an error message if the exact same python interpreter build cannot be found on the computer, with a message prompt to select one of the available python interpreters. --- doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index 5c3d65854a5..63faf782d34 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -926,9 +926,9 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3.8.5 64-bit", + "display_name": "Python 3", "language": "python", - "name": "python38564bit31b7b39ba35a4507a8c61be7a2214877" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -940,7 +940,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.5" + "version": "3.6.9" }, "latex_envs": { "LaTeX_envs_menu_present": true, From 4e7ee290d193601617c385977137e36b58d49c37 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 30 Sep 2020 12:31:09 +0200 Subject: [PATCH 43/46] Remove LaTeX plugin metadata --- .../01-lennard_jones/01-lennard_jones.ipynb | 18 ------------------ 1 file changed, 18 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index 63faf782d34..3446a2ca9cd 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -941,24 +941,6 @@ "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.9" - }, - "latex_envs": { - "LaTeX_envs_menu_present": true, - "autoclose": false, - "autocomplete": true, - "bibliofile": "biblio.bib", - "cite_by": "apalike", - "current_citInitial": 1, - "eqLabelWithNumbers": true, - "eqNumInitial": 1, - "hotkeys": { - "equation": "Ctrl-E", - "itemize": "Ctrl-I" - }, - "labels_anchors": false, - "latex_user_defs": false, - "report_style_numbering": false, - "user_envs_cfg": false } }, "nbformat": 4, From 6426cbff60ed4bd50088f3199f8fa187f1617147 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 30 Sep 2020 12:35:39 +0200 Subject: [PATCH 44/46] Fix overlapping legends in matplotlib plot --- .../01-lennard_jones/01-lennard_jones.ipynb | 23 ++++++++++--------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index 3446a2ca9cd..0a5b18bc2f8 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -699,23 +699,24 @@ "# Plot total energy and instantaneous temperature\n", "import matplotlib.pyplot as plt\n", "\n", - "fig1,ax = plt.subplots(num=None, figsize=(10, 6),\n", + "fig1,ax1 = plt.subplots(num=None, figsize=(10, 6),\n", " dpi=80, facecolor='w', edgecolor='k')\n", "fig1.set_tight_layout(False)\n", "\n", "t = np.arange(0,SAMPLING_ITERATIONS*SAMPLING_INTERVAL,SAMPLING_INTERVAL)\n", - "ax.plot(t,e_total, lw=2, label='e_total')\n", - "ax.set_xlabel('time step', fontsize=20)\n", - "ax.set_ylabel('total energy', fontsize=20)\n", - "ax.legend(fontsize=16)\n", - "ax.tick_params(axis='both', labelsize=16)\n", - "ax2 = ax.twinx()\n", - "ax2.plot(t,e_kin/ (1.5 * N_PART), color='r', ls='--', lw=2,\n", + "ln1 = ax1.plot(t, e_total, lw=2, label='total energy')\n", + "ax1.set_xlabel('time step', fontsize=20)\n", + "ax1.set_ylabel('Energy', fontsize=20)\n", + "ax1.tick_params(axis='both', labelsize=16)\n", + "ax2 = ax1.twinx()\n", + "ln2 = ax2.plot(t, e_kin/ (1.5 * N_PART), color='r', ls='--', lw=2,\n", " label='instantaneous temperature')\n", "ax2.axhline(TEMPERATURE, ls='-.', lw=3, color='k')\n", - "ax2.set_ylabel('instantaneous temperature', fontsize=20)\n", - "ax2.legend(fontsize=16)\n", + "ax2.set_ylabel('Temperature', fontsize=20)\n", "ax2.tick_params(axis='both', labelsize=16)\n", + "lns = ln1 + ln2\n", + "labs = [l.get_label() for l in lns]\n", + "ax1.legend(lns, labs, fontsize=16)\n", "print(\"average simulation temperature: {:1.3f}\".format(np.round(e_kin.mean()/ (1.5 * N_PART),3)),\n", " \"\\t(which differs by {:1.3f}% from TEMPERATURE)\".format(np.round(((e_kin.mean()/ (1.5 * N_PART) - TEMPERATURE) / TEMPERATURE)*100,3)))\n", "```" @@ -740,7 +741,7 @@ "solution2": "hidden" }, "source": [ - "``` python\n", + "```python\n", "# autocorrelation using numpy.correlate\n", "def autocor(x):\n", " mean=x.mean()\n", From d6b4d97d83386c3d44387b055ccaf2ddf93f2963 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 30 Sep 2020 12:36:19 +0200 Subject: [PATCH 45/46] Fix handling of matplotlib in hidden cells --- doc/tutorials/html_runner.py | 7 ++++--- testsuite/scripts/tutorials/test_html_runner.py | 7 +++++-- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/doc/tutorials/html_runner.py b/doc/tutorials/html_runner.py index f067d35e0ca..cbaf3b95403 100644 --- a/doc/tutorials/html_runner.py +++ b/doc/tutorials/html_runner.py @@ -196,15 +196,16 @@ def execute_notebook(nb, src, cell_separator): for filepath in args.scripts: add_cell_from_script(nb, filepath) +# convert solution cells to code cells +if args.exercise2: + convert_exercise2_to_code(nb) + # disable plot interactivity disable_plot_interactivity(nb) # guard against a jupyter bug involving matplotlib split_matplotlib_cells(nb) -if args.exercise2: - convert_exercise2_to_code(nb) - if args.substitutions or args.execute: # substitute global variables cell_separator = '\n##{}\n'.format(uuid.uuid4().hex) diff --git a/testsuite/scripts/tutorials/test_html_runner.py b/testsuite/scripts/tutorials/test_html_runner.py index 89b625046cb..8d8221675e9 100644 --- a/testsuite/scripts/tutorials/test_html_runner.py +++ b/testsuite/scripts/tutorials/test_html_runner.py @@ -157,7 +157,7 @@ def test_exercise2_plugin(self): cell_md['metadata']['solution2_first'] = True cell_md['metadata']['solution2'] = 'hidden' nb['cells'].append(cell_md) - code = '```python\n2\nglobal_var = 5\n```' + code = '```python\n2\nimport matplotlib.pyplot\nglobal_var = 5\n```' cell_md = nbformat.v4.new_markdown_cell(source=code) cell_md['metadata']['solution2'] = 'hidden' nb['cells'].append(cell_md) @@ -195,7 +195,10 @@ def test_exercise2_plugin(self): self.assertEqual(cell['source'], 'Question 2') cell = next(cells) self.assertEqual(cell['cell_type'], 'code') - self.assertEqual(cell['source'], '2\nglobal_var = 20') + self.assertEqual(cell['source'], '2\nimport matplotlib.pyplot') + cell = next(cells) + self.assertEqual(cell['cell_type'], 'code') + self.assertEqual(cell['source'], 'global_var = 20') cell = next(cells) self.assertEqual(cell['cell_type'], 'code') self.assertEqual(cell['source'], '3') From b5dd583c49436dbdd2874e69a4a7f0c71d53ea5e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 30 Sep 2020 12:52:30 +0200 Subject: [PATCH 46/46] Add plot legend --- doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index 0a5b18bc2f8..008e1bccd74 100644 --- a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb +++ b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb @@ -880,8 +880,9 @@ "\n", "plt.plot(r, theo_rdf, '-', color=\"blue\", linewidth=2, alpha=1, ls='--', label=r'theoretical')\n", "\n", - "plt.xlabel('$r [\\sigma]$', fontsize=20)\n", + "plt.xlabel('$r\\, [\\sigma]$', fontsize=20)\n", "plt.ylabel('$g(r)$', fontsize=20)\n", + "plt.legend()\n", "plt.gca().tick_params(axis='both', labelsize=16)\n", "```" ]