[SD] thermalization test

libs/stokesian_dynamics/sd.hpp

@@ -1015,7 +1015,9 @@ struct thermalizer {
         uint4 rnd_ints = curand_Philox4x32_10(make_uint4(offset >> 32, seed >> 32, index >> 32, index),
                                               make_uint2(offset, seed));
         T rnd = _curand_uniform_double_hq(rnd_ints.w, rnd_ints.x);
-        return 12 * std::sqrt(T{2.0}) * kT * (rnd - 0.5);
+        return 12 * kT * std::sqrt(T{12.0}) * (rnd - 0.5);
+        // 12 * kT * time_step is the desired variance
+        // the rest is a random number with unit variance and zero mean
     }
 };
 

samples/dancing.py

@@ -20,7 +20,7 @@
 gravity = constraints.Gravity(g=[0, -1, 0])
 system.constraints.add(gravity)
 
-intsteps = 13000
+intsteps = int(13000 / system.time_step)
 pos = np.empty([intsteps, 3 * len(system.part)])
 for i in range(intsteps):
     system.integrator.run(1)

src/core/stokesian_dynamics/sd_interface.cpp

@@ -138,15 +138,15 @@ void propagate_vel_pos_sd() {
 
 #if defined(BLAS) && defined(LAPACK)
     case CPU:
-      v_sd = sd_cpu(x_host, f_host, a_host, n_part, sd_viscosity, sd_kT, offset,
-                    sd_seed, sd_flags);
+      v_sd = sd_cpu(x_host, f_host, a_host, n_part, sd_viscosity,
+                    sd_kT / time_step, offset, sd_seed, sd_flags);
       break;
 #endif
 
 #ifdef CUDA
     case GPU:
-      v_sd = sd_gpu(x_host, f_host, a_host, n_part, sd_viscosity, sd_kT, offset,
-                    sd_seed, sd_flags);
+      v_sd = sd_gpu(x_host, f_host, a_host, n_part, sd_viscosity,
+                    sd_kT / time_step, offset, sd_seed, sd_flags);
       break;
 #endif
 

testsuite/python/stokesian_dynamics.py

@@ -1,15 +1,17 @@
 import espressomd
 from espressomd import constraints
 import numpy as np
-import matplotlib.pyplot as plt
 import unittest as ut
 import unittest_decorators as utx
 from tests_common import abspath
+import scipy.optimize
+
+s = espressomd.System(box_l=[1.0, 1.0, 1.0])
 
 
 @utx.skipIfMissingFeatures(["STOKESIAN_DYNAMICS"])
 class StokesianDynamicsTest(ut.TestCase):
-    system = espressomd.System(box_l=[1.0, 1.0, 1.0])
+    system = s
 
     # Digitized reference data of Figure 5b from
     #
@@ -20,6 +22,7 @@ class StokesianDynamicsTest(ut.TestCase):
     def setUp(self):
         self.system.part.clear()
         self.system.box_l = [10] * 3
+        self.system.periodicity = [0, 0, 0]
 
         self.system.time_step = 1.0
         self.system.cell_system.skin = 0.4
@@ -39,21 +42,21 @@ def setUp(self):
         self.system.constraints.add(gravity)
 
     def test(self):
-        intsteps = 8000
-        self.pos = np.empty([intsteps + 1, 3 * len(self.system.part)])
+        intsteps = int(8000 / self.system.time_step)
+        pos = np.empty([intsteps + 1, 3 * len(self.system.part)])
 
-        self.pos[0, :] = self.system.part[:].pos.flatten()
+        pos[0, :] = self.system.part[:].pos.flatten()
         for i in range(intsteps):
             self.system.integrator.run(1)
             for n, p in enumerate(self.system.part):
-                self.pos[i + 1, 3 * n:3 * n + 3] = p.pos
+                pos[i + 1, 3 * n:3 * n + 3] = p.pos
 
         for i in range(self.data.shape[0]):
             for n in range(3):
                 x_ref = self.data[i, 2 * n]
                 y_ref = self.data[i, 2 * n + 1]
-                x = self.pos[:, 3 * n]
-                y = self.pos[:, 3 * n + 1]
+                x = pos[:, 3 * n]
+                y = pos[:, 3 * n + 1]
 
                 if y_ref < -555:
                     continue
@@ -63,5 +66,65 @@ def test(self):
                 self.assertLess(dist, 0.5)
 
 
+@utx.skipIfMissingFeatures(["STOKESIAN_DYNAMICS"])
+class StokesianDiffusionTest(ut.TestCase):
+    system = s
+
+    kT = 1e-4
+    R = 1.0
+    eta = 1.0
+
+    def setUp(self):
+        self.system.part.clear()
+        self.system.box_l = [10] * 3
+        self.system.periodicity = [0, 0, 0]
+
+        self.system.time_step = 1.0
+        self.system.cell_system.skin = 0.4
+
+        self.system.part.add(pos=[0, 0, 0], rotation=[1, 1, 1])
+
+        from espressomd.thermostat import flags
+        self.system.thermostat.set_sd(viscosity=self.eta,
+                                      device="cpu",
+                                      radii={0: self.R},
+                                      kT=self.kT,
+                                      seed=0,
+                                      flags=flags.SELF_MOBILITY | flags.PAIR_MOBILITY | flags.FTS)
+        self.system.integrator.set_sd()
+
+    def test(self):
+        intsteps = int(1000000 / self.system.time_step)
+        pos = np.empty([intsteps + 1, 3])
+        orientation = np.empty((intsteps + 1, 3))
+
+        pos[0, :] = self.system.part[0].pos
+        orientation[0, :] = self.system.part[0].director
+        for i in range(intsteps):
+            self.system.integrator.run(1)
+            pos[i + 1, :] = self.system.part[0].pos
+            orientation[i + 1, :] = self.system.part[0].director
+
+        t = np.arange(0, intsteps + 1)
+        msd = np.linalg.norm(pos - pos[0, :], axis=1)**2
+        costheta = np.dot(orientation[:, :], orientation[0, :])
+
+        D_expected = self.kT / (6 * np.pi * self.eta * self.R)
+        fit = scipy.optimize.curve_fit(lambda t, D: 6 * D * t, t, msd)
+        D_measured = fit[0][0]
+        self.assertAlmostEqual(D_expected, D_measured, delta=D_expected * 0.2)
+
+        Dr_expected = self.kT / (8 * np.pi * self.eta * self.R**3)
+        tr_expected = int(1 / (2 * Dr_expected))
+        fit = scipy.optimize.curve_fit(lambda t, Dr: np.exp(-2 * Dr * t),
+                                       t[:2 * tr_expected],
+                                       costheta[:2 * tr_expected])
+        Dr_measured = fit[0][0]
+        self.assertAlmostEqual(
+            Dr_expected,
+            Dr_measured,
+            delta=Dr_expected * 1)
+
+
 if __name__ == '__main__':
     ut.main()