Line data Source code
1 : /*
2 : * Copyright (c) 2016 Sebastian Weber, Henri Menke, Johannes Block. 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 : #ifndef SYSTEMTWO_H
21 : #define SYSTEMTWO_H
22 :
23 : #include "State.hpp"
24 : #include "SystemBase.hpp"
25 : #include "SystemOne.hpp"
26 :
27 : #include "EigenCompat.hpp"
28 : #include <Eigen/SparseCore>
29 : #include <boost/math/special_functions/binomial.hpp>
30 : #include <cereal/types/array.hpp>
31 : #include <cereal/types/base_class.hpp>
32 : #include <cereal/types/polymorphic.hpp>
33 : #include <cereal/types/string.hpp>
34 : #include <cereal/types/unordered_map.hpp>
35 :
36 : #include <cmath>
37 : #include <set>
38 : #include <type_traits>
39 :
40 : template <typename Scalar_>
41 : class SystemTwo : public SystemBase<Scalar_, StateTwo> {
42 : public:
43 : using Scalar = Scalar_;
44 : SystemTwo(const SystemOne<Scalar_> &b1, const SystemOne<Scalar_> &b2,
45 : MatrixElementCache &cache);
46 : SystemTwo(const SystemOne<Scalar_> &b1, const SystemOne<Scalar_> &b2, MatrixElementCache &cache,
47 : bool memory_saving);
48 :
49 : const std::array<std::string, 2> &getSpecies();
50 : std::vector<StateOne> getStatesFirst();
51 : std::vector<StateOne> getStatesSecond();
52 : void enableGreenTensor(bool GTboolean);
53 : void setSurfaceDistance(double d);
54 : void setAngle(double a);
55 : void setDistance(double d);
56 : void setDistanceVector(std::array<double, 3> d);
57 : void setOrder(double o);
58 :
59 : void setConservedParityUnderPermutation(parity_t parity);
60 : void setConservedParityUnderInversion(parity_t parity);
61 : void setConservedParityUnderReflection(parity_t parity);
62 : void setConservedMomentaUnderRotation(const std::set<int> &momenta);
63 : void setOneAtomBasisvectors(const std::vector<std::array<size_t, 2>> &indices);
64 :
65 : protected:
66 : void initializeBasis() override;
67 : void initializeInteraction() override;
68 : void addInteraction() override;
69 : void transformInteraction(const Eigen::SparseMatrix<Scalar_> &transformator) override;
70 : void deleteInteraction() override;
71 : Eigen::SparseMatrix<Scalar_> rotateStates(const std::vector<size_t> &states_indices,
72 : double alpha, double beta, double gamma) override;
73 : Eigen::SparseMatrix<Scalar_> buildStaterotator(double alpha, double beta,
74 : double gamma) override;
75 : void incorporate(SystemBase<Scalar_, StateTwo> &system) override;
76 : void onStatesChange() override;
77 :
78 : private:
79 : std::array<std::string, 2> species;
80 : SystemOne<Scalar_> system1; // is needed in the initializeBasis method and afterwards deleted
81 : SystemOne<Scalar_> system2; // is needed in the initializeBasis method and afterwards deleted
82 :
83 : std::unordered_map<int, Eigen::SparseMatrix<Scalar_>> interaction_angulardipole;
84 : std::unordered_map<int, Eigen::SparseMatrix<Scalar_>> interaction_multipole;
85 : std::unordered_map<int, Eigen::SparseMatrix<Scalar_>> interaction_greentensor_dd;
86 : std::unordered_map<int, Eigen::SparseMatrix<Scalar_>> interaction_greentensor_dq;
87 : std::unordered_map<int, Eigen::SparseMatrix<Scalar_>> interaction_greentensor_qd;
88 :
89 : double minimal_le_roy_radius;
90 : double distance;
91 : double distance_x; // NOLINT
92 : double distance_y; // NOLINT
93 : double distance_z;
94 : bool GTbool; // NOLINT
95 : double surface_distance;
96 : unsigned int ordermax; // NOLINT
97 :
98 : parity_t sym_permutation; // NOLINT
99 : parity_t sym_inversion; // NOLINT
100 : parity_t sym_reflection; // NOLINT
101 : std::set<int> sym_rotation;
102 :
103 : std::unordered_map<int, double> angle_terms;
104 : std::unordered_map<int, double> greentensor_terms_dd;
105 : std::unordered_map<int, double> greentensor_terms_dq;
106 : std::unordered_map<int, double> greentensor_terms_qd;
107 : std::vector<std::array<size_t, 2>> one_atom_basisvectors_indices;
108 :
109 : ////////////////////////////////////////////////////////////////////
110 : /// Utility methods ////////////////////////////////////////////////
111 : ////////////////////////////////////////////////////////////////////
112 :
113 : void checkDistance(const double &distance);
114 :
115 : void addBasisvectors(const StateTwo &state, const size_t &col_new, const Scalar_ &value_new,
116 : std::vector<Eigen::Triplet<Scalar_>> &basisvectors_triplets,
117 : std::vector<double> &sqnorm_list);
118 :
119 : template <typename T>
120 326677 : void addTriplet(std::vector<Eigen::Triplet<T>> &triplets, size_t r_idx, size_t c_idx, T val) {
121 326677 : triplets.emplace_back(r_idx, c_idx, val);
122 326677 : }
123 :
124 : template <class T>
125 307 : void addRotated(const StateTwo &state, const size_t &idx,
126 : std::vector<Eigen::Triplet<T>> &triplets, WignerD &wigner, const double &alpha,
127 : const double &beta, const double &gamma) {
128 : // Check whether the angles are compatible to the used data type
129 307 : double tolerance = 1e-16;
130 38 : if (std::is_same<T, double>::value &&
131 38 : std::abs(std::remainder(alpha, 2 * M_PI)) > tolerance) {
132 0 : throw std::runtime_error(
133 : "If the Euler angle alpha is not a multiple of 2 pi, the Wigner D matrix element "
134 : "is complex and cannot be converted to double.");
135 : }
136 38 : if (std::is_same<T, double>::value &&
137 38 : std::abs(std::remainder(gamma, 2 * M_PI)) > tolerance) {
138 0 : throw std::runtime_error(
139 : "If the Euler angle gamma is not a multiple of 2 pi, the Wigner D matrix element "
140 : "is complex and cannot be converted to double.");
141 : }
142 :
143 : // Add rotated triplet entries
144 614 : std::vector<T> val2_vector;
145 307 : val2_vector.reserve(2 * state.getSecondState().getJ() + 1);
146 :
147 1787 : for (float m2 = -state.getSecondState().getJ(); m2 <= state.getSecondState().getJ(); ++m2) {
148 1480 : val2_vector.push_back(
149 2884 : utils::convert<T>(wigner(state.getSecondState().getJ(), m2,
150 1480 : state.getSecondState().getM(), -gamma, -beta, -alpha)));
151 : }
152 :
153 1787 : for (float m1 = -state.getFirstState().getJ(); m1 <= state.getFirstState().getJ(); ++m1) {
154 : auto val1 =
155 1480 : utils::convert<T>(wigner(state.getFirstState().getJ(), m1,
156 1480 : state.getFirstState().getM(), -gamma, -beta, -alpha));
157 :
158 8352 : for (float m2 = -state.getSecondState().getJ(); m2 <= state.getSecondState().getJ();
159 : ++m2) {
160 13744 : StateTwo newstate(state.getSpecies(), state.getN(), state.getL(), state.getJ(),
161 : {{m1, m2}});
162 6872 : auto state_iter = this->states.template get<1>().find(newstate);
163 :
164 6872 : if (state_iter != this->states.template get<1>().end()) {
165 6872 : auto val = val1 * val2_vector[m2 + state.getSecondState().getJ()];
166 6872 : triplets.push_back(Eigen::Triplet<T>(state_iter->idx, idx, val));
167 : } else {
168 0 : std::cerr << "Warning: Incomplete rotation because the basis is lacking some "
169 : "Zeeman levels."
170 0 : << std::endl;
171 : }
172 : }
173 : }
174 307 : }
175 :
176 : bool isRefelectionAndRotationCompatible();
177 :
178 : ////////////////////////////////////////////////////////////////////
179 : /// Method for serialization ///////////////////////////////////////
180 : ////////////////////////////////////////////////////////////////////
181 :
182 : friend class cereal::access;
183 : SystemTwo();
184 :
185 : template <class Archive>
186 5 : void serialize(Archive &ar, unsigned int /* version */) {
187 5 : ar &cereal::make_nvp("base_class", cereal::base_class<SystemBase<Scalar_, StateTwo>>(this));
188 5 : ar &CEREAL_NVP(species) & CEREAL_NVP(system1) & CEREAL_NVP(system2);
189 5 : ar &CEREAL_NVP(distance) & CEREAL_NVP(distance_x) & CEREAL_NVP(distance_y) &
190 10 : CEREAL_NVP(distance_z) & CEREAL_NVP(surface_distance) & CEREAL_NVP(ordermax);
191 5 : ar &CEREAL_NVP(sym_permutation) & CEREAL_NVP(sym_inversion) & CEREAL_NVP(sym_reflection) &
192 5 : CEREAL_NVP(sym_rotation);
193 5 : ar &CEREAL_NVP(angle_terms) & CEREAL_NVP(greentensor_terms_dd) &
194 10 : CEREAL_NVP(greentensor_terms_dq) & CEREAL_NVP(greentensor_terms_qd);
195 5 : ar &CEREAL_NVP(interaction_angulardipole) & CEREAL_NVP(interaction_multipole) &
196 10 : CEREAL_NVP(interaction_greentensor_dd) & CEREAL_NVP(interaction_greentensor_dq) &
197 5 : CEREAL_NVP(interaction_greentensor_qd);
198 5 : }
199 : };
200 :
201 : extern template class SystemTwo<std::complex<double>>;
202 : extern template class SystemTwo<double>;
203 :
204 : #ifndef SWIG
205 64 : CEREAL_REGISTER_TYPE(SystemTwo<std::complex<double>>) // NOLINT
206 64 : CEREAL_REGISTER_TYPE(SystemTwo<double>) // NOLINT
207 : #endif
208 :
209 : #endif
|
