GCC Code Coverage Report

Directory: ./
File: src/potentials.cpp
Date: 2026-04-08 15:34:21
Exec Total Coverage
Lines: 512 565 90.6%
Functions: 72 82 87.8%
Branches: 310 561 55.3%
Line Branch Exec Source
1 #include <cmath>
2 #include <iostream>
3 #include <iterator>
4 #include <string>
5
6 #include "particle.hpp"
7 #include "potentials.hpp"
8 #include "system.hpp"
9 #include "types.hpp"
10
11 // --- Factory ---
12
13 #define RegisterExternal(A) \
14 { \
15 #A, [](System *s, const YAML::Node &c) -> External * { return new A(s, c); } \
16 }
17
18 std::map<std::string, External *(*)(System *s, const YAML::Node &c)> PotentialFactory::externalMap =
19 {
20
1/2
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
 20 RegisterExternal(TestExternal),
21
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 RegisterExternal(Const),
22 ✗ RegisterExternal(PiecewiseLinear),
23
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 RegisterExternal(WallsHard),
24
1/2
✓ Branch 2 taken 28 times.
✗ Branch 3 not taken.
 28 RegisterExternal(WallLJ),
25
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 RegisterExternal(MovingInterface),
26
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 RegisterExternal(FixedLJ),
27 ✗ RegisterExternal(FixedHard),
28
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 RegisterExternal(HarmonicTrapPlanar),
29 ✗ RegisterExternal(Sine),
30
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 RegisterExternal(ExternalTable),
31 };
32
33 #undef RegisterExternal
34
35 #define RegisterInternal(A) \
36 { \
37 #A, [](System *s, const YAML::Node &c, bool newton3) -> Internal * { \
38 return new A(s, c, newton3); \
39 } \
40 }
41
42 std::map<std::string, Internal *(*)(System * s, const YAML::Node &c, bool newton3)>
43 PotentialFactory::internalMap = {
44
1/2
✓ Branch 2 taken 53 times.
✗ Branch 3 not taken.
 106 RegisterInternal(TestInternal), RegisterInternal(Ideal),
45
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
 6 RegisterInternal(Hard), RegisterInternal(SquareShoulder),
46
1/2
✓ Branch 2 taken 26 times.
✗ Branch 3 not taken.
 52 RegisterInternal(LJ), RegisterInternal(WCA),
47
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 4 RegisterInternal(Yukawa), RegisterInternal(Stockmayer),
48
1/2
✓ Branch 2 taken 22 times.
✗ Branch 3 not taken.
 44 RegisterInternal(GayBerne), RegisterInternal(SW),
49
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 RegisterInternal(SWOptimized),
50 };
51
52 #undef RegisterInternal
53
54 54 External *PotentialFactory::createExternal(System *s, const YAML::Node &config) {
55
3/6
✓ Branch 2 taken 54 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 54 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 54 times.
✗ Branch 8 not taken.
 108 std::string externalName = config.begin()->first.as<std::string>();
56
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 const YAML::Node &c = config[externalName];
57
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 if (auto it = externalMap.find(externalName); it != externalMap.end()) {
58
2/4
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 54 times.
✗ Branch 5 not taken.
 54 std::clog << "\n\ttype: " << externalName;
59
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 return it->second(s, c);
60 } else {
61 ✗ throw std::invalid_argument("No external interaction of type " + externalName + " found");
62 }
63 108 }
64
65 109 std::vector<External *> PotentialFactory::createExternals(System *s, const YAML::Node &config) {
66
2/2
✓ Branch 0 taken 56 times.
✓ Branch 1 taken 53 times.
 109 std::vector<External *> externals;
67
2/2
✓ Branch 0 taken 56 times.
✓ Branch 1 taken 53 times.
 109 if (!config) {
68 56 return {};
69 }
70
1/2
✓ Branch 1 taken 53 times.
✗ Branch 2 not taken.
 53 std::clog << "\nExternal interaction:";
71
3/4
✓ Branch 1 taken 53 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 52 times.
 53 if (config.IsSequence()) {
72
5/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✓ Branch 11 taken 1 times.
 5 for (auto &c : config) {
73
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 externals.push_back(createExternal(s, c));
74
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 4 }
75 } else {
76
2/4
✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 52 times.
✗ Branch 5 not taken.
 52 externals.push_back(createExternal(s, config));
77 }
78 53 return externals;
79 109 }
80
81 108 Internal *PotentialFactory::createInternal(System *s, const YAML::Node &config, bool newton3) {
82
3/6
✓ Branch 2 taken 108 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 108 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 108 times.
✗ Branch 8 not taken.
 216 std::string internalName = config.begin()->first.as<std::string>();
83
1/2
✓ Branch 1 taken 108 times.
✗ Branch 2 not taken.
 108 const YAML::Node &c = config[internalName];
84
1/2
✓ Branch 1 taken 108 times.
✗ Branch 2 not taken.
 108 if (auto it = internalMap.find(internalName); it != internalMap.end()) {
85 108 std::clog << "\nInternal interaction:"
86
3/6
✓ Branch 1 taken 108 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 108 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 108 times.
✗ Branch 8 not taken.
 108 << "\n\ttype: " << internalName;
87
1/2
✓ Branch 1 taken 108 times.
✗ Branch 2 not taken.
 108 return it->second(s, c, newton3);
88 } else {
89 ✗ throw std::invalid_argument("No internal interaction of type " + internalName + " found");
90 }
91 216 }
92
93 // --- External ---
94
95 54 External::External(int props) : props(props) {}
96
97 52 bool External::doesE() { return props & DOES_E; }
98
99 52 bool External::doesF() { return props & DOES_F; }
100
101 // This just counts the interactions that would be calculated
102 20 TestExternal::TestExternal(System *s, const YAML::Node &c) : External(DOES_E | DOES_F) {}
103
104 2976 double TestExternal::operator()(Particle &p) {
105 #ifdef FORCE_DEBUG
106 for (int d = 0; d < p.r.size(); ++d) {
107 std::clog << p.r[d] << " ";
108 }
109 std::clog << std::endl;
110 #endif
111 2976 return 1;
112 }
113
114 1 Const::Const(System *s, const YAML::Node &c)
115
5/10
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
 3 : External(DOES_F), d(YAML::parse<int>(c, "d")), mag(YAML::parse<double>(c, "mag")) {
116
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 std::clog << "\n\tmag: " << mag;
117 1 }
118
119 1 double Const::operator()(Particle &p) {
120 1 p.Fext[d] += mag;
121 1 return 0;
122 }
123
124 ✗ PiecewiseLinear::PiecewiseLinear(System *s, const YAML::Node &c)
125 ✗ : External(DOES_E | DOES_F), d(YAML::parse<int>(c, "d")), pos1(YAML::parse<double>(c, "pos1")),
126 ✗ pos2(YAML::parse<double>(c, "pos2")), posDeltaInv(1. / (pos2 - pos1)),
127 ✗ E1(YAML::parse<double>(c, "E1")), E2(YAML::parse<double>(c, "E2")),
128 ✗ F((E1 - E2) * posDeltaInv) {
129 ✗ std::clog << "\n\tpos1: " << pos1;
130 ✗ std::clog << "\n\tpos2: " << pos2;
131 ✗ std::clog << "\n\tE1: " << E1;
132 ✗ std::clog << "\n\tE2: " << E2;
133 }
134
135 ✗ double PiecewiseLinear::operator()(Particle &p) {
136 ✗ Real pos = p.r[d];
137 ✗ if (pos >= pos1 && pos <= pos2) {
138 ✗ p.Fext[d] += F;
139 ✗ Real alpha = (pos2 - pos) * posDeltaInv;
140 ✗ return alpha * E1 + (1 - alpha) * E2;
141 }
142 return 0;
143 }
144
145 1 WallsHard::WallsHard(System *s, const YAML::Node &c)
146
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 2 : External(DOES_E), L(s->L), d(YAML::parse<int>(c, "d")),
147
4/8
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
 3 width(YAML::parse<double>(c, "width")) {
148
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 std::clog << "\n\twidth: " << width;
149 1 }
150
151 1 double WallsHard::operator()(Particle &p) {
152
1/3
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 5 not taken.
 1 if (p.r[d] < width || p.r[d] > L[d] - width) {
153 1 return INFINITY;
154 }
155 return 0.;
156 }
157
158 28 WallLJ::WallLJ(System *s, const YAML::Node &c)
159
2/4
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 28 times.
✗ Branch 6 not taken.
 56 : External(DOES_E | DOES_F), L(s->L), d(YAML::parse<int>(c, "d")),
