GCC Code Coverage Report

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