ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
oif_local_forces.hpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2012-2022 The ESPResSo project
3 *
4 * This file is part of ESPResSo.
5 *
6 * ESPResSo is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * ESPResSo is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20#pragma once
21
22/** \file
23 * Routines to calculate the OIF local forces for a particle quadruple
24 * (two neighboring triangles with common edge).
25 * See @cite dupin07a.
26 */
27
28#include "config/config.hpp"
29
30#include "BoxGeometry.hpp"
31#include "Particle.hpp"
32
33#include <utils/Vector.hpp>
35
36#include <cmath>
37#include <tuple>
38
39/** Parameters for OIF local forces
40 *
41 * Characterize the deformation of two triangles sharing an edge.
42 */
44 /** Equilibrium bond length of triangle edges */
45 double r0;
46 /** Non-linear stretching coefficient of triangle edges */
47 double ks;
48 /** Linear stretching coefficient of triangle edges */
49 double kslin;
50 /** Equilibrium angle between the two triangles */
51 double phi0;
52 /** Bending coefficient for the angle between the two triangles */
53 double kb;
54 /** Equilibrium surface of the first triangle */
55 double A01;
56 /** Equilibrium surface of the second triangle */
57 double A02;
58 /** Stretching coefficient of a triangle surface */
59 double kal;
60 /** Viscous coefficient of the triangle vertices */
61 double kvisc;
62
63 double cutoff() const { return 0.; }
64
65 static constexpr int num = 3;
66
67 OifLocalForcesBond(double r0, double ks, double kslin, double phi0, double kb,
68 double A01, double A02, double kal, double kvisc) {
69 this->phi0 = phi0;
70 this->kb = kb;
71 this->r0 = r0;
72 this->ks = ks;
73 this->kslin = kslin;
74 this->A01 = A01;
75 this->A02 = A02;
76 this->kal = kal;
77 this->kvisc = kvisc;
78 }
79
80 std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
81 calc_forces(BoxGeometry const &box_geo, Particle const &p2,
82 Particle const &p1, Particle const &p3, Particle const &p4) const;
83};
84
85/** Compute the OIF local forces.
86 * See @cite dupin07a, @cite jancigova16a.
87 * @param box_geo Box geometry.
88 * @param p2 Particle of triangle 1.
89 * @param p1 , p3 Particles common to triangle 1 and triangle 2.
90 * @param p4 Particle of triangle 2.
91 * @return forces on @p p1, @p p2, @p p3, @p p4
92 */
96 Particle const &p1, Particle const &p3,
97 Particle const &p4) const {
98
99 // first-fold-then-the-same approach
100 auto const fp2 = box_geo.unfolded_position(p2.pos(), p2.image_box());
101 auto const fp1 = fp2 + box_geo.get_mi_vector(p1.pos(), fp2);
102 auto const fp3 = fp2 + box_geo.get_mi_vector(p3.pos(), fp2);
103 auto const fp4 = fp2 + box_geo.get_mi_vector(p4.pos(), fp2);
104
105 Utils::Vector3d force1{}, force2{}, force3{}, force4{};
106
107 // surface strain constraint
110 auto const dx = fp2 - fp3;
111 auto const len = dx.norm();
112 auto const dr = len - r0;
113 auto fac = 0.;
114 // linear stretching
116 fac -= kslin * dr; // no normalization
117 }
118 // non-linear stretching
120 /** For non-linear stretching, see eq. (19) in @cite dupin07a */
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; // no normalization
125 }
126 auto const f = (fac / len) * dx;
127 force2 += f;
128 force3 -= f;
129 }
130
131 // viscous force
132 if (kvisc > TINY_OIF_ELASTICITY_COEFFICIENT) { // to be implemented....
133 auto const dx = fp2 - fp3;
134 auto const len2 = dx.norm2();
135 auto const v_ij = p3.v() - p2.v();
136
137 // Variant A
138 // Here the force is in the direction of relative velocity btw points
139
140 // Code:
141 // force2 += kvisc * v;
142 // force3 -= kvisc * v;
143
144 // Variant B
145 // Here the force is the projection of relative velocity btw points onto
146 // line btw the points
147
148 // denote p vector between p2 and p3
149 // denote v the velocity difference between the points p2 and p3
150 // denote alpha the angle between p and v
151 // denote x the projected v onto p
152 // cos alpha = |x|/|v|
153 // cos alpha = (v,p)/(|v||p|)
154 // together we get |x|=(v,p)/|p|
155 // also, x is along p, so x = |x|.p/|p|
156 // so x = p/|p| . (v,p)/|p|
157 // altogether x = p . (v,p)/(|p|)^2
158 // |p|^2 is stored in len2
159
160 // Code:
161 auto const fac = kvisc * (dx * v_ij) / len2;
162 auto const f = fac * dx;
163 force2 += f;
164 force3 -= f;
165 }
166
167 /* bending
168 forceT1 is restoring force for triangle p1,p2,p3 and forceT2 restoring
169 force for triangle p2,p3,p4 p1 += forceT1; p2 -= 0.5*forceT1+0.5*forceT2;
170 p3 -= 0.5*forceT1+0.5*forceT2; p4 += forceT2; */
172 auto const phi = Utils::angle_btw_triangles(fp1, fp2, fp3, fp4);
173 auto const aa = (phi - phi0); // no renormalization by phi0, to be
174 // consistent with Krueger and Fedosov
175 auto const fac = kb * aa;
176
177 auto const Nc = Utils::get_n_triangle(
178 fp1, fp2,
179 fp3); // returns (fp2 - fp1)x(fp3 - fp1), thus Nc = (A - C)x(B - C)
180 auto const Nd = Utils::get_n_triangle(
181 fp4, fp3,
182 fp2); // returns (fp3 - fp4)x(fp2 - fp4), thus Nd = (B - D)x(A - D)
183
184 auto const BminA = fp3 - fp2;
185
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();
194
195 force1 -= fac * BminA.norm() / Nc.norm2() * Nc; // Fc
196 force2 += fac * (factorFaNc * Nc + factorFaNd * Nd); // Fa
197 force3 += fac * (factorFbNc * Nc + factorFbNd * Nd); // Fb
198 force4 -= fac * BminA.norm() / Nd.norm2() * Nd; // Fc
199 }
200
201 /* local area
202 * for both triangles, only 1/3 of calculated forces are added, because each
203 * triangle will enter this calculation 3 times (one time per edge).
204 * Proportional distribution of forces, implemented according to
205 * @cite jancigova16a.
206 */
208
209 auto handle_triangle = [](double kal, double A0, Utils::Vector3d const &fp1,
210 Utils::Vector3d const &fp2,
211 Utils::Vector3d const &fp3,
212 Utils::Vector3d &force1, Utils::Vector3d &force2,
213 Utils::Vector3d &force3) {
214 auto const h = (1. / 3.) * (fp1 + fp2 + fp3);
215 auto const A = Utils::area_triangle(fp1, fp2, fp3);
216 auto const t = sqrt(A / A0) - 1.0;
217
218 auto const m1 = h - fp1;
219 auto const m2 = h - fp2;
220 auto const m3 = h - fp3;
221
222 auto const m1_length2 = m1.norm2();
223 auto const m2_length2 = m2.norm2();
224 auto const m3_length2 = m3.norm2();
225
226 auto const fac =
227 kal * A0 * (2 * t + t * t) / (m1_length2 + m2_length2 + m3_length2);
228
229 // local area force for p1
230 force1 += (fac / 3.0) * m1;
231 force2 += (fac / 3.0) * m2;
232 force3 += (fac / 3.0) * m3;
233 };
234
235 handle_triangle(kal, A01, fp1, fp2, fp3, force1, force2, force3);
236 handle_triangle(kal, A02, fp2, fp3, fp4, force2, force3, force4);
237 }
238 return std::make_tuple(force2, force1, force3, force4);
239}
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.
T norm() const
Definition Vector.hpp:138
This file contains the defaults for ESPResSo.
#define TINY_OIF_ELASTICITY_COEFFICIENT
Tiny oif elasticity cutoff.
Definition config.hpp:83
VectorXd< 3 > Vector3d
Definition Vector.hpp:164
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.
static constexpr int num
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.
Definition Particle.hpp:395
auto const & v() const
Definition Particle.hpp:433
auto const & image_box() const
Definition Particle.hpp:444
auto const & pos() const
Definition Particle.hpp:431