ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
gay_berne.hpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2022 The ESPResSo project
3 * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
4 * Max-Planck-Institute for Polymer Research, Theory Group
5 *
6 * This file is part of ESPResSo.
7 *
8 * ESPResSo is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * ESPResSo is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with this program. If not, see <http://www.gnu.org/licenses/>.
20 */
21#ifndef CORE_NB_IA_GB_HPP
22#define CORE_NB_IA_GB_HPP
23
24/** \file
25 * Routines to calculate the Gay-Berne potential between particle pairs.
26 *
27 * Please note that contrary to the original paper @cite gay81a, here the GB
28 * prefactor @f$ \varepsilon_0 @f$ is not raised to the power of @f$ \nu @f$,
29 * in agreement with the convention used in the GB literature.
30 *
31 * Implementation in \ref gay_berne.cpp.
32 */
33
34#include "config/config.hpp"
35
36#ifdef GAY_BERNE
37
38#include "Particle.hpp"
40
41#include <utils/Vector.hpp>
43#include <utils/math/sqr.hpp>
44#include <utils/quaternion.hpp>
45
46#include <cmath>
47
48/** Calculate Gay-Berne force and torques */
51 IA_parameters const &ia_params,
52 Utils::Vector3d const &d, double dist) {
53 using Utils::int_pow;
54 using Utils::sqr;
55
56 if (dist >= ia_params.gay_berne.cut) {
57 return {};
58 }
59
60 auto const ui = Utils::convert_quaternion_to_director(qi);
61 auto const uj = Utils::convert_quaternion_to_director(qj);
62 auto const e0 = ia_params.gay_berne.eps;
63 auto const s0 = ia_params.gay_berne.sig;
64 auto const chi1 = ia_params.gay_berne.chi1;
65 auto const chi2 = ia_params.gay_berne.chi2;
66 auto const mu = ia_params.gay_berne.mu;
67 auto const nu = ia_params.gay_berne.nu;
68 auto const r = Utils::Vector3d(d).normalize();
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);
84
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);
88
89 auto const s1 = 1. / sqrt(1. - 0.5 * Brhi1);
90 auto const s = s0 * s1;
91 auto Koef1 = mu / e2;
92 auto Koef2 = int_pow<3>(s1) * 0.5;
93
94 auto const X = s0 / (dist - s + s0);
95 auto const Xcut = s0 / (ia_params.gay_berne.cut - s + s0);
96
97 auto const X6 = int_pow<6>(X);
98 auto const Xcut6 = int_pow<6>(Xcut);
99
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);
104
105 /*-------- Here we calculate derivatives -----------------------------*/
106
107 auto const dU_dr = e *
108 (Koef1 * Brhi2 * (Brack - BrackCut) -
109 Koef2 * Brhi1 * (Bra12 - Bra12Cut) - Bra12 * dist / s0) /
110 sqr(dist);
111
112 Koef1 *= chi2 / sqr(dist);
113 Koef2 *= chi1 / sqr(dist);
114
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));
119 auto const dU_dc =
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)));
123
124 // torque = u_i x G
125 auto const G2 = -dU_da * d - dU_dc * uj;
126
127 ParticleForce pf{-dU_dr * d - dU_da * ui - dU_db * uj,
128 vector_product(ui, G2)};
129
130 return pf;
131}
132
133/** Calculate Gay-Berne energy */
136 IA_parameters const &ia_params,
137 Utils::Vector3d const &d, double dist) {
138 using Utils::int_pow;
139 using Utils::sqr;
140
141 if (dist >= ia_params.gay_berne.cut) {
142 return 0.0;
143 }
144
145 auto const ui = Utils::convert_quaternion_to_director(qi);
146 auto const uj = Utils::convert_quaternion_to_director(qj);
147 auto const e0 = ia_params.gay_berne.eps;
148 auto const s0 = ia_params.gay_berne.sig;
149 auto const chi1 = ia_params.gay_berne.chi1;
150 auto const chi2 = ia_params.gay_berne.chi2;
151 auto const mu = ia_params.gay_berne.mu;
152 auto const nu = ia_params.gay_berne.nu;
153 auto const r = Utils::Vector3d(d).normalize();
154
155 auto const uij = ui * uj;
156 auto const rui = r * ui;
157 auto const ruj = r * uj;
158
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);
163
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;
167
168 auto const s1 = 1. / std::sqrt(1. - 0.5 * chi1 * (o1 + o2));
169 auto const s = s0 * s1;
170
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));
174 };
175
176 return E(r_eff(dist)) - E(r_eff(ia_params.gay_berne.cut));
177}
178
179#endif // GAY_BERNE
180#endif
Vector implementation and trait types for boost qvm interoperability.
This file contains the defaults for ESPResSo.
__device__ void vector_product(float const *a, float const *b, float *out)
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.
Definition gay_berne.hpp:49
VectorXd< 3 > Vector3d
Definition Vector.hpp:164
Vector< T, 3 > convert_quaternion_to_director(Quaternion< T > const &quat)
Convert quaternion to director.
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Definition sqr.hpp:28
DEVICE_QUALIFIER constexpr T int_pow(T x)
Calculate integer powers.
Definition int_pow.hpp:61
Various procedures concerning interactions between particles.
Quaternion implementation and trait types for boost qvm interoperability.
Parameters for non-bonded interactions.
GayBerne_Parameters gay_berne
Force information on a particle.
Definition Particle.hpp:290
Quaternion representation.