Implement setter for particle director. (#4053)

Fixes #1470

src/core/particle_data.cpp

@@ -733,13 +733,13 @@ void set_particle_rotational_inertia(int part, double *rinertia) {
 #else
 constexpr Utils::Vector3d ParticleProperties::rinertia;
 #endif
+
 #ifdef ROTATION
 void set_particle_rotation(int part, int rot) {
   mpi_update_particle_property<uint8_t, &ParticleProperties::rotation>(part,
                                                                        rot);
 }
-#endif
-#ifdef ROTATION
+
 void rotate_particle(int part, const Utils::Vector3d &axis, double angle) {
   mpi_send_update_message(part, UpdateOrientation{axis, angle});
 }
@@ -759,7 +759,6 @@ void set_particle_dip(int part, double const *const dip) {
   set_particle_dipm(part, dipm);
   set_particle_quat(part, quat.data());
 }
-
 #endif
 
 #ifdef VIRTUAL_SITES
@@ -844,6 +843,12 @@ void set_particle_quat(int part, double *quat) {
       part, Utils::Quaternion<double>{quat[0], quat[1], quat[2], quat[3]});
 }
 
+void set_particle_director(int part, const Utils::Vector3d &director) {
+  Utils::Quaternion<double> quat =
+      convert_director_to_quaternion(director.normalized());
+  set_particle_quat(part, quat.data());
+}
+
 void set_particle_omega_lab(int part, const Utils::Vector3d &omega_lab) {
   auto const &particle = get_particle_data(part);
 
@@ -1200,7 +1205,6 @@ void pointer_to_omega_body(Particle const *p, double const *&res) {
 void pointer_to_quat(Particle const *p, double const *&res) {
   res = p->r.quat.data();
 }
-
 #endif
 
 void pointer_to_q(Particle const *p, double const *&res) { res = &(p->p.q); }

src/core/particle_data.hpp

@@ -194,6 +194,13 @@ void set_particle_mol_id(int part, int mid);
  */
 void set_particle_quat(int part, double *quat);
 
+/** Call only on the master node: set particle orientation using director.
+ *  The particle director defines the z-axis in the body-fixed frame.
+ *  @param part the particle.
+ *  @param director its new director vector (will be normalized if necessary)
+ */
+void set_particle_director(int part, const Utils::Vector3d &director);
+
 /** Call only on the master node: set particle angular velocity from lab frame.
  *  @param part the particle.
  *  @param omega_lab its new angular velocity.
@@ -211,7 +218,6 @@ void set_particle_omega_body(int part, const Utils::Vector3d &omega);
  *  @param torque_lab its new torque.
  */
 void set_particle_torque_lab(int part, const Utils::Vector3d &torque_lab);
-
 #endif
 
 #ifdef DIPOLES
@@ -366,7 +372,6 @@ void pointer_to_omega_body(Particle const *p, double const *&res);
 inline Utils::Vector3d get_torque_body(const Particle &p) { return p.f.torque; }
 
 void pointer_to_quat(Particle const *p, double const *&res);
-
 #endif
 
 void pointer_to_q(Particle const *p, double const *&res);

src/python/espressomd/particle_data.pxd

@@ -113,6 +113,7 @@ cdef extern from "particle_data.hpp":
     IF ROTATION:
         void set_particle_quat(int part, double quat[4])
         void pointer_to_quat(const particle * p, const double * & res)
+        void set_particle_director(int part, Vector3d director)
         void set_particle_omega_lab(int part, Vector3d omega)
         void set_particle_omega_body(int part, Vector3d omega)
         void set_particle_torque_lab(int part, Vector3d torque)

src/python/espressomd/particle_data.pyx

@@ -411,20 +411,35 @@ cdef class ParticleHandle:
                 pointer_to_quat(self.particle_data, x)
                 return array_locked([x[0], x[1], x[2], x[3]])
 
-        # Director (z-axis in body fixed frame)
         property director:
             """
-            Director.
+            The particle director.
+
+            The ``director`` defines the the z-axis in the body-fixed frame.
+            If particle rotations happen, the director, i.e., the body-fixed
+            coordinate system co-rotates. Properties such as the angular
+            velocity :attr:`espressomd.particle_data.ParticleHandle.omega_body`
+            are evaluated in this body-fixed coordinate system.
+            When using particle dipoles, the dipole moment is co-aligned with
+            the particle director. Setting the director thus modifies the
+            dipole moment orientation (:attr:`espressomd.particle_data.ParticleHandle.dip`)
+            and vice verca.
+            See also :ref:`Rotational degrees of freedom and particle anisotropy`.
+
+            director : (3,) array_like of :obj:`float`
 
             .. note::
-               Setting the director is not implemented.
                This needs the feature ``ROTATION``.
 
             """
 
-            def __set__(self, _q):
-                raise AttributeError(
-                    "Setting the director is not implemented in the C++-core of ESPResSo.")
+            def __set__(self, _d):
+                cdef Vector3d myd
+                check_type_or_throw_except(
+                    _d, 3, float, "Particle director has to be 3 floats.")
+                for i in range(3):
+                    myd[i] = _d[i]
+                set_particle_director(self._id, myd)
 
             def __get__(self):
                 self.update_particle_data()

testsuite/python/particle.py

@@ -120,6 +120,21 @@ def func(self):
         test_quat = generateTestForVectorProperty(
             "quat", np.array([0.5, 0.5, 0.5, 0.5]))
 
+        def test_director(self):
+            """
+            Test `director`. When set, it should get normalized.
+
+            """
+            sample_vector = np.array([0.5, -0.4, 1.3])
+            sample_vector_normalized = sample_vector / \
+                np.linalg.norm(sample_vector)
+
+            setattr(self.system.part[self.pid], "director", sample_vector)
+            np.testing.assert_allclose(
+                np.array(getattr(self.system.part[self.pid], "director")),
+                sample_vector_normalized
+            )
+
         if espressomd.has_features(["LANGEVIN_PER_PARTICLE"]):
             if espressomd.has_features(["PARTICLE_ANISOTROPY"]):
                 test_gamma = generateTestForVectorProperty(
@@ -147,8 +162,6 @@ def test_gamma_rot_single(self):
             else:
                 test_gamma_rot = generateTestForScalarProperty(
                     "gamma_rot", 14.23)
-    # test_director=generateTestForVectorProperty("director",
-    # np.array([0.5,0.4,0.3]))
 
     if espressomd.has_features(["ELECTROSTATICS"]):
         test_charge = generateTestForScalarProperty("q", -19.7)