From 4f41ee8a9238637b94f11b5325ff55591c6830f7 Mon Sep 17 00:00:00 2001 From: RudolfWeeber Date: Mon, 7 Jan 2019 09:32:21 +0100 Subject: [PATCH] Gay Berne potential (#2424) * core: gb refactor * Work on gb test * Testsuite: Gay Berne pair potential (no forces/torques) * Core: Apply corrections suggested in #2229 to forces/torques * GAY_BERNE depends on EXPERIMENTAL_FEATURES * Testsuite: Gay-Berne, avoid dipole hack to set orientation * Formatting * Testsuite: Fix gb energy test (ia deactivation) * Testsuite: Gay-Berne on 3 cores and silent --- src/core/nonbonded_interactions/gb.hpp | 67 ++++++++++----------- src/features.def | 2 +- testsuite/python/interactions_non-bonded.py | 32 +++++++++- testsuite/python/tests_common.py | 26 ++++++++ 4 files changed, 91 insertions(+), 36 deletions(-) diff --git a/src/core/nonbonded_interactions/gb.hpp b/src/core/nonbonded_interactions/gb.hpp index 5461dc27aa0..66dd6fcaeb5 100644 --- a/src/core/nonbonded_interactions/gb.hpp +++ b/src/core/nonbonded_interactions/gb.hpp @@ -30,6 +30,8 @@ #include "nonbonded_interaction_data.hpp" #include "particle_data.hpp" #include "utils.hpp" +#include "utils/math/int_pow.hpp" +#include "utils/math/sqr.hpp" #ifdef GAY_BERNE @@ -80,8 +82,8 @@ add_gb_pair_force(const Particle *const p1, const Particle *const p2, Koef1 = ia_params->GB_mu / E2; Koef2 = Sigma * Sigma * Sigma * 0.5; - X = 1 / (dist - Sigma + ia_params->GB_sig); - Xcut = 1 / (ia_params->GB_cut - Sigma + ia_params->GB_sig); + X = ia_params->GB_sig / (dist - Sigma + ia_params->GB_sig); + Xcut = ia_params->GB_sig / (ia_params->GB_cut - Sigma + ia_params->GB_sig); if (X < 1.25) { /* 1.25 corresponds to the interparticle penetration of 0.2 units of length. @@ -163,47 +165,44 @@ add_gb_pair_force(const Particle *const p1, const Particle *const p2, inline double gb_pair_energy(const Particle *p1, const Particle *p2, const IA_parameters *ia_params, const double d[3], double dist) { + using Utils::int_pow; + using Utils::sqr; + if (!(dist < ia_params->GB_cut)) return 0.0; - double a, b, c, X, Xcut, Brack, BrackCut, u1x, u1y, u1z, u2x, u2y, u2z, E, E1, - E2, Sigma, Plus1, Minus1, Plus2, Minus2; + auto const e0 = ia_params->GB_eps; + auto const s0 = ia_params->GB_sig; + auto const chi1 = ia_params->GB_chi1; + auto const chi2 = ia_params->GB_chi2; + auto const mu = ia_params->GB_mu; + auto const nu = ia_params->GB_nu; + auto const r = Vector3d({d[0], d[1], d[2]}).normalize(); - u1x = p1->r.calc_director()[0]; - u1y = p1->r.calc_director()[1]; - u1z = p1->r.calc_director()[2]; - u2x = p2->r.calc_director()[0]; - u2y = p2->r.calc_director()[1]; - u2z = p2->r.calc_director()[2]; + auto const ui = p1->r.calc_director(); + auto const uj = p2->r.calc_director(); + auto const uij = ui * uj; + auto const rui = r * ui; + auto const ruj = r * uj; - a = d[0] * u1x + d[1] * u1y + d[2] * u1z; - b = d[0] * u2x + d[1] * u2y + d[2] * u2z; - c = u1x * u2x + u1y * u2y + u1z * u2z; + auto const e1 = std::pow(1. - sqr(chi1 * uij), -0.5 * nu); - Plus1 = (a + b) / (1 + ia_params->GB_chi1 * c); - Plus2 = (a + b) / (1 + ia_params->GB_chi2 * c); - Minus1 = (a - b) / (1 - ia_params->GB_chi1 * c); - Minus2 = (a - b) / (1 - ia_params->GB_chi2 * c); - E1 = 1 / sqrt(1 - ia_params->GB_chi1 * ia_params->GB_chi1 * c * c); - E2 = 1 - 0.5 * (ia_params->GB_chi2 / dist / dist) * - (Plus2 * (a + b) + Minus2 * (a - b)); - E = 4 * ia_params->GB_eps * pow(E1, ia_params->GB_nu) * - pow(E2, ia_params->GB_mu); - Sigma = - ia_params->GB_sig / sqrt(1 - 0.5 * (ia_params->GB_chi1 / dist / dist) * - (Plus1 * (a + b) + Minus1 * (a - b))); + auto const t1 = sqr(rui + ruj) / (1 + chi2 * uij); + auto const t2 = sqr(rui - ruj) / (1 - chi2 * uij); + auto const e2 = std::pow(1. - 0.5 * chi2 * (t1 + t2), mu); + + auto const e = e0 * e1 * e2; - X = 1 / (dist - Sigma + ia_params->GB_sig); - Xcut = 1 / (ia_params->GB_cut - Sigma + ia_params->GB_sig); + auto const s = s0 / std::sqrt(1. - 0.5 * chi1 * + (sqr(rui + ruj) / (1 + chi1 * uij) + + sqr(rui - ruj) / (1 - chi1 * uij))); - Brack = X * X * X; - BrackCut = Xcut * Xcut * Xcut; - Brack = Brack * Brack; - BrackCut = BrackCut * BrackCut; - Brack = Brack * (Brack - 1); - BrackCut = BrackCut * (BrackCut - 1); + auto r_eff = [=](double r) { return (r - s + s0) / s0; }; + auto E = [=](double r) { + return 4. * e * (int_pow<12>(1. / r) - int_pow<6>(1. / r)); + }; - return E * (Brack - BrackCut); + return E(r_eff(dist)) - E(r_eff(ia_params->GB_cut)); } #endif diff --git a/src/features.def b/src/features.def index a2fb9fe91a1..9e6e3f7e180 100644 --- a/src/features.def +++ b/src/features.def @@ -92,7 +92,7 @@ LJCOS LJCOS2 LJGEN_SOFTCORE implies LENNARD_JONES_GENERIC COS2 -GAY_BERNE +GAY_BERNE depends EXPERIMENTAL_FEATURES SMOOTH_STEP HERTZIAN GAUSSIAN diff --git a/testsuite/python/interactions_non-bonded.py b/testsuite/python/interactions_non-bonded.py index 37a0e0fe8b8..12a2f8865e8 100644 --- a/testsuite/python/interactions_non-bonded.py +++ b/testsuite/python/interactions_non-bonded.py @@ -36,7 +36,7 @@ class InteractionsNonBondedTest(ut.TestCase): def setUp(self): self.system.box_l = [self.box_l] * 3 - self.system.cell_system.skin = 0.4 + self.system.cell_system.skin = 0. self.system.time_step = .1 self.system.part.add(id=0, pos=self.start_pos, type=0) @@ -585,6 +585,36 @@ def test_gaussian(self): self.system.non_bonded_inter[0, 0].gaussian.set_params(eps=0.) + @ut.skipIf(not espressomd.has_features("GAY_BERNE"), "skipped for lack of Gay Berne") + def test_gb(self): + k_1 = 1.2 + k_2 = 2.4 + mu = 2. + nu = 5. + sigma_0 = 1.2 + epsilon_0 = 0.8 + cut = 3.3 + + self.system.part[:].pos = ((1, 2, 3), (2.2, 2.1, 2.9)) + self.system.non_bonded_inter[0, 0].gay_berne.set_params( + sig=sigma_0, cut=cut, eps=epsilon_0, k1=k_1, k2=k_2, mu=mu, nu=nu) + p1 = self.system.part[0] + p2 = self.system.part[1] + p1.rotate(axis=(1, 2, 3), angle=0.3) + p1.rotate(axis=(1, -2, -4), angle=1.2) + + r = self.system.distance_vec(p1, p2) + + r_cut = r * cut / numpy.linalg.norm(r) + self.assertAlmostEqual(self.system.analysis.energy()["non_bonded"], + tests_common.gay_berne_potential(r, p1.director, p2.director, epsilon_0, sigma_0, mu, nu, k_1, k_2) - + tests_common.gay_berne_potential(r_cut, p1.director, p2.director, epsilon_0, sigma_0, mu, nu, k_1, k_2), delta=1E-14) + + self.system.integrator.run(0) + self.system.non_bonded_inter[0, 0].gay_berne.set_params( + sig=sigma_0, cut=0, eps=0, k1=k_1, k2=k_2, mu=mu, nu=nu) + self.system.integrator.run(0) + self.assertEqual(self.system.analysis.energy()["non_bonded"], 0.0) if __name__ == '__main__': print("Features: ", espressomd.features()) diff --git a/testsuite/python/tests_common.py b/testsuite/python/tests_common.py index c6a3058756a..0cf3e47cfca 100644 --- a/testsuite/python/tests_common.py +++ b/testsuite/python/tests_common.py @@ -583,6 +583,32 @@ def gaussian_force(r, eps, sig, cutoff): return f +def gay_berne_potential(r_ij, u_i, u_j, epsilon_0, sigma_0, mu, nu, k_1, k_2): + r_normed = r_ij / np.linalg.norm(r_ij) + r_u_i = np.dot(r_normed, u_i) + r_u_j = np.dot(r_normed, u_j) + u_i_u_j = np.dot(u_i, u_j) + + chi = (k_1**2 - 1.) / (k_1**2 + 1.) + chi_d = (k_2**(1. / mu) - 1) / (k_2**(1. / mu) + 1) + + sigma = sigma_0 \ + / np.sqrt( + (1 - 0.5 * chi * ( + (r_u_i + r_u_j)**2 / (1 + chi * u_i_u_j) + + (r_u_i - r_u_j)**2 / (1 - chi * u_i_u_j)))) + + epsilon = epsilon_0 *\ + (1 - chi**2 * u_i_u_j**2)**(-nu / 2.) *\ + (1 - chi_d / 2. * ( + (r_u_i + r_u_j)**2 / (1 + chi_d * u_i_u_j) + + (r_u_i - r_u_j)**2 / (1 - chi_d * u_i_u_j)))**mu + + rr = np.linalg.norm((np.linalg.norm(r_ij) - sigma + sigma_0) / sigma_0) + + return 4. * epsilon * (rr**-12 - rr**-6) + + class DynamicDict(dict): def __getitem__(self, key):