Parametric weight function for DPD....

doc/sphinx/inter_non-bonded.rst

@@ -529,12 +529,13 @@ The parameters can be set via::
     system.non_bonded_inter[type1, type2].dpd.set_params(**kwargs)
 
 This command defines an interaction between particles of the types ``type1`` and ``type2``
-that contains velocity-dependent friction and noise according 
-to a temperature set by :py:meth:`espressomd.thermostat.Thermostat.set_dpd()`. 
+that contains velocity-dependent friction and noise according
+to a temperature set by :py:meth:`espressomd.thermostat.Thermostat.set_dpd()`.
 The parameters for the interaction are
 
 * ``gamma``
 * ``weight_function``
+* ``k``
 * ``r_cut``
 * ``trans_gamma``
 * ``trans_weight_function``
@@ -556,9 +557,9 @@ for the dissipative force and
 .. math:: \vec{F}_{ij}^R = \sqrt{2 k_B T \gamma_\parallel w_\parallel (r_{ij}) }  \eta_{ij}(t) \hat{r}_{ij}
 
 for the random force. This introduces the friction coefficient :math:`\gamma_\parallel` (parameter ``gamma``) and the weight function
-:math:`w_\parallel`. The thermal energy :math:`k_B T` is not set by the interaction, 
-but by the DPD thermostat (:py:meth:`espressomd.thermostat.Thermostat.set_dpd()`) 
-to be equal for all particles. The weight function can be specified via the ``weight_function`` switch. 
+:math:`w_\parallel`. The thermal energy :math:`k_B T` is not set by the interaction,
+but by the DPD thermostat (:py:meth:`espressomd.thermostat.Thermostat.set_dpd()`)
+to be equal for all particles. The weight function can be specified via the ``weight_function`` switch.
 The possible values for ``weight_function`` are 0 and 1, corresponding to the
 order of :math:`w_\parallel`:
 
@@ -567,22 +568,26 @@ order of :math:`w_\parallel`:
    w_\parallel (r_{ij}) =  \left\{
    \begin{array}{clcr}
                 1                      & , \; \text{weight_function} = 0 \\
-                {( 1 - \frac{r_{ij}}{r^\text{cut}_\parallel}} )^2 & , \; \text{weight_function} = 1
+                {( 1 - (\frac{r_{ij}}{r^\text{cut}_\parallel})^k} )^2 & , \; \text{weight_function} = 1
       \end{array}
       \right.
-      
+
 Both weight functions are set to zero for :math:`r_{ij}>r^\text{cut}_\parallel` (parameter ``r_cut``).
 
+In case the ``weight_function`` 1 is selected the parameter ``k`` can be chosen. :math:`k = 1` is the
+default and recovers the linear ramp. :math:`k > 1` enhances the dissipative nature of the interaction
+and thus yields higher Schmidt numbers :cite:`yaghoubi2015a`.
+
 The random force has the properties
 
 .. math:: <\eta_{ij}(t)> = 0 , <\eta_{ij}^\alpha(t)\eta_{kl}^\beta(t')> = \delta_{\alpha\beta} \delta_{ik}\delta_{jl}\delta(t-t')
 
-and is numerically discretized to a random number :math:`\overline{\eta}` for each spatial 
+and is numerically discretized to a random number :math:`\overline{\eta}` for each spatial
 component for each particle pair drawn from a uniform distribution
 with properties
 
 .. math:: <\overline{\eta}> = 0 , <\overline{\eta}\overline{\eta}> = 1/dt
-    
+
 For the perpendicular part, the dissipative and random force are calculated analogously
 
 .. math:: \vec{F}_{ij}^{D} = -\gamma_\bot w^D (r_{ij}) (I-\hat{r}_{ij}\otimes\hat{r}_{ij}) \cdot \vec{v}_{ij}

doc/sphinx/zrefs.bib

@@ -1064,3 +1064,14 @@ @Article{hofling11a
   doi = {10.1039/C0SM00718H},
   owner = {rajarshi}
 }
+
+@article{yaghoubi2015a,
+  title={New modified weight function for the dissipative force in the DPD method to increase the Schmidt number},
+  author={Yaghoubi, S and Shirani, E and Pishevar, AR and Afshar, Yaser},
+  journal={EPL (Europhysics Letters)},
+  volume={110},
+  number={2},
+  pages={24002},
+  year={2015},
+  publisher={IOP Publishing}
+}

src/core/dpd.cpp

@@ -52,14 +52,14 @@ Vector3d dpd_noise(uint32_t pid1, uint32_t pid2) {
       dpd.rng_get(), (pid1 < pid2) ? pid2 : pid1, (pid1 < pid2) ? pid1 : pid2);
 }
 
-int dpd_set_params(int part_type_a, int part_type_b, double gamma, double r_c,
-                   int wf, double tgamma, double tr_c, int twf) {
+int dpd_set_params(int part_type_a, int part_type_b, double gamma, double k,
+                   double r_c, int wf, double tgamma, double tr_c, int twf) {
   IA_parameters &ia_params = *get_ia_param_safe(part_type_a, part_type_b);
 
   ia_params.dpd_radial = DPDParameters{
-      gamma, r_c, wf, sqrt(24.0 * temperature * gamma / time_step)};
+      gamma, k, r_c, wf, sqrt(24.0 * temperature * gamma / time_step)};
   ia_params.dpd_trans = DPDParameters{
-      tgamma, tr_c, twf, sqrt(24.0 * temperature * tgamma / time_step)};
+      tgamma, k, tr_c, twf, sqrt(24.0 * temperature * tgamma / time_step)};
 
   /* broadcast interaction parameters */
   mpi_bcast_ia_params(part_type_a, part_type_b);
