Merge pull request #2433 from jngrad/fix-tutorials-samples

Fix tutorials and samples
espressomd · Jan 10, 2019 · 84491e5 · 84491e5
2 parents 070ca8f + 92ad6de
commit 84491e5
Show file tree

Hide file tree

Showing 12 changed files with 205 additions and 154 deletions.
diff --git a/doc/tutorials/02-charged_system/scripts/nacl.py b/doc/tutorials/02-charged_system/scripts/nacl.py
@@ -61,11 +61,11 @@
 
 # Place particles
 for i in range(int(n_ionpairs)):
-    system.part.add(id=len(system.part), type=types["Anion"], pos=numpy.random.random(
-        3) * box_l, q=charges["Anion"])
+    system.part.add(id=len(system.part), type=types["Anion"],
+                    pos=numpy.random.random(3) * box_l, q=charges["Anion"])
 for i in range(int(n_ionpairs)):
-    system.part.add(id=len(system.part), type=types["Cation"], pos=numpy.random.random(
-        3) * box_l, q=charges["Cation"])
+    system.part.add(id=len(system.part), type=types["Cation"],
+                    pos=numpy.random.random(3) * box_l, q=charges["Cation"])
 
 
 def combination_rule_epsilon(rule, eps1, eps2):
@@ -103,8 +103,8 @@ def combination_rule_sigma(rule, sig1, sig2):
 
 while min_dist < max_sigma:
     # Warmup Helper: Cap max. force, increase slowly for overlapping particles
-    min_dist = system.analysis.min_dist([types["Anion"], types["Cation"]], [
-        types["Anion"], types["Cation"]])
+    min_dist = system.analysis.min_dist([types["Anion"], types["Cation"]],
+                                        [types["Anion"], types["Cation"]])
     cap += min_dist
 # print min_dist, cap
     system.force_cap = cap
@@ -121,15 +121,10 @@ def combination_rule_sigma(rule, sig1, sig2):
 print("\n--->Temperature Equilibration")
 system.time = 0.0
 for i in range(int(num_steps_equilibration / 100)):
-    temp_measured = system.analysis.energy(
-    )['kinetic'] / ((3.0 / 2.0) * n_part)
-    print(
-        "t={0:.1f}, E_total={1:.2f}, E_coulomb={2:.2f}, T_cur={3:.4f}".format(system.time,
-                                                                              system.analysis.energy()[
-                                                                                  'total'],
-                                                                              system.analysis.energy()[
-                                                                                  'coulomb'],
-                                                                              temp_measured))
+    temp_measured = system.analysis.energy()['kinetic'] / ((3 / 2.0) * n_part)
+    print("t={0:.1f}, E_total={1:.2f}, E_coulomb={2:.2f}, T_cur={3:.4f}".format(
+        system.time, system.analysis.energy()['total'],
+        system.analysis.energy()['coulomb'], temp_measured))
     system.integrator.run(100)
 
 print("\n--->Integration")
@@ -138,13 +133,9 @@ def combination_rule_sigma(rule, sig1, sig2):
 for i in range(num_configs):
     temp_measured.append(system.analysis.energy()[
                          'kinetic'] / ((3.0 / 2.0) * n_part))
-    print(
-        "t={0:.1f}, E_total={1:.2f}, E_coulomb={2:.2f}, T_cur={3:.4f}".format(system.time,
-                                                                              system.analysis.energy()[
-                                                                                  'total'],
-                                                                              system.analysis.energy()[
-                                                                                  'coulomb'],
-                                                                              temp_measured[-1]))
+    print("t={0:.1f}, E_total={1:.2f}, E_coulomb={2:.2f}, T_cur={3:.4f}".format(
+        system.time, system.analysis.energy()['total'],
+        system.analysis.energy()['coulomb'], temp_measured[-1]))
     system.integrator.run(integ_steps_per_config)
 
     # Internally append particle configuration

diff --git a/doc/tutorials/02-charged_system/scripts/nacl_units.py b/doc/tutorials/02-charged_system/scripts/nacl_units.py
@@ -18,6 +18,7 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
+import espressomd
 from espressomd import electrostatics, assert_features
 import numpy
 
@@ -72,11 +73,19 @@
         for k in range(n_ppside):
             p = numpy.array([i, j, k]) * l
             if q < 0:
-                system.part.add(id=len(
-                    system.part), type=types["Cl"], pos=p, q=charges["Cl"], mass=masses["Cl"])
+                system.part.add(
+                    id=len(system.part),
+                    type=types["Cl"],
+                    pos=p,
+                    q=charges["Cl"],
+                    mass=masses["Cl"])
             else:
-                system.part.add(id=len(
-                    system.part), type=types["Na"], pos=p, q=charges["Na"], mass=masses["Na"])
+                system.part.add(
+                    id=len(system.part),
+                    type=types["Na"],
+                    pos=p,
+                    q=charges["Na"],
+                    mass=masses["Na"])
 
             q *= -1
         q *= -1
@@ -112,35 +121,25 @@ def combination_rule_sigma(rule, sig1, sig2):
 print("\n--->Tuning Electrostatics")
 # p3m = electrostatics.P3M(bjerrum_length=l_bjerrum, accuracy=1e-2,
 # mesh=[84,84,84], cao=6)
-p3m = electrostatics.P3M(bjerrum_length=l_bjerrum, accuracy=1e-2)
+p3m = electrostatics.P3M(prefactor=l_bjerrum * temp, accuracy=1e-2)
 system.actors.add(p3m)
 
 print("\n--->Temperature Equilibration")
 system.time = 0.0
 for i in range(int(num_steps_equilibration / 100)):
     energy = system.analysis.energy()
     temp_measured = energy['kinetic'] / ((3.0 / 2.0) * n_part)
-    print(
-        "t={0:.1f}, E_total={1:.2f}, E_coulomb={2:.2f}, T_cur={3:.4f}".format(system.time,
-                                                                              energy[
-                                                                              'total'],
-                                                                              energy[
-                                                                              'coulomb'],
-                                                                              temp_measured))
+    print("t={0:.1f}, E_total={1:.2f}, E_coulomb={2:.2f}, T_cur={3:.4f}".format(
+        system.time, energy['total'], energy['coulomb'], temp_measured))
     system.integrator.run(100)
 
 print("\n--->Integration")
 system.time = 0.0
 for i in range(num_configs):
     energy = system.analysis.energy()
     temp_measured = energy['kinetic'] / ((3.0 / 2.0) * n_part)
-    print(
-        "t={0:.1f}, E_total={1:.2f}, E_coulomb={2:.2f}, T_cur={3:.4f}".format(system.time,
-                                                                              energy[
-                                                                              'total'],
-                                                                              energy[
-                                                                              'coulomb'],
-                                                                              temp_measured))
+    print("t={0:.1f}, E_total={1:.2f}, E_coulomb={2:.2f}, T_cur={3:.4f}".format(
+        system.time, energy['total'], energy['coulomb'], temp_measured))
     system.integrator.run(integ_steps_per_config)
 
     # Internally append particle configuration

diff --git a/doc/tutorials/02-charged_system/scripts/nacl_units_confined.py b/doc/tutorials/02-charged_system/scripts/nacl_units_confined.py
@@ -18,11 +18,12 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
+import espressomd
 from espressomd import electrostatics, electrostatic_extensions, assert_features
 from espressomd.shapes import Wall
 import numpy
 
-assert_features(["ELECTROSTATICS", "CONSTRAINTS", "MASS", "LENNARD_JONES"])
+assert_features(["ELECTROSTATICS", "MASS", "LENNARD_JONES"])
 
 system = espressomd.System(box_l=[1.0, 1.0, 1.0])
 system.seed = system.cell_system.get_state()['n_nodes'] * [1234]
@@ -73,22 +74,30 @@
 system.thermostat.set_langevin(kT=temp, gamma=gamma)
 
 # Walls
-system.constraints.add(shape=Wall(
-    dist=0, normal=[0, 0, 1]), particle_type=types["Electrode"])
-system.constraints.add(shape=Wall(
-    dist=-box_z, normal=[0, 0, -1]), particle_type=types["Electrode"])
+system.constraints.add(shape=Wall(dist=0, normal=[0, 0, 1]),
+                       particle_type=types["Electrode"])
+system.constraints.add(shape=Wall(dist=-box_z, normal=[0, 0, -1]),
+                       particle_type=types["Electrode"])
 
 # Place particles
 for i in range(int(n_ionpairs)):
     p = numpy.random.random(3) * box_l
     p[2] += lj_sigmas["Electrode"]
