From c2c2a424a87843ef3b958a3a4232356777067568 Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Mon, 7 Oct 2019 15:24:33 +0200 Subject: [PATCH 01/10] added test for second charged system tutorial --- testsuite/scripts/tutorials/test_02-charged_system-2.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/testsuite/scripts/tutorials/test_02-charged_system-2.py b/testsuite/scripts/tutorials/test_02-charged_system-2.py index e0cbca913ce..ceed7a6ca14 100644 --- a/testsuite/scripts/tutorials/test_02-charged_system-2.py +++ b/testsuite/scripts/tutorials/test_02-charged_system-2.py @@ -17,6 +17,7 @@ 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", @@ -26,6 +27,14 @@ @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__": From 0f4ec356a856520234863ab444e706cdcc8923e9 Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Mon, 7 Oct 2019 15:27:49 +0200 Subject: [PATCH 02/10] fixed typo --- doc/sphinx/analysis.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/sphinx/analysis.rst b/doc/sphinx/analysis.rst index e7c90cff07e..5dcf88f44d0 100644 --- a/doc/sphinx/analysis.rst +++ b/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: From e55071cfc2cb20fd0c791be169bccb77b9defee4 Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Tue, 8 Oct 2019 13:59:27 +0200 Subject: [PATCH 03/10] added mean value check for charge test --- testsuite/scripts/tutorials/test_02-charged_system-2.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/testsuite/scripts/tutorials/test_02-charged_system-2.py b/testsuite/scripts/tutorials/test_02-charged_system-2.py index ceed7a6ca14..d51fe221cc4 100644 --- a/testsuite/scripts/tutorials/test_02-charged_system-2.py +++ b/testsuite/scripts/tutorials/test_02-charged_system-2.py @@ -35,6 +35,8 @@ def test_distribution(self): self.assertGreater(tutorial.res[pos[0], 0], tutorial.wall_margin) self.assertLess(tutorial.res[pos[-1], 0], tutorial.box_z - tutorial.wall_margin) + self.assertAlmostEqual(np.mean(tutorial.res[pos, i]), 0.0148, delta=1e-3) + if __name__ == "__main__": From ac47eef1e24c4ffb439797225411b29e4de901ee Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Tue, 8 Oct 2019 14:10:23 +0200 Subject: [PATCH 04/10] changed test criteria --- .../scripts/tutorials/test_02-charged_system-2.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/testsuite/scripts/tutorials/test_02-charged_system-2.py b/testsuite/scripts/tutorials/test_02-charged_system-2.py index d51fe221cc4..f20ac9d28f6 100644 --- a/testsuite/scripts/tutorials/test_02-charged_system-2.py +++ b/testsuite/scripts/tutorials/test_02-charged_system-2.py @@ -29,13 +29,19 @@ class Tutorial(ut.TestCase): system = tutorial.system def test_distribution(self): - """checks if the particle distribution is within the box""" + """ + checks if the particle distribution is within the box + and the mean of the distribution matches to a reference value + """ + avg = 0 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) - self.assertAlmostEqual(np.mean(tutorial.res[pos, i]), 0.0148, delta=1e-3) + avg += np.mean(tutorial.res[pos, i]) + avg *= 0.5 + self.assertAlmostEqual(avg, 0.0148, delta=5e-4) From 4bc24b0bf1aeb4424e043f516069f735b6758cf5 Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Tue, 8 Oct 2019 14:16:09 +0200 Subject: [PATCH 05/10] updated docstring --- .../scripts/tutorials/test_02-charged_system-2.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/testsuite/scripts/tutorials/test_02-charged_system-2.py b/testsuite/scripts/tutorials/test_02-charged_system-2.py index f20ac9d28f6..59cefa7a8f0 100644 --- a/testsuite/scripts/tutorials/test_02-charged_system-2.py +++ b/testsuite/scripts/tutorials/test_02-charged_system-2.py @@ -30,16 +30,16 @@ class Tutorial(ut.TestCase): def test_distribution(self): """ - checks if the particle distribution is within the box - and the mean of the distribution matches to a reference value + checks if the particle distribution is within the box and + the mean of the distribution matches with a reference value """ avg = 0 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) + 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) - avg += np.mean(tutorial.res[pos, i]) + avg += np.mean(tutorial.res[pos, i]) avg *= 0.5 self.assertAlmostEqual(avg, 0.0148, delta=5e-4) From e9edf01cc001373a33dd3fad1357f78461824e39 Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Tue, 8 Oct 2019 14:21:40 +0200 Subject: [PATCH 06/10] formatting --- testsuite/scripts/tutorials/test_02-charged_system-2.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/testsuite/scripts/tutorials/test_02-charged_system-2.py b/testsuite/scripts/tutorials/test_02-charged_system-2.py index 59cefa7a8f0..e4f32b821ff 100644 --- a/testsuite/scripts/tutorials/test_02-charged_system-2.py +++ b/testsuite/scripts/tutorials/test_02-charged_system-2.py @@ -27,7 +27,7 @@ @skipIfMissingFeatures class Tutorial(ut.TestCase): system = tutorial.system - + def test_distribution(self): """ checks if the particle distribution is within the box and @@ -37,12 +37,12 @@ def test_distribution(self): 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) - + self.assertLess(tutorial.res[pos[-1], 0], + tutorial.box_z - tutorial.wall_margin) + avg += np.mean(tutorial.res[pos, i]) avg *= 0.5 self.assertAlmostEqual(avg, 0.0148, delta=5e-4) - if __name__ == "__main__": From 64de065b186c1fbf7a7681e2de4baa61b992ee55 Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Fri, 11 Oct 2019 09:51:36 +0200 Subject: [PATCH 07/10] removed useless mean density check --- testsuite/scripts/tutorials/test_02-charged_system-2.py | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/testsuite/scripts/tutorials/test_02-charged_system-2.py b/testsuite/scripts/tutorials/test_02-charged_system-2.py index e4f32b821ff..b72184b85c1 100644 --- a/testsuite/scripts/tutorials/test_02-charged_system-2.py +++ b/testsuite/scripts/tutorials/test_02-charged_system-2.py @@ -30,20 +30,14 @@ class Tutorial(ut.TestCase): def test_distribution(self): """ - checks if the particle distribution is within the box and - the mean of the distribution matches with a reference value + checks if the particle distribution is within the box """ - avg = 0 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) - avg += np.mean(tutorial.res[pos, i]) - avg *= 0.5 - self.assertAlmostEqual(avg, 0.0148, delta=5e-4) - if __name__ == "__main__": ut.main() From ec1e7853572e4de63c58c4463ce2f7208b0ca7af Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Fri, 11 Oct 2019 13:39:48 +0200 Subject: [PATCH 08/10] changed LJ to WCA in nacl tutorial script --- .../02-charged_system/scripts/nacl.py | 27 ++++++++----------- 1 file changed, 11 insertions(+), 16 deletions(-) diff --git a/doc/tutorials/02-charged_system/scripts/nacl.py b/doc/tutorials/02-charged_system/scripts/nacl.py index 899ca7b38d7..0a1c38645b2 100644 --- a/doc/tutorials/02-charged_system/scripts/nacl.py +++ b/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) From ec2e064e79c82b01b7ab27a891f8197825e947bb Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Fri, 11 Oct 2019 13:46:05 +0200 Subject: [PATCH 09/10] replaced LJ with WCA in charged_system-1 tutorial --- .../02-charged_system-1.ipynb | 36 ++++++++----------- 1 file changed, 15 insertions(+), 21 deletions(-) diff --git a/doc/tutorials/02-charged_system/02-charged_system-1.ipynb b/doc/tutorials/02-charged_system/02-charged_system-1.ipynb index 9b3dc777b46..a9bb5e41d6d 100644 --- a/doc/tutorials/02-charged_system/02-charged_system-1.ipynb +++ b/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 system.non_bonded_inter[\\*,\\*] and set the \n", - "pre-calculated LJ parameters epsilon, sigma and cutoff. With shift=\"auto\",\n", - "we shift the interaction potential to the cutoff so that $U_\\mathrm{LJ}(r_\\mathrm{cutoff})=0$." + "pre-calculated WCA parameters epsilon and sigma." ] }, { @@ -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 /doc/tutorials/02-charged_system/scripts/nacl_units.py (or nacl_units_vis.py with visualization)." From 129029a882e29927ceb01d0ff23d2e1e11821430 Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Fri, 11 Oct 2019 14:40:33 +0200 Subject: [PATCH 10/10] added proper equil.time to tutorial --- testsuite/scripts/tutorials/test_02-charged_system-2.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testsuite/scripts/tutorials/test_02-charged_system-2.py b/testsuite/scripts/tutorials/test_02-charged_system-2.py index b72184b85c1..44d1e036d9f 100644 --- a/testsuite/scripts/tutorials/test_02-charged_system-2.py +++ b/testsuite/scripts/tutorials/test_02-charged_system-2.py @@ -21,7 +21,7 @@ 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