160
4/8
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 28 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 28 times.
✗ Branch 12 not taken.
 56 cutoff(YAML::parse<double>(c, "cutoff")), pos([this, c]() {
161 28 Vec result = Vec::Zero();
162
1/2
✓ Branch 2 taken 28 times.
✗ Branch 3 not taken.
 56 result[d] = YAML::parse<double>(c, "pos", 0.);
163 28 return result;
164
1/2
✓ Branch 2 taken 28 times.
✗ Branch 3 not taken.
 28 }()) {
165
2/4
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
 28 std::clog << "\n\td: " << d;
166
2/4
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
 28 std::clog << "\n\tcutoff: " << cutoff;
167
2/4
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 28 times.
✗ Branch 6 not taken.
 28 std::clog << "\n\tpos: " << pos[d];
168 28 }
169
170 5965 double WallLJ::operator()(Particle &p) {
171 5965 Real dist = -r_ij(L, p.r, pos)[d];
172
2/2
✓ Branch 0 taken 2829 times.
✓ Branch 1 taken 3136 times.
 5965 if (std::abs(dist) >= cutoff) {
173 return 0;
174 }
175 2829 Real distsqinv = 1. / (dist * dist);
176 2829 Real dist6inv = distsqinv * distsqinv * distsqinv;
177 2829 p.Fext[d] += 24. * dist6inv * (2. * dist6inv - 1.) * distsqinv * dist;
178 2829 return 4. * dist6inv * (dist6inv - 1.);
179 }
180
181 1 MovingInterface::MovingInterface(System *s, const YAML::Node &c)
182
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 2 : External(DOES_E | DOES_F), t(s->t), d(YAML::parse<int>(c, "d")),
183
5/10
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
✓ Branch 15 taken 1 times.
✗ Branch 16 not taken.
 3 B(YAML::parse<double>(c, "B")), kappa(YAML::parse<double>(c, "kappa")),
184
7/14
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 1 times.
✗ Branch 18 not taken.
✓ Branch 21 taken 1 times.
✗ Branch 22 not taken.
✓ Branch 24 taken 1 times.
✗ Branch 25 not taken.
 6 v(YAML::parse<double>(c, "v")), posStart(YAML::parse<double>(c, "posStart", s->L[d])) {
185
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 std::clog << "\n\td: " << d;
186
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 std::clog << "\n\tB: " << B;
187
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 std::clog << "\n\tkappa: " << kappa;
188
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 std::clog << "\n\tv: " << v;
189
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 std::clog << "\n\tposStart: " << posStart;
190 1 }
191
192 2 double MovingInterface::operator()(Particle &p) {
193 2 double R = 0.5 * p.sigma;
194 2 double pos = posStart - v * t;
195 2 double leftexp = exp(-kappa * (p.r[d] - R));
196 2 double rightexp = exp(kappa * (p.r[d] + R - pos));
197 2 p.Fext[d] += B * kappa * (leftexp - rightexp);
198 2 return B * (leftexp + rightexp);
199 }
200
201 1 FixedLJ::FixedLJ(System *s, const YAML::Node &c)
202
3/6
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
 2 : External(DOES_E | DOES_F), L(s->L), pos(L / 2.), cutoff(YAML::parse<double>(c, "cutoff")) {
203
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 std::clog << "\n\tcutoff: " << cutoff;
204 1 }
205
206 1 double FixedLJ::operator()(Particle &p) {
207 1 Vec rVec = r_ij(L, p.r, pos);
208 1 Real rsq = rVec.squaredNorm();
209
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (rsq >= cutoff * cutoff) {
210 return 0.;
211 }
212 1 Real rsqinv = 1. / rsq;
213 1 Real r6inv = rsqinv * rsqinv * rsqinv;
214 1 p.Fext += -24. * r6inv * (2. * r6inv - 1.) * rsqinv * rVec;
215 1 return 4 * r6inv * (r6inv - 1.);
216 }
217
218 ✗ FixedHard::FixedHard(System *s, const YAML::Node &c)
219 ✗ : External(DOES_E), L(s->L), pos(L / 2.), sigma(YAML::parse<double>(c, "sigma")) {
220 ✗ std::clog << "\n\tsigma: " << sigma;
221 }
222
223 ✗ double FixedHard::operator()(Particle &p) {
224 ✗ Vec rVec = r_ij(L, p.r, pos);
225 ✗ Real rsq = rVec.squaredNorm();
226 ✗ if (rsq >= sigma * sigma) {
227 ✗ return 0.;
228 }
229 return INFINITY;
230 }
231
232
233 1 HarmonicTrapPlanar::HarmonicTrapPlanar(System *s, const YAML::Node &c)
234
4/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
 3 : External(DOES_E | DOES_F), d(YAML::parse<int>(c, "d")), A(YAML::parse<double>(c, "A")),
