Rewrite particle serialization mechanism

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/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

@@ -249,25 +249,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))
@@ -281,17 +280,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")
@@ -301,9 +289,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)
@@ -372,6 +357,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)
 

testsuite/python/test_checkpoint.py

@@ -157,12 +157,87 @@ def test_lees_edwards(self):
         self.assertAlmostEqual(protocol.time_0, 0.2, delta=1e-10)
         self.assertAlmostEqual(protocol.shear_velocity, 1.2, delta=1e-10)
 
-    def test_part(self):
-        p1, p2 = system.part.by_ids([0, 1])
+    def test_particle_properties(self):
+        p1, p2, p3, p4, p8 = system.part.by_ids([0, 1, 3, 4, 8])
         np.testing.assert_allclose(np.copy(p1.pos), np.array([1.0, 2.0, 3.0]))
         np.testing.assert_allclose(np.copy(p2.pos), np.array([1.0, 1.0, 2.0]))
         np.testing.assert_allclose(np.copy(p1.f), particle_force0)
         np.testing.assert_allclose(np.copy(p2.f), particle_force1)
+        self.assertEqual(p1.type, 0)
+        self.assertEqual(p2.type, 0)
+        self.assertEqual(p3.type, 1)
+        self.assertEqual(p4.type, 1)
+        np.testing.assert_allclose(np.copy(p3.v), [0., 0., 0.])
+        np.testing.assert_allclose(np.copy(p4.v), [-1., 2., -4.])
+        np.testing.assert_allclose(p8.lees_edwards_offset, 0.2)
+        np.testing.assert_equal(np.copy(p1.image_box), [0, 0, 0])
+        np.testing.assert_equal(np.copy(p8.image_box),
+                                np.copy(system.periodicity).astype(int))
+        self.assertEqual(p8.lees_edwards_flag, 0)
+        if espressomd.has_features('MASS'):
+            np.testing.assert_allclose(p3.mass, 1.5)
+            np.testing.assert_allclose(p4.mass, 1.)
+        if espressomd.has_features('ROTATIONAL_INERTIA'):
+            np.testing.assert_allclose(np.copy(p3.rinertia), [2., 3., 4.])
+            np.testing.assert_allclose(np.copy(p4.rinertia), [1., 1., 1.])
+        if espressomd.has_features('ELECTROSTATICS'):
+            np.testing.assert_allclose(p1.q, 1.)
+            if espressomd.has_features(['MASS', 'ROTATION']):
+                # check Drude particles
+                p5 = system.part.by_id(5)
+                np.testing.assert_allclose(p2.q, +0.118, atol=1e-3)
+                np.testing.assert_allclose(p5.q, -1.118, atol=1e-3)
+                np.testing.assert_allclose(p2.mass, 0.4)
+                np.testing.assert_allclose(p5.mass, 0.6)
+            else:
+                np.testing.assert_allclose(p2.q, -1.)
+                np.testing.assert_allclose(p2.mass, 1.)
+            np.testing.assert_allclose(p3.q, 0.)
+            np.testing.assert_allclose(p4.q, 0.)
+        if espressomd.has_features('DIPOLES'):
+            np.testing.assert_allclose(np.copy(p1.dip), [1.3, 2.1, -6.])
+        if espressomd.has_features('ROTATION'):
+            np.testing.assert_allclose(np.copy(p3.quat), [1., 2., 3., 4.])
+            np.testing.assert_allclose(np.copy(p4.director) * np.sqrt(14.),
+                                       [3., 2., 1.])
+            np.testing.assert_allclose(np.copy(p3.omega_body), [0., 0., 0.])
+            np.testing.assert_allclose(np.copy(p4.omega_body), [0.3, 0.5, 0.7])
+            np.testing.assert_equal(np.copy(p3.rotation), [True, False, True])
+        if espressomd.has_features('EXTERNAL_FORCES'):
+            np.testing.assert_equal(np.copy(p3.fix), [False, True, False])
+            np.testing.assert_allclose(np.copy(p3.ext_force), [-0.6, 0.1, 0.2])
+        if espressomd.has_features(['EXTERNAL_FORCES', 'ROTATION']):
+            np.testing.assert_allclose(np.copy(p3.ext_torque), [0.3, 0.5, 0.7])
+        if espressomd.has_features('ROTATIONAL_INERTIA'):
+            np.testing.assert_allclose(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.])
+            np.testing.assert_allclose(p4.gamma, gamma)
+            if espressomd.has_features('ROTATION'):
+                np.testing.assert_allclose(p3.gamma_rot, 2. * gamma)
+        if espressomd.has_features('ENGINE'):
+            self.assertEqual(p3.swimming, {"f_swim": 0.03, "mode": "N/A",
+                                           "v_swim": 0., "dipole_length": 0.})
+        if espressomd.has_features('ENGINE') and has_lb_mode:
+            self.assertEqual(p4.swimming, {"v_swim": 0.02, "mode": "puller",
+                                           "f_swim": 0., "dipole_length": 1.})
+        if espressomd.has_features('LB_ELECTROHYDRODYNAMICS') and has_lb_mode:
+            np.testing.assert_allclose(np.copy(p8.mu_E), [-0.1, 0.2, -0.3])
+        if espressomd.has_features('VIRTUAL_SITES_RELATIVE'):
+            from scipy.spatial.transform import Rotation as R
+            q_ind = ([1, 2, 3, 0],)  # convert from scalar-first to scalar-last
+            vs_id, vs_dist, vs_quat = p2.vs_relative
+            d = p2.pos - p1.pos
+            theta = np.arccos(d[2] / np.linalg.norm(d))
+            assert abs(theta - 3. * np.pi / 4.) < 1e-8
+            q = np.array([0., 0., np.sin(theta / 2.), -np.cos(theta / 2.)])
+            r = R.from_quat(p1.quat[q_ind]) * R.from_quat(vs_quat[q_ind])
+            self.assertEqual(vs_id, p1.id)
+            np.testing.assert_allclose(vs_dist, np.sqrt(2.))
+            np.testing.assert_allclose(q[q_ind], r.as_quat(), atol=1e-10)
+            np.testing.assert_allclose(np.copy(p2.vs_quat), [1., 0., 0., 0.])
 
     def test_part_slice(self):
         np.testing.assert_allclose(np.copy(p_slice.id), [4, 1])
