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 : #ifndef SYSTEMONE_H
21 : #define SYSTEMONE_H
22 :
23 : #include "State.hpp"
24 : #include "Symmetry.hpp"
25 : #include "SystemBase.hpp"
26 : #include "WignerD.hpp"
27 : #include "utils.hpp"
28 :
29 : #include <cereal/types/array.hpp>
30 : #include <cereal/types/base_class.hpp>
31 : #include <cereal/types/polymorphic.hpp>
32 : #include <cereal/types/unordered_map.hpp>
33 :
34 : #include <array>
35 : #include <cmath>
36 : #include <set>
37 : #include <unordered_map>
38 :
39 : template <typename Scalar_>
40 : class SystemTwo;
41 :
42 : template <typename Scalar_>
43 : class SystemOne : public SystemBase<Scalar_, StateOne> {
44 : public:
45 : using Scalar = Scalar_;
46 : SystemOne(std::string species, MatrixElementCache &cache);
47 : SystemOne(std::string species, MatrixElementCache &cache, bool memory_saving);
48 :
49 : const std::string &getSpecies() const;
50 : void setEfield(std::array<double, 3> field);
51 : void setBfield(std::array<double, 3> field);
52 : void setEfield(std::array<double, 3> field, std::array<double, 3> to_z_axis,
53 : std::array<double, 3> to_y_axis);
54 : void setBfield(std::array<double, 3> field, std::array<double, 3> to_z_axis,
55 : std::array<double, 3> to_y_axis);
56 : void setEfield(std::array<double, 3> field, double alpha, double beta, double gamma);
57 : void setBfield(std::array<double, 3> field, double alpha, double beta, double gamma);
58 : void enableDiamagnetism(bool enable);
59 : void setIonCharge(int c);
60 : void setRydIonOrder(unsigned int o);
61 : void setRydIonDistance(double d);
62 : void setConservedParityUnderReflection(parity_t parity);
63 : void setConservedMomentaUnderRotation(const std::set<float> &momenta);
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_, StateOne> &system) override;
76 :
77 : private:
78 : std::array<double, 3> efield, bfield;
79 : std::unordered_map<int, Scalar_> efield_spherical, bfield_spherical;
80 : bool diamagnetism;
81 : std::unordered_map<std::array<int, 2>, Scalar_, utils::hash<std::array<int, 2>>>
82 : diamagnetism_terms;
83 : int charge;
84 : unsigned int ordermax;
85 : double distance;
86 : std::string species;
87 :
88 : std::unordered_map<int, Eigen::SparseMatrix<Scalar_>> interaction_efield;
89 : std::unordered_map<int, Eigen::SparseMatrix<Scalar_>> interaction_bfield;
90 : std::unordered_map<std::array<int, 2>, Eigen::SparseMatrix<Scalar_>,
91 : utils::hash<std::array<int, 2>>>
92 : interaction_diamagnetism;
93 : std::unordered_map<int, Eigen::SparseMatrix<Scalar_>> interaction_multipole;
94 : parity_t sym_reflection;
95 : std::set<float> sym_rotation;
96 :
97 : ////////////////////////////////////////////////////////////////////
98 : /// Utility methods ////////////////////////////////////////////////
99 : ////////////////////////////////////////////////////////////////////
100 :
101 : void addSymmetrizedBasisvectors(const StateOne &state, size_t &idx, const double &energy,
102 : std::vector<Eigen::Triplet<Scalar_>> &basisvectors_triplets,
103 : std::vector<Eigen::Triplet<Scalar_>> &hamiltonian_triplets,
104 : parity_t &sym_reflection_local);
105 :
106 : void addBasisvectors(const StateOne &state, const size_t &idx, const Scalar_ &value,
107 : std::vector<Eigen::Triplet<Scalar_>> &basisvectors_triplets);
108 :
109 : void changeToSphericalbasis(std::array<double, 3> field,
110 : std::unordered_map<int, double> &field_spherical);
111 : void changeToSphericalbasis(std::array<double, 3> field,
112 : std::unordered_map<int, std::complex<double>> &field_spherical);
113 : void addTriplet(std::vector<Eigen::Triplet<Scalar_>> &triplets, size_t r_idx, size_t c_idx,
114 : Scalar_ val);
115 : void rotateVector(std::array<double, 3> &field, std::array<double, 3> &to_z_axis,
116 : std::array<double, 3> &to_y_axis);
117 : void rotateVector(std::array<double, 3> &field, double alpha, double beta, double gamma);
118 :
119 : template <class T>
120 275 : void addRotated(const StateOne &state, const size_t &idx,
121 : std::vector<Eigen::Triplet<T>> &triplets, WignerD &wigner, const double &alpha,
122 : const double &beta, const double &gamma) {
123 : // Check whether the angles are compatible to the used data type
124 275 : double tolerance = 1e-16;
125 0 : if (std::is_same<T, double>::value &&
126 0 : std::abs(std::remainder(alpha, 2 * M_PI)) > tolerance) {
127 0 : throw std::runtime_error(
128 : "If the Euler angle alpha is not a multiple of 2 pi, the Wigner D matrix element "
129 : "is complex and cannot be converted to double.");
130 : }
131 0 : if (std::is_same<T, double>::value &&
132 0 : std::abs(std::remainder(gamma, 2 * M_PI)) > tolerance) {
133 0 : throw std::runtime_error(
134 : "If the Euler angle gamma is not a multiple of 2 pi, the Wigner D matrix element "
135 : "is complex and cannot be converted to double.");
136 : }
137 :
138 : // Add rotated triplet entries
139 1847 : for (float m = -state.getJ(); m <= state.getJ(); ++m) {
140 3144 : StateOne newstate(state.getSpecies(), state.getN(), state.getL(), state.getJ(), m);
141 1572 : auto state_iter = this->states.template get<1>().find(newstate);
142 :
143 1572 : if (state_iter != this->states.template get<1>().end()) {
144 : // We calculate the matrix element <m|d_mM|state.getM()>|m>. Note that we must
145 : // invert the angles since they are given for a passive rotation, i.e. rotation of
146 : // the coordinate system, and the WignerD matrix is for an active rotation of spins.
147 0 : auto val =
148 1572 : utils::convert<T>(wigner(state.getJ(), m, state.getM(), -gamma, -beta, -alpha));
149 1572 : triplets.push_back(Eigen::Triplet<T>(state_iter->idx, idx, val));
150 : } else {
151 0 : std::cerr << "Warning: Incomplete rotation because the basis is lacking some "
152 : "Zeeman levels."
153 0 : << std::endl;
154 : }
155 : }
156 275 : }
157 :
158 : bool isRefelectionAndRotationCompatible();
159 :
160 : ////////////////////////////////////////////////////////////////////
161 : /// Method for serialization ///////////////////////////////////////
162 : ////////////////////////////////////////////////////////////////////
163 :
164 : friend class cereal::access;
165 : friend class SystemTwo<Scalar_>;
166 : SystemOne();
167 :
168 : template <class Archive>
169 15 : void serialize(Archive &ar, unsigned int /* version */) {
170 15 : ar &cereal::make_nvp("base_class", cereal::base_class<SystemBase<Scalar_, StateOne>>(this));
171 15 : ar &CEREAL_NVP(species);
172 15 : ar &CEREAL_NVP(efield) & CEREAL_NVP(bfield) & CEREAL_NVP(diamagnetism) &
173 30 : CEREAL_NVP(sym_reflection) & CEREAL_NVP(sym_rotation);
174 15 : ar &CEREAL_NVP(efield_spherical) & CEREAL_NVP(bfield_spherical) &
175 15 : CEREAL_NVP(diamagnetism_terms);
176 15 : ar &CEREAL_NVP(interaction_efield) & CEREAL_NVP(interaction_bfield) &
177 15 : CEREAL_NVP(interaction_diamagnetism);
178 15 : }
179 : };
180 :
181 : extern template class SystemOne<std::complex<double>>;
182 : extern template class SystemOne<double>;
183 :
184 : #ifndef SWIG
185 64 : CEREAL_REGISTER_TYPE(SystemOne<std::complex<double>>) // NOLINT
186 64 : CEREAL_REGISTER_TYPE(SystemOne<double>) // NOLINT
187 : #endif
188 :
189 : #endif
|