235
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 2 pos(s->L[d] / 2.) {
236
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 std::clog << "\n\td: " << d;
237
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 std::clog << "\n\tA: " << A;
238 1 }
239
240 1 double HarmonicTrapPlanar::operator()(Particle &p) {
241 1 Real dist = p.r[d] - pos;
242 1 Real tmp = A * dist;
243 1 p.Fext[d] += -2. * tmp;
244 1 return tmp * dist;
245 }
246
247 ✗ Sine::Sine(System *s, const YAML::Node &c)
248 ✗ : External(DOES_E | DOES_F), dMod(YAML::parse<int>(c, "dMod")),
249 ✗ dForce(YAML::parse<int>(c, "dForce")), L(s->L[dMod]), A(YAML::parse<double>(c, "A")),
250 ✗ k(YAML::parse<double>(c, "k", 1)), phase(YAML::parse<double>(c, "phase", 0)),
251 ✗ Linv2pik(2 * M_PI * k / L), ALby2pik(A * L / (2 * M_PI * k)) {
252 ✗ std::clog << "\n\tdMod: " << dMod;
253 ✗ std::clog << "\n\tdForce: " << dForce;
254 ✗ std::clog << "\n\tk: " << k;
255 ✗ std::clog << "\n\tA: " << A;
256 ✗ std::clog << "\n\tphase: " << phase;
257 }
258
259 ✗ double Sine::operator()(Particle &p) {
260 ✗ Real arg = p.r[dMod] * Linv2pik + phase;
261 ✗ Real sin = std::sin(arg);
262 ✗ Real cos = std::cos(arg);
263 ✗ p.Fext[dForce] += A * sin;
264 ✗ return dForce == dMod ? ALby2pik * cos : 0;
265 }
266
267 1 ExternalTable::ExternalTable(System *s, const YAML::Node &c)
268
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 : External(DOES_F), bins(YAML::parse<ArrayI>(c, "bins")),
269
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 2 dr(s->L.array() / bins.cast<double>()) {
270 1 int linsize = 1;
271
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.
 4 for (int i = bins.size() - 1; i >= 0; --i) {
272 3 this->stride[i] = linsize;
273 3 linsize *= bins[i];
274 }
275
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 std::vector<double> Fx, Fy, Fz;
276
4/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
 3 Fx = YAML::parse<std::vector<double>>(c, "fextx", std::vector<double>(bins.prod(), 0.));
277
4/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
 3 Fy = YAML::parse<std::vector<double>>(c, "fexty", std::vector<double>(bins.prod(), 0.));
278
2/4
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
 1 if ((int)Fx.size() != bins.prod() || Fx.size() != Fy.size()) {
279 ✗ throw std::invalid_argument("Wrong number of entries in ExternalTable");
280 }
281 1 if constexpr (DIM == 3) {
282
4/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
 3 Fz = YAML::parse<std::vector<double>>(c, "fextz", std::vector<double>(bins.prod(), 0.));
283
2/4
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 if (Fx.size() != Fz.size() || Fy.size() != Fz.size()) {
284 ✗ throw std::invalid_argument("Wrong number of entries in ExternalTable");
285 }
286 }
287 1 Fext.clear();
288
1/2
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 1 Fext.resize(Fx.size(), Vec::Constant(0.));
289
2/2
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 1 times.
 25 for (size_t i = 0; i < Fx.size(); ++i) {
290 24 Fext[i][0] = Fx[i];
291 24 Fext[i][1] = Fy[i];
292 24 if constexpr (DIM == 3) {
293 24 Fext[i][2] = Fz[i];
294 }
295 }
296 1 }
297
298 1 double ExternalTable::operator()(Particle &p) {
299 1 p.Fext += Fext[index(p.r)];
300 1 return 0.;
301 }
302
303 1 int ExternalTable::index(const Vec &r) {
304 1 return ((r.array() / dr.array()).cast<int>() * stride).sum();
305 }
306
307 // --- Internal ---
308
309 108 Internal::Internal(int props, double range, bool newton3, System *s)
310 108 : props(props), range(range), newton3(newton3), L(s->L) {}
311
312 // When implementing a new force, pay attention to the following aspects:
313 // - Register the force in the corresponding map (see top of file)
314 // - Pass correct properties to the base class constructor (DOES_E, DOES_F, THREEBODY and range)
315 // - You can assume that only distinct particles will be given to operator()
316
317 // This just counts the interactions that would be calculated
318 18 TestInternal::TestInternal(System *s, const YAML::Node &c, bool newton3)
319 : Internal(DOES_E | DOES_F | THREEBODY, YAML::parse<double>(c, "range", INFINITY), newton3, s),
320
4/6
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 18 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 8 times.
✓ Branch 8 taken 10 times.
 36 doThreebody(YAML::parse<bool>(c, "doThreebody", true)),
