doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 1: Lennard-Jones Liquid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of Contents\n",
    "1. [Introduction](#Introduction)\n",
    "2. [Background](#Background)\n",
    "3. [The Lennard-Jones Potential](#The-Lennard-Jones-Potential)\n",
    "4. [Units](#Units)\n",
    "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",
    "    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",
    "    1. [Binary Lennard-Jones Liquid](#Binary-Lennard-Jones-Liquid)\n",
    "8. [References](#References)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "Welcome to the basic ESPResSo tutorial!\n",
    "\n",
    "In this tutorial, you will learn, how to use the ESPResSo package for your \n",
    "research. We will cover the basics of ESPResSo, i.e., how to set up and modify a physical system, how to run a simulation, and how to load, save and analyze the produced simulation data.\n",
    "\n",
    "More advanced features and algorithms available in the ESPResSo package are \n",
    "described in additional tutorials."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Background\n",
    "\n",
    "Today's research on Soft Condensed Matter has brought the needs for having a flexible, extensible, reliable, and efficient (parallel) molecular simulation package. For this reason ESPResSo (Extensible Simulation Package for Research on Soft Matter Systems) [1] has been developed at the Max Planck Institute for Polymer Research, Mainz, and at the Institute for Computational Physics at the University of Stuttgart in  the group of Prof. Dr. Christian Holm [2,3]. The ESPResSo package is probably the most flexible and extensible simulation package in the market. It is specifically developed for coarse-grained molecular dynamics (MD) simulation of polyelectrolytes but is not necessarily limited to this. For example, it could also be used to simulate granular media. ESPResSo has been nominated for the Heinz-Billing-Preis for Scientific Computing in 2003 [4]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Lennard-Jones Potential\n",
    "\n",
    "A pair  of  neutral  atoms  or  molecules  is  subject  to  two  distinct  forces  in  the  limit  of large separation and small separation:  an attractive force at long ranges (van der Waals force, or dispersion force) and a repulsive force at short ranges (the result of overlapping electron  orbitals,  referred  to  as  Pauli  repulsion  from  the Pauli  exclusion principle). The Lennard-Jones potential (also  referred  to  as  the  L-J potential, 6-12 potential  or, less commonly, 12-6 potential) is a simple mathematical model that represents this behavior. It  was  proposed  in  1924  by  John  Lennard-Jones. The  L-J  potential  is  of  the  form\n",
    "\n",
    "\\begin{equation}\n",
    "V(r) = 4\\epsilon \\left[ \\left( \\dfrac{\\sigma}{r} \\right)^{12} - \\left( \\dfrac{\\sigma}{r} \\right)^{6} \\right]\n",
    "\\end{equation}\n",
    "\n",
    "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",
    "\n",
    "<figure>\n",
    "<img src='figures/lennard-jones-potential.png' alt='missing' style='width: 600px;'/>\n",
    "<center>\n",
    "<figcaption>Figure 1: Lennard-Jones potential</figcaption>\n",
    "</center>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Units\n",
    "\n",
    "Novice users must understand that Espresso has no fixed unit system. The unit \n",
    "system is set by the user. Conventionally, reduced units are employed, in other \n",
    "words LJ units."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## First steps\n",
    "\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",
    "\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",
    "The scripting interface (Python) is used to setup the system (particles, boundary conditions, interactions etc.), control the simulation, run analysis, and store and load results. The user has at hand the full reliability and functionality of the scripting language. For instance, it is possible to use the SciPy package for analysis and PyPlot for plotting.\n",
    "With a certain overhead in efficiency, it can also be bused to reject/accept new configurations in combined MD/MC schemes. In principle, any parameter which is accessible from the scripting level can be changed at any moment of runtime. In this way methods like thermodynamic integration become readily accessible.\n",
    "\n",
    "_Note: This tutorial assumes that you already have a working ESPResSo\n",
    "installation on your system. If this is not the case, please consult the first chapters of the user's guide for installation instructions._\n",
    "\n",
    "Python simulation scripts can be run conveniently:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import espressomd\n",
    "print(espressomd.features())\n",
    "required_features = [\"LENNARD_JONES\"]\n",
    "espressomd.assert_features(required_features)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview of a simulation script\n",
    "\n",
    "Typically, a simulation script consists of the following parts:\n",
    "\n",
    "* System setup (box geometry, thermodynamic ensemble, integrator parameters)\n",
    "* Placing the particles\n",
    "* Setup of interactions between particles\n",
    "* Warm up (bringing the system into a state suitable for measurements)\n",
    "* Integration loop (propagate the system in time and record measurements)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### System setup\n",
    "\n",
    "The functionality of ESPResSo for python is provided via a python module called <tt>espressomd</tt>. At the beginning of the simulation script, it has to be imported."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Importing other relevant python modules\n",
    "import numpy as np\n",
    "# System parameters\n",
    "n_part = 100\n",
    "density = 0.5\n",
    "\n",
    "box_l = np.power(n_part / density, 1.0 / 3.0) * np.ones(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step would be to create an instance of the System class and to seed espresso. 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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "system = espressomd.System(box_l=box_l)\n",
    "system.seed = 42"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can be used to manipulate the crucial system parameters like the time step and the size of the simulation box (<tt>time_step</tt>, and <tt>box_l</tt>)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "skin = 0.4\n",
    "time_step = 0.01\n",
    "eq_tstep = 0.001\n",
    "temperature = 0.728\n",
    "\n",
    "system.time_step = time_step\n",
    "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "system.thermostat.turn_off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "system.thermostat.set_langevin(kT=temperature, gamma=1.0, seed=42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use a Langevin thermostat (NVT or NPT ensemble) with temperature set to 1.0 and damping coefficient to 0.5."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Placing and accessing particles\n",
    "\n",
    "Particles in the simulation can be added and accessed via the <tt>part</tt> property of the System class. Individual  particles  are  referred  to  by  an  integer  id, e.g., <tt>system.part[0]</tt>. If <tt>id</tt> 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' <tt>type</tt> attribute. If it is not specified, zero is implied."
   ]
  },
  {
   "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(system.part[0].pos)\n",
    "\n",
    "# Iterate over the first five particles for the purpose of demonstration.\n",
    "# For accessing all particles, do not splice system.part\n",
    "for p in system.part[:5]:\n",
    "    print(p.pos)\n",
    "    print(p.v)\n",
    "\n",
    "# Obtain all particle positions\n",
    "cur_pos = system.part[:].pos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many objects in ESPResSo have a string representation, and thus can be displayed via python's <tt>print</tt> function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(system.part[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting up non-bonded interactions\n",
    "\n",
    "Non-bonded interactions act between all particles of a given combination of particle types. In this tutorial, we use the Lennard-Jones non-bonded interaction. The interaction of two particles of type 0 can be setup as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Warmup\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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Integrating equations of motion and taking measurements\n",
    "\n",
    "Once warmup is done, the force capping is switched off by setting it to zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "system.force_cap = 0"
   ]
  },
  {
   "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",
    "\n",
    "The potential and kinetic energies can be monitored using the analysis method <tt>system.analysis.energy()</tt>. <tt>kinetic_temperature</tt> here refers to the measured temperature obtained from kinetic energy and the number of degrees of freedom in the system. It should fluctuate around the preset temperature of the thermostat.\n",
    "\n",
    "The particles' mean square displacement,\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathrm{msd}(t) =\\langle (x(t_0+t) -x(t_0))^2\\rangle,\n",
    "\\end{equation}\n",
    "\n",
    "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Integration parameters\n",
    "sampling_interval = 100\n",
    "sampling_iterations = 100\n",
    "\n",
    "from espressomd.observables import ParticlePositions\n",
    "from espressomd.accumulators import Correlator\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=10000 * 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 = 70\n",
    "r_min = 0.0\n",
    "r_max = system.box_l[0] / 2.0\n",
    "\n",
    "avg_rdf = np.zeros((r_bins,))\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 radial distribution function\n",
    "    r, rdf = system.analysis.rdf(rdf_type=\"rdf\", type_list_a=[0], type_list_b=[0],\n",
    "                                 r_min=r_min, r_max=r_max, r_bins=r_bins)\n",
    "    avg_rdf += rdf / sampling_iterations\n",
    "\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()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now use the plotting library <tt>matplotlib</tt> available in Python to visualize the measurements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.ion()"
   ]
  },
  {
   "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, avg_rdf, '-', color=\"#A60628\", linewidth=2, alpha=1)\n",
    "plt.xlabel('$r$', fontsize=20)\n",
    "plt.ylabel('$g(r)$', fontsize=20)\n",
    "plt.show()"
   ]
  },
  {
   "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('Time', fontsize=20)\n",
    "plt.ylabel('Temperature', fontsize=20)\n",
    "plt.legend(fontsize=16, loc=0)\n",
    "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, we may compute the actual temperature of the system via $k_B T= 2/(3N) \\langle E_\\text{kin}\\rangle$. The temperature is fixed and does not fluctuate in the NVT ensemble! The instantaneous temperature is calculated via $2/(3N) E_\\text{kin}$ (without ensemble averaging), but it is not the temperature of the system."
   ]
  },
  {
   "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",
    "plt.plot(msd[:, 0], msd[:, 2] + msd[:, 3] + msd[:, 4],\n",
    "         'o-', color=\"#348ABD\", linewidth=2, alpha=1)\n",
    "plt.xlabel('Time', fontsize=20)\n",
    "plt.ylabel('Mean squared displacement', fontsize=20)\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "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\n",
    "standard_error_total_energy = np.sqrt(etotal.var()) / np.sqrt(sampling_iterations)\n",
    "print(standard_error_total_energy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 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 <tt>type=1</tt>. 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 <tt>system.non_bonded_inter[i,j].lennard_jones</tt> where <tt>{i,j}</tt> can be <tt>{0,0}</tt>, <tt>{1,1}</tt>, and <tt>{0,1}</tt>. Use the same Lennard-Jones parameters for interactions within a component, but use a different <tt>lj_cut_mixed</tt> 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 <tt>type_list</tt> arguments of the <tt>system.analysis.rdf()</tt> 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 <tt>lj_cut_mixed</tt>. Explain these differences."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "[1] <a href=\"http://espressomd.org\">http://espressomd.org</a>  \n",
    "[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."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}