Skip to content

Commit

Permalink
Gay Berne potential (#2424)
Browse files Browse the repository at this point in the history
* 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
  • Loading branch information
RudolfWeeber authored and fweik committed Jan 7, 2019
1 parent 9f9a125 commit 4f41ee8
Showing 4 changed files with 91 additions and 36 deletions.
67 changes: 33 additions & 34 deletions src/core/nonbonded_interactions/gb.hpp
Original file line number Diff line number Diff line change
@@ -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
2 changes: 1 addition & 1 deletion src/features.def
Original file line number Diff line number Diff line change
@@ -92,7 +92,7 @@ LJCOS
LJCOS2
LJGEN_SOFTCORE implies LENNARD_JONES_GENERIC
COS2
GAY_BERNE
GAY_BERNE depends EXPERIMENTAL_FEATURES
SMOOTH_STEP
HERTZIAN
GAUSSIAN
32 changes: 31 additions & 1 deletion testsuite/python/interactions_non-bonded.py
Original file line number Diff line number Diff line change
@@ -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())
26 changes: 26 additions & 0 deletions testsuite/python/tests_common.py
Original file line number Diff line number Diff line change
@@ -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):

0 comments on commit 4f41ee8

Please sign in to comment.