-    system.part.add(id=len(system.part),
-                    type=types["Cl"], pos=p, q=charges["Cl"], mass=masses["Cl"])
+    system.part.add(
+        id=len(system.part),
+        type=types["Cl"],
+        pos=p,
+        q=charges["Cl"],
+        mass=masses["Cl"])
 for i in range(int(n_ionpairs)):
     p = numpy.random.random(3) * box_l
     p[2] += lj_sigmas["Electrode"]
-    system.part.add(id=len(system.part),
-                    type=types["Na"], pos=p, q=charges["Na"], mass=masses["Na"])
+    system.part.add(
+        id=len(system.part),
+        type=types["Na"],
+        pos=p,
+        q=charges["Na"],
+        mass=masses["Na"])
 
 # Lennard-Jones interactions parameters
 
@@ -107,7 +116,8 @@ def combination_rule_sigma(rule, sig1, sig2):
         return ValueError("No combination rule defined")
 
 
-for s in [["Cl", "Na"], ["Cl", "Cl"], ["Na", "Na"], ["Na", "Electrode"], ["Cl", "Electrode"]]:
+for s in [["Cl", "Na"], ["Cl", "Cl"], ["Na", "Na"],
+          ["Na", "Electrode"], ["Cl", "Electrode"]]:
     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]])
@@ -127,7 +137,7 @@ def combination_rule_sigma(rule, sig1, sig2):
 print("After Minimization: E_total=", energy['total'])
 
 print("\n--->Tuning Electrostatics")
-p3m = electrostatics.P3M(bjerrum_length=l_bjerrum, accuracy=1e-2)
+p3m = electrostatics.P3M(prefactor=l_bjerrum * temp, accuracy=1e-2)
 system.actors.add(p3m)
 elc = electrostatic_extensions.ELC(gap_size=elc_gap, maxPWerror=1e-3)
 system.actors.add(elc)
@@ -140,13 +150,8 @@ def combination_rule_sigma(rule, sig1, sig2):
 for i in range(int(num_steps_equilibration / 100)):
     energy = system.analysis.energy()
     temp_measured = energy['kinetic'] / ((3.0 / 2.0) * n_part)
-    print(
-        "t={0:.1f}, E_total={1:.2f}, E_coulomb={2:.2f}, T_cur={3:.4f}".format(system.time,
-                                                                              energy[
-                                                                              'total'],
-                                                                              energy[
-                                                                              'coulomb'],
-                                                                              temp_measured))
+    print("t={0:.1f}, E_total={1:.2f}, E_coulomb={2:.2f}, T_cur={3:.4f}".format(
+        system.time, energy['total'], energy['coulomb'], temp_measured))
     system.integrator.run(100)
 
 
@@ -158,15 +163,10 @@ def combination_rule_sigma(rule, sig1, sig2):
 cnt = 0
 
 for i in range(num_configs):
-    temp_measured = system.analysis.energy(
-    )['kinetic'] / ((3.0 / 2.0) * n_part)
-    print(
-        "t={0:.1f}, E_total={1:.2f}, E_coulomb={2:.2f}, T_cur={3:.4f}".format(system.time,
-                                                                              system.analysis.energy()[
-                                                                              'total'],
-                                                                              system.analysis.energy()[
-                                                                              'coulomb'],
-                                                                              temp_measured))
+    temp_measured = system.analysis.energy()['kinetic'] / ((3 / 2.0) * n_part)
+    print("t={0:.1f}, E_total={1:.2f}, E_coulomb={2:.2f}, T_cur={3:.4f}".format(
+        system.time, system.analysis.energy()['total'],
+        system.analysis.energy()['coulomb'], temp_measured))
     system.integrator.run(integ_steps_per_config)
 
     for p in system.part:

diff --git a/doc/tutorials/02-charged_system/scripts/nacl_units_confined_vis.py b/doc/tutorials/02-charged_system/scripts/nacl_units_confined_vis.py
@@ -18,14 +18,15 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
+import espressomd
 from espressomd import assert_features, electrostatics, electrostatic_extensions
 from espressomd.shapes import Wall
 from espressomd.visualization_opengl import *
 import numpy
 from threading import Thread
 from time import sleep
 
-assert_features(["ELECTROSTATICS", "CONSTRAINTS", "MASS", "LENNARD_JONES"])
+assert_features(["ELECTROSTATICS", "MASS", "LENNARD_JONES"])
 
 system = espressomd.System(box_l=[1.0, 1.0, 1.0])
 system.seed = system.cell_system.get_state()['n_nodes'] * [1234]