@@ -361,12 +436,11 @@ def test_bonded_inter(self):
         reference = {'r_0': 0.0, 'k': 1.0, 'r_cut': 0.0}
         self.assertEqual(partcl_1.bonds[0][0].params, reference)
         self.assertEqual(system.bonded_inter[0].params, reference)
-        if 'THERM.LB' not in modes:
-            # all thermalized bonds should be identical
-            reference = {**thermalized_bond_params, 'seed': 3}
-            self.assertEqual(partcl_1.bonds[1][0].params, reference)
-            self.assertEqual(system.bonded_inter[1].params, reference)
-            self.assertEqual(thermalized_bond2.params, reference)
+        # all thermalized bonds should be identical
+        reference = {**therm_params, 'seed': 3}
+        self.assertEqual(partcl_1.bonds[1][0].params, reference)
+        self.assertEqual(system.bonded_inter[1].params, reference)
+        self.assertEqual(therm_bond2.params, reference)
         # immersed boundary bonds
         self.assertEqual(
             ibm_volcons_bond.params, {'softID': 15, 'kappaV': 0.01})
@@ -394,7 +468,6 @@ def test_bond_breakage_specs(self):
             delta=1e-10)
         self.assertEqual(break_spec.action_type, cpt_spec.action_type)
 
-    @ut.skipIf('THERM.LB' in modes, 'LB thermostat in modes')
     @utx.skipIfMissingFeatures(['ELECTROSTATICS', 'MASS', 'ROTATION'])
     def test_drude_helpers(self):
         drude_type = 10
@@ -467,9 +540,8 @@ def test_h5md(self):
             def predicate(cur, key):
                 np.testing.assert_allclose(cur[key][0], cur[key][1],
                                            err_msg=f"mismatch for '{key}'")
-            # note: cannot compare forces since they are NaN in frame #1
-            for key in ('position', 'image', 'velocity',
-                        'id', 'species', 'mass', 'charge'):
+            for key in ('id', 'position', 'image', 'velocity',
+                        'species', 'mass', 'charge', 'force'):
                 predicate(cur, f'particles/atoms/{key}/value')
             for key in ('offset', 'direction', 'normal'):
                 predicate(cur, f'particles/atoms/lees_edwards/{key}/value')