From 14afa09639c879f11ae42254b6a590387328b262 Mon Sep 17 00:00:00 2001 From: Sebastian Bindgen Date: Mon, 9 Mar 2020 19:55:55 +0100 Subject: [PATCH 1/2] Exponent for DPD weight function --- src/core/dpd.cpp | 14 +++++----- src/core/dpd.hpp | 5 ++-- src/python/espressomd/interactions.pxd | 3 ++- src/python/espressomd/interactions.pyx | 7 ++++- testsuite/python/dpd.py | 37 ++++++++++++++++++++++++++ 5 files changed, 55 insertions(+), 11 deletions(-) diff --git a/src/core/dpd.cpp b/src/core/dpd.cpp index 1de4dc45475..87846abacfd 100644 --- a/src/core/dpd.cpp +++ b/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 ¶ms, 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; diff --git a/src/core/dpd.hpp b/src/core/dpd.hpp index a65433659d3..39f119c25d4 100644 --- a/src/core/dpd.hpp +++ b/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); diff --git a/src/python/espressomd/interactions.pxd b/src/python/espressomd/interactions.pxd index aa960cfdedc..316d5a52cea 100644 --- a/src/python/espressomd/interactions.pxd +++ b/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: diff --git a/src/python/espressomd/interactions.pyx b/src/python/espressomd/interactions.pyx index ce01b4c5bf4..e6ecf20e850 100644 --- a/src/python/espressomd/interactions.pyx +++ b/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): diff --git a/testsuite/python/dpd.py b/testsuite/python/dpd.py index 2e255abeffe..f95f22acc31 100644 --- a/testsuite/python/dpd.py +++ b/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 From 56362fcf0c1eca72832636520fb598ff8f54b364 Mon Sep 17 00:00:00 2001 From: Sebastian Bindgen Date: Tue, 10 Mar 2020 13:49:06 +0100 Subject: [PATCH 2/2] Documention for k in DPD interaction --- doc/sphinx/inter_non-bonded.rst | 23 ++++++++++++++--------- doc/sphinx/zrefs.bib | 11 +++++++++++ 2 files changed, 25 insertions(+), 9 deletions(-) diff --git a/doc/sphinx/inter_non-bonded.rst b/doc/sphinx/inter_non-bonded.rst index e2cf4bfb9b5..4624ff558fa 100644 --- a/doc/sphinx/inter_non-bonded.rst +++ b/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} diff --git a/doc/sphinx/zrefs.bib b/doc/sphinx/zrefs.bib index 82d4d9b7ce2..6aacd1f93f5 100644 --- a/doc/sphinx/zrefs.bib +++ b/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} +}