Update tutorials and improve tutorial testing framework

doc/sphinx/electrostatics.rst

@@ -221,13 +221,13 @@ sets up parallel metallic plates and activates ICC::
     # Left electrode (normal [0,0,1])
     for xi in xrange(nicc):
         for yi in xrange(nicc):
-            system.part.add(pos=[l * xi, l * yi, 0], q=-0.0001, fix=[1, 1, 1], type=icc_type)
+            system.part.add(pos=[l * xi, l * yi, 0], q=-0.0001, fix=3*[True], type=icc_type)
     iccNormals.extend([0, 0, 1] * nicc_per_electrode)
 
     # Right electrode (normal [0,0,-1])
     for xi in xrange(nicc):
         for yi in xrange(nicc):
-            system.part.add(pos=[l * xi, l * yi, box_l], q=0.0001, fix=[1, 1, 1], type=icc_type)
+            system.part.add(pos=[l * xi, l * yi, box_l], q=0.0001, fix=3*[True], type=icc_type)
     iccNormals.extend([0, 0, -1] * nicc_per_electrode)
 
     # Common area, sigma and metallic epsilon

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

@@ -27,8 +27,9 @@
     "    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",
-    "    "
+    "8. [Comparison with the literature](#Comparison-with-the-literature)\n",
+    "9. [References](#References)\n",
+    ""
    ]
   },
   {
@@ -52,7 +53,7 @@
    "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]."
+    "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) <a href='#[1]'>[1]</a> 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 <a href='#[2]'>[2,3]</a>. 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 <a href='#[4]'>[4]</a>."
    ]
   },
   {
@@ -97,7 +98,7 @@
    "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",
+    "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 <a href='#[1]'>[1,4,2,3]</a>.\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",
@@ -120,7 +121,6 @@
    "outputs": [],
    "source": [
     "import espressomd\n",
-    "print(espressomd.features())\n",
     "required_features = [\"LENNARD_JONES\"]\n",
     "espressomd.assert_features(required_features)"
    ]
@@ -397,7 +397,7 @@
    "source": [
     "# Integration parameters\n",
     "sampling_interval = 100\n",
-    "sampling_iterations = 100\n",
+    "sampling_iterations = 200\n",
     "\n",
     "from espressomd.observables import ParticlePositions\n",
     "from espressomd.accumulators import Correlator\n",
@@ -412,7 +412,7 @@
     "system.auto_update_accumulators.add(msd_corr)\n",
     "\n",
     "# Set parameters for the radial distribution function\n",
-    "r_bins = 70\n",
+    "r_bins = 100\n",
     "r_min = 0.0\n",
     "r_max = system.box_l[0] / 2.0\n",
     "\n",
@@ -459,6 +459,13 @@
     "plt.ion()"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Plot the experimental radial distribution."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -468,11 +475,18 @@
     "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 $[\\sigma]$', fontsize=20)\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,
@@ -600,16 +614,74 @@
     "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": [
+    "## Comparison with the literature\n",
+    "\n",
+    "Empirical radial distribution functions have been determined for pure fluids <a href='#[5]'>[5]</a>, mixtures <a href='#[6]'>[6]</a> and confined fluids <a href='#[7]'>[7]</a>. We will compare our distribution $g(r)$ to the theoretical distribution $g(r^*, \\rho^*, T^*)$ of a pure fluid <a href='#[5]'>[5]</a>."
+   ]
+  },
+  {
+   "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, avg_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": {},
    "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."
+    "<a id='[1]'></a>[1] <a href=\"http://espressomd.org\">http://espressomd.org</a>  \n",
+    "<a id='[2]'></a>[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",
+    "<a id='[3]'></a>[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",
+    "<a id='[4]'></a>[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",
+    "<a id='[5]'></a>[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. <small>DOI:</small><a href=\"https://doi.org/10.1016/j.chemphys.2004.09.027\">10.1016/j.chemphys.2004.09.027</a>  \n",
+    "<a id='[6]'></a>[6] Matteoli. A simple expression for radial distribution functions of pure fluids and mixtures. *The Journal of Chemical Physics*, 103(11):4672, 1995. <small>DOI:</small><a href=\"https://doi.org/10.1063/1.470654\">10.1063/1.470654</a>  \n",
+    "<a id='[7]'></a>[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. <small>DOI:</small><a href=\"https://doi.org/10.1039/C5RA16151G\">10.1039/C5RA16151G</a>  "
    ]
   }
  ],

doc/tutorials/01-lennard_jones/CMakeLists.txt