@@ -71,26 +72,41 @@
 system.thermostat.set_langevin(kT=temp, gamma=gamma)
 
 # Visualizer
-visualizer = openGLLive(system, camera_position=[-3 * box_l, box_l * 0.5, box_l * 0.5], camera_right=[
-                        0, 0, 1], drag_force=5 * 298, background_color=[1, 1, 1], light_pos=[30, 30, 30], ext_force_arrows_scale=[0.0001], ext_force_arrows=False)
+visualizer = openGLLive(
+    system,
+    camera_position=[-3 * box_l, box_l * 0.5, box_l * 0.5],
+    camera_right=[0, 0, 1],
+    drag_force=5 * 298,
+    background_color=[1, 1, 1],
+    light_pos=[30, 30, 30],
+    ext_force_arrows_type_scale=[0.0001],
+    ext_force_arrows=False)
 
 # Walls
-system.constraints.add(shape=Wall(
-    dist=0, normal=[0, 0, 1]), particle_type=types["Electrode"])
-system.constraints.add(shape=Wall(
-    dist=-box_z, normal=[0, 0, -1]), particle_type=types["Electrode"])
+system.constraints.add(shape=Wall(dist=0, normal=[0, 0, 1]),
+                       particle_type=types["Electrode"])
+system.constraints.add(shape=Wall(dist=-box_z, normal=[0, 0, -1]),
+                       particle_type=types["Electrode"])
 
 # Place particles
 for i in range(int(n_ionpairs)):
     p = numpy.random.random(3) * box_l
     p[2] += lj_sigmas["Electrode"]
-    system.part.add(id=len(system.part),
-                    type=types["Cl"], pos=p, q=charges["Cl"], mass=masses["Cl"])
+    system.part.add(
+        id=len(system.part),
+        type=types["Cl"],
+        pos=p,
+        q=charges["Cl"],
+        mass=masses["Cl"])
 for i in range(int(n_ionpairs)):
     p = numpy.random.random(3) * box_l
     p[2] += lj_sigmas["Electrode"]
-    system.part.add(id=len(system.part),
-                    type=types["Na"], pos=p, q=charges["Na"], mass=masses["Na"])
+    system.part.add(
+        id=len(system.part),
+        type=types["Na"],
+        pos=p,
+        q=charges["Na"],
+        mass=masses["Na"])
 
 # Lennard-Jones interactions parameters
 
@@ -109,7 +125,8 @@ def combination_rule_sigma(rule, sig1, sig2):
         return ValueError("No combination rule defined")
 
 
-for s in [["Cl", "Na"], ["Cl", "Cl"], ["Na", "Na"], ["Na", "Electrode"], ["Cl", "Electrode"]]:
+for s in [["Cl", "Na"], ["Cl", "Cl"], ["Na", "Na"],
+          ["Na", "Electrode"], ["Cl", "Electrode"]]:
     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]])
@@ -119,12 +136,12 @@ def combination_rule_sigma(rule, sig1, sig2):
     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.minimize_energy.init(
-    f_max=10, gamma=10, max_steps=2000, max_displacement=0.1)
+system.minimize_energy.init(f_max=10, gamma=10, max_steps=2000,
+                            max_displacement=0.1)
 system.minimize_energy.minimize()
 
 print("\n--->Tuning Electrostatics")
-p3m = electrostatics.P3M(bjerrum_length=l_bjerrum, accuracy=1e-2)
+p3m = electrostatics.P3M(prefactor=l_bjerrum * temp, accuracy=1e-2)
 system.actors.add(p3m)
 elc = electrostatic_extensions.ELC(gap_size=elc_gap, maxPWerror=1e-3)
 system.actors.add(elc)
@@ -147,10 +164,10 @@ def decreaseElectricField():
 
 
 # Register buttons
-visualizer.keyboardManager.registerButton(KeyboardButtonEvent(
-    'u', KeyboardFireEvent.Hold, increaseElectricField))
-visualizer.keyboardManager.registerButton(KeyboardButtonEvent(
-    'j', KeyboardFireEvent.Hold, decreaseElectricField))
+visualizer.keyboardManager.register_button(
+    KeyboardButtonEvent('u', KeyboardFireEvent.Hold, increaseElectricField))
+visualizer.keyboardManager.register_button(
+    KeyboardButtonEvent('j', KeyboardFireEvent.Hold, decreaseElectricField))
 
 
 def main():