321
5/8
✓ Branch 2 taken 18 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 18 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 18 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 8 times.
✓ Branch 11 taken 10 times.
 54 rangesq(std::isfinite(range) ? range * range : INFINITY) {
322
1/2
✓ Branch 0 taken 18 times.
✗ Branch 1 not taken.
 18 if (doThreebody) {
323 18 range *= 2;
324 }
325 18 }
326
327 589515 double TestInternal::operator()(Particle &p1, Particle &p2) {
328 589515 Vec rVec = r_ij(L, p1.r, p2.r);
329 589515 Real rsq = rVec.squaredNorm();
330
2/2
✓ Branch 0 taken 562830 times.
✓ Branch 1 taken 26685 times.
 589515 if (rsq >= rangesq) {
331 562830 return 0;
332 }
333 #ifdef FORCE_DEBUG
334 for (int d = 0; d < p1.r.size(); ++d) {
335 std::clog << p1.r[d] << " ";
336 }
337 std::clog << "with ";
338 for (int d = 0; d < p2.r.size(); ++d) {
339 std::clog << p2.r[d] << " ";
340 }
341 std::clog << std::endl;
342 #endif
343 return 1;
344 }
345
346 106703600 double TestInternal::operator()(Particle &p1, Particle &p2, Particle &p3) {
347
1/2
✓ Branch 0 taken 106703600 times.
✗ Branch 1 not taken.
 106703600 if (!doThreebody) {
348 return 0;
349 }
350 106703600 Vec r12Vec = r_ij(L, p1.r, p2.r);
351 106703600 Vec r13Vec = r_ij(L, p1.r, p3.r);
352 106703600 Real r12sq = r12Vec.squaredNorm();
353 106703600 Real r13sq = r13Vec.squaredNorm();
354
4/4
✓ Branch 0 taken 3802166 times.
✓ Branch 1 taken 102901434 times.
✓ Branch 2 taken 3609706 times.
✓ Branch 3 taken 192460 times.
 106703600 if (r12sq > rangesq || r13sq > rangesq) {
355 106511140 return 0.;
356 }
357 #ifdef DEBUG
358 for (int d = 0; d < p1.r.size(); ++d) {
359 std::clog << p1.r[d] << " ";
360 }
361 std::clog << "with ";
362 for (int d = 0; d < p2.r.size(); ++d) {
363 std::clog << p2.r[d] << " ";
364 }
365 std::clog << "with ";
366 for (int d = 0; d < p3.r.size(); ++d) {
367 std::clog << p3.r[d] << " ";
368 }
369 std::clog << std::endl;
370 #endif
371 return 1;
372 }
373
374 35 Ideal::Ideal(System *s, const YAML::Node &c, bool newton3)
375 35 : Internal(DOES_E | DOES_F, 1., newton3, s) {}
376
377 2420 double Ideal::operator()(Particle &p1, Particle &p2) { return 0; }
378
379 2 Hard::Hard(System *s, const YAML::Node &c, bool newton3) : Internal(DOES_E, 1., newton3, s) {}
380
381 2 double Hard::operator()(Particle &p1, Particle &p2) {
382 2 Vec rVec = r_ij(L, p1.r, p2.r);
383 2 Real rsq = rVec.squaredNorm();
384 2 Real rmin = 0.5 * (p1.sigma + p2.sigma);
385
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 if (rsq >= rmin * rmin) {
386 1 return 0.;
387 }
388 return INFINITY;
389 }
390
391 1 SquareShoulder::SquareShoulder(System *s, const YAML::Node &c, bool newton3)
392
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 : Internal(DOES_E, YAML::parse<double>(c, "cutoff", 2.0), newton3, s), rangesq(range * range),
393
5/10
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
✓ Branch 15 taken 1 times.
✗ Branch 16 not taken.
 3 epsilon(YAML::parse<double>(c, "epsilon", 1.0)) {
394
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 std::clog << "\n\tcutoff: " << range;
395
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 std::clog << "\n\tepsilon: " << epsilon;
396 1 }
397
398 ✗ double SquareShoulder::operator()(Particle &p1, Particle &p2) {
399 ✗ Vec rVec = r_ij(L, p1.r, p2.r);
400 ✗ Real rsq = rVec.squaredNorm();
401 ✗ Real rmin = 0.5 * (p1.sigma + p2.sigma);
402 ✗ if (rsq <= rmin * rmin) {
403 return INFINITY;
404 }
405 ✗ if (rsq <= rangesq) {
406 ✗ return epsilon;
407 }
408 return 0;
409 }
410
411 24 LJ::LJ(System *s, const YAML::Node &c, bool newton3)
412 : Internal(DOES_E, YAML::parse<double>(c, "cutoff", 2.5), newton3, s),
413
3/6
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 24 times.
✗ Branch 9 not taken.
 72 rangesq(range * range), Eshift(YAML::parse<bool>(c, "shift", false)
414
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 13 times.
 24 ? 4. * std::pow(range, -6) * (std::pow(range, -6) - 1.)
415
3/6
✓ Branch 2 taken 24 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 24 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 24 times.
✗ Branch 9 not taken.
 48 : 0.) {
416
2/4
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
 24 std::clog << "\n\tcutoff: " << range;
417
2/4
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
 24 std::clog << "\n\tEshift: " << Eshift;
418
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 13 times.
 24 if (std::abs(Eshift) > 0.) {
419 11 props ^= DOES_F;
420 }
421 24 }
422
423 425625 double LJ::operator()(Particle &p1, Particle &p2) {
424 425625 Vec rVec = r_ij(L, p1.r, p2.r);
425 425625 Real rsq = rVec.squaredNorm();
426 425625 Real sigmamix = 0.5 * (p1.sigma + p2.sigma);
427 425625 double sigmamixsq = sigmamix * sigmamix;
428
2/2
✓ Branch 0 taken 42763 times.
✓ Branch 1 taken 382862 times.
 425625 if (rsq >= sigmamixsq * rangesq) {
429 return 0.;
430 }
431 42763 Real rsqinv = 1 / rsq;
432 42763 Real srsqinv = sigmamixsq * rsqinv;
433 42763 Real sr6inv = srsqinv * srsqinv * srsqinv;
434 42763 Vec FVec = -24. * sr6inv * (2. * sr6inv - 1.) * rsqinv * rVec;
435 42763 p1.Fint += FVec;
436
2/2
✓ Branch 0 taken 8843 times.
✓ Branch 1 taken 33920 times.
 42763 if (newton3) {
437 8843 p2.Fint -= FVec;
438 }
439 42763 return 4. * sr6inv * (sr6inv - 1.) - Eshift;
440 }
441
442 2 WCA::WCA(System *s, const YAML::Node &c, bool newton3)
443 2 : Internal(DOES_E | DOES_F, std::pow(2., 1. / 6.), newton3, s), rangesq(range * range) {}
444
445 1858 double WCA::operator()(Particle &p1, Particle &p2) {
446 1858 Vec rVec = r_ij(L, p1.r, p2.r);
447 1858 Real rsq = rVec.squaredNorm();
448 1858 Real sigmamix = 0.5 * (p1.sigma + p2.sigma);
449 1858 double sigmamixsq = sigmamix * sigmamix;
450
2/2
✓ Branch 0 taken 576 times.
✓ Branch 1 taken 1282 times.
 1858 if (rsq >= sigmamixsq * rangesq) {
451 return 0.;
452 }
453 576 Real rsqinv = 1 / rsq;
454 576 Real srsqinv = sigmamixsq * rsqinv;
455 576 Real sr6inv = srsqinv * srsqinv * srsqinv;
456 576 Vec FVec = -24. * sr6inv * (2. * sr6inv - 1.) * rsqinv * rVec;
457 576 p1.Fint += FVec;
458
1/2
✓ Branch 0 taken 576 times.
✗ Branch 1 not taken.
 576 if (newton3) {
459 576 p2.Fint -= FVec;
460 }
461 576 return 4. * sr6inv * (sr6inv - 1.) + 1.;
462 }
463
464 1 Yukawa::Yukawa(System *s, const YAML::Node &c, bool newton3)
465 1 : Internal(DOES_E | DOES_F, s->L.minCoeff(), newton3, s),
466
2/4
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
 2 kappa(YAML::parse<double>(c, "kappa", 1.)) {}