@@ -3,4 +3,4 @@ configure_file(01-lennard_jones.ipynb ${CMAKE_CURRENT_BINARY_DIR}/01-lennard_jon
 
 add_custom_target(tutorials_01)
 
-nb_export(TUTORIAL tutorials_01 FILE "01-lennard_jones.ipynb" HTML_RUN VAR_SUBST "sampling_iterations=4000")
+nb_export(TUTORIAL tutorials_01 FILE "01-lennard_jones.ipynb" HTML_RUN VAR_SUBST)

doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part1.ipynb

@@ -54,8 +54,11 @@
     "as it is much (100x) faster than the CPU code.\n",
     "\n",
     "For the tutorial you will have to compile in the following features:\n",
-    "* <tt>LB_BOUNDARIES_GPU</tt>,\n",
-    "* <tt>LENNARD_JONES</tt>."
+    "```c++\n",
+    "#define LB_BOUNDARIES_GPU\n",
+    "#define LENNARD_JONES\n",
+    "```\n",
+    "Please uncomment them in the <tt>myconfig.hpp</tt> and compile **ESPResSo** using this <tt>myconfig.hpp</tt>."
    ]
   },
   {
@@ -325,7 +328,7 @@
     "import espressomd.shapes\n",
     "\n",
     "wall = espressomd.lbboundaries.LBBoundary(shape=espressomd.shapes.Wall(normal=[1, 0, 0], dist=1),\n",
-    "                               velocity=[0, 0, 0.01])\n",
+    "                                          velocity=[0, 0, 0.01])\n",
     "system.lbboundaries.add(wall)\n",
     "```\n",
     "\n",

doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part4.ipynb

@@ -80,7 +80,7 @@
     "system.actors.add(lbf)\n",
     "\n",
     "logging.info(\"Setup LB boundaries.\")\n",
-    "WALL_OFFSET = 2.0 * AGRID\n",
+    "WALL_OFFSET = AGRID\n",
     "top_wall = espressomd.lbboundaries.LBBoundary(shape=espressomd.shapes.Wall(normal=[1, 0, 0], dist=WALL_OFFSET))\n",
     "bottom_wall = espressomd.lbboundaries.LBBoundary(shape=espressomd.shapes.Wall(normal=[-1, 0, 0], dist=-(BOX_L - WALL_OFFSET)))\n",
     "\n",

doc/tutorials/04-lattice_boltzmann/CMakeLists.txt

@@ -7,9 +7,24 @@ configure_file(scripts/04-lattice_boltzmann_part4_solution.py ${CMAKE_CURRENT_BI
 configure_file(figures/latticeboltzmann-momentumexchange.png ${CMAKE_CURRENT_BINARY_DIR}/figures/latticeboltzmann-momentumexchange.png COPYONLY)
 configure_file(figures/latticeboltzmann-grid.png ${CMAKE_CURRENT_BINARY_DIR}/figures/latticeboltzmann-grid.png COPYONLY)
 
+add_custom_command(
+      OUTPUT "${CMAKE_CURRENT_BINARY_DIR}/scripts/04-lattice_boltzmann_part4_solution_cut.py"
+      DEPENDS "${CMAKE_CURRENT_BINARY_DIR}/scripts/04-lattice_boltzmann_part4_solution.py"
+      COMMAND
+      "cp"
+      "${CMAKE_CURRENT_BINARY_DIR}/scripts/04-lattice_boltzmann_part4_solution.py"
+      "${CMAKE_CURRENT_BINARY_DIR}/scripts/04-lattice_boltzmann_part4_solution_cut.py"
+      COMMAND
+      "sed"
+      "-n"
+      "-i"
+      "'/# Extract fluid velocity along the x-axis/,$$p;1i\\import matplotlib.pyplot as plt'"
+      "${CMAKE_CURRENT_BINARY_DIR}/scripts/04-lattice_boltzmann_part4_solution_cut.py"
+      )
+
 add_custom_target(tutorials_04)
 
 nb_export(TUTORIAL tutorials_04_1 FILE "04-lattice_boltzmann_part1.ipynb" HTML_RUN)
 nb_export(TUTORIAL tutorials_04_2 FILE "04-lattice_boltzmann_part2.ipynb" HTML_RUN)
 nb_export(TUTORIAL tutorials_04_3 FILE "04-lattice_boltzmann_part3.ipynb" HTML_RUN)
-nb_export(TUTORIAL tutorials_04_4 FILE "04-lattice_boltzmann_part4.ipynb" HTML_RUN)
+nb_export(TUTORIAL tutorials_04_4 FILE "04-lattice_boltzmann_part4.ipynb" HTML_RUN ADD_SCRIPTS "${CMAKE_CURRENT_BINARY_DIR}/scripts/04-lattice_boltzmann_part4_solution_cut.py")

doc/tutorials/04-lattice_boltzmann/scripts/04-lattice_boltzmann_part4_solution.py

@@ -69,19 +69,18 @@ def poiseuille_flow(x, force_density, dynamic_viscosity, height):
         (height**2.0 / 4.0 - x**2.0)
 
 
+# Note that the LB viscosity is not the dynamic viscosity but the
+# kinematic viscosity (mu=LB_viscosity * density)
 x_values = np.linspace(0.0, BOX_L, lbf.shape[0])
 HEIGHT = BOX_L - 2.0 * AGRID
+# analytical curve
+y_values = poiseuille_flow(x_values - (HEIGHT / 2.0 + AGRID), FORCE_DENSITY[1],
+                           VISCOSITY * DENSITY, HEIGHT)
+# velocity is zero outside the walls
+y_values[np.nonzero(x_values < WALL_OFFSET)] = 0.0
+y_values[np.nonzero(x_values > BOX_L - WALL_OFFSET)] = 0.0
 
-# Note that the LB viscosity is not the dynamic viscosity but the
-# kinematic viscosity (mu=LB_viscosity * density)
-plt.plot(
-    x_values,
-    poiseuille_flow(x_values - (HEIGHT / 2.0 + AGRID),
-                    FORCE_DENSITY[1],
-                    VISCOSITY * DENSITY,
-                    HEIGHT),
-    'o-',
-    label='analytical')
+plt.plot(x_values, y_values, 'o-', label='analytical')
 plt.plot(fluid_velocities[:, 0], fluid_velocities[:, 1], label='simulation')
 plt.legend()
 plt.show()