GCC Code Coverage Report

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