Commit 4f41ee8a authored by RudolfWeeber's avatar RudolfWeeber Committed by Florian Weik

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
parent 9f9a1255
Pipeline #4800 passed with stages
in 72 minutes and 14 seconds
......@@ -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
......
......@@ -92,7 +92,7 @@ LJCOS
LJCOS2
LJGEN_SOFTCORE implies LENNARD_JONES_GENERIC
COS2
GAY_BERNE
GAY_BERNE depends EXPERIMENTAL_FEATURES
SMOOTH_STEP
HERTZIAN
GAUSSIAN
......
......@@ -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 @@ class InteractionsNonBondedTest(ut.TestCase):
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())
......
......@@ -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):
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment