Line data Source code
1 : /*
2 : * Copyright (c) 2016 Sebastian Weber, Henri Menke. All rights reserved.
3 : *
4 : * This file is part of the pairinteraction library.
5 : *
6 : * The pairinteraction library is free software: you can redistribute it and/or modify
7 : * it under the terms of the GNU Lesser General Public License as published by
8 : * the Free Software Foundation, either version 3 of the License, or
9 : * (at your option) any later version.
10 : *
11 : * The pairinteraction library is distributed in the hope that it will be useful,
12 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : * GNU Lesser General Public License for more details.
15 : *
16 : * You should have received a copy of the GNU Lesser General Public License
17 : * along with the pairinteraction library. If not, see <http://www.gnu.org/licenses/>.
18 : */
19 :
20 : #include "PerturbativeInteraction.hpp"
21 : #include "Constants.hpp"
22 :
23 : #include "EigenCompat.hpp"
24 : #include <Eigen/Core>
25 :
26 : #include "Constants.hpp"
27 :
28 : #include <unordered_set>
29 :
30 2 : PerturbativeInteraction::PerturbativeInteraction(MatrixElementCache &cache) : cache(cache) {
31 2 : initializeAngleTerms(0);
32 2 : }
33 :
34 3 : PerturbativeInteraction::PerturbativeInteraction(double angle, MatrixElementCache &cache)
35 3 : : cache(cache) {
36 3 : initializeAngleTerms(angle);
37 3 : }
38 :
39 2 : double PerturbativeInteraction::getC6(const StateTwo &state, double deltaN) {
40 2 : double C6 = 0;
41 :
42 18 : for (int n0 = state.getN(0) - deltaN; n0 <= state.getN(0) + deltaN; ++n0) {
43 :
44 146 : for (int n1 = state.getN(1) - deltaN; n1 <= state.getN(1) + deltaN; ++n1) {
45 :
46 260 : std::vector<int> array_l0;
47 130 : if (state.getL(0) > 0) {
48 81 : array_l0.push_back(state.getL(0) - 1);
49 : }
50 130 : if (state.getL(0) < n0 - 1) {
51 130 : array_l0.push_back(state.getL(0) + 1);
52 : }
53 341 : for (int l0 : array_l0) {
54 :
55 422 : std::vector<int> array_l1;
56 211 : if (state.getL(1) > 0) {
57 162 : array_l1.push_back(state.getL(1) - 1);
58 : }
59 211 : if (state.getL(1) < n1 - 1) {
60 211 : array_l1.push_back(state.getL(1) + 1);
61 : }
62 584 : for (int l1 : array_l1) {
63 :
64 746 : std::set<float> array_j0;
65 373 : if (std::abs(std::abs(l0 - state.getS(0)) - state.getJ(0)) < 2) {
66 373 : array_j0.insert(std::abs(l0 - state.getS(0)));
67 : }
68 373 : if (std::abs(l0 + state.getS(0) - state.getJ(0)) < 2) {
69 211 : array_j0.insert(l0 + state.getS(0));
70 : }
71 373 : if (array_j0.size() == 2) {
72 211 : for (float j0 = std::abs(l0 - state.getS(0)) + 1; j0 < l0 + state.getS(0);
73 : ++j0) {
74 0 : array_j0.insert(j0);
75 : }
76 : }
77 957 : for (float j0 : array_j0) {
78 :
79 1168 : std::set<float> array_j1;
80 584 : if (std::abs(std::abs(l1 - state.getS(1)) - state.getJ(1)) < 2) {
81 584 : array_j1.insert(std::abs(l1 - state.getS(1)));
82 : }
83 584 : if (std::abs(l1 + state.getS(1) - state.getJ(1)) < 2) {
84 341 : array_j1.insert(l1 + state.getS(1));
85 : }
86 584 : if (array_j1.size() == 2) {
87 341 : for (float j1 = std::abs(l1 - state.getS(1)) + 1;
88 341 : j1 < l1 + state.getS(1); ++j1) {
89 0 : array_j1.insert(j1);
90 : }
91 : }
92 1509 : for (float j1 : array_j1) {
93 :
94 3700 : for (auto &q : array_q) {
95 : // q = final.m - initial.m
96 :
97 2775 : float m0 = state.getM(0) + q[0];
98 2775 : if (std::abs(m0) > j0) {
99 1573 : continue;
100 : }
101 :
102 1948 : float m1 = state.getM(1) + q[1];
103 1948 : if (std::abs(m1) > j1) {
104 746 : continue;
105 : }
106 :
107 : StateTwo state_virtual =
108 0 : StateTwo(state.getSpecies(), {{n0, n1}}, {{l0, l1}}, {{j0, j1}},
109 1202 : {{m0, m1}});
110 :
111 : double energydiff =
112 1202 : state.getEnergy(cache) - state_virtual.getEnergy(cache);
113 :
114 1202 : C6 +=
115 1202 : std::pow(
116 1202 : coulombs_constant *
117 1202 : array_angle_term[3 * (q[0] + 1) + (q[1] + 1)] *
118 1202 : cache.getElectricDipole(state_virtual.getFirstState(),
119 : state.getFirstState()) *
120 1202 : cache.getElectricDipole(state_virtual.getSecondState(),
121 : state.getSecondState()),
122 1202 : 2) /
123 : energydiff; // getMultipole(final, inital, 1)
124 : }
125 : }
126 : }
127 : }
128 : }
129 : }
130 : }
131 :
132 2 : return C6;
133 : }
134 :
135 3 : Eigen::MatrixX<double> PerturbativeInteraction::getC6(const std::vector<StateTwo> &states,
136 : double deltaN) {
137 6 : Eigen::MatrixX<double> C6_matrix = Eigen::MatrixX<double>::Zero(states.size(), states.size());
138 :
139 3 : std::unordered_set<StateTwo> set_states(states.begin(), states.end());
140 :
141 15 : for (size_t idx_row = 0; idx_row < states.size(); ++idx_row) {
142 12 : auto &state_row = states[idx_row];
143 :
144 42 : for (size_t idx_col = idx_row; idx_col < states.size(); ++idx_col) {
145 30 : auto &state_col = states[idx_col];
146 :
147 30 : if (state_row.getS(0) != state_col.getS(0) || state_row.getS(1) != state_col.getS(1)) {
148 0 : continue;
149 : }
150 30 : auto s = state_col.getS();
151 :
152 60 : if (state_row.getSpecies(0) != state_col.getSpecies(0) ||
153 30 : state_row.getSpecies(1) != state_col.getSpecies(1)) {
154 0 : continue;
155 : }
156 30 : auto species = state_col.getSpecies();
157 :
158 30 : int n0_min = std::min(state_row.getN(0), state_col.getN(0));
159 30 : int n0_max = std::max(state_row.getN(0), state_col.getN(0));
160 30 : int n1_min = std::min(state_row.getN(1), state_col.getN(1));
161 30 : int n1_max = std::max(state_row.getN(1), state_col.getN(1));
162 :
163 30 : double C6 = 0;
164 :
165 200 : for (int n0 = n0_min - deltaN; n0 <= n0_max + deltaN; ++n0) {
166 :
167 1160 : for (int n1 = n1_min - deltaN; n1 <= n1_max + deltaN; ++n1) {
168 :
169 1980 : std::vector<int> array_l0;
170 990 : if (state_row.getL(0) == state_col.getL(0)) {
171 890 : if (state_col.getL(0) < n0 - 1) {
172 890 : array_l0.push_back(state_col.getL(0) + 1);
173 : }
174 890 : if (state_col.getL(0) > 0) {
175 75 : array_l0.push_back(state_col.getL(0) - 1);
176 : }
177 100 : } else if (state_row.getL(0) + 2 == state_col.getL(0) &&
178 0 : state_col.getL(0) < n0 - 1) {
179 0 : array_l0.push_back(state_col.getL(0) + 1);
180 100 : } else if (state_row.getL(0) - 2 == state_col.getL(0) &&
181 0 : state_col.getL(0) > 0) {
182 0 : array_l0.push_back(state_col.getL(0) - 1);
183 : }
184 1955 : for (int l0 : array_l0) {
185 :
186 1930 : std::vector<int> array_l1;
187 965 : if (state_row.getL(1) == state_col.getL(1)) {
188 965 : if (state_col.getL(1) < n1 - 1) {
189 965 : array_l1.push_back(state_col.getL(1) + 1);
190 : }
191 965 : if (state_col.getL(1) > 0) {
192 75 : array_l1.push_back(state_col.getL(1) - 1);
193 : }
194 0 : } else if (state_row.getL(1) + 2 == state_col.getL(1) &&
195 0 : state_col.getL(1) < n1 - 1) {
196 0 : array_l1.push_back(state_col.getL(1) + 1);
197 0 : } else if (state_row.getL(1) - 2 == state_col.getL(1) &&
198 0 : state_col.getL(1) > 0) {
199 0 : array_l1.push_back(state_col.getL(1) - 1);
200 : }
201 2005 : for (int l1 : array_l1) {
202 :
203 2080 : std::set<float> array_j0;
204 2080 : if (std::abs(std::abs(l0 - s[0]) - state_col.getJ(0)) < 2 &&
205 1040 : std::abs(std::abs(l0 - s[0]) - state_row.getJ(0)) < 2) {
206 1040 : array_j0.insert(std::abs(l0 - s[0]));
207 : }
208 2005 : if (std::abs(l0 + s[0] - state_col.getJ(0)) < 2 &&
209 965 : std::abs(l0 + s[0] - state_row.getJ(0)) < 2) {
210 965 : array_j0.insert(l0 + s[0]);
211 : }
212 2970 : for (float j0 : array_j0) {
213 :
214 3860 : std::set<float> array_j1;
215 3860 : if (std::abs(std::abs(l1 - s[1]) - state_col.getJ(1)) < 2 &&
216 1930 : std::abs(std::abs(l1 - s[1]) - state_row.getJ(1)) < 2) {
217 1930 : array_j1.insert(std::abs(l1 - s[1]));
218 : }
219 3710 : if (std::abs(l1 + s[1] - state_col.getJ(1)) < 2 &&
220 1780 : std::abs(l1 + s[1] - state_row.getJ(1)) < 2) {
221 1780 : array_j1.insert(l1 + s[1]);
222 : }
223 5490 : for (float j1 : array_j1) {
224 :
225 32000 : for (auto &idx : array_q) {
226 28440 : int q0_forth = idx[0]; // q = final.m - initial.m
227 28440 : int q1_forth = idx[1];
228 :
229 28440 : float m0 = state_col.getM(0) + q0_forth;
230 28440 : if (std::abs(m0) > j0) {
231 11464 : continue;
232 : }
233 :
234 23700 : float m1 = state_col.getM(1) + q1_forth;
235 23700 : if (std::abs(m1) > j1) {
236 3850 : continue;
237 : }
238 :
239 19850 : int q0_back = state_row.getM(0) - m0;
240 19850 : int q1_back = state_row.getM(1) - m1;
241 19850 : if (std::abs(q0_back) > 1 || std::abs(q1_back) > 1) {
242 2862 : continue;
243 : }
244 16988 : if (array_q.size() == 3 && q0_back + q1_back != 0) {
245 0 : continue;
246 : }
247 :
248 : StateTwo state_virtual =
249 : StateTwo(species, {{n0, n1}}, {{l0, l1}}, {{j0, j1}},
250 16988 : {{m0, m1}});
251 16988 : if (set_states.find(state_virtual) != set_states.end()) {
252 12 : continue;
253 : }
254 :
255 16976 : double energydiff_row = state_row.getEnergy(cache) -
256 16976 : state_virtual.getEnergy(cache);
257 :
258 16976 : double energydiff_col = state_col.getEnergy(cache) -
259 16976 : state_virtual.getEnergy(cache);
260 :
261 33952 : C6 += coulombs_constant *
262 16976 : array_angle_term[3 * (q0_back + 1) + (q1_back + 1)] *
263 16976 : cache.getElectricDipole(state_row.getFirstState(),
264 16976 : state_virtual.getFirstState()) *
265 16976 : cache.getElectricDipole(
266 : state_row.getSecondState(),
267 16976 : state_virtual.getSecondState()) *
268 16976 : coulombs_constant *
269 16976 : array_angle_term[3 * (q0_forth + 1) + (q1_forth + 1)] *
270 16976 : cache.getElectricDipole(state_virtual.getFirstState(),
271 16976 : state_col.getFirstState()) *
272 16976 : cache.getElectricDipole(state_virtual.getSecondState(),
273 16976 : state_col.getSecondState()) *
274 16976 : 0.5 *
275 16976 : (1 / energydiff_row +
276 16976 : 1 / energydiff_col); // getMultipole(final, inital, 1)
277 : }
278 : }
279 : }
280 : }
281 : }
282 : }
283 : }
284 :
285 30 : C6_matrix(idx_row, idx_col) = C6;
286 : }
287 : }
288 :
289 9 : return C6_matrix.selfadjointView<Eigen::Upper>();
290 : }
291 :
292 2 : Eigen::MatrixX<double> PerturbativeInteraction::getC3(const std::vector<StateTwo> &states) {
293 2 : Eigen::MatrixX<double> C3_matrix = Eigen::MatrixX<double>::Zero(states.size(), states.size());
294 :
295 10 : for (size_t idx_row = 0; idx_row < states.size(); ++idx_row) {
296 8 : auto &state_row = states[idx_row];
297 :
298 20 : for (size_t idx_col = idx_row + 1; idx_col < states.size(); ++idx_col) {
299 12 : auto &state_col = states[idx_col];
300 :
301 12 : int q0 = state_row.getM(0) - state_col.getM(0);
302 12 : int q1 = state_row.getM(1) - state_col.getM(1);
303 :
304 12 : if (array_q.size() == 3 && q0 + q1 != 0) {
305 0 : continue;
306 : }
307 12 : if (std::abs(q0) > 1 || std::abs(q1) > 1) {
308 0 : continue;
309 : }
310 :
311 12 : C3_matrix(idx_row, idx_col) = coulombs_constant *
312 12 : array_angle_term[3 * (q0 + 1) + (q1 + 1)] *
313 12 : cache.getElectricDipole(state_row.getFirstState(), state_col.getFirstState()) *
314 12 : cache.getElectricDipole(state_row.getSecondState(), state_col.getSecondState());
315 : }
316 : }
317 :
318 6 : return C3_matrix.selfadjointView<Eigen::Upper>();
319 : }
320 :
321 2 : Eigen::MatrixX<double> PerturbativeInteraction::getEnergy(const std::vector<StateTwo> &states) {
322 : Eigen::MatrixX<double> energies_matrix =
323 2 : Eigen::MatrixX<double>::Zero(states.size(), states.size());
324 :
325 10 : for (size_t idx = 0; idx < states.size(); ++idx) {
326 8 : auto &state = states[idx];
327 :
328 8 : double energy = state.getEnergy(cache);
329 8 : energies_matrix(idx, idx) = energy;
330 : }
331 :
332 2 : return energies_matrix;
333 : }
334 :
335 5 : void PerturbativeInteraction::initializeAngleTerms(double angle) {
336 5 : if (angle == 0) {
337 3 : array_q.reserve(3);
338 :
339 3 : array_q.push_back({{0, 0}});
340 3 : array_angle_term[3 * (0 + 1) + (0 + 1)] = -2;
341 :
342 3 : array_q.push_back({{1, -1}});
343 3 : array_angle_term[3 * (1 + 1) + (-1 + 1)] = -1;
344 :
345 3 : array_q.push_back({{-1, 1}});
346 3 : array_angle_term[3 * (-1 + 1) + (1 + 1)] = -1;
347 : } else {
348 2 : array_q.reserve(9);
349 :
350 2 : array_q.push_back({{0, 0}});
351 2 : array_angle_term[3 * (0 + 1) + (0 + 1)] = 1. - 3. * std::pow(std::cos(angle), 2);
352 :
353 2 : array_q.push_back({{1, -1}});
354 2 : array_angle_term[3 * (1 + 1) + (-1 + 1)] = -1 + 1.5 * std::pow(std::sin(angle), 2);
355 :
356 2 : array_q.push_back({{-1, 1}});
357 2 : array_angle_term[3 * (-1 + 1) + (1 + 1)] = -1 + 1.5 * std::pow(std::sin(angle), 2);
358 :
359 2 : array_q.push_back({{1, 1}});
360 2 : array_angle_term[3 * (1 + 1) + (1 + 1)] = -1.5 * std::pow(std::sin(angle), 2);
361 :
362 2 : array_q.push_back({{-1, -1}});
363 2 : array_angle_term[3 * (-1 + 1) + (-1 + 1)] = -1.5 * std::pow(std::sin(angle), 2);
364 :
365 2 : array_q.push_back({{0, 1}});
366 4 : array_angle_term[3 * (0 + 1) + (1 + 1)] =
367 2 : 3. / std::sqrt(2) * std::sin(angle) * std::cos(angle);
368 :
369 2 : array_q.push_back({{1, 0}});
370 4 : array_angle_term[3 * (1 + 1) + (0 + 1)] =
371 2 : 3. / std::sqrt(2) * std::sin(angle) * std::cos(angle);
372 :
373 2 : array_q.push_back({{0, -1}});
374 4 : array_angle_term[3 * (0 + 1) + (-1 + 1)] =
375 2 : -3. / std::sqrt(2) * std::sin(angle) * std::cos(angle);
376 :
377 2 : array_q.push_back({{-1, 0}});
378 2 : array_angle_term[3 * (-1 + 1) + (0 + 1)] =
379 2 : -3. / std::sqrt(2) * std::sin(angle) * std::cos(angle);
380 : }
381 5 : }
|
