69 auto const dui = d * ui;
70 auto const duj = d * uj;
71 auto const rui = r * ui;
72 auto const ruj = r * uj;
73 auto const uij = ui * uj;
74 auto const oo1 = (dui + duj) / (1 + chi1 * uij);
75 auto const oo2 = (dui - duj) / (1 - chi1 * uij);
76 auto const tt1 = (dui + duj) / (1 + chi2 * uij);
77 auto const tt2 = (dui - duj) / (1 - chi2 * uij);
78 auto const o1 = sqr(rui + ruj) / (1 + chi1 * uij);
79 auto const o2 = sqr(rui - ruj) / (1 - chi1 * uij);
80 auto const t1 = sqr(rui + ruj) / (1 + chi2 * uij);
81 auto const t2 = sqr(rui - ruj) / (1 - chi2 * uij);
82 auto const Brhi1 = chi1 * (o1 + o2);
83 auto const Brhi2 = chi2 * (t1 + t2);
85 auto const e1 = 1. / (1. - sqr(chi1 * uij));
86 auto const e2 = 1. - 0.5 * Brhi2;
87 auto const e = 4 * e0 * pow(e1, 0.5 * nu) * pow(e2, mu);
89 auto const s1 = 1. / sqrt(1. - 0.5 * Brhi1);
90 auto const s = s0 * s1;
92 auto Koef2 = int_pow<3>(s1) * 0.5;
94 auto const X = s0 / (dist - s + s0);
95 auto const Xcut = s0 / (ia_params.
gay_berne.
cut - s + s0);
97 auto const X6 = int_pow<6>(X);
98 auto const Xcut6 = int_pow<6>(Xcut);
100 auto const Bra12 = 6 * X6 * X * (2 * X6 - 1);
101 auto const Bra12Cut = 6 * Xcut6 * Xcut * (2 * Xcut6 - 1);
102 auto const Brack = X6 * (X6 - 1);
103 auto const BrackCut = Xcut6 * (Xcut6 - 1);
107 auto const dU_dr = e *
108 (Koef1 * Brhi2 * (Brack - BrackCut) -
109 Koef2 * Brhi1 * (Bra12 - Bra12Cut) - Bra12 * dist / s0) /
112 Koef1 *= chi2 / sqr(dist);
113 Koef2 *= chi1 / sqr(dist);
115 auto const dU_da = e * (Koef1 * (tt1 + tt2) * (BrackCut - Brack) +
116 Koef2 * (oo1 + oo2) * (Bra12 - Bra12Cut));
117 auto const dU_db = e * (Koef1 * (tt2 - tt1) * (Brack - BrackCut) +
118 Koef2 * (oo1 - oo2) * (Bra12 - Bra12Cut));
120 e * ((Brack - BrackCut) * (nu * e1 * sqr(chi1) * uij +
121 0.5 * Koef1 * chi2 * (sqr(tt1) - sqr(tt2))) -
122 (Bra12 - Bra12Cut) * 0.5 * Koef2 * chi1 * (sqr(oo1) - sqr(oo2)));
125 auto const G2 = -dU_da * d - dU_dc * uj;
155 auto const uij = ui * uj;
156 auto const rui = r * ui;
157 auto const ruj = r * uj;
159 auto const o1 = sqr(rui + ruj) / (1 + chi1 * uij);
160 auto const o2 = sqr(rui - ruj) / (1 - chi1 * uij);
161 auto const t1 = sqr(rui + ruj) / (1 + chi2 * uij);
162 auto const t2 = sqr(rui - ruj) / (1 - chi2 * uij);
164 auto const e1 = std::pow(1. - sqr(chi1 * uij), -0.5 * nu);
165 auto const e2 = std::pow(1. - 0.5 * chi2 * (t1 + t2), mu);
166 auto const e = e0 * e1 * e2;
168 auto const s1 = 1. / std::sqrt(1. - 0.5 * chi1 * (o1 + o2));
169 auto const s = s0 * s1;
171 auto r_eff = [=](
double r) {
return (r - s + s0) / s0; };
172 auto E = [=](
double r) {
173 return 4. * e * (int_pow<12>(1. / r) - int_pow<6>(1. / r));
176 return E(r_eff(dist)) - E(r_eff(ia_params.
gay_berne.
cut));
double gb_pair_energy(Utils::Quaternion< double > const &qi, Utils::Quaternion< double > const &qj, IA_parameters const &ia_params, Utils::Vector3d const &d, double dist)
Calculate Gay-Berne energy.
ParticleForce gb_pair_force(Utils::Quaternion< double > const &qi, Utils::Quaternion< double > const &qj, IA_parameters const &ia_params, Utils::Vector3d const &d, double dist)
Calculate Gay-Berne force and torques.
Force information on a particle.