467
468 1500 double Yukawa::operator()(Particle &p1, Particle &p2) {
469 1500 Vec rVec = r_ij(L, p1.r, p2.r);
470 1500 Real r = rVec.norm();
471
1/2
✓ Branch 0 taken 1500 times.
✗ Branch 1 not taken.
 1500 if (r >= range) {
472 return 0.;
473 }
474 1500 Real exp = std::exp(kappa * (1. - r));
475 1500 Vec FVec = -exp * (r + 1.) / (r * r * r) * rVec;
476 1500 p1.Fint += FVec;
477
1/2
✓ Branch 0 taken 1500 times.
✗ Branch 1 not taken.
 1500 if (newton3) {
478 1500 p2.Fint -= FVec;
479 }
480 1500 return exp / r;
481 }
482
483 1 Stockmayer::Stockmayer(System *s, const YAML::Node &c, bool newton3)
484
3/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
 2 : Internal(DOES_E | DOES_F, INFINITY, newton3, s), mu(YAML::parse<double>(c, "mu", 1.)),
485
3/6
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
 3 musq(mu * mu), epsilon6(YAML::parse<double>(c, "epsilon6", 0.)) {}
486
487 1512 double Stockmayer::operator()(Particle &p1, Particle &p2) {
488 1512 Vec rVec = r_ij(L, p1.r, p2.r);
489 1512 Real r = rVec.norm();
490 1512 Real rinv = 1. / r;
491 1512 Vec er = -rinv * rVec;
492 1512 Real rsqinv = rinv * rinv;
493 1512 Real r3inv = rsqinv * rinv;
494 1512 Real r4inv = rsqinv * rsqinv;
495 1512 Real r6inv = r4inv * rsqinv;
496 1512 Vec Flj = -24. * r6inv * (2. * r6inv - epsilon6) * rsqinv * rVec;
497 1512 Real Elj = 4. * r6inv * (r6inv - epsilon6);
498 1512 Real u1dotu2 = p1.u.dot(p2.u);
499 1512 Real erdotu1 = er.dot(p1.u);
500 1512 Real erdotu2 = er.dot(p2.u);
501 1512 Vec u1perp = p1.u - erdotu1 * er;
502 1512 Vec u2perp = p2.u - erdotu2 * er;
503 1512 Real u1dotu2minus3erdotu1erdotu2 = u1dotu2 - 3. * erdotu1 * erdotu2;
504 1512 Vec Fdd =
505 1512 3. * musq * r4inv * (u1dotu2minus3erdotu1erdotu2 * er + erdotu2 * u1perp + erdotu1 * u2perp);
506 1512 Real Edd = musq * r3inv * u1dotu2minus3erdotu1erdotu2;
507 1512 Vec FVec = Flj + Fdd;
508 1512 p1.Fint += FVec;
509
1/2
✓ Branch 0 taken 1512 times.
✗ Branch 1 not taken.
 1512 if (newton3) {
510 1512 p2.Fint -= FVec;
511 }
512 1512 return Elj + Edd;
513 }
514
515 2 GayBerne::GayBerne(System *s, const YAML::Node &c, bool newton3)
516
2/4
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
 4 : Internal(DOES_E | DOES_F, YAML::parse<double>(c, "kappa", 4.) + 1., newton3, s),
