Merge branch 'python' into improve_testing

src/script_interface/particle_data/ParticleHandle.cpp

@@ -32,6 +32,8 @@
 
 #include <utils/Vector.hpp>
 
+#include <boost/format.hpp>
+
 #include <algorithm>
 #include <cmath>
 #include <memory>
@@ -458,32 +460,78 @@ Variant ParticleHandle::do_call_method(std::string const &name,
   return {};
 }
 
+#ifdef ROTATION
+static auto const contradicting_arguments_quat = std::vector<
+    std::array<std::string, 3>>{{
+    {{"dip", "dipm",
+      "Setting 'dip' is sufficient as it defines the scalar dipole moment."}},
+    {{"quat", "director",
+      "Setting 'quat' is sufficient as it defines the director."}},
+    {{"dip", "quat",
+      "Setting 'dip' would overwrite 'quat'. Set 'quat' and 'dipm' instead."}},
+    {{"dip", "director",
+      "Setting 'dip' would overwrite 'director'. Set 'director' and "
+      "'dipm' instead."}},
+}};
+#endif // ROTATION
+
 void ParticleHandle::do_construct(VariantMap const &params) {
-  m_pid = (params.count("id")) ? get_value<int>(params, "id")
-                               : get_maximal_particle_id() + 1;
+  auto const n_extra_args = params.size() - params.count("id");
+  auto const has_param = [&params](std::string const key) {
+    return params.count(key) == 1;
+  };
+  m_pid = (has_param("id")) ? get_value<int>(params, "id")
+                            : get_maximal_particle_id() + 1;
 
   // create a new particle if extra arguments were passed
-  if ((params.size() - params.count("id")) > 0) {
+  if (n_extra_args > 0) {
     if (particle_exists(m_pid)) {
       throw std::invalid_argument("Particle " + std::to_string(m_pid) +
                                   " already exists");
     }
+#ifdef ROTATION
+    // if we are not constructing a particle from a checkpoint file,
+    // check the quaternion is not accidentally set twice by the user
+    if (not has_param("__cpt_sentinel")) {
+      auto formatter =
+          boost::format("Contradicting particle attributes: '%s' and '%s'. %s");
+      for (auto const &[prop1, prop2, reason] : contradicting_arguments_quat) {
+        if (has_param(prop1) and has_param(prop2)) {
+          auto const err_msg = boost::str(formatter % prop1 % prop2 % reason);
+          throw std::invalid_argument(err_msg);
+        }
+      }
+    }
+#endif // ROTATION
 
     // create a default-constructed particle
     auto const pos = get_value<Utils::Vector3d>(params, "pos");
     mpi_make_new_particle(m_pid, pos);
 
     // set particle properties (filter out read-only and deferred properties)
     std::vector<std::string> skip = {
-        "id",        "pos",  "dipm",  "pos_folded", "lees_edwards_flag",
-        "image_box", "node", "bonds", "exclusions",
+        "pos_folded", "pos", "quat", "director",  "id",    "lees_edwards_flag",
+        "exclusions", "dip", "node", "image_box", "bonds", "__cpt_sentinel",
     };
+#ifdef ROTATION
+    // multiple parameters can potentially set the quaternion, but only one
+    // can be allowed to; these conditionals are required to handle a reload
+    // from a checkpoint file, where all properties exist (avoids accidentally
+    // overwriting the quaternion by the default-constructed dipole moment)
+    if (has_param("quat")) {
+      do_set_parameter("quat", params.at("quat"));
+    } else if (has_param("director")) {
+      do_set_parameter("director", params.at("director"));
+    } else if (has_param("dip")) {
+      do_set_parameter("dip", params.at("dip"));
+    }
+#endif // ROTATION
     for (auto const &kv : params) {
       if (std::find(skip.begin(), skip.end(), kv.first) == skip.end()) {
         do_set_parameter(kv.first, kv.second);
       }
     }
-    if (params.count("type") == 0) {
+    if (not has_param("type")) {
       do_set_parameter("type", 0);
     }
   }

src/script_interface/particle_data/ParticleList.cpp

@@ -74,6 +74,8 @@ std::string ParticleList::get_internal_state() const {
     state.params.emplace_back(
         std::pair<std::string, PackedVariant>{"exclusions", pack(exclusions)});
 #endif // EXCLUSIONS
+    state.params.emplace_back(
+        std::pair<std::string, PackedVariant>{"__cpt_sentinel", pack(None{})});
     return Utils::pack(state);
   });
 

testsuite/python/observable_profile.py

@@ -51,7 +51,7 @@ def test_density_profile(self):
         obs_bin_centers = density_profile.bin_centers()
         np_hist, np_edges = tests_common.get_histogram(
             np.copy(self.system.part.all().pos), self.kwargs, 'cartesian',
-            normed=True)
+            density=True)
         np_hist *= len(self.system.part)
         np.testing.assert_array_almost_equal(obs_data, np_hist)
         for i in range(3):

testsuite/python/observable_profileLB.py

@@ -91,7 +91,7 @@ def test_velocity_profile(self):
         obs_edges = obs.call_method("edges")
         _, np_edges = tests_common.get_histogram(
             np.zeros([1, 3]), LB_VELOCITY_PROFILE_PARAMS, 'cartesian',
-            normed=True)
+            density=True)
         for i in range(3):
             np.testing.assert_array_almost_equal(obs_edges[i], np_edges[i])
         for x in range(obs_data.shape[0]):

