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 | 425586 | double LJ::operator()(Particle &p1, Particle &p2) { | |
385 | 425586 | Vec rVec = r_ij(L, p1.r, p2.r); | |
386 | 425586 | Real rsq = rVec.squaredNorm(); | |
387 | 425586 | Real sigmamix = 0.5 * (p1.sigma + p2.sigma); | |
388 | 425586 | double sigmamixsq = sigmamix * sigmamix; | |
389 |
2/2✓ Branch 0 taken 42776 times.
✓ Branch 1 taken 382810 times.
|
425586 | if (rsq >= sigmamixsq * rangesq) { |
390 | return 0.; | ||
391 | } | ||
392 | 42776 | Real rsqinv = 1 / rsq; | |
393 | 42776 | Real srsqinv = sigmamixsq * rsqinv; | |
394 | 42776 | Real sr6inv = srsqinv * srsqinv * srsqinv; | |
395 | 42776 | Vec FVec = -24. * sr6inv * (2. * sr6inv - 1.) * rsqinv * rVec; | |
396 | 42776 | p1.Fint += FVec; | |
397 |
2/2✓ Branch 0 taken 8843 times.
✓ Branch 1 taken 33933 times.
|
42776 | if (newton3) { |
398 | 8843 | p2.Fint -= FVec; | |
399 | } | ||
400 | 42776 | 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 | 2040 | double WCA::operator()(Particle &p1, Particle &p2) { | |
407 | 2040 | Vec rVec = r_ij(L, p1.r, p2.r); | |
408 | 2040 | Real rsq = rVec.squaredNorm(); | |
409 | 2040 | Real sigmamix = 0.5 * (p1.sigma + p2.sigma); | |
410 | 2040 | double sigmamixsq = sigmamix * sigmamix; | |
411 |
2/2✓ Branch 0 taken 572 times.
✓ Branch 1 taken 1468 times.
|
2040 | if (rsq >= sigmamixsq * rangesq) { |
412 | return 0.; | ||
413 | } | ||
414 | 572 | Real rsqinv = 1 / rsq; | |
415 | 572 | Real srsqinv = sigmamixsq * rsqinv; | |
416 | 572 | Real sr6inv = srsqinv * srsqinv * srsqinv; | |
417 | 572 | Vec FVec = -24. * sr6inv * (2. * sr6inv - 1.) * rsqinv * rVec; | |
418 | 572 | p1.Fint += FVec; | |
419 |
1/2✓ Branch 0 taken 572 times.
✗ Branch 1 not taken.
|
572 | if (newton3) { |
420 | 572 | p2.Fint -= FVec; | |
421 | } | ||
422 | 572 | 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 | 1527 | double Stockmayer::operator()(Particle &p1, Particle &p2) { | |
449 | 1527 | Vec rVec = r_ij(L, p1.r, p2.r); | |
450 | 1527 | Real r = rVec.norm(); | |
451 | 1527 | Real rinv = 1. / r; | |
452 | 1527 | Vec er = -rinv * rVec; | |
453 | 1527 | Real rsqinv = rinv * rinv; | |
454 | 1527 | Real r3inv = rsqinv * rinv; | |
455 | 1527 | Real r4inv = rsqinv * rsqinv; | |
456 | 1527 | Real r6inv = r4inv * rsqinv; | |
457 | 1527 | Vec Flj = -24. * r6inv * (2. * r6inv - epsilon6) * rsqinv * rVec; | |
458 | 1527 | Real Elj = 4. * r6inv * (r6inv - epsilon6); | |
459 | 1527 | Real u1dotu2 = p1.u.dot(p2.u); | |
460 | 1527 | Real erdotu1 = er.dot(p1.u); | |
461 | 1527 | Real erdotu2 = er.dot(p2.u); | |
462 | 1527 | Vec u1perp = p1.u - erdotu1 * er; | |
463 | 1527 | Vec u2perp = p2.u - erdotu2 * er; | |
464 | 1527 | Real u1dotu2minus3erdotu1erdotu2 = u1dotu2 - 3. * erdotu1 * erdotu2; | |
465 | 1527 | Vec Fdd = | |
466 | 1527 | 3. * musq * r4inv * (u1dotu2minus3erdotu1erdotu2 * er + erdotu2 * u1perp + erdotu1 * u2perp); | |
467 | 1527 | Real Edd = musq * r3inv * u1dotu2minus3erdotu1erdotu2; | |
468 | 1527 | Vec FVec = Flj + Fdd; | |
469 | 1527 | p1.Fint += FVec; | |
470 |
1/2✓ Branch 0 taken 1527 times.
✗ Branch 1 not taken.
|
1527 | if (newton3) { |
471 | 1527 | p2.Fint -= FVec; | |
472 | } | ||
473 | 1527 | 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 | 1641 | double GayBerne::operator()(Particle &p1, Particle &p2) { | |
488 | 1641 | if constexpr (DIM == 3) { | |
489 | 1641 | Vec rVec = r_ij(L, p1.r, p2.r); | |
490 | 1641 | Real rsq = rVec.squaredNorm(); | |
491 |
2/2✓ Branch 0 taken 1491 times.
✓ Branch 1 taken 150 times.
|
1641 | if (rsq >= rangesq) { |
492 | return 0.; | ||
493 | } | ||
494 | 1491 | Real r = std::sqrt(rsq); | |
495 | 1491 | Real rinv = 1. / r; | |
496 | 1491 | Vec er = -rVec * rinv; | |
497 | 1491 | Vec uplus = p1.u + p2.u; | |
498 | 1491 | Vec uminus = p1.u - p2.u; | |
499 | 1491 | Real u1dotu2 = p1.u.dot(p2.u); | |
500 | 1491 | Real erdotuplus = er.dot(uplus); | |
501 | 1491 | Real erdotuminus = er.dot(uminus); | |
502 | 1491 | Real erdotuplussq = erdotuplus * erdotuplus; | |
503 | 1491 | Real erdotuminussq = erdotuminus * erdotuminus; | |
504 | 1491 | Real chiu1dotu2 = chi * u1dotu2; | |
505 | 1491 | Real chiprimeu1dotu2 = chiprime * u1dotu2; | |
506 | 1491 | Real onepluschiu1dotu2inv = 1. / (1. + chiu1dotu2); | |
507 | 1491 | Real oneminuschiu1dotu2inv = 1. / (1. - chiu1dotu2); | |
508 | 1491 | Real onepluschiprimeu1dotu2inv = 1. / (1. + chiprimeu1dotu2); | |
509 | 1491 | Real oneminuschiprimeu1dotu2inv = 1. / (1. - chiprimeu1dotu2); | |
510 | 1491 | Real xi = chi * (erdotuplussq * onepluschiu1dotu2inv + erdotuminussq * oneminuschiu1dotu2inv); | |
511 | 1491 | Real xiprime = chiprime * (erdotuplussq * onepluschiprimeu1dotu2inv + | |
512 | 1491 | erdotuminussq * oneminuschiprimeu1dotu2inv); | |
513 | 1491 | Real sigma = 1. / std::sqrt(1. - 0.5 * xi); | |
514 | 1491 | Real epsilon = 1. / std::sqrt(1. - chiu1dotu2 * chiu1dotu2); | |
515 | 1491 | Real epsilonprime = 1. - 0.5 * xiprime; | |
516 | 1491 | Real Rinv = 1. / (r - sigma + 1.); | |
517 |
2/2✓ Branch 0 taken 1479 times.
✓ Branch 1 taken 12 times.
|
1491 | if (Rinv < 0) { |
518 | return INFINITY; | ||
519 | } | ||
520 | 1479 | Real R3inv = Rinv * Rinv * Rinv; | |
521 | 1479 | Real R6inv = R3inv * R3inv; | |
522 | 1479 | Real R12inv = R6inv * R6inv; | |
523 | 1479 | Real Rcutinv = 1. / (range - sigma + 1.); | |
524 | 1479 | Real Rcut3inv = Rcutinv * Rcutinv * Rcutinv; | |
525 | 1479 | Real Rcut6inv = Rcut3inv * Rcut3inv; | |
526 | 1479 | Real Rcut12inv = Rcut6inv * Rcut6inv; | |
527 | 1479 | Real R12invminusR6inv = R12inv - R6inv; | |
528 | 1479 | Real Rcut12invminusRcut6inv = Rcut12inv - Rcut6inv; | |
529 | 1479 | Real withcutR12invminusR6inv = R12invminusR6inv - Rcut12invminusRcut6inv; | |
530 | 1479 | Real twoR13invminusR7inv = (R12inv + R12invminusR6inv) * Rinv; | |
531 | 1479 | Real twoRcut13invminusRcut7inv = (Rcut12inv + Rcut12invminusRcut6inv) * Rcutinv; | |
532 | 1479 | Real withcuttwoR13invminusR7inv = twoR13invminusR7inv - twoRcut13invminusRcut7inv; | |
533 | 1479 | Real erdotfracplus = erdotuplus * onepluschiu1dotu2inv; | |
534 | 1479 | Real erdotfracminus = erdotuminus * oneminuschiu1dotu2inv; | |
535 | 1479 | Real erdotfracprimeplus = erdotuplus * onepluschiprimeu1dotu2inv; | |
536 | 1479 | Real erdotfracprimeminus = erdotuminus * oneminuschiprimeu1dotu2inv; | |
537 | 1479 | Real gamma1 = chi * (erdotfracplus + erdotfracminus); | |
538 | 1479 | Real gamma2 = chi * (erdotfracplus - erdotfracminus); | |
539 | 1479 | Real gammaprime1 = chiprime * (erdotfracprimeplus + erdotfracprimeminus); | |
540 | 1479 | Real gammaprime2 = chiprime * (erdotfracprimeplus - erdotfracprimeminus); | |
541 | 1479 | Real pre = 3. * epsilonprime * withcuttwoR13invminusR7inv * sigma * sigma * sigma; | |
542 | 1479 | Real preprime = -2. * withcutR12invminusR6inv; | |
543 | 1479 | Real fourepsilonepsilonprime = 4. * epsilon * epsilonprime; | |
544 | 1479 | Real fourepsilonepsilonprimesq = fourepsilonepsilonprime * epsilonprime; | |
545 | 1479 | Real drU = -6. * fourepsilonepsilonprimesq * twoR13invminusR7inv; | |
546 | 1479 | Real derdotu1U = fourepsilonepsilonprime * (pre * gamma1 + preprime * gammaprime1); | |
547 | 1479 | Real derdotu2U = fourepsilonepsilonprime * (pre * gamma2 + preprime * gammaprime2); | |
548 | 1479 | Vec u1perp = p1.u - p1.u.dot(er) * er; | |
549 | 1479 | Vec u2perp = p2.u - p2.u.dot(er) * er; | |
550 | 1479 | Vec FVec = -drU * er - rinv * (derdotu1U * u1perp + derdotu2U * u2perp); | |
551 | 1479 | Real iotaprime = chiprime * chiprime * | |
552 | 1479 | (erdotfracprimeplus * onepluschiprimeu1dotu2inv - | |
553 | 1479 | erdotfracprimeminus * oneminuschiprimeu1dotu2inv); | |
554 | 1479 | Real du1dotu2U = fourepsilonepsilonprime * | |
555 | 1479 | (chi * chiu1dotu2 * epsilon * epsilon * epsilonprime + iotaprime) * | |
556 | withcutR12invminusR6inv; | ||
557 | 1479 | Vec ercrossu1 = er.cross(p1.u); | |
558 | 1479 | Vec ercrossu2 = er.cross(p2.u); | |
559 | 1479 | Vec u1crossu2 = p1.u.cross(p2.u); | |
560 | 1479 | Vec Tau1Vec = derdotu1U * ercrossu1 - du1dotu2U * u1crossu2; | |
561 | 1479 | Vec Tau2Vec = derdotu2U * ercrossu2 + du1dotu2U * u1crossu2; | |
562 | 1479 | p1.Fint += FVec; | |
563 | 1479 | p1.Tauint += Tau1Vec; | |
564 |
1/2✓ Branch 0 taken 1479 times.
✗ Branch 1 not taken.
|
1479 | if (newton3) { |
565 | 1479 | p2.Fint -= FVec; | |
566 | 1479 | p2.Tauint += Tau2Vec; | |
567 | } | ||
568 | 1479 | 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 | 591459 | double SW::operator()(Particle &p1, Particle &p2) { | |
591 | 591459 | Vec rVec = r_ij(L, p1.r, p2.r); | |
592 | 591459 | Real rsq = rVec.squaredNorm(); | |
593 |
2/2✓ Branch 0 taken 25674 times.
✓ Branch 1 taken 565785 times.
|
591459 | if (rsq >= asq) { |
594 | return 0.; | ||
595 | } | ||
596 | 25674 | Real r = std::sqrt(rsq); | |
597 | 25674 | Real rsqinv = 1. / rsq; | |
598 | 25674 | Real r4inv = rsqinv * rsqinv; | |
599 | 25674 | Real rainv = 1. / (r - a); | |
600 | 25674 | Real exp = std::exp(rainv); | |
601 | 25674 | Vec FVec = -A * (4. * B * r4inv + (B * r4inv - 1.) * rainv * rainv * r) * exp * rsqinv * rVec; | |
602 | 25674 | p1.Fint += FVec; | |
603 |
2/2✓ Branch 0 taken 5663 times.
✓ Branch 1 taken 20011 times.
|
25674 | if (newton3) { |
604 | 5663 | p2.Fint -= FVec; | |
605 | } | ||
606 | 25674 | return A * (B * r4inv - 1.) * exp; | |
607 | } | ||
608 | |||
609 | 106705496 | double SW::operator()(Particle &p1, Particle &p2, Particle &p3) { | |
610 | 106705496 | double E = 0.; | |
611 |
2/2✓ Branch 0 taken 40015830 times.
✓ Branch 1 taken 66689666 times.
|
106705496 | if (newton3) { |
612 | 40015830 | 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 | 106705496 | return E; | |
619 | } | ||
620 | |||
621 | 240084828 | double SW::threebody(Particle &p1, Particle &p2, Particle &p3, int update) { | |
622 | 240084828 | Vec r12Vec = r_ij(L, p1.r, p2.r); | |
623 | 240084828 | Vec r13Vec = r_ij(L, p1.r, p3.r); | |
624 | 240084828 | Real r12sq = r12Vec.squaredNorm(); | |
625 | 240084828 | Real r13sq = r13Vec.squaredNorm(); | |
626 |
4/4✓ Branch 0 taken 9062135 times.
✓ Branch 1 taken 231022693 times.
✓ Branch 2 taken 407853 times.
✓ Branch 3 taken 8654282 times.
|
240084828 | if (r12sq >= asq || r13sq >= asq) { |
627 | return 0.; | ||
628 | } | ||
629 | 407853 | Real r12 = std::sqrt(r12sq); | |
630 | 407853 | Real r13 = std::sqrt(r13sq); | |
631 | 407853 | Real r12ainv = 1. / (r12 - a); | |
632 | 407853 | Real r13ainv = 1. / (r13 - a); | |
633 | 407853 | Real gr12ainv = gamma * r12ainv; | |
634 | 407853 | Real gr13ainv = gamma * r13ainv; | |
635 | 407853 | Real r12r13inv = 1. / (r12 * r13); | |
636 | 407853 | Real costheta = r12Vec.dot(r13Vec) * r12r13inv; | |
637 | 407853 | Real cosdelta = costheta - costheta0; | |
638 | 407853 | Real exp1213 = std::exp(gr12ainv + gr13ainv); | |
639 | 407853 | Real E = lambda * cosdelta * cosdelta * exp1213; | |
640 | 407853 | Real Frad12 = E * gr12ainv * r12ainv / r12; | |
641 | 407853 | Real Frad13 = E * gr13ainv * r13ainv / r13; | |
642 | 407853 | Real Fang = 2. * lambda * cosdelta * exp1213; | |
643 | 407853 | Real Fang1213 = r12r13inv * Fang; | |
644 | 407853 | Real cosFang = costheta * Fang; | |
645 | 407853 | Vec F2Vec = r12Vec * (Frad12 + cosFang / r12sq) - r13Vec * Fang1213; | |
646 | 407853 | 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 68104 times.
|
407853 | 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 | 68104 | default: | |
658 | 68104 | p1.Fint -= F2Vec + F3Vec; | |
659 | 68104 | p2.Fint += F2Vec; | |
660 | 68104 | 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 | 647 | void SWOptimized::reset() { | |
682 | 647 | int N = particles.size(); | |
683 | 647 | h.resize(N); | |
684 | 647 | s.resize(N); | |
685 | 647 | T.resize(N); | |
686 |
2/2✓ Branch 1 taken 1941 times.
✓ Branch 2 taken 647 times.
|
2588 | for (int i = 0; i < N; ++i) { |
687 | 1941 | h[i] = 0; | |
688 | } | ||
689 |
2/2✓ Branch 0 taken 1941 times.
✓ Branch 1 taken 647 times.
|
2588 | for (int i = 0; i < N; ++i) { |
690 | 1941 | s[i].setZero(); | |
691 | } | ||
692 |
2/2✓ Branch 0 taken 1941 times.
✓ Branch 1 taken 647 times.
|
2588 | for (int i = 0; i < N; ++i) { |
693 | 1941 | T[i].setZero(); | |
694 | } | ||
695 | 647 | } | |
696 | |||
697 | 1941 | double SWOptimized::operator()(Particle &p) { | |
698 | 1941 | int i = std::distance(&particles.front(), &p); | |
699 | 1941 | return 0.5 * lambda * | |
700 | 1941 | (costheta0 * costheta0 * h[i] * h[i] - 2 * costheta0 * s[i].squaredNorm() + | |
701 | 1941 | (T[i] * T[i]).trace()); | |
702 | } | ||
703 | |||
704 | 3346 | double SWOptimized::operator()(Particle &p1, Particle &p2) { | |
705 |
2/2✓ Branch 0 taken 1673 times.
✓ Branch 1 taken 1673 times.
|
3346 | if (gather) { |
706 | 1673 | return gatherPhase(p1, p2); | |
707 | } | ||
708 | 1673 | return accumulationPhase(p1, p2); | |
709 | } | ||
710 | |||
711 | 1673 | double SWOptimized::accumulationPhase(Particle &p1, Particle &p2) { | |
712 | 1673 | Vec rVec = r_ij(L, p1.r, p2.r); | |
713 | 1673 | Real rsq = rVec.squaredNorm(); | |
714 |
2/2✓ Branch 0 taken 651 times.
✓ Branch 1 taken 1022 times.
|
1673 | if (rsq >= asq) { |
715 | return 0.; | ||
716 | } | ||
717 | 651 | int i = std::distance(&particles.front(), &p1); | |
718 | 651 | int j = std::distance(&particles.front(), &p2); | |
719 | // Two-body | ||
720 | 651 | Real r = std::sqrt(rsq); | |
721 | 651 | Vec rVecUnit = rVec / r; | |
722 | 651 | Real rsqinv = 1. / rsq; | |
723 | 651 | Real r4inv = rsqinv * rsqinv; | |
724 | 651 | Real rainv = 1. / (r - a); | |
725 | 651 | Real rainvsq = rainv * rainv; | |
726 | 651 | Real exp = std::exp(rainv); | |
727 | 651 | double ETwobody = A * (B * r4inv - 1.) * exp; | |
728 | 651 | Vec FVecTwobody = -A * (4. * B * r4inv + (B * r4inv - 1.) * rainvsq * r) * exp * rsqinv * rVec; | |
729 | // Correction from three-body | ||
730 | 651 | Real gij = std::exp(gamma * rainv); | |
731 | 651 | Real gijsq = gij * gij; | |
732 | 651 | double oneminuscostheta0sq = (1 - costheta0) * (1 - costheta0); | |
733 | 651 | double ECorr = lambda * gijsq * oneminuscostheta0sq; | |
734 | 651 | Vec FVecCorr = 2 * lambda * oneminuscostheta0sq * gijsq * gamma * rainvsq * rVecUnit; | |
735 | // Update accumulators for rest of three-body | ||
736 | 651 | double hi = gij; | |
737 | 651 | Vec si = gij * rVecUnit; | |
738 | 651 | Eigen::Matrix<double, DIM, DIM> Ti = gij * (rVecUnit * rVecUnit.transpose()); | |
739 | 651 | h[i] += hi; | |
740 | 651 | s[i] += si; | |
741 | 651 | T[i] += Ti; | |
742 |
1/2✓ Branch 0 taken 651 times.
✗ Branch 1 not taken.
|
651 | if (newton3) { |
743 | 651 | h[j] += hi; | |
744 | 651 | s[j] -= si; // rVecUnit is flipped -> change of sign | |
745 | 651 | 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 | 651 | Vec FVec = FVecTwobody + FVecCorr; | |
749 | 651 | p1.Fint += FVec; | |
750 |
1/2✓ Branch 0 taken 651 times.
✗ Branch 1 not taken.
|
651 | if (newton3) { |
751 | 651 | p2.Fint -= FVec; | |
752 | } | ||
753 | 651 | return ETwobody - ECorr; | |
754 | } | ||
755 | |||
756 | 1673 | double SWOptimized::gatherPhase(Particle &p1, Particle &p2) { | |
757 | 1673 | Vec rVec = r_ij(L, p1.r, p2.r); | |
758 | 1673 | Real rsq = rVec.squaredNorm(); | |
759 |
2/2✓ Branch 0 taken 651 times.
✓ Branch 1 taken 1022 times.
|
1673 | if (rsq >= asq) { |
760 | return 0.; | ||
761 | } | ||
762 | 651 | int i = std::distance(&particles.front(), &p1); | |
763 | 651 | int j = std::distance(&particles.front(), &p2); | |
764 | 651 | Real r = std::sqrt(rsq); | |
765 | 651 | Real rinv = 1 / r; | |
766 | 651 | Real rsqinv = rinv * rinv; | |
767 | 651 | Real rainv = 1. / (r - a); | |
768 | 651 | Real rainvsq = rainv * rainv; | |
769 | 651 | Real gij = std::exp(gamma * rainv); | |
770 | 651 | Real dgijdrij = -gamma * rainvsq * gij; | |
771 | 651 | Real cijcji = lambda * (costheta0 * costheta0 * dgijdrij * (h[i] + h[j]) * rinv - | |
772 | 651 | 2 * costheta0 * (dgijdrij - gij * rinv) * rVec.dot(s[i] - s[j]) * rsqinv + | |
773 | 1302 | (dgijdrij - 2 * gij * rinv) * | |
774 | 651 | (rVec.transpose() * (T[i] + T[j]) * rVec).value() * rsqinv * rinv); | |
775 | 651 | Vec FVec = cijcji * rVec - 2 * lambda * costheta0 * gij * rinv * (s[i] - s[j]) + | |
776 | 651 | 2 * lambda * gij * rsqinv * ((T[i] + T[j]) * rVec); | |
777 | 651 | p1.Fint += FVec; | |
778 |
1/2✓ Branch 0 taken 651 times.
✗ Branch 1 not taken.
|
651 | if (newton3) { |
779 | 651 | p2.Fint -= FVec; | |
780 | } | ||
781 | return 0.; | ||
782 | } | ||
783 |