65 static constexpr int num = 3;
80 std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
110 auto const dx = fp2 - fp3;
111 auto const len = dx.
norm();
112 auto const dr = len -
r0;
121 auto const lambda = len /
r0;
122 auto const spring = (std::pow(lambda, 0.5) + std::pow(lambda, -2.5)) /
123 (lambda + std::pow(lambda, -3.));
124 fac -=
ks * spring * dr;
126 auto const f = (fac / len) * dx;
133 auto const dx = fp2 - fp3;
134 auto const len2 = dx.norm2();
135 auto const v_ij = p3.
v() - p2.
v();
161 auto const fac =
kvisc * (dx * v_ij) / len2;
162 auto const f = fac * dx;
173 auto const aa = (phi -
phi0);
175 auto const fac =
kb * aa;
184 auto const BminA = fp3 - fp2;
186 auto const factorFaNc =
187 (fp2 - fp3) * (fp1 - fp3) / BminA.norm() / Nc.norm2();
188 auto const factorFaNd =
189 (fp2 - fp3) * (fp4 - fp3) / BminA.norm() / Nd.norm2();
190 auto const factorFbNc =
191 (fp2 - fp3) * (fp2 - fp1) / BminA.norm() / Nc.norm2();
192 auto const factorFbNd =
193 (fp2 - fp3) * (fp2 - fp4) / BminA.norm() / Nd.norm2();
195 force1 -= fac * BminA.norm() / Nc.norm2() * Nc;
196 force2 += fac * (factorFaNc * Nc + factorFaNd * Nd);
197 force3 += fac * (factorFbNc * Nc + factorFbNd * Nd);
198 force4 -= fac * BminA.norm() / Nd.norm2() * Nd;
214 auto const h = (1. / 3.) * (fp1 + fp2 + fp3);
216 auto const t = sqrt(A / A0) - 1.0;
218 auto const m1 = h - fp1;
219 auto const m2 = h - fp2;
220 auto const m3 = h - fp3;
222 auto const m1_length2 = m1.norm2();
223 auto const m2_length2 = m2.norm2();
224 auto const m3_length2 = m3.norm2();
227 kal * A0 * (2 * t + t * t) / (m1_length2 + m2_length2 + m3_length2);
230 force1 += (fac / 3.0) * m1;
231 force2 += (fac / 3.0) * m2;
232 force3 += (fac / 3.0) * m3;
235 handle_triangle(
kal,
A01, fp1, fp2, fp3, force1, force2, force3);
236 handle_triangle(
kal,
A02, fp2, fp3, fp4, force2, force3, force4);
238 return std::make_tuple(force2, force1, force3, force4);
Vector implementation and trait types for boost qvm interoperability.
auto unfolded_position(Utils::Vector3d const &pos, Utils::Vector3i const &image_box) const
Unfold particle coordinates to image box.
Utils::Vector< T, 3 > get_mi_vector(const Utils::Vector< T, 3 > &a, const Utils::Vector< T, 3 > &b) const
Get the minimum-image vector between two coordinates.
This file contains the defaults for ESPResSo.
#define TINY_OIF_ELASTICITY_COEFFICIENT
Tiny oif elasticity cutoff.
double angle_btw_triangles(const Vector3d &P1, const Vector3d &P2, const Vector3d &P3, const Vector3d &P4)
Computes the angle between two triangles in 3D space.
double area_triangle(const Vector3d &P1, const Vector3d &P2, const Vector3d &P3)
Computes the area of triangle between vectors P1,P2,P3, by computing the cross product P1P2 x P1P3 an...
Vector3d get_n_triangle(const Vector3d &P1, const Vector3d &P2, const Vector3d &P3)
Computes the normal vector of a triangle.
Parameters for OIF local forces.
double A02
Equilibrium surface of the second triangle.
OifLocalForcesBond(double r0, double ks, double kslin, double phi0, double kb, double A01, double A02, double kal, double kvisc)
double phi0
Equilibrium angle between the two triangles.
double r0
Equilibrium bond length of triangle edges.
double kslin
Linear stretching coefficient of triangle edges.
double kb
Bending coefficient for the angle between the two triangles.
double A01
Equilibrium surface of the first triangle.
double kvisc
Viscous coefficient of the triangle vertices.
double ks
Non-linear stretching coefficient of triangle edges.
double kal
Stretching coefficient of a triangle surface.
std::tuple< Utils::Vector3d, Utils::Vector3d, Utils::Vector3d, Utils::Vector3d > calc_forces(BoxGeometry const &box_geo, Particle const &p2, Particle const &p1, Particle const &p3, Particle const &p4) const
Compute the OIF local forces.
Struct holding all information for one particle.
auto const & image_box() const