testsuite/python/particle.py

@@ -228,16 +228,23 @@ def test_invalid_exclusions(self):
             p1.add_exclusion(pid2)
             p1.add_exclusion(pid2)
 
-    @utx.skipIfMissingFeatures("DIPOLES")
-    def test_contradicting_properties_dip_dipm(self):
-        with self.assertRaises(ValueError):
-            self.system.part.add(pos=[0, 0, 0], dip=[1, 1, 1], dipm=1.0)
-
-    @utx.skipIfMissingFeatures(["DIPOLES", "ROTATION"])
-    def test_contradicting_properties_dip_quat(self):
-        with self.assertRaises(ValueError):
-            self.system.part.add(pos=[0, 0, 0], dip=[1, 1, 1],
-                                 quat=[1.0, 1.0, 1.0, 1.0])
+    @utx.skipIfMissingFeatures(["ROTATION"])
+    def test_contradicting_properties_quat(self):
+        invalid_combinations = [
+            {'quat': [1., 1., 1., 1.], 'director': [1., 1., 1.]},
+        ]
+        if espressomd.has_features(["DIPOLES"]):
+            invalid_combinations += [
+                {'dip': [1., 1., 1.], 'dipm': 1.},
+                {'dip': [1., 1., 1.], 'quat': [1., 1., 1., 1.]},
+                {'dip': [1., 1., 1.], 'director': [1., 1., 1.]},
+            ]
+        # check all methods that can instantiate particles
+        for kwargs in invalid_combinations:
+            for make_new_particle in (
+                    self.system.part.add, espressomd.particle_data.ParticleHandle):
+                with self.assertRaises(ValueError):
+                    make_new_particle(pos=[0., 0., 0.], **kwargs)
 
     @utx.skipIfMissingFeatures(["ROTATION"])
     def test_invalid_quat(self):

testsuite/python/save_checkpoint.py

@@ -251,25 +251,24 @@
 harmonic_bond = espressomd.interactions.HarmonicBond(r_0=0.0, k=1.0)
 system.bonded_inter.add(harmonic_bond)
 p2.add_bond((harmonic_bond, p1))
-if 'THERM.LB' not in modes:
-    # create 3 thermalized bonds that will overwrite each other's seed
-    thermalized_bond_params = dict(temp_com=0.1, temp_distance=0.2,
-                                   gamma_com=0.3, gamma_distance=0.5, r_cut=2.)
-    thermalized_bond1 = espressomd.interactions.ThermalizedBond(
-        seed=1, **thermalized_bond_params)
-    thermalized_bond2 = espressomd.interactions.ThermalizedBond(
-        seed=2, **thermalized_bond_params)
-    thermalized_bond3 = espressomd.interactions.ThermalizedBond(
-        seed=3, **thermalized_bond_params)
-    system.bonded_inter.add(thermalized_bond1)
-    p2.add_bond((thermalized_bond1, p1))
-    checkpoint.register("thermalized_bond2")
-    checkpoint.register("thermalized_bond_params")
-    if espressomd.has_features(['ELECTROSTATICS', 'MASS', 'ROTATION']):
-        dh = espressomd.drude_helpers.DrudeHelpers()
-        dh.add_drude_particle_to_core(system, harmonic_bond, thermalized_bond1,
-                                      p2, 10, 1., 4.6, 0.8, 2.)
-        checkpoint.register("dh")
+# create 3 thermalized bonds that will overwrite each other's seed
+therm_params = dict(temp_com=0.1, temp_distance=0.2, gamma_com=0.3,
+                    gamma_distance=0.5, r_cut=2.)
+therm_bond1 = espressomd.interactions.ThermalizedBond(seed=1, **therm_params)
+therm_bond2 = espressomd.interactions.ThermalizedBond(seed=2, **therm_params)
+therm_bond3 = espressomd.interactions.ThermalizedBond(seed=3, **therm_params)
+system.bonded_inter.add(therm_bond1)
+p2.add_bond((therm_bond1, p1))
+checkpoint.register("therm_bond2")
+checkpoint.register("therm_params")
+# create Drude particles
+if espressomd.has_features(['ELECTROSTATICS', 'MASS', 'ROTATION']):
+    dh = espressomd.drude_helpers.DrudeHelpers()
+    dh.add_drude_particle_to_core(
+        system=system, harmonic_bond=harmonic_bond,
+        thermalized_bond=therm_bond1, p_core=p2, type_drude=10,
+        alpha=1., mass_drude=0.6, coulomb_prefactor=0.8, thole_damping=2.)
+    checkpoint.register("dh")
 strong_harmonic_bond = espressomd.interactions.HarmonicBond(r_0=0.0, k=5e5)
 system.bonded_inter.add(strong_harmonic_bond)
 p4.add_bond((strong_harmonic_bond, p3))
