Skip to content

Commit

Permalink
Merge #3091
Browse files Browse the repository at this point in the history 
3091: Tests and Bugfixes for the GB potential r=KaiSzuttor a=Kischy

Fixes #2922

Description of changes:
 - Added tests for the gb potential for the force and torque
 - Fixed bugs in gay_berne.hpp
 - Removed condition in gay_berne.hpp

PR Checklist
------------
 - [x] Tests?
   - [ ] Interface
   - [x] Core 
 - [ ] Docs?


Co-authored-by: Christian Kischy Haege <[email protected]>
Co-authored-by: Jean-Noël Grad <[email protected]>
Co-authored-by: Christian Haege <[email protected]>
Co-authored-by: Kai Szuttor <[email protected]>
  • Loading branch information
@bors @Kischy
@jngrad @KaiSzuttor
5 people authored Aug 23, 2019
2 parents 6ef8d36 + 975bc21 commit 852b1a0
Show file tree
Hide file tree
Showing 2 changed files with 217 additions and 110 deletions.
170 changes: 84 additions & 86 deletions src/core/nonbonded_interactions/gay_berne.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,86 +52,84 @@ inline void add_gb_pair_force(Particle const *const p1,
return;
}

auto const u1 = p1->r.calc_director();
auto const u2 = p2->r.calc_director();
auto const a = d * u1;
auto const b = d * u2;
auto const c = u1 * u2;
auto const E1 = 1 / sqrt(1 - ia_params->gay_berne.chi1 *
ia_params->gay_berne.chi1 * c * c);
auto const Plus1 = (a + b) / (1 + ia_params->gay_berne.chi1 * c);
auto const Plus2 = (a + b) / (1 + ia_params->gay_berne.chi2 * c);
auto const Minus1 = (a - b) / (1 - ia_params->gay_berne.chi1 * c);
auto const Minus2 = (a - b) / (1 - ia_params->gay_berne.chi2 * c);
auto const Brhi2 = (ia_params->gay_berne.chi2 / dist / dist) *
(Plus2 * (a + b) + Minus2 * (a - b));
auto const E2 = 1 - 0.5 * Brhi2;
auto const E = 4 * ia_params->gay_berne.eps *
pow(E1, ia_params->gay_berne.nu) *
pow(E2, ia_params->gay_berne.mu);
auto const Brhi1 = (ia_params->gay_berne.chi1 / dist / dist) *
(Plus1 * (a + b) + Minus1 * (a - b));
auto const Sigma = ia_params->gay_berne.sig / sqrt(1 - 0.5 * Brhi1);
auto Koef1 = ia_params->gay_berne.mu / E2;
auto Koef2 = int_pow<3>(Sigma) * 0.5;

auto const X =
ia_params->gay_berne.sig / (dist - Sigma + ia_params->gay_berne.sig);
auto const Xcut =
ia_params->gay_berne.sig /
(ia_params->gay_berne.cut - Sigma + ia_params->gay_berne.sig);

