Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Tests and Bugfixes for the GB potential #3091

Merged
merged 20 commits into from
Aug 23, 2019
Merged

Tests and Bugfixes for the GB potential #3091

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
20 commits
Select commit Hold shift + click to select a range
efe22bd
Issue #2922: Refactoring energy tests gb potential
Kischy Aug 21, 2019
d5bc3f1
Issue #2922: Added tests for gb forces and corrected gay_berne.hpp. R…
Kischy Aug 21, 2019
100362e
Issue #2922: Added gb torque tests. Some refactoring
Kischy Aug 21, 2019
7631e45
Issue #2922: Some refactoring
Kischy Aug 21, 2019
7986e5b
Issue #2922: Updated style to meet code formatting rules
Kischy Aug 21, 2019
d79e9c1
#2922: Changed misleading variable name
Kischy Aug 21, 2019
4fd0c68
Formatting
jngrad Aug 21, 2019
e421226
Refactor Gay-Berne force/energy kernels
jngrad Aug 22, 2019
d0f9a5c
#2922: Changed delta value
Kischy Aug 22, 2019
4ee19cd
Resolved small merge conflict
Kischy Aug 22, 2019
20e904c
Merge fixes of gay berne potential
Kischy Aug 22, 2019
26889e0
#2922:
Kischy Aug 22, 2019
a910bab
Merge branch 'python' into Bugfix_Tests_Gay_Berne_Issue_2922
Kischy Aug 22, 2019
e60d7e7
#2922:
Kischy Aug 22, 2019
4e3d366
#2922: Applied styling rules
Kischy Aug 22, 2019
377b216
Rename gradient function
jngrad Aug 23, 2019
9cb0b7c
testsuite: improve readability of gradient calculation.
KaiSzuttor Aug 23, 2019
27b283d
Merge branch 'Bugfix_Tests_Gay_Berne_Issue_2922' of github.com:Kischy…
KaiSzuttor Aug 23, 2019
9a0afdd
numpy as np.
KaiSzuttor Aug 23, 2019
975bc21
testsuite: docstring.
KaiSzuttor Aug 23, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
170 changes: 84 additions & 86 deletions src/core/nonbonded_interactions/gay_berne.hpp
View file
Open in desktop
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