@@ -91,17 +91,17 @@ void dpd_update_params(double pref_scale) {
   }
 }
 
-static double weight(int type, double r_cut, double r) {
+static double weight(int type, double r_cut, double k, double r) {
   if (type == 0) {
     return 1.;
   }
-  return 1. - r / r_cut;
+  return 1. - pow((r / r_cut), k);
 }
 
 Vector3d dpd_pair_force(DPDParameters const &params, Vector3d const &v,
                         double dist, Vector3d const &noise) {
   if (dist < params.cutoff) {
-    auto const omega = weight(params.wf, params.cutoff, dist);
+    auto const omega = weight(params.wf, params.cutoff, params.k, dist);
     auto const omega2 = Utils::sqr(omega);
 
     auto const f_d = params.gamma * omega2 * v;

src/core/dpd.hpp

@@ -37,13 +37,14 @@ struct IA_parameters;
 
 struct DPDParameters {
   double gamma = 0.;
+  double k = 1.;
   double cutoff = -1.;
   int wf = 0;
   double pref = 0.0;
 };
 
-int dpd_set_params(int part_type_a, int part_type_b, double gamma, double r_c,
-                   int wf, double tgamma, double tr_c, int twf);
+int dpd_set_params(int part_type_a, int part_type_b, double gamma, double k,
+                   double r_c, int wf, double tgamma, double tr_c, int twf);
 void dpd_init();
 void dpd_update_params(double pref2_scale);
 

src/python/espressomd/interactions.pxd

@@ -40,6 +40,7 @@ cdef extern from "TabulatedPotential.hpp":
 cdef extern from "dpd.hpp":
     cdef struct DPDParameters:
         double gamma
+        double k
         double cutoff
         int wf
         double pref
@@ -282,7 +283,7 @@ IF GAUSSIAN:
 IF DPD:
     cdef extern from "dpd.hpp":
         int dpd_set_params(int part_type_a, int part_type_b,
-                           double gamma, double r_c, int wf,
+                           double gamma, double k, double r_c, int wf,
                            double tgamma, double tr_c, int twf)
 
 IF HAT:

src/python/espressomd/interactions.pyx

@@ -859,6 +859,7 @@ IF DPD:
             return {
                 "weight_function": ia_params.dpd_radial.wf,
                 "gamma": ia_params.dpd_radial.gamma,
+                "k": ia_params.dpd_radial.k,
                 "r_cut": ia_params.dpd_radial.cutoff,
                 "trans_weight_function": ia_params.dpd_trans.wf,
                 "trans_gamma": ia_params.dpd_trans.gamma,
@@ -879,6 +880,8 @@ IF DPD:
                 either 0 (constant) or 1 (linear)
             gamma : :obj:`float`
                 Friction coefficient of the parallel part
+            k : :obj:`float`
+                Exponent in the modified weight function
             r_cut : :obj:`float`
                 Cutoff of the parallel part
             trans_weight_function : :obj:`int`, \{0, 1\}
@@ -896,6 +899,7 @@ IF DPD:
             if dpd_set_params(self._part_types[0],
                               self._part_types[1],
                               self._params["gamma"],
+                              self._params["k"],
                               self._params["r_cut"],
                               self._params["weight_function"],
                               self._params["trans_gamma"],
@@ -907,6 +911,7 @@ IF DPD:
             return {
                 "weight_function": 0,
                 "gamma": 0.0,
+                "k": 1.0,
                 "r_cut": -1.0,
                 "trans_weight_function": 0,
                 "trans_gamma": 0.0,
@@ -916,7 +921,7 @@ IF DPD:
             return "DPD"
 
         def valid_keys(self):
-            return {"weight_function", "gamma", "r_cut",
+            return {"weight_function", "gamma", "k", "r_cut",
                     "trans_weight_function", "trans_gamma", "trans_r_cut"}
 
         def required_keys(self):

testsuite/python/dpd.py

@@ -279,6 +279,43 @@ def calc_omega(dist, r_cut):
         np.testing.assert_array_equal(
             np.copy(s.part[0].f), -np.copy(s.part[1].f))
 
+    def test_parabolic_weight_function(self):
+        s = self.s
+        kT = 0.
+        gamma = 1.42
+        kappa = 7.8
+        r_cut = 1.2
+        s.thermostat.set_dpd(kT=kT, seed=42)
+        s.non_bonded_inter[0, 0].dpd.set_params(
+            weight_function=1, gamma=gamma, k=kappa, r_cut=r_cut,
+            trans_weight_function=1, trans_gamma=0.0, trans_r_cut=0.0)
+
+        def calc_omega(dist, r_cut):
+            return (1. - (dist / r_cut) ** kappa)
+
+        s.part.add(id=0, pos=[5, 5, 5], type=0, v=[0, 0, 0])
+        v = np.array([.5, 0., 0.])
+        s.part.add(id=1, pos=[3, 5, 5], type=0, v=v)
+
+        # Outside of both cutoffs, forces should be 0
+        for f in s.part[:].f:
+            np.testing.assert_array_equal(f, [0., 0., 0.])
+
+        # Place the particle at different positions to test the parabolic
+        # weight function
+        for dist in np.arange(0.1, 1.2, 50):
+
+            s.part[1].pos = [5. + dist, 5., 5.]
+            s.integrator.run(0)
+            omega = calc_omega(dist, r_cut)**2
+
+            # The particle is moved along the x-direction. Hence, we are
+            # testing the x element.
+            np.testing.assert_allclose(
+                np.copy(s.part[0].f), omega * gamma * v, rtol=0, atol=1e-11)
+            np.testing.assert_array_equal(
+                np.copy(s.part[0].f), -np.copy(s.part[1].f))
+
     def test_ghosts_have_v(self):
         s = self.s