if (X < 1.25) { /* 1.25 corresponds to the interparticle penetration of 0.2
units of length.
If they are not that close, the GB forces and torques are
calculated */

auto const X6 = int_pow<6>(X);
auto const Xcut6 = int_pow<6>(Xcut);

auto const Bra12 = 6 * X6 * X * (2 * X6 - 1);
auto const Bra12Cut = 6 * Xcut6 * Xcut * (2 * Xcut6 - 1);
auto const Brack = X6 * (X6 - 1);
auto const BrackCut = Xcut6 * (Xcut6 - 1);

/*-------- Here we calculate derivatives -----------------------------*/

auto const dU_dr = E *
(Koef1 * Brhi2 * (Brack - BrackCut) -
Koef2 * Brhi1 * (Bra12 - Bra12Cut) - Bra12 * dist) /
sqr(dist);
Koef1 *= ia_params->gay_berne.chi2 / sqr(dist);
Koef2 *= ia_params->gay_berne.chi1 / sqr(dist);
auto const dU_da = E * (Koef1 * (Minus2 + Plus2) * (BrackCut - Brack) +
Koef2 * (Plus1 + Minus1) * (Bra12 - Bra12Cut));
auto const dU_db = E * (Koef1 * (Minus2 - Plus2) * (Brack - BrackCut) +
Koef2 * (Plus1 - Minus1) * (Bra12 - Bra12Cut));
auto const dU_dc =
E * ((Brack - BrackCut) * (ia_params->gay_berne.nu *
sqr(E1 * ia_params->gay_berne.chi1) * c +
0.5 * Koef1 * ia_params->gay_berne.chi2 *
(sqr(Plus2) - sqr(Minus2))) -
(Bra12 - Bra12Cut) * 0.5 * Koef2 * ia_params->gay_berne.chi1 *
(sqr(Plus1) - sqr(Minus1)));

/*--------------------------------------------------------------------*/

force -= dU_dr * d + dU_da * u1 + dU_db * u2;

if (torque1 != nullptr) {
/* calculate torque: torque = u_1 x G */
auto const G2 = -dU_da * d - dU_dc * u2;
*torque1 += vector_product(u1, G2);

if (torque2 != nullptr) {
/* calculate torque: torque = u_2 x G */
auto const G1 = -dU_db * d - dU_dc * u1;
*torque2 += vector_product(u2, G1);
}
auto const e0 = ia_params->gay_berne.eps;
auto const s0 = ia_params->gay_berne.sig;
auto const chi1 = ia_params->gay_berne.chi1;
auto const chi2 = ia_params->gay_berne.chi2;
auto const mu = ia_params->gay_berne.mu;
auto const nu = ia_params->gay_berne.nu;
auto const ui = p1->r.calc_director();
auto const uj = p2->r.calc_director();
auto const r = Utils::Vector3d(d).normalize();
auto const dui = d * ui;
auto const duj = d * uj;
auto const rui = r * ui;
auto const ruj = r * uj;
auto const uij = ui * uj;
auto const oo1 = (dui + duj) / (1 + chi1 * uij);
auto const oo2 = (dui - duj) / (1 - chi1 * uij);
auto const tt1 = (dui + duj) / (1 + chi2 * uij);
auto const tt2 = (dui - duj) / (1 - chi2 * uij);
auto const o1 = sqr(rui + ruj) / (1 + chi1 * uij);
auto const o2 = sqr(rui - ruj) / (1 - chi1 * uij);
auto const t1 = sqr(rui + ruj) / (1 + chi2 * uij);
auto const t2 = sqr(rui - ruj) / (1 - chi2 * uij);
auto const Brhi1 = chi1 * (o1 + o2);
auto const Brhi2 = chi2 * (t1 + t2);

auto const e1 = 1. / (1. - sqr(chi1 * uij));
auto const e2 = 1. - 0.5 * Brhi2;
auto const e = 4 * e0 * pow(e1, 0.5 * nu) * pow(e2, mu);

auto const s1 = 1. / sqrt(1. - 0.5 * Brhi1);
auto const s = s0 * s1;
auto Koef1 = mu / e2;
auto Koef2 = int_pow<3>(s1) * 0.5;

auto const X = s0 / (dist - s + s0);
auto const Xcut = s0 / (ia_params->gay_berne.cut - s + s0);

auto const X6 = int_pow<6>(X);
auto const Xcut6 = int_pow<6>(Xcut);

auto const Bra12 = 6 * X6 * X * (2 * X6 - 1);
auto const Bra12Cut = 6 * Xcut6 * Xcut * (2 * Xcut6 - 1);
auto const Brack = X6 * (X6 - 1);
auto const BrackCut = Xcut6 * (Xcut6 - 1);

/*-------- Here we calculate derivatives -----------------------------*/

auto const dU_dr = e *
(Koef1 * Brhi2 * (Brack - BrackCut) -
Koef2 * Brhi1 * (Bra12 - Bra12Cut) - Bra12 * dist / s0) /
sqr(dist);

Koef1 *= chi2 / sqr(dist);
Koef2 *= chi1 / sqr(dist);

auto const dU_da = e * (Koef1 * (tt1 + tt2) * (BrackCut - Brack) +
Koef2 * (oo1 + oo2) * (Bra12 - Bra12Cut));
auto const dU_db = e * (Koef1 * (tt2 - tt1) * (Brack - BrackCut) +
Koef2 * (oo1 - oo2) * (Bra12 - Bra12Cut));
auto const dU_dc =
e * ((Brack - BrackCut) * (nu * e1 * sqr(chi1) * uij +
0.5 * Koef1 * chi2 * (sqr(tt1) - sqr(tt2))) -
(Bra12 - Bra12Cut) * 0.5 * Koef2 * chi1 * (sqr(oo1) - sqr(oo2)));

/*--------------------------------------------------------------------*/

force -= dU_dr * d + dU_da * ui + dU_db * uj;

if (torque1 != nullptr) {
/* calculate torque: torque = u_1 x G */
auto const G2 = -dU_da * d - dU_dc * uj;
*torque1 += vector_product(ui, G2);

if (torque2 != nullptr) {
/* calculate torque: torque = u_2 x G */
auto const G1 = -dU_db * d - dU_dc * ui;
*torque2 += vector_product(uj, G1);
}
} else { /* the particles are too close to each other */
Koef1 = 100;
force += Koef1 * d;
}
}

Expand All @@ -151,25 +149,25 @@ inline double gb_pair_energy(Particle const *const p1, Particle const *const p2,
auto const chi2 = ia_params->gay_berne.chi2;
auto const mu = ia_params->gay_berne.mu;
auto const nu = ia_params->gay_berne.nu;
auto const r = Utils::Vector3d({d[0], d[1], d[2]}).normalize();
auto const r = Utils::Vector3d(d).normalize();

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;

auto const e1 = std::pow(1. - sqr(chi1 * uij), -0.5 * nu);

auto const o1 = sqr(rui + ruj) / (1 + chi1 * uij);
auto const o2 = sqr(rui - ruj) / (1 - chi1 * uij);
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 e1 = std::pow(1. - sqr(chi1 * uij), -0.5 * nu);
auto const e2 = std::pow(1. - 0.5 * chi2 * (t1 + t2), mu);
auto const e = e0 * e1 * e2;

auto const s = s0 / std::sqrt(1. - 0.5 * chi1 *
(sqr(rui + ruj) / (1 + chi1 * uij) +
sqr(rui - ruj) / (1 - chi1 * uij)));
auto const s1 = 1. / std::sqrt(1. - 0.5 * chi1 * (o1 + o2));
auto const s = s0 * s1;

auto r_eff = [=](double r) { return (r - s + s0) / s0; };
auto E = [=](double r) {
Expand Down
Loading

0 comments on commit 852b1a0

Please sign in to comment.