517 2 rangesq(range * range), kappa(range - 1.),
518
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 4 kappaprime(YAML::parse<double>(c, "kappaprime", 5.)),
519 2 chi((kappa * kappa - 1.) / (kappa * kappa + 1.)),
520
2/4
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
 6 chiprime((std::sqrt(kappaprime) - 1.) / (std::sqrt(kappaprime) + 1.)) {
521 2 if (DIM != 3) {
522 throw std::invalid_argument("Gay-Berne interaction is only implemented for DIM == 3.");
523 }
524 2 }
525
526 1650 double GayBerne::operator()(Particle &p1, Particle &p2) {
527 #if DIM == 3
528 1650 Vec rVec = r_ij(L, p1.r, p2.r);
529 1650 Real rsq = rVec.squaredNorm();
530
2/2
✓ Branch 0 taken 1500 times.
✓ Branch 1 taken 150 times.
 1650 if (rsq >= rangesq) {
531 return 0.;
532 }
533 1500 Real r = std::sqrt(rsq);
534 1500 Real rinv = 1. / r;
535 1500 Vec er = -rVec * rinv;
536 1500 Vec uplus = p1.u + p2.u;
537 1500 Vec uminus = p1.u - p2.u;
538 1500 Real u1dotu2 = p1.u.dot(p2.u);
539 1500 Real erdotuplus = er.dot(uplus);
540 1500 Real erdotuminus = er.dot(uminus);
541 1500 Real erdotuplussq = erdotuplus * erdotuplus;
542 1500 Real erdotuminussq = erdotuminus * erdotuminus;
543 1500 Real chiu1dotu2 = chi * u1dotu2;
544 1500 Real chiprimeu1dotu2 = chiprime * u1dotu2;
545 1500 Real onepluschiu1dotu2inv = 1. / (1. + chiu1dotu2);
546 1500 Real oneminuschiu1dotu2inv = 1. / (1. - chiu1dotu2);
547 1500 Real onepluschiprimeu1dotu2inv = 1. / (1. + chiprimeu1dotu2);
548 1500 Real oneminuschiprimeu1dotu2inv = 1. / (1. - chiprimeu1dotu2);
549 1500 Real xi = chi * (erdotuplussq * onepluschiu1dotu2inv + erdotuminussq * oneminuschiu1dotu2inv);
550 1500 Real xiprime = chiprime * (erdotuplussq * onepluschiprimeu1dotu2inv +
551 1500 erdotuminussq * oneminuschiprimeu1dotu2inv);
552 1500 Real sigma = 1. / std::sqrt(1. - 0.5 * xi);
553 1500 Real epsilon = 1. / std::sqrt(1. - chiu1dotu2 * chiu1dotu2);
554 1500 Real epsilonprime = 1. - 0.5 * xiprime;
555 1500 Real Rinv = 1. / (r - sigma + 1.);
556
2/2
✓ Branch 0 taken 1484 times.
✓ Branch 1 taken 16 times.
 1500 if (Rinv < 0) {
557 return INFINITY;
558 }
559 1484 Real R3inv = Rinv * Rinv * Rinv;
560 1484 Real R6inv = R3inv * R3inv;
561 1484 Real R12inv = R6inv * R6inv;
562 1484 Real Rcutinv = 1. / (range - sigma + 1.);
563 1484 Real Rcut3inv = Rcutinv * Rcutinv * Rcutinv;
564 1484 Real Rcut6inv = Rcut3inv * Rcut3inv;
565 1484 Real Rcut12inv = Rcut6inv * Rcut6inv;
566 1484 Real R12invminusR6inv = R12inv - R6inv;
567 1484 Real Rcut12invminusRcut6inv = Rcut12inv - Rcut6inv;
568 1484 Real withcutR12invminusR6inv = R12invminusR6inv - Rcut12invminusRcut6inv;
569 1484 Real twoR13invminusR7inv = (R12inv + R12invminusR6inv) * Rinv;
570 1484 Real twoRcut13invminusRcut7inv = (Rcut12inv + Rcut12invminusRcut6inv) * Rcutinv;
571 1484 Real withcuttwoR13invminusR7inv = twoR13invminusR7inv - twoRcut13invminusRcut7inv;
572 1484 Real erdotfracplus = erdotuplus * onepluschiu1dotu2inv;
573 1484 Real erdotfracminus = erdotuminus * oneminuschiu1dotu2inv;
574 1484 Real erdotfracprimeplus = erdotuplus * onepluschiprimeu1dotu2inv;
575 1484 Real erdotfracprimeminus = erdotuminus * oneminuschiprimeu1dotu2inv;
576 1484 Real gamma1 = chi * (erdotfracplus + erdotfracminus);
577 1484 Real gamma2 = chi * (erdotfracplus - erdotfracminus);
578 1484 Real gammaprime1 = chiprime * (erdotfracprimeplus + erdotfracprimeminus);
579 1484 Real gammaprime2 = chiprime * (erdotfracprimeplus - erdotfracprimeminus);
580 1484 Real pre = 3. * epsilonprime * withcuttwoR13invminusR7inv * sigma * sigma * sigma;
581 1484 Real preprime = -2. * withcutR12invminusR6inv;
582 1484 Real fourepsilonepsilonprime = 4. * epsilon * epsilonprime;
583 1484 Real fourepsilonepsilonprimesq = fourepsilonepsilonprime * epsilonprime;
584 1484 Real drU = -6. * fourepsilonepsilonprimesq * twoR13invminusR7inv;
585 1484 Real derdotu1U = fourepsilonepsilonprime * (pre * gamma1 + preprime * gammaprime1);
586 1484 Real derdotu2U = fourepsilonepsilonprime * (pre * gamma2 + preprime * gammaprime2);
587 1484 Vec u1perp = p1.u - p1.u.dot(er) * er;
588 1484 Vec u2perp = p2.u - p2.u.dot(er) * er;
589 1484 Vec FVec = -drU * er - rinv * (derdotu1U * u1perp + derdotu2U * u2perp);
590 1484 Real iotaprime = chiprime * chiprime *
591 1484 (erdotfracprimeplus * onepluschiprimeu1dotu2inv -
592 1484 erdotfracprimeminus * oneminuschiprimeu1dotu2inv);
593 1484 Real du1dotu2U = fourepsilonepsilonprime *
594 1484 (chi * chiu1dotu2 * epsilon * epsilon * epsilonprime + iotaprime) *
595 withcutR12invminusR6inv;
596 1484 Vec ercrossu1 = er.cross(p1.u);
597 1484 Vec ercrossu2 = er.cross(p2.u);
598 1484 Vec u1crossu2 = p1.u.cross(p2.u);
599 1484 Vec Tau1Vec = derdotu1U * ercrossu1 - du1dotu2U * u1crossu2;
600 1484 Vec Tau2Vec = derdotu2U * ercrossu2 + du1dotu2U * u1crossu2;
601 1484 p1.Fint += FVec;
602 1484 p1.Tauint += Tau1Vec;
603
1/2
✓ Branch 0 taken 1484 times.
✗ Branch 1 not taken.
 1484 if (newton3) {
604 1484 p2.Fint -= FVec;
605 1484 p2.Tauint += Tau2Vec;
606 }
607 1484 return fourepsilonepsilonprimesq * withcutR12invminusR6inv;
608 #else
609 return 0;
610 #endif
611 }
612
613 20 SW::SW(System *s, const YAML::Node &c, bool newton3)
614
2/4
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 20 times.
✗ Branch 6 not taken.
 40 : Internal(DOES_E | DOES_F | THREEBODY, 2. * YAML::parse<double>(c, "a", 1.8), newton3, s),
615
2/4
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 20 times.
✗ Branch 6 not taken.
 40 a(0.5 * range), A(YAML::parse<double>(c, "A", 7.04955627)),
616
5/10
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 20 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 20 times.
✗ Branch 12 not taken.
✓ Branch 15 taken 20 times.
✗ Branch 16 not taken.
 60 B(YAML::parse<double>(c, "B", 0.6022245584)), lambda(YAML::parse<double>(c, "lambda", 23.15)),
617
3/6
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 20 times.
✗ Branch 9 not taken.
 40 gamma(YAML::parse<double>(c, "gamma", 1.2)),
