Rewrite particle serialization mechanism (espressomd#4637)

Fixes espressomd#4633

Release 4.2.0 introduced a regression that causes checkpoint files to overwrite the particle quaternion/director by a unit vector pointing along the z direction, when the `DIPOLES` feature is part of the myconfig file. This leads to incorrect trajectories when reloading a simulation from a checkpoint file, if the particle director plays a role in the simulation (ex: relative virtual sites, Gay-Berne potential, anisotropic particles, active particles, etc.). Since the default myconfig file contains `DIPOLES`, most ESPResSo users are affected.

Description of changes:
- write new checkpointing logic to avoid overwriting the particle director with the particle dipole moment
- add checks to verify particle properties are correctly reloaded from a checkpoint file
- fix regressions in the checkpointing tests

src/python/espressomd/particle_data.pyx

@@ -1434,7 +1434,13 @@ cdef class ParticleHandle:
         if "id" in new_properties:
             raise Exception("Cannot change particle id.")
 
+        for k in ("quat", "director", "dip"):
+            if k in new_properties:
+                setattr(self, k, new_properties[k])
+                break
         for k, v in new_properties.items():
+            if k in ("quat", "director", "dip"):
+                continue
             setattr(self, k, v)
 
     IF ROTATION:
@@ -1795,13 +1801,18 @@ cdef class ParticleList:
         # Prevent setting of contradicting attributes
         IF DIPOLES:
             if 'dip' in p_dict and 'dipm' in p_dict:
-                raise ValueError("Contradicting attributes: dip and dipm. Setting \
-dip is sufficient as the length of the vector defines the scalar dipole moment.")
+                raise ValueError("Contradicting attributes: 'dip' and 'dipm'. Setting \
+'dip' is sufficient as the length of the vector defines the scalar dipole moment.")
             IF ROTATION:
-                if 'dip' in p_dict and 'quat' in p_dict:
-                    raise ValueError("Contradicting attributes: dip and quat. \
-Setting dip overwrites the rotation of the particle around the dipole axis. \
-Set quat and scalar dipole moment (dipm) instead.")
+                for key in ('quat', 'director'):
+                    if 'dip' in p_dict and key in p_dict:
+                        raise ValueError(f"Contradicting attributes: 'dip' and '{key}'. \
+Setting 'dip' overwrites the rotation of the particle around the dipole axis. \
+Set '{key}' and 'dipm' instead.")
+        IF ROTATION:
+            if 'director' in p_dict and 'quat' in p_dict:
+                raise ValueError("Contradicting attributes: 'director' and 'quat'. \
+Setting 'quat' is sufficient as it defines the director.")
 
         # The ParticleList can not be used yet, as the particle
         # doesn't yet exist. Hence, the setting of position has to be

testsuite/python/particle.py

@@ -229,16 +229,20 @@ 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.]},
+            ]
+        for kwargs in invalid_combinations:
+            with self.assertRaises(ValueError):
+                self.system.part.add(pos=[0., 0., 0.], **kwargs)
 
     @utx.skipIfMissingFeatures(["ROTATION"])
     def test_invalid_quat(self):

testsuite/python/save_checkpoint.py

@@ -230,22 +230,35 @@
         epsilon=1.2, sigma=1.7, cutoff=2.0, shift=0.1)
     system.non_bonded_inter[1, 17].lennard_jones.set_params(
         epsilon=1.2e5, sigma=1.7, cutoff=2.0, shift=0.1)
+if espressomd.has_features(['DPD']):
+    dpd_params = {"weight_function": 1, "gamma": 2., "trans_r_cut": 2., "k": 2.,
+                  "trans_weight_function": 0, "trans_gamma": 1., "r_cut": 2.}
+    dpd_ia = espressomd.interactions.DPDInteraction(**dpd_params)
+    checkpoint.register("dpd_ia")
+    checkpoint.register("dpd_params")
 
 # bonded interactions
 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:
-    thermalized_bond = espressomd.interactions.ThermalizedBond(
-        temp_com=0.0, gamma_com=0.0, temp_distance=0.2, gamma_distance=0.5,
-        r_cut=2, seed=51)
-    system.bonded_inter.add(thermalized_bond)
-    p2.add_bond((thermalized_bond, p1))
-    if espressomd.has_features(['ELECTROSTATICS', 'MASS', 'ROTATION']):
-        dh = espressomd.drude_helpers.DrudeHelpers()
-        dh.add_drude_particle_to_core(system, harmonic_bond, thermalized_bond,
-                                      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))
@@ -259,17 +272,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")
@@ -278,9 +280,6 @@
 checkpoint.register("ibm_tribend_bond")
 checkpoint.register("ibm_triel_bond")
 checkpoint.register("break_spec")
-if espressomd.has_features("H5MD"):
-    checkpoint.register("h5")
-    checkpoint.register("h5_units")
 
 # calculate forces
 system.integrator.run(0)
@@ -349,6 +348,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

@@ -154,12 +154,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_bonded_interactions_serialization(self):
         '''
@@ -335,24 +410,24 @@ def test_non_bonded_inter(self):
         self.assertEqual(params1, reference1)
         self.assertEqual(params2, reference2)
 
+    @utx.skipIfMissingFeatures('DPD')
+    def test_non_bonded_inter_dpd(self):
+        self.assertEqual(dpd_ia.get_params(), dpd_params)
+
     def test_bonded_inter(self):
         # check the ObjectHandle was correctly initialized (including MPI)
         bond_ids = system.bonded_inter.call_method('get_bond_ids')
         self.assertEqual(len(bond_ids), len(system.bonded_inter))
         # check bonded interactions
         partcl_1 = system.part.by_id(1)
-        state = partcl_1.bonds[0][0].params
         reference = {'r_0': 0.0, 'k': 1.0, 'r_cut': 0.0}
-        self.assertEqual(state, reference)
-        state = partcl_1.bonds[0][0].params
-        self.assertEqual(state, reference)
-        if 'THERM.LB' not in modes:
-            state = partcl_1.bonds[1][0].params
-            reference = {'temp_com': 0., 'gamma_com': 0., 'temp_distance': 0.2,
-                         'gamma_distance': 0.5, 'r_cut': 2.0, 'seed': 51}
-            self.assertEqual(state, reference)
-            state = partcl_1.bonds[1][0].params
-            self.assertEqual(state, reference)
+        self.assertEqual(partcl_1.bonds[0][0].params, reference)
+        self.assertEqual(system.bonded_inter[0].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})
@@ -380,7 +455,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
@@ -453,9 +527,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')