diff --git a/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb b/doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb index e2c79277544..008e1bccd74 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" ] }, { @@ -70,7 +68,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", @@ -99,8 +97,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", @@ -155,10 +152,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)" ] @@ -170,13 +175,49 @@ "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", + "\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 documentation](http://espressomd.org/html/doc/espressomd.html#module-espressomd.system)." + ] + }, + { + "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))" ] }, { @@ -195,9 +236,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" ] @@ -206,29 +244,40 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Choosing the thermodynamic ensemble, thermostat\n", + "### Placing and accessing particles\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", + "Particles in the simulation can be added and accessed via the part property of the System class. Individual particles are referred to by an integer id, e.g., system.part[0]. If id is unspecified, an unused particle id is automatically assigned. It is also possible to use common python iterators and slicing operations to add or access several particles at once.\n", "\n", - "The NVE ensemble is simulated without a thermostat. A previously enabled thermostat can be switched off as follows:" + "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": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, "source": [ - "system.thermostat.turn_off()" + "\n", + "\n", + "* Create N_PART particles at random positions.\n", + "\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." ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "solution2": "hidden" + }, "source": [ - "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:" + "```python\n", + "# Add particles to the simulation box at random positions\n", + "system.part.add(type=[0]*N_PART, pos=np.random.random((N_PART,3)) * system.box_l)\n", + "```" ] }, { @@ -236,26 +285,23 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "system.thermostat.set_langevin(kT=TEMPERATURE, gamma=GAMMA, seed=42)" - ] + "source": [] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "Use a Langevin thermostat (NVT or NPT ensemble) with temperature set to `temperature` and damping coefficient to `GAMMA`." + "# Test that now we have indeed N_PART particles in the system\n", + "assert(len(system.part) == N_PART)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Placing and accessing particles\n", - "\n", - "Particles in the simulation can be added and accessed via the part property of the System class. Individual particles are referred to by an integer id, e.g., system.part[0]. If id is unspecified, an unused particle id is automatically assigned. It is also possible to use common python iterators and slicing operations to add or access several particles at once.\n", - "\n", - "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." + "The particle properties can be accessed using standard numpy slicing syntax:" ] }, { @@ -264,10 +310,6 @@ "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", @@ -312,22 +354,25 @@ "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\n", - "LJ_CAP = 0.5\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" + "LJ_CUT = 2.5 * LJ_SIG" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Warmup\n", + "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", - "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:" + "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:" ] }, { @@ -336,27 +381,35 @@ "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" + "assert((BOX_L-2*SKIN > LJ_CUT).all())" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, "source": [ - "### Integrating equations of motion and taking measurements\n", + "\n", "\n", - "Once warmup is done, the force capping is switched off by setting it to zero." + "* 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)." + ] + }, + { + "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=0)\n", + "```" ] }, { @@ -364,27 +417,20 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "system.force_cap = 0" - ] + "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "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", + "### Energy minimization\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", + "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", - "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." + "**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." ] }, { @@ -393,78 +439,70 @@ "metadata": {}, "outputs": [], "source": [ - "# Integration parameters\n", - "sampling_interval = 100\n", - "sampling_iterations = 200\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(1, sampling_iterations + 1):\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", - "\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()" + "F_TOL = 1e-2\n", + "DAMPING = 30\n", + "MAX_STEPS = 10000\n", + "MAX_DISPLACEMENT = 0.01 * LJ_SIG\n", + "EM_STEP = 10" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, "source": [ - "We now use the plotting library matplotlib available in Python to visualize the measurements." + "**Exercise:**\n", + "\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", + " 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, i.e. call system.integrator.run(0) before the force array can be accessed." ] }, { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "cell_type": "markdown", + "metadata": { + "code_folding": [], + "solution2": "hidden" + }, "source": [ - "import matplotlib.pyplot as plt\n", - "plt.ion()" + "```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", + "max_force = np.max(np.linalg.norm(system.part[:].f, axis=1))\n", + "\n", + "i = 0\n", + "while i < MAX_STEPS//EM_STEP:\n", + " prev_max_force = max_force\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", + " 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", + "```" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, - "source": [ - "Plot the experimental radial distribution." - ] + "outputs": [], + "source": [] }, { "cell_type": "code", @@ -472,19 +510,27 @@ "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()" + "# check that after the exercise the total energy is negative\n", + "assert(system.analysis.energy()['total'] < 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Now plot the instantaneous temperature." + "### 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:" ] }, { @@ -493,33 +539,39 @@ "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()" + "# 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" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, "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." + "**Exercise:**\n", + "\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", + "For details see the [online documentation](http://espressomd.org/html/doc/espressomd.html#module-espressomd.thermostat)." ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "solution2": "hidden" + }, "source": [ - "The correlator output is stored in the array `msd` and has the following shape:" + "```python\n", + "system.integrator.set_vv()\n", + "system.thermostat.set_langevin(kT=TEMPERATURE, gamma=GAMMA, seed=42)\n", + "```" ] }, { @@ -527,146 +579,333 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "print(msd.shape)" - ] + "source": [] }, { "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..." + "### Integrating equations of motion and taking measurements\n", + "\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", + "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." ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "code_folding": [] + }, "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()" + "# Integration parameters\n", + "SAMPLING_INTERVAL = 20\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": {}, + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, "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", + "**Exercise:**\n", + "\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.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", + " 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", + " 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." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "code_folding": [], + "solution2": "hidden" + }, + "source": [ + "```python\n", + "e_total = np.zeros(SAMPLING_ITERATIONS)\n", + "e_kin = np.zeros(SAMPLING_ITERATIONS)\n", + "time = np.zeros(SAMPLING_ITERATIONS)\n", "\n", - "\\begin{equation}\n", - " SE = \\sqrt{\\frac{\\sigma^2}{N}},\n", - "\\end{equation}\n", + "from espressomd.observables import RDF\n", + "from espressomd.accumulators import Correlator, MeanVarianceCalculator\n", "\n", - "where $\\sigma^2$ is the variance\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", - "\\begin{equation}\n", - " \\sigma^2 = \\left\\langle x^2 - \\langle x\\rangle^2 \\right\\rangle\n", - "\\end{equation}" + "for i in range(SAMPLING_ITERATIONS):\n", + " time[i] = system.time\n", + " print(i + 1, '/', SAMPLING_ITERATIONS, end='\\r')\n", + " system.integrator.run(SAMPLING_INTERVAL)\n", + " # Measure energies\n", + " energies = system.analysis.energy()\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", + "r = rdf_obs.bin_centers()\n", + "```" ] }, { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "cell_type": "markdown", + "metadata": { + "code_folding": [], + "solution2": "hidden" + }, "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)" + "```python\n", + "# Plot total energy and instantaneous temperature\n", + "import matplotlib.pyplot as plt\n", + "\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", + "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('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", + "```" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "solution2": "hidden" + }, "source": [ - "## Exercises" + "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": {}, + "metadata": { + "code_folding": [], + "solution2": "hidden" + }, "source": [ - "### Binary Lennard-Jones Liquid\n", + "```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", - "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", + "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", - "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", - "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." + "from scipy.optimize import curve_fit\n", + "fitfn = lambda x,a,t: a*np.exp(-x/t)\n", + "popt,pcov = curve_fit(fitfn, time/(system.time_step*SAMPLING_INTERVAL), acf)\n", + "tcor = popt[1]\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_ylim(bottom=1e-3)\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(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", + "```" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "solution2": "hidden" + }, "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]." + "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, + "cell_type": "markdown", "metadata": { - "scrolled": true + "solution2": "hidden" }, - "outputs": [], "source": [ + "```python\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", + "the value -5.38 given in [6]\".format(*np.round(np.array([sim_energy,dsim_energy]),3)))\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "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]." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "code_folding": [], + "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=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 + 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", - "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", - "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()" + "\n", + "# plot fitted expression (=theoretical curve)\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**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, 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.legend()\n", + "plt.gca().tick_params(axis='both', labelsize=16)\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Futher Exercises" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Binary Lennard-Jones Liquid\n", + "\n", + "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^{\\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." ] }, { @@ -679,9 +918,11 @@ "[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 " ] } ], @@ -701,7 +942,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.8" + "version": "3.6.9" } }, "nbformat": 4, diff --git a/doc/tutorials/01-lennard_jones/NotesForTutor.md b/doc/tutorials/01-lennard_jones/NotesForTutor.md new file mode 100644 index 00000000000..f84d39db9e6 --- /dev/null +++ b/doc/tutorials/01-lennard_jones/NotesForTutor.md @@ -0,0 +1,44 @@ +# Hints for Tutors: Tutorial 01 - Lennard Jones + +## Learning objectives (physics) + +After the tutorial, students should be able to + +* explain + * 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! + * 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 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 +* 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. +* 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 - 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. +* Statistically independent measurements for t > correlation time. 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_01-lennard_jones.py b/testsuite/scripts/tutorials/test_01-lennard_jones.py index ad7e060efd9..62e99085dc7 100644 --- a/testsuite/scripts/tutorials/test_01-lennard_jones.py +++ b/testsuite/scripts/tutorials/test_01-lennard_jones.py @@ -20,20 +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_global_variables(self): - self.assertLess(tutorial.standard_error_total_energy, 2.5) + def test_rdf(self): + np.testing.assert_allclose(tutorial.rdf, tutorial.theo_rdf, + rtol=0.0, atol=0.1) - def test_rdf_curve(self): - self.assertGreater( - np.corrcoef(tutorial.rdf, tutorial.theo_rdf)[1, 0], 0.985) + def test_potential_energy(self): + # 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__": 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')