Added test criteria for the charged_system-2 tutorial

doc/sphinx/analysis.rst

@@ -669,7 +669,7 @@ which the observable should take into account.
 
 The current value of an observable can be obtained using its calculate()-method::
 
-    print(par_pos.calculate())
+    print(part_pos.calculate())
 
 .. _Available observables:
 

doc/tutorials/02-charged_system/02-charged_system-1.ipynb

@@ -19,7 +19,7 @@
     "#define EXTERNAL_FORCES\n",
     "#define MASS\n",
     "#define ELECTROSTATICS\n",
-    "#define LENNARD_JONES\n",
+    "#define WCA\n",
     "```"
    ]
   },
@@ -46,7 +46,7 @@
     "plt.ion()\n",
     "\n",
     "# Print enabled features\n",
-    "required_features = [\"EXTERNAL_FORCES\", \"MASS\", \"ELECTROSTATICS\", \"LENNARD_JONES\"]\n",
+    "required_features = [\"EXTERNAL_FORCES\", \"MASS\", \"ELECTROSTATICS\", \"WCA\"]\n",
     "espressomd.assert_features(required_features)\n",
     "print(espressomd.features())\n",
     "\n",
@@ -67,12 +67,8 @@
     "types       = {\"Anion\":          0, \"Cation\": 1}\n",
     "numbers     = {\"Anion\": n_ionpairs, \"Cation\": n_ionpairs}\n",
     "charges     = {\"Anion\":       -1.0, \"Cation\": 1.0}\n",
