Add noise term in dpd pressure calculation

src/core/dpd.cpp

@@ -51,7 +51,7 @@
  *  3. Two particle IDs (order-independent, decorrelates particles, gets rid of
  *     seed-per-node)
  */
-static Utils::Vector3d dpd_noise(DPDThermostat const &dpd, int pid1, int pid2) {
+Utils::Vector3d dpd_noise(DPDThermostat const &dpd, int pid1, int pid2) {
   auto const pref = (pid1 < pid2) ? 1.0 : -1.0;
   return pref * Random::noise_uniform<RNGSalt::SALT_DPD>(
                     dpd.rng_counter(), dpd.rng_seed(),
@@ -104,20 +104,25 @@ static auto dpd_viscous_stress_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](
-                                     Particle const &p1, Particle const &p2,
-                                     Distance const &d) {
+  cell_structure.non_bonded_loop([&stress, &box_geo, &nonbonded_ias,
+                                  &dpd](Particle const &p1, Particle const &p2,
+                                        Distance const &d) {
     auto const v21 =
         box_geo.velocity_difference(p1.pos(), p2.pos(), p1.v(), p2.v());
 
     auto const &ia_params = nonbonded_ias.get_ia_param(p1.type(), p2.type());
+    auto const noise_vec =
+        (ia_params.dpd.radial.pref > 0.0 || ia_params.dpd.trans.pref > 0.0)
+            ? dpd_noise(dpd, p1.id(), p2.id())
+            : Utils::Vector3d{};
     auto const dist = std::sqrt(d.dist2);
 
-    auto const f_r = dpd_pair_force(ia_params.dpd.radial, v21, dist, {});
-    auto const f_t = dpd_pair_force(ia_params.dpd.trans, v21, dist, {});
+    auto const f_r = dpd_pair_force(ia_params.dpd.radial, v21, dist, noise_vec);
+    auto const f_t = dpd_pair_force(ia_params.dpd.trans, v21, dist, noise_vec);
 
     /* Projection operator to radial direction */
     auto const P = tensor_product(d.vec21 / d.dist2, d.vec21);

src/core/dpd.hpp

@@ -83,4 +83,13 @@ Utils::Vector9d dpd_stress(System::System &system,
  */
 Utils::Vector9d dpd_pressure_local(System::System &system);
 
+/** 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::Vector3d dpd_noise(DPDThermostat const &dpd, int pid1, int pid2);
+
 #endif // ESPRESSO_DPD

src/core/pressure.cpp

@@ -107,6 +107,9 @@ std::shared_ptr<Observable_stat> System::calculate_pressure() {
                                get_ptr(coulomb_force_kernel),
                                get_ptr(coulomb_pressure_kernel),
                                *box_geo,
+#ifdef ESPRESSO_DPD
+                               thermostat->dpd.get(),
+#endif
                                cell_structure->get_unique_particles(),
                                local_pressure,
                                layout,

src/core/pressure_cabana.hpp

@@ -102,6 +102,9 @@ struct PressureKernel {
   Coulomb::ShortRangeForceKernel::kernel_type const *coulomb_f_kernel;
   Coulomb::ShortRangePressureKernel::kernel_type const *coulomb_p_kernel;
   BoxGeometry const &box_geo;
+#ifdef ESPRESSO_DPD
+  DPDThermostat const *dpd;
+#endif
   std::vector<Particle *> const &unique_particles;
   Kokkos::View<double **, Kokkos::LayoutRight> local_pressure;
   PressureBinLayout layout;
@@ -117,6 +120,9 @@ struct PressureKernel {
       Coulomb::ShortRangeForceKernel::kernel_type const *coulomb_f_kernel_,
       Coulomb::ShortRangePressureKernel::kernel_type const *coulomb_p_kernel_,
       BoxGeometry const &box_geo_,
+#ifdef ESPRESSO_DPD
+      DPDThermostat const *dpd_,
+#endif
       std::vector<Particle *> const &unique_particles_,
       Kokkos::View<double **, Kokkos::LayoutRight> const &local_pressure_,
       PressureBinLayout layout_, CellStructure::AoSoA_pack const &aosoa_,
@@ -125,9 +131,13 @@ struct PressureKernel {
       : bonded_ias(bonded_ias_), nonbonded_ias(nonbonded_ias_),
         coulomb(coulomb_), coulomb_f_kernel(coulomb_f_kernel_),
         coulomb_p_kernel(coulomb_p_kernel_), box_geo(box_geo_),
+#ifdef ESPRESSO_DPD
+        dpd(dpd_),
+#endif
         unique_particles(unique_particles_), local_pressure(local_pressure_),
         layout(layout_), aosoa(aosoa_), mol_id_view(std::move(mol_id_view_)),
-        system_max_cutoff(system_max_cutoff_), thermo_switch(thermo_switch_) {}
+        system_max_cutoff(system_max_cutoff_), thermo_switch(thermo_switch_) {
+  }
 
   KOKKOS_INLINE_FUNCTION
   void operator()(std::size_t i, std::size_t j) const {
@@ -179,13 +189,21 @@ struct PressureKernel {
 
 #ifdef ESPRESSO_DPD
         if (dpd_active(ia_params, thermo_switch)) {
+          auto const pid1 = aosoa.id(i);
+          auto const pid2 = aosoa.id(j);
+          auto const noise_vec = (ia_params.dpd.radial.pref > 0.0 ||
+                                  ia_params.dpd.trans.pref > 0.0)
+                                     ? dpd_noise(*dpd, pid1, pid2)
+                                     : Utils::Vector3d{};
           auto const vel1 = aosoa.get_vector_at(aosoa.velocity, i);
           auto const vel2 = aosoa.get_vector_at(aosoa.velocity, j);
           auto const v21 = box_geo.velocity_difference(pos1, pos2, vel1, vel2);
           auto const dist2 = d.norm2();
           // f_r/f_t: dissipative force from radial/transverse DPD channel
-          auto const f_r = dpd_pair_force(ia_params.dpd.radial, v21, dist, {});
-          auto const f_t = dpd_pair_force(ia_params.dpd.trans, v21, dist, {});
+          auto const f_r =
+              dpd_pair_force(ia_params.dpd.radial, v21, dist, noise_vec);
+          auto const f_t =
+              dpd_pair_force(ia_params.dpd.trans, v21, dist, noise_vec);
           auto const P = Utils::tensor_product(d / dist2, d);
           auto const f_d = P * (f_r - f_t) + f_t;
           auto const s = Utils::flatten(Utils::tensor_product(d, f_d));

testsuite/python/dpd.py

@@ -456,33 +456,84 @@ def calc_stress(dist, vel_diff):
         partcls = system.part.add(pos=pos)
         system.integrator.run(10)
 
-        for kT in [0., 2.]:
-            system.thermostat.set_dpd(kT=kT, seed=3)
-            # run 1 integration step to get velocities
-            partcls.v = np.zeros((n_part, 3))
-            system.integrator.run(steps=1)
+        # 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()
+        pairs = system.part.pairs()
 
-            stress = np.zeros([3, 3])
+        stress = 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)
+        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)
 
-            stress /= system.volume()
+        stress /= system.volume()
 
-            dpd_stress = system.analysis.dpd_stress()
+        dpd_stress = system.analysis.dpd_stress()
 
-            dpd_obs = espressomd.observables.DPDStress()
-            obs_stress = dpd_obs.calculate()
-            pressure = system.analysis.pressure_tensor()["dpd"]
+        dpd_obs = espressomd.observables.DPDStress()
+        obs_stress = 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_stress), stress)
+        np.testing.assert_array_almost_equal(np.copy(obs_stress), stress)
+        np.testing.assert_array_almost_equal(np.copy(pressure), -stress)
+
+    @utx.skipIfMissingFeatures("EXTERNAL_FORCES")
+    def test_dpd_stress_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 nonzero and distinct.
+        """
+        system = self.system
+        kT, gamma_r, gamma_t, r_cut = 2.0, 1.5, 0.7, 1.5
+        dt = system.time_step
+
+        system.thermostat.set_dpd(kT=kT, seed=42)
+        system.non_bonded_inter[0, 0].dpd.set_params(
+            weight_function=0, gamma=gamma_r, r_cut=r_cut,
+            trans_weight_function=0, trans_gamma=gamma_t, trans_r_cut=r_cut)
+
+        # Both particles are fixed with zero velocity so v21 == 0 exactly and
+        # the geometry never changes; each run(1) only advances the RNG counter.
+        d = np.array([0.5, 0.7, 0.9])  # |d| = sqrt(1.55) ~ 1.245 < r_cut
+        pos0 = np.array([5., 5., 5.])
+        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)
+        # with omega = 1 (weight_function=0), A/B the radial/trans amplitudes,
+        # and noise components iid with variance s2 = 1/12.
+        V = system.volume()
+        s2 = 1. / 12.
+        A = np.sqrt(24. * kT * gamma_r / dt)  # radial amplitude (omega=1)
+        B = np.sqrt(24. * kT * gamma_t / dt)  # trans amplitude  (omega=1)
+        dhat = d / np.linalg.norm(d)
+        sum_Mjk2 = (A**2 - B**2) * dhat**2 + B**2  # length-3 over j
+        var_analytic = np.outer(d**2, sum_Mjk2) * \
+            (s2 / V**2)  # (3, 3), all > 0
+
+        N = 5000
+        samples = np.empty((N, 3, 3))
+        for n in range(N):
+            system.integrator.run(1)
+            samples[n] = system.analysis.dpd_stress()
+
+        mean = samples.mean(axis=0)
+        var = samples.var(axis=0)
+        # (1) all 9 means ~ 0  (~6 sigma of the sample mean)
+        atol_mean = 6. * np.sqrt(var_analytic / N)
+        np.testing.assert_array_less(np.abs(mean), atol_mean)
+        # (2) all 9 variances ~ analytic
+        np.testing.assert_allclose(var, var_analytic, rtol=0.05)
 
     @utx.skipIfMissingFeatures("MASS")
     def test_momentum_conservation(self):

testsuite/samples/test_dpd.py

@@ -29,7 +29,7 @@ class Sample(ut.TestCase):
     def test_final_pressure(self):
         # test simulated value to +/- 10% accuracy
         self.assertLess(abs(sample.p_avg - 0.23), 0.02)
-        self.assertLess(abs(sample.p_std - 0.01), 0.01)
+        self.assertLess(abs(sample.p_std - 0.05), 0.02)
 
 
 if __name__ == "__main__":