Fix pressure tensor sign and ordering conventions

doc/sphinx/analysis.rst

@@ -523,7 +523,7 @@ documentation for all available observables in :mod:`espressomd.observables`.
 
    - :class:`~espressomd.observables.PressureTensor`: Total pressure tensor (see :ref:`Pressure Tensor`)
 
-   - :class:`~espressomd.observables.DPDStress`
+   - :class:`~espressomd.observables.DPDPressure`
 
 
 .. _Accumulators:

src/core/dpd.cpp

@@ -100,15 +100,15 @@ Utils::Vector3d dpd_pair_force(
   return force;
 }
 
-static auto dpd_viscous_stress_local(System::System &system) {
+static auto dpd_pressure_local(System::System &system) {
   auto const &box_geo = *system.box_geo;
   auto const &nonbonded_ias = *system.nonbonded_ias;
   auto &cell_structure = *system.cell_structure;
   auto &dpd = *system.thermostat->dpd;
   system.on_observable_calc();
 
-  Utils::Matrix<double, 3, 3> stress{};
-  cell_structure.non_bonded_loop([&stress, &box_geo, &nonbonded_ias,
+  Utils::Matrix<double, 3, 3> pressure{};
+  cell_structure.non_bonded_loop([&pressure, &box_geo, &nonbonded_ias,
                                   &dpd](Particle const &p1, Particle const &p2,
                                         Distance const &d) {
     auto const v21 =
@@ -130,40 +130,35 @@ static auto dpd_viscous_stress_local(System::System &system) {
      * doing only one matrix-vector multiplication */
     auto const f = P * (f_r - f_t) + f_t;
 
-    stress += tensor_product(d.vec21, f);
+    pressure += tensor_product(d.vec21, f);
   });
 
-  return stress;
-}
-
-Utils::Vector9d dpd_pressure_local(System::System &system) {
-  auto const local_stress = dpd_viscous_stress_local(system);
-  return -Utils::flatten(local_stress);
+  return pressure;
 }
 
 /**
- * @brief Viscous stress tensor of the DPD interaction.
+ * @brief Pressure tensor contribution of the DPD interaction.
  *
- * This calculates the total viscous stress contribution of the
- * DPD interaction. It contains only the dissipative contributions
- * of the interaction without noise. It's calculated as the
- * sum over all pair virials as
+ * This calculates the total contribution of the DPD interaction to
+ * the pressure tensor, i.e. the dissipative and random (noise) forces.
+ * It's calculated as the sum over all pair virials as
  * \f[
- *    \sigma^{\nu\mu} = V^{-1}\sum_i \sum_{j < i} r_{i,j}^{\nu} (- \gamma_{i,j}
- * v_{i,j})^{\mu} \f] where \f$\gamma_{i,j}\f$ is the (in general tensor valued)
- * DPD friction coefficient for particles i and j, \f$v_{i,j}\f$, \f$r_{i,j}\f$
- * are their relative velocity and distance and \f$V\f$ is the box volume.
+ *    P^{\nu\mu} = V^{-1}\sum_i \sum_{j < i} r_{i,j}^{\nu} F_{i,j}^{\mu}
+ * \f]
+ * where \f$F_{i,j}\f$ is the DPD force (dissipative plus noise) exerted
+ * by particle j on particle i, \f$r_{i,j}\f$ is their distance vector
+ * and \f$V\f$ is the box volume.
  *
- * @return Stress tensor contribution.
+ * @return Pressure tensor contribution.
  */
-Utils::Vector9d dpd_stress(System::System &system,
-                           boost::mpi::communicator const &comm) {
-  auto const local_stress = dpd_viscous_stress_local(system);
-  std::remove_const_t<decltype(local_stress)> global_stress{};
+Utils::Vector9d dpd_pressure(System::System &system,
+                             boost::mpi::communicator const &comm) {
+  auto const local_pressure = dpd_pressure_local(system);
+  std::remove_const_t<decltype(local_pressure)> global_pressure{};
 
-  boost::mpi::reduce(comm, local_stress, global_stress, std::plus<>(), 0);
+  boost::mpi::reduce(comm, local_pressure, global_pressure, std::plus<>(), 0);
 
-  return Utils::flatten(global_stress) / system.box_geo->volume();
+  return Utils::flatten(global_pressure) / system.box_geo->volume();
 }
 
 #endif // ESPRESSO_DPD

src/core/dpd.hpp

@@ -74,14 +74,17 @@ Utils::Vector3d dpd_pair_force(
     BoxGeometry const &box_geo, IA_parameters const &ia_params,
     Utils::Vector3d const &d, double dist, double dist2);
 
-Utils::Vector9d dpd_stress(System::System &system,
-                           boost::mpi::communicator const &comm);
+Utils::Vector9d dpd_pressure(System::System &system,
+                             boost::mpi::communicator const &comm);
 
-/**
- * @brief Local contribution to the pressure tensor.
- * Needs to be rescaled by the box volume.
+/** Return a random uniform 3D vector with the Philox thermostat.
+ *  Random numbers depend on
+ *  1. dpd_rng_counter (initialized by seed) which is increased on integration
+ *  2. Salt (decorrelates different counters)
+ *  3. Two particle IDs (order-independent, decorrelates particles, gets rid of
+ *     seed-per-node)
  */
-Utils::Vector9d dpd_pressure_local(System::System &system);
+Utils::Vector3d dpd_noise(DPDThermostat const &dpd, int pid1, int pid2);
 
 /** Return a random uniform 3D vector with the Philox thermostat.
  *  Random numbers depend on

src/core/electrostatics/coulomb_inline.hpp

@@ -92,7 +92,7 @@ struct ShortRangePressureKernel {
     if constexpr (traits::has_pressure<T>::value) {
       return kernel_type{
           [&actor = *ptr](double q1q2, Utils::Vector3d const &d, double dist) {
-            return Utils::tensor_product(actor.pair_force(q1q2, d, dist), d);
+            return Utils::tensor_product(d, actor.pair_force(q1q2, d, dist));
           }};
     }
     return std::nullopt;

src/core/observables/DPDStress.hpp → src/core/observables/DPDPressure.hpp

@@ -31,12 +31,12 @@
 
 namespace Observables {
 
-class DPDStress : public Observable {
+class DPDPressure : public Observable {
 public:
   std::vector<std::size_t> shape() const override { return {3, 3}; }
   std::vector<double>
   operator()(boost::mpi::communicator const &comm) const override {
-    return dpd_stress(System::get_system(), comm);
+    return dpd_pressure(System::get_system(), comm);
   }
 };
 

src/core/pressure_cabana.hpp

@@ -208,7 +208,7 @@ struct PressureKernel {
           auto const f_d = P * (f_r - f_t) + f_t;
           auto const s = Utils::flatten(Utils::tensor_product(d, f_d));
           for (std::size_t k = 0; k < 9; ++k)
-            local_pressure(tid, layout.tensor_offset(layout.dpd_idx(), k)) -=
+            local_pressure(tid, layout.tensor_offset(layout.dpd_idx(), k)) +=
                 s[k];
         }
 #endif

src/core/pressure_inline.hpp

@@ -60,7 +60,7 @@ calc_bonded_virial_pressure_tensor(
   );
   std::optional<Utils::Matrix<double, 3, 3>> pressure{std::nullopt};
   if (pair_force) {
-    pressure = Utils::tensor_product(*pair_force, dx);
+    pressure = Utils::tensor_product(dx, *pair_force);
   }
   return pressure;
 }
@@ -100,8 +100,8 @@ calc_bonded_three_body_pressure_tensor(Bonded_IA_Parameters const &iaparams,
     if (result) {
       Utils::Vector3d force2, force3;
       std::tie(std::ignore, force2, force3) = result.value();
-      return Utils::tensor_product(force2, dx21) +
-             Utils::tensor_product(force3, dx31);
+      return Utils::tensor_product(dx21, force2) +
+             Utils::tensor_product(dx31, force3);
     }
   } else {
     runtimeWarningMsg() << "Unsupported bond type " +
@@ -134,9 +134,9 @@ calc_bonded_four_body_pressure_tensor(Bonded_IA_Parameters const &iaparams,
       Utils::Vector3d force2, force3, force4;
       std::tie(std::ignore, force2, force3, force4) = result.value();
 
-      return -Utils::tensor_product(force2, v12) +
-             Utils::tensor_product(force3, v23) +
-             Utils::tensor_product(force4, v23 + v34);
+      return -Utils::tensor_product(v12, force2) +
+             Utils::tensor_product(v23, force3) +
+             Utils::tensor_product(v23 + v34, force4);
     }
   } else {
     runtimeWarningMsg() << "Unsupported bond type " +

src/core/virtual_sites/relative.cpp

@@ -99,7 +99,7 @@ static auto constraint_stress(Particle const &p_ref, Particle const &p_vs) {
   /* The constraint force is minus the force on the particle, make it force
    * free. The counter force is translated by the connection vector to the
    * real particle, hence the virial stress is */
-  return tensor_product(-p_vs.force(), connection_vector(p_ref, p_vs));
+  return tensor_product(connection_vector(p_ref, p_vs), -p_vs.force());
 }
 
 static bool is_vs_relative_trans(Particle const &p) {

src/python/espressomd/analyze.py

@@ -552,19 +552,19 @@ def particle_bond_energy(self, particle, bond):
         return self.call_method("particle_bond_energy", pid=particle.id,
                                 bond_id=interaction._bond_id, partners=partners)
 
-    def dpd_stress(self):
+    def dpd_pressure(self):
         """
-        Calculate the total viscous stress contribution of the DPD interaction.
-        It contains only the dissipative contributions, without noise:
-        :math:`\\sigma^{\\nu\\mu} = V^{-1}\\sum_i \\sum_{j < i}
-        r_{i,j}^{\\nu} (- \\gamma_{i,j} v_{i,j})^{\\mu}`
-        where :math:`\\gamma_{i,j}` is the (in general tensor-valued) DPD
-        friction coefficient for particles i and j, :math:`v_{i,j}`,
-        :math:`r_{i,j}` are their relative velocity and distance,
-        and :math:`V` is the box volume.
+        Calculate the DPD interaction contribution to the pressure tensor.
+        It contains the dissipative and random (noise) contributions:
+        :math:`P^{\\nu\\mu} = V^{-1}\\sum_i \\sum_{j < i}
+        r_{i,j}^{\\nu} F_{i,j}^{\\mu}`
+        where :math:`F_{i,j}` is the DPD force (dissipative plus noise)
+        exerted by particle j on particle i, :math:`r_{i,j}` is their
+        distance vector, and :math:`V` is the box volume. This equals the
+        ``"dpd"`` contribution of :meth:`pressure_tensor`.
         """
         assert_features("DPD")
-        return np.reshape(self.call_method("dpd_stress"), (3, 3))
+        return np.reshape(self.call_method("dpd_pressure"), (3, 3))
 
     def gyration_tensor(self, p_type=None):
         """

src/python/espressomd/observables.py

@@ -829,10 +829,10 @@ class PressureTensor(Observable):
 
 
 @script_interface_register
-class DPDStress(Observable):
+class DPDPressure(Observable):
 
     """Calculates the non-equilibrium contribution of the DPD interaction
-    to the stress tensor.
+    to the pressure tensor.
 
     Parameters
     ----------
@@ -848,7 +848,7 @@ class DPDStress(Observable):
         (3, 3) :obj:`ndarray` of :obj:`float`
 
     """
-    _so_name = "Observables::DPDStress"
+    _so_name = "Observables::DPDPressure"
 
 
 @script_interface_register

src/script_interface/analysis/Analysis.cpp

@@ -149,8 +149,8 @@ Variant Analysis::do_call_method(std::string const &name,
     return make_unordered_map_of_variants(dict);
   }
 #ifdef ESPRESSO_DPD
-  if (name == "dpd_stress") {
-    auto const result = dpd_stress(get_system(), context()->get_comm());
+  if (name == "dpd_pressure") {
+    auto const result = dpd_pressure(get_system(), context()->get_comm());
     return result.as_vector();
   }
 #endif // ESPRESSO_DPD

src/script_interface/observables/ParamlessObservable.hpp

@@ -26,7 +26,7 @@
 #include "script_interface/ScriptInterface.hpp"
 
 #include "Observable.hpp"
-#include "core/observables/DPDStress.hpp"
+#include "core/observables/DPDPressure.hpp"
 #include "core/observables/EnergyObservable.hpp"
 #include "core/observables/LBFluidPressureTensor.hpp"
 #include "core/observables/Observable.hpp"
@@ -64,7 +64,7 @@ NEW_PARAMLESS_OBSERVABLE(Pressure)
 NEW_PARAMLESS_OBSERVABLE(PressureTensor)
 NEW_PARAMLESS_OBSERVABLE(LBFluidPressureTensor)
 #ifdef ESPRESSO_DPD
-NEW_PARAMLESS_OBSERVABLE(DPDStress)
+NEW_PARAMLESS_OBSERVABLE(DPDPressure)
 #endif
 
 } /* namespace Observables */

src/script_interface/observables/initialize.cpp

@@ -151,7 +151,7 @@ void initialize(Utils::Factory<ObjectHandle> *om) {
   REGISTER_PAIRWISE_DISTANCES(PairwiseDistances);
 
 #ifdef ESPRESSO_DPD
-  REGISTER(DPDStress)
+  REGISTER(DPDPressure)
 #endif
   REGISTER(LBFluidPressureTensor);
   REGISTER_CYLPID_PROFILE_OBS(

testsuite/python/dpd.py

@@ -398,7 +398,7 @@ def test_zz_constraint(self):
         np.testing.assert_array_almost_equal(np.copy(p3.f), 0.)
         system.constraints.remove(constraint)
 
-    def test_dpd_stress(self):
+    def test_dpd_pressure(self):
 
         def calc_omega(dist):
             return (1. / dist - 1. / r_cut) ** 2.0
@@ -435,11 +435,11 @@ def diss_force_2(dist, vel_diff):
 
             return f
 
-        def calc_stress(dist, vel_diff):
+        def calc_pressure(dist, vel_diff):
             force_pair = diss_force_1(dist, vel_diff) +\
                 diss_force_2(dist, vel_diff)
-            stress_pair = np.outer(dist, force_pair)
-            return stress_pair
+            pressure_pair = np.outer(dist, force_pair)
+            return pressure_pair
 
         n_part = 200
         r_cut = 1.0
@@ -453,44 +453,43 @@ def calc_stress(dist, vel_diff):
             trans_weight_function=1, trans_gamma=gamma / 2.0, trans_r_cut=r_cut)
 
         pos = system.box_l * np.random.random((n_part, 3))
-        partcls = system.part.add(pos=pos)
+        vel = np.random.random((n_part, 3))
+        system.part.add(pos=pos, v=vel)
         system.integrator.run(10)
 
         # test non-thermalized case
         system.thermostat.set_dpd(kT=0, seed=3)
-        # run 1 integration step to get velocities
-        partcls.v = np.zeros((n_part, 3))
-        system.integrator.run(steps=1)
 
         pairs = system.part.pairs()
 
-        stress = np.zeros([3, 3])
+        pressure_ref = np.zeros([3, 3])
 
         for pair in pairs:
             dist = system.distance_vec(pair[0], pair[1])
             if np.linalg.norm(dist) < r_cut:
                 vel_diff = pair[1].v - pair[0].v
-                stress += calc_stress(dist, vel_diff)
+                pressure_ref += calc_pressure(dist, vel_diff)
 
-        stress /= system.volume()
+        pressure_ref /= system.volume()
 
-        dpd_stress = system.analysis.dpd_stress()
+        dpd_pressure = system.analysis.dpd_pressure()
 
-        dpd_obs = espressomd.observables.DPDStress()
-        obs_stress = dpd_obs.calculate()
+        dpd_obs = espressomd.observables.DPDPressure()
+        obs_pressure = dpd_obs.calculate()
         pressure = system.analysis.pressure_tensor()["dpd"]
 
-        np.testing.assert_array_almost_equal(np.copy(dpd_stress), stress)
-        np.testing.assert_array_almost_equal(np.copy(obs_stress), stress)
-        np.testing.assert_array_almost_equal(np.copy(pressure), -stress)
+        np.testing.assert_array_almost_equal(
+            np.copy(dpd_pressure), pressure_ref)
+        np.testing.assert_array_almost_equal(
+            np.copy(obs_pressure), pressure_ref)
+        np.testing.assert_array_almost_equal(np.copy(pressure), pressure_ref)
 
     @utx.skipIfMissingFeatures("EXTERNAL_FORCES")
-    def test_dpd_stress_noise_statistics(self):
+    def test_dpd_pressure_noise_statistics(self):
         """
-        Thermalized DPD stress: for a fixed pair with zero relative velocity
-        the stress is pure noise. Check its mean (=0) and per-component variance
-        against the analytic fluctuation-dissipation result. A generic off-axis
-        separation makes all 9 stress components non-zero and distinct.
+        Thermalized DPD pressure: for a fixed pair with zero relative velocity
+        the pressure is pure noise. Check its mean (=0) and per-component
+        variance against the analytic fluctuation-dissipation result.
         """
         system = self.system
         kT, gamma_r, gamma_t, r_cut = 2.0, 1.5, 0.7, 1.5
@@ -508,9 +507,9 @@ def test_dpd_stress_noise_statistics(self):
         system.part.add(pos=pos0, v=[0., 0., 0.], fix=[True, True, True])
         system.part.add(pos=pos0 + d, v=[0., 0., 0.], fix=[True, True, True])
 
-        # analytic mean (0) and variance of the single-pair noise stress:
-        #   stress_ij = d_i * (M @ noise)_j / V,   M = a*P + b*I
-        #   Var(stress_ij) = (d_i^2 / V^2) * s2 * ((A^2 - B^2)*dhat_j^2 + B^2)
+        # analytic mean (0) and variance of the single-pair noise pressure:
+        #   pressure_ij = d_i * (M @ noise)_j / V,   M = a*P + b*I
+        #   Var(pressure_ij) = (d_i^2 / V^2) * s2 * ((A^2 - B^2)*dhat_j^2 + B^2)
         # with omega = 1 (weight_function=0), A/B the radial/trans amplitudes,
         # and noise components iid with variance s2 = 1/12.
         V = system.volume()
@@ -525,7 +524,7 @@ def test_dpd_stress_noise_statistics(self):
         samples = np.empty((N, 3, 3))
         for n in range(N):
             system.integrator.run(1)
-            samples[n] = system.analysis.dpd_stress()
+            samples[n] = system.analysis.dpd_pressure()
 
         mean = samples.mean(axis=0)
         var = samples.var(axis=0)