@@ -283,17 +282,6 @@
     breakage_length=5., action_type="delete_bond")
 system.bond_breakage[strong_harmonic_bond._bond_id] = break_spec
 
-# h5md output
-if espressomd.has_features("H5MD"):
-    h5_units = espressomd.io.writer.h5md.UnitSystem(
-        time="ps", mass="u", length="m", charge="e")
-    h5 = espressomd.io.writer.h5md.H5md(
-        file_path=str(path_cpt_root / "test.h5"),
-        unit_system=h5_units)
-    h5.write()
-    h5.flush()
-    h5.close()
-
 checkpoint.register("system")
 checkpoint.register("acc_mean_variance")
 checkpoint.register("acc_time_series")
@@ -303,9 +291,6 @@
 checkpoint.register("ibm_triel_bond")
 checkpoint.register("break_spec")
 checkpoint.register("p_slice")
-if espressomd.has_features("H5MD"):
-    checkpoint.register("h5")
-    checkpoint.register("h5_units")
 
 # calculate forces
 system.integrator.run(0)
@@ -374,6 +359,51 @@
     lbf_cpt_path = path_cpt_root / "lb.cpt"
     lbf.save_checkpoint(str(lbf_cpt_path), lbf_cpt_mode)
 
+# set various properties
+p8 = system.part.add(id=8, pos=[2.0] * 3 + system.box_l)
+p8.lees_edwards_offset = 0.2
+p4.v = [-1., 2., -4.]
+if espressomd.has_features('MASS'):
+    p3.mass = 1.5
+if espressomd.has_features('ROTATION'):
+    p3.quat = [1., 2., 3., 4.]
+    p4.director = [3., 2., 1.]
+    p4.omega_body = [0.3, 0.5, 0.7]
+    p3.rotation = [True, False, True]
+if espressomd.has_features('EXTERNAL_FORCES'):
+    p3.fix = [False, True, False]
+    p3.ext_force = [-0.6, 0.1, 0.2]
+if espressomd.has_features(['EXTERNAL_FORCES', 'ROTATION']):
+    p3.ext_torque = [0.3, 0.5, 0.7]
+if espressomd.has_features('ROTATIONAL_INERTIA'):
+    p3.rinertia = [2., 3., 4.]
+if espressomd.has_features('THERMOSTAT_PER_PARTICLE'):
+    gamma = 2.
+    if espressomd.has_features('PARTICLE_ANISOTROPY'):
+        gamma = np.array([2., 3., 4.])
+    p4.gamma = gamma
+    if espressomd.has_features('ROTATION'):
+        p3.gamma_rot = 2. * gamma
+if espressomd.has_features('ENGINE'):
+    p3.swimming = {"f_swim": 0.03}
+if espressomd.has_features('ENGINE') and lbf_actor:
+    p4.swimming = {"v_swim": 0.02, "mode": "puller", "dipole_length": 1.}
+if espressomd.has_features('LB_ELECTROHYDRODYNAMICS') and lbf_actor:
+    p8.mu_E = [-0.1, 0.2, -0.3]
+
+# h5md output
+if espressomd.has_features("H5MD"):
+    h5_units = espressomd.io.writer.h5md.UnitSystem(
+        time="ps", mass="u", length="m", charge="e")
+    h5 = espressomd.io.writer.h5md.H5md(
+        file_path=str(path_cpt_root / "test.h5"),
+        unit_system=h5_units)
+    h5.write()
+    h5.flush()
+    h5.close()
+    checkpoint.register("h5")
+    checkpoint.register("h5_units")
+
 # save checkpoint file
 checkpoint.save(0)