-    "lj_sigmas   = {\"Anion\":        1.0, \"Cation\": 1.0}\n",
-    "lj_epsilons = {\"Anion\":        1.0, \"Cation\": 1.0}\n",
-    "\n",
-    "WCA_cut = 2.**(1. / 6.)\n",
-    "lj_cuts     = {\"Anion\":  WCA_cut * lj_sigmas[\"Anion\"],\n",
-    "               \"Cation\": WCA_cut * lj_sigmas[\"Cation\"]}"
+    "wca_sigmas   = {\"Anion\":        1.0, \"Cation\": 1.0}\n",
+    "wca_epsilons = {\"Anion\":        1.0, \"Cation\": 1.0}"
    ]
   },
   {
@@ -139,12 +135,11 @@
    "metadata": {},
    "source": [
     "Before we can really start the simulation, we have to specify the\n",
-    "interactions between our particles. We already defined the Lennard-Jones parameters at the beginning,\n",
+    "interactions between our particles. We already defined the WCA parameters at the beginning,\n",
     "what is left is to specify the combination rule and to iterate over particle type pairs. For simplicity, \n",
     "we implement only the *Lorentz-Berthelot* rules. \n",
     "We pass our interaction pair to <tt>system.non_bonded_inter[\\*,\\*]</tt> and set the \n",
-    "pre-calculated LJ parameters <tt>epsilon</tt>, <tt>sigma</tt> and <tt>cutoff</tt>. With <tt>shift=\"auto\"</tt>,\n",
-    "we shift the interaction potential to the cutoff so that $U_\\mathrm{LJ}(r_\\mathrm{cutoff})=0$."
+    "pre-calculated WCA parameters <tt>epsilon</tt> and <tt>sigma</tt>."
    ]
   },
   {
@@ -165,14 +160,13 @@
     "    else:\n",
     "        return ValueError(\"No combination rule defined\")\n",
     "\n",
-    "# Lennard-Jones interactions parameters\n",
+    "# WCA interactions parameters\n",
     "for s in [[\"Anion\", \"Cation\"], [\"Anion\", \"Anion\"], [\"Cation\", \"Cation\"]]:\n",
-    "    lj_sig = combination_rule_sigma(\"Berthelot\", lj_sigmas[s[0]], lj_sigmas[s[1]])\n",
-    "    lj_cut = combination_rule_sigma(\"Berthelot\", lj_cuts[s[0]], lj_cuts[s[1]])\n",
-    "    lj_eps = combination_rule_epsilon(\"Lorentz\", lj_epsilons[s[0]], lj_epsilons[s[1]])\n",
+    "    wca_sig = combination_rule_sigma(\"Berthelot\", wca_sigmas[s[0]], wca_sigmas[s[1]])\n",
+    "    wca_eps = combination_rule_epsilon(\"Lorentz\", wca_epsilons[s[0]], wca_epsilons[s[1]])\n",
     "\n",
-    "    system.non_bonded_inter[types[s[0]], types[s[1]]].lennard_jones.set_params(\n",
-    "        epsilon=lj_eps, sigma=lj_sig, cutoff=lj_cut, shift=\"auto\")"
+    "    system.non_bonded_inter[types[s[0]], types[s[1]]].wca.set_params(\n",
+    "        epsilon=wca_eps, sigma=wca_sig)"
    ]
   },
   {
@@ -182,7 +176,7 @@
     "## 3 Equilibration\n",
     "\n",
     "With randomly positioned particles, we most likely have huge overlap and the strong repulsion will\n",
-    "cause the simulation to crash. The next step in our script therefore is a suitable LJ equilibration.\n",
+    "cause the simulation to crash. The next step in our script therefore is a suitable WCA equilibration.\n",
     "This is known to be a tricky part of a simulation and several approaches exist to reduce the particle overlap.\n",
     "Here, we use the steepest descent algorithm and cap the maximal particle displacement per integration step\n",
     "to 1% of sigma.\n",
@@ -197,8 +191,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Lennard-Jones Equilibration\n",
-    "max_sigma = max(lj_sigmas.values())\n",
+    "# WCA Equilibration\n",
+    "max_sigma = max(wca_sigmas.values())\n",
     "min_dist = 0.0\n",
     "system.minimize_energy.init(f_max=0, gamma=10.0, max_steps=10,\n",
     "                            max_displacement=max_sigma * 0.01)\n",
@@ -423,7 +417,7 @@
     "\n",
     "To make your life easier, don't try to equilibrate randomly positioned particles,\n",
     "but set them up in a crystal structure close to equilibrium. If you do it right,\n",
-    "you don't even need the Lennard-Jones equilibration. \n",
+    "you don't even need the WCA equilibration. \n",
     "To speed things up, don't go further than 1000 particles and use a P$^3$M accuracy of $10^{-2}$.\n",
     "Your RDF should look like the plot in figure 2. If you get stuck,\n",
     "you can look at the solution script <tt>/doc/tutorials/02-charged_system/scripts/nacl_units.py</tt> (or <tt>nacl_units_vis.py</tt> with visualization)."

doc/tutorials/02-charged_system/scripts/nacl.py

@@ -22,7 +22,7 @@
 from espressomd import assert_features, electrostatics
 import numpy
 
-assert_features(["ELECTROSTATICS", "LENNARD_JONES"])
+assert_features(["ELECTROSTATICS", "WCA"])
 
 print("\n--->Setup system")
 
@@ -43,12 +43,8 @@
 types = {"Anion": 0, "Cation": 1}
 numbers = {"Anion": n_ionpairs, "Cation": n_ionpairs}
 charges = {"Anion": -1.0, "Cation": 1.0}
-lj_sigmas = {"Anion": 1.0, "Cation": 1.0}
-lj_epsilons = {"Anion": 1.0, "Cation": 1.0}
-
-WCA_cut = 2.**(1. / 6.)
-lj_cuts = {"Anion": WCA_cut * lj_sigmas["Anion"],
-           "Cation": WCA_cut * lj_sigmas["Cation"]}
+wca_sigmas = {"Anion": 1.0, "Cation": 1.0}
+wca_epsilons = {"Anion": 1.0, "Cation": 1.0}
 
 # Setup System
 box_l = (n_part / density)**(1. / 3.)
@@ -83,18 +79,17 @@ def combination_rule_sigma(rule, sig1, sig2):
 
 # Lennard-Jones interactions parameters
 for s in [["Anion", "Cation"], ["Anion", "Anion"], ["Cation", "Cation"]]:
-    lj_sig = combination_rule_sigma(
-        "Berthelot", lj_sigmas[s[0]], lj_sigmas[s[1]])
-    lj_cut = combination_rule_sigma("Berthelot", lj_cuts[s[0]], lj_cuts[s[1]])
-    lj_eps = combination_rule_epsilon(
-        "Lorentz", lj_epsilons[s[0]], lj_epsilons[s[1]])
+    wca_sig = combination_rule_sigma(
+        "Berthelot", wca_sigmas[s[0]], wca_sigmas[s[1]])
+    wca_eps = combination_rule_epsilon(
+        "Lorentz", wca_epsilons[s[0]], wca_epsilons[s[1]])
 
-    system.non_bonded_inter[types[s[0]], types[s[1]]].lennard_jones.set_params(
-        epsilon=lj_eps, sigma=lj_sig, cutoff=lj_cut, shift="auto")
+    system.non_bonded_inter[types[s[0]], types[s[1]]].wca.set_params(
+        epsilon=wca_eps, sigma=wca_sig)
 
 
-print("\n--->Lennard-Jones Equilibration")
-max_sigma = max(lj_sigmas.values())
+print("\n--->WCA Equilibration")
+max_sigma = max(wca_sigmas.values())
 min_dist = 0.0
 system.minimize_energy.init(f_max=0, gamma=10, max_steps=10,
                             max_displacement=max_sigma * 0.01)

testsuite/scripts/tutorials/test_02-charged_system-2.py

@@ -17,16 +17,27 @@
 
 import unittest as ut
 import importlib_wrapper
+import numpy as np
 
 tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import(
     "@TUTORIALS_DIR@/02-charged_system/02-charged_system-2.py",
-    num_steps_equilibration=60, num_configs=5, integ_steps_per_config=60)
+    num_steps_equilibration=200, num_configs=5, integ_steps_per_config=60)
 
 
 @skipIfMissingFeatures
 class Tutorial(ut.TestCase):
     system = tutorial.system
 
+    def test_distribution(self):
+        """
+        checks if the particle distribution is within the box
+        """
+        for i in range(1, 3):
+            pos = np.flatnonzero(tutorial.res[:, i] > 0)
+            self.assertGreater(tutorial.res[pos[0], 0], tutorial.wall_margin)
+            self.assertLess(tutorial.res[pos[-1], 0],
+                            tutorial.box_z - tutorial.wall_margin)
+
 
 if __name__ == "__main__":
     ut.main()