618
3/6
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 20 times.
 20 costheta0(c["theta0"] ? std::cos(c["theta0"].as<double>() * M_PI / 180.) : -1. / 3.),
619
2/4
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 20 times.
✗ Branch 6 not taken.
 60 asq(a * a) {
620
2/4
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
 20 std::clog << "\n\ta: " << a;
621
2/4
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
 20 std::clog << "\n\tA: " << A;
622
2/4
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
 20 std::clog << "\n\tB: " << B;
623
2/4
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
 20 std::clog << "\n\tlambda: " << lambda;
624
2/4
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
 20 std::clog << "\n\tgamma: " << gamma;
625
4/8
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 20 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 20 times.
✗ Branch 11 not taken.
 20 std::clog << "\n\tcostheta0: " << costheta0 << " (theta0: " << std::acos(costheta0) * 180. / M_PI
626
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 << ")";
627 20 }
628
629 591456 double SW::operator()(Particle &p1, Particle &p2) {
630 591456 Vec rVec = r_ij(L, p1.r, p2.r);
631 591456 Real rsq = rVec.squaredNorm();
632
2/2
✓ Branch 0 taken 25749 times.
✓ Branch 1 taken 565707 times.
 591456 if (rsq >= asq) {
633 return 0.;
634 }
635 25749 Real r = std::sqrt(rsq);
636 25749 Real rsqinv = 1. / rsq;
637 25749 Real r4inv = rsqinv * rsqinv;
638 25749 Real rainv = 1. / (r - a);
639 25749 Real exp = std::exp(rainv);
640 25749 Vec FVec = -A * (4. * B * r4inv + (B * r4inv - 1.) * rainv * rainv * r) * exp * rsqinv * rVec;
641 25749 p1.Fint += FVec;
642
2/2
✓ Branch 0 taken 5738 times.
✓ Branch 1 taken 20011 times.
 25749 if (newton3) {
643 5738 p2.Fint -= FVec;
644 }
645 25749 return A * (B * r4inv - 1.) * exp;
646 }
647
648 106705493 double SW::operator()(Particle &p1, Particle &p2, Particle &p3) {
649 106705493 double E = 0.;
650
2/2
✓ Branch 0 taken 40015827 times.
✓ Branch 1 taken 66689666 times.
 106705493 if (newton3) {
651 40015827 E = threebody(p1, p2, p3);
652 } else {
653 66689666 E = threebody(p1, p2, p3, 1);
654 66689666 threebody(p2, p1, p3, 2);
655 66689666 threebody(p3, p2, p1, 3);
656 }
657 106705493 return E;
658 }
659
660 240084825 double SW::threebody(Particle &p1, Particle &p2, Particle &p3, int update) {
661 240084825 Vec r12Vec = r_ij(L, p1.r, p2.r);
662 240084825 Vec r13Vec = r_ij(L, p1.r, p3.r);
663 240084825 Real r12sq = r12Vec.squaredNorm();
664 240084825 Real r13sq = r13Vec.squaredNorm();
665
4/4
✓ Branch 0 taken 9062195 times.
✓ Branch 1 taken 231022630 times.
✓ Branch 2 taken 407957 times.
✓ Branch 3 taken 8654238 times.
 240084825 if (r12sq >= asq || r13sq >= asq) {
666 return 0.;
667 }
668 407957 Real r12 = std::sqrt(r12sq);
669 407957 Real r13 = std::sqrt(r13sq);
670 407957 Real r12ainv = 1. / (r12 - a);
671 407957 Real r13ainv = 1. / (r13 - a);
672 407957 Real gr12ainv = gamma * r12ainv;
673 407957 Real gr13ainv = gamma * r13ainv;
674 407957 Real r12r13inv = 1. / (r12 * r13);
675 407957 Real costheta = r12Vec.dot(r13Vec) * r12r13inv;
676 407957 Real cosdelta = costheta - costheta0;
677 407957 Real exp1213 = std::exp(gr12ainv + gr13ainv);
678 407957 Real E = lambda * cosdelta * cosdelta * exp1213;
679 407957 Real Frad12 = E * gr12ainv * r12ainv / r12;
680 407957 Real Frad13 = E * gr13ainv * r13ainv / r13;
681 407957 Real Fang = 2. * lambda * cosdelta * exp1213;
682 407957 Real Fang1213 = r12r13inv * Fang;
683 407957 Real cosFang = costheta * Fang;
684 407957 Vec F2Vec = r12Vec * (Frad12 + cosFang / r12sq) - r13Vec * Fang1213;
685 407957 Vec F3Vec = r13Vec * (Frad13 + cosFang / r13sq) - r12Vec * Fang1213;
686
4/4
✓ Branch 0 taken 101069 times.
✓ Branch 1 taken 101904 times.
✓ Branch 2 taken 136776 times.
✓ Branch 3 taken 68208 times.
 407957 switch (update) {
687 101069 case 1:
688 101069 p1.Fint -= F2Vec + F3Vec;
689 101069 break;
690 101904 case 2:
691 101904 p2.Fint += F2Vec;
692 break;
693 136776 case 3:
694 136776 p3.Fint += F3Vec;
695 break;
696 68208 default:
697 68208 p1.Fint -= F2Vec + F3Vec;
698 68208 p2.Fint += F2Vec;
699 68208 p3.Fint += F3Vec;
700 }
701 return E;
702 }
703
704 2 SWOptimized::SWOptimized(System *s, const YAML::Node &c, bool newton3)
705 : Internal(DOES_E | DOES_F | NEEDS_GATHER, YAML::parse<double>(c, "a", 1.8), newton3, s),
706
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
 4 particles(s->particles), a(range), A(YAML::parse<double>(c, "A", 7.04955627)),
707
5/10
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 2 times.
✗ Branch 12 not taken.
✓ Branch 15 taken 2 times.
✗ Branch 16 not taken.
 6 B(YAML::parse<double>(c, "B", 0.6022245584)), lambda(YAML::parse<double>(c, "lambda", 23.15)),
708
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
 4 gamma(YAML::parse<double>(c, "gamma", 1.2)),
709
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2 times.
 2 costheta0(c["theta0"] ? std::cos(c["theta0"].as<double>() * M_PI / 180.) : -1. / 3.),
710
3/6
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
 4 asq(a * a) {
711
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 std::clog << "\n\ta: " << a;
712
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 std::clog << "\n\tA: " << A;
713
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 std::clog << "\n\tB: " << B;
714
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 std::clog << "\n\tlambda: " << lambda;
715
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 std::clog << "\n\tgamma: " << gamma;
716
4/8
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
 2 std::clog << "\n\tcostheta0: " << costheta0 << " (theta0: " << std::acos(costheta0) * 180. / M_PI
717
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 << ")";
718 2 }
719
720 631 void SWOptimized::reset() {
721 631 int N = particles.size();
722 631 h.resize(N);
723 631 s.resize(N);
724 631 T.resize(N);
725
2/2
✓ Branch 1 taken 1893 times.
✓ Branch 2 taken 631 times.
 2524 for (int i = 0; i < N; ++i) {
726 1893 h[i] = 0;
727 }
728
2/2
✓ Branch 0 taken 1893 times.
✓ Branch 1 taken 631 times.
 2524 for (int i = 0; i < N; ++i) {
729 1893 s[i].setZero();
730 }
731
2/2
✓ Branch 0 taken 1893 times.
✓ Branch 1 taken 631 times.
 2524 for (int i = 0; i < N; ++i) {
732 1893 T[i].setZero();
733 }
734 631 }
735
736 1893 double SWOptimized::operator()(Particle &p) {
737 1893 int i = std::distance(&particles.front(), &p);
738 1893 return 0.5 * lambda *
739 1893 (costheta0 * costheta0 * h[i] * h[i] - 2 * costheta0 * s[i].squaredNorm() +
740 1893 (T[i] * T[i]).trace());
741 }
742
743 3304 double SWOptimized::operator()(Particle &p1, Particle &p2) {
744
2/2
✓ Branch 0 taken 1652 times.
✓ Branch 1 taken 1652 times.
 3304 if (gather) {
745 1652 return gatherPhase(p1, p2);
746 }
747 1652 return accumulationPhase(p1, p2);
748 }
749
750 1652 double SWOptimized::accumulationPhase(Particle &p1, Particle &p2) {
751 1652 Vec rVec = r_ij(L, p1.r, p2.r);
752 1652 Real rsq = rVec.squaredNorm();
753
2/2
✓ Branch 0 taken 613 times.
✓ Branch 1 taken 1039 times.
 1652 if (rsq >= asq) {
754 return 0.;
755 }
756 613 int i = std::distance(&particles.front(), &p1);
757 613 int j = std::distance(&particles.front(), &p2);
758 // Two-body
759 613 Real r = std::sqrt(rsq);
760 613 Vec rVecUnit = rVec / r;
761 613 Real rsqinv = 1. / rsq;
762 613 Real r4inv = rsqinv * rsqinv;
763 613 Real rainv = 1. / (r - a);
764 613 Real rainvsq = rainv * rainv;
765 613 Real exp = std::exp(rainv);
766 613 double ETwobody = A * (B * r4inv - 1.) * exp;
767 613 Vec FVecTwobody = -A * (4. * B * r4inv + (B * r4inv - 1.) * rainvsq * r) * exp * rsqinv * rVec;
768 // Correction from three-body
769 613 Real gij = std::exp(gamma * rainv);
770 613 Real gijsq = gij * gij;
771 613 double oneminuscostheta0sq = (1 - costheta0) * (1 - costheta0);
772 613 double ECorr = lambda * gijsq * oneminuscostheta0sq;
773 613 Vec FVecCorr = 2 * lambda * oneminuscostheta0sq * gijsq * gamma * rainvsq * rVecUnit;
774 // Update accumulators for rest of three-body
775 613 double hi = gij;
776 613 Vec si = gij * rVecUnit;
777 613 Eigen::Matrix<double, DIM, DIM> Ti = gij * (rVecUnit * rVecUnit.transpose());
778 613 h[i] += hi;
779 613 s[i] += si;
780 613 T[i] += Ti;
781
1/2
✓ Branch 0 taken 613 times.
✗ Branch 1 not taken.
 613 if (newton3) {
782 613 h[j] += hi;
783 613 s[j] -= si; // rVecUnit is flipped -> change of sign
784 613 T[j] += Ti; // quadratic in rVecUnit -> no change of sign
785 }
786 // Update particle with two-body and correction right away (they do not depend on accumulators)
787 613 Vec FVec = FVecTwobody + FVecCorr;
788 613 p1.Fint += FVec;
789
1/2
✓ Branch 0 taken 613 times.
✗ Branch 1 not taken.
 613 if (newton3) {
790 613 p2.Fint -= FVec;
791 }
792 613 return ETwobody - ECorr;
793 }
794
795 1652 double SWOptimized::gatherPhase(Particle &p1, Particle &p2) {
796 1652 Vec rVec = r_ij(L, p1.r, p2.r);
797 1652 Real rsq = rVec.squaredNorm();
798
2/2
✓ Branch 0 taken 613 times.
✓ Branch 1 taken 1039 times.
 1652 if (rsq >= asq) {
799 return 0.;
800 }
801 613 int i = std::distance(&particles.front(), &p1);
802 613 int j = std::distance(&particles.front(), &p2);
803 613 Real r = std::sqrt(rsq);
804 613 Real rinv = 1 / r;
805 613 Real rsqinv = rinv * rinv;
806 613 Real rainv = 1. / (r - a);
807 613 Real rainvsq = rainv * rainv;
808 613 Real gij = std::exp(gamma * rainv);
809 613 Real dgijdrij = -gamma * rainvsq * gij;
810 613 Real cijcji = lambda * (costheta0 * costheta0 * dgijdrij * (h[i] + h[j]) * rinv -
811 613 2 * costheta0 * (dgijdrij - gij * rinv) * rVec.dot(s[i] - s[j]) * rsqinv +
812 1226 (dgijdrij - 2 * gij * rinv) *
813 613 (rVec.transpose() * (T[i] + T[j]) * rVec).value() * rsqinv * rinv);
814 613 Vec FVec = cijcji * rVec - 2 * lambda * costheta0 * gij * rinv * (s[i] - s[j]) +
815 613 2 * lambda * gij * rsqinv * ((T[i] + T[j]) * rVec);
816 613 p1.Fint += FVec;
817
1/2
✓ Branch 0 taken 613 times.
✗ Branch 1 not taken.
 613 if (newton3) {
818 613 p2.Fint -= FVec;
819 }
820 return 0.;
821 }
822