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 SYSTEMBASE_H
21 : #define SYSTEMBASE_H
22 :
23 : #include "State.hpp"
24 : #include "serialization_eigen.hpp"
25 : #include "utils.hpp"
26 :
27 : #include <boost/multi_index/hashed_index.hpp>
28 : #include <boost/multi_index/member.hpp>
29 : #include <boost/multi_index/random_access_index.hpp>
30 : #include <boost/multi_index_container.hpp>
31 : #include <cereal/types/complex.hpp>
32 : #include <cereal/types/set.hpp>
33 : #include <cereal/types/vector.hpp>
34 :
35 : #include <algorithm>
36 : #include <complex>
37 : #include <functional>
38 : #include <set>
39 : #include <vector>
40 :
41 : class MatrixElementCache;
42 :
43 : template <class T>
44 : class enumerated_state {
45 : public:
46 27699 : enumerated_state(size_t idx, T state) : idx(idx), state(std::move(state)) {}
47 1102 : enumerated_state()
48 1102 : : state() { // TODO remove and use
49 : // http://www.boost.org/doc/libs/1_46_1/libs/serialization/doc/serialization.html#constructors
50 : // instead
51 1102 : }
52 : size_t idx{0};
53 : T state;
54 :
55 : private:
56 : ////////////////////////////////////////////////////////////////////
57 : /// Method for serialization ///////////////////////////////////////
58 : ////////////////////////////////////////////////////////////////////
59 :
60 : friend class cereal::access;
61 :
62 : template <class Archive>
63 1102 : void serialize(Archive &ar, unsigned int /* version */) {
64 1102 : ar &CEREAL_NVP(idx) & CEREAL_NVP(state);
65 1102 : }
66 : };
67 :
68 : #ifndef SWIG
69 :
70 : namespace utils {
71 : template <class T>
72 : struct hash<enumerated_state<T>> {
73 48 : size_t operator()(enumerated_state<T> const &s) const { return std::hash<T>{}(s.state); }
74 : };
75 : } // namespace utils
76 :
77 : #endif
78 :
79 : template <class T>
80 : struct states_set {
81 : typedef typename boost::multi_index_container<
82 : enumerated_state<T>,
83 : boost::multi_index::indexed_by<
84 : boost::multi_index::random_access<>,
85 : boost::multi_index::hashed_unique<
86 : boost::multi_index::member<enumerated_state<T>, T, &enumerated_state<T>::state>,
87 : std::hash<T>>>>
88 : type;
89 : };
90 :
91 : namespace cereal {
92 :
93 : template <class Archive, typename Value, typename IndexSpecifierList, typename Allocator>
94 0 : void save(
95 : Archive &ar,
96 : boost::multi_index::multi_index_container<Value, IndexSpecifierList, Allocator> const &vector) {
97 0 : ar(make_size_tag(static_cast<size_type>(vector.size())));
98 0 : for (auto &&v : vector) {
99 0 : ar(v);
100 : }
101 0 : }
102 :
103 : template <class Archive, typename Value, typename IndexSpecifierList, typename Allocator>
104 20 : void load(Archive &ar,
105 : boost::multi_index::multi_index_container<Value, IndexSpecifierList, Allocator> &vector) {
106 : size_type size;
107 20 : ar(make_size_tag(size));
108 20 : vector.reserve(static_cast<size_type>(size));
109 1122 : for (size_type i = 0; i < size; ++i) {
110 2204 : Value v;
111 1102 : ar(v);
112 1102 : vector.push_back(v);
113 : }
114 20 : }
115 :
116 : } // namespace cereal
117 :
118 : template <class Scalar, class State>
119 : class SystemBase {
120 : public:
121 385 : virtual ~SystemBase() = default;
122 :
123 : void setMinimalNorm(const double &threshold);
124 :
125 : ////////////////////////////////////////////////////////////////////
126 : /// Methods to restrict the number of states inside the basis //////
127 : ////////////////////////////////////////////////////////////////////
128 :
129 : void restrictEnergy(double e_min, double e_max);
130 :
131 : void restrictN(int n_min, int n_max);
132 :
133 : void restrictN(std::set<int> n);
134 :
135 : void restrictL(int l_min, int l_max);
136 :
137 : void restrictL(std::set<int> l);
138 :
139 : void restrictJ(float j_min, float j_max);
140 :
141 : void restrictJ(std::set<float> j);
142 :
143 : void restrictM(float m_min, float m_max);
144 :
145 : void restrictM(std::set<float> m);
146 :
147 : ////////////////////////////////////////////////////////////////////
148 : /// Method for adding user-defined states //////////////////////////
149 : ////////////////////////////////////////////////////////////////////
150 :
151 : void addStates(const State &s);
152 :
153 : void addStates(const std::set<State> &s);
154 :
155 : // TODO make it possible to just use added states, i.e. use no restrictions on quantum numbers
156 : // and energies
157 :
158 : ////////////////////////////////////////////////////////////////////
159 : /// Methods to get overlaps ////////////////////////////////////////
160 : ////////////////////////////////////////////////////////////////////
161 :
162 : Eigen::VectorX<double> getOverlap(const State &generalizedstate);
163 :
164 : Eigen::VectorX<double> getOverlap(const std::vector<State> &generalizedstates);
165 :
166 : Eigen::VectorX<double> getOverlap(const size_t &state_index);
167 :
168 : Eigen::VectorX<double> getOverlap(const std::vector<size_t> &states_indices);
169 :
170 : Eigen::VectorX<double> getOverlap(const State &generalizedstate,
171 : std::array<double, 3> to_z_axis,
172 : std::array<double, 3> to_y_axis);
173 :
174 : Eigen::VectorX<double> getOverlap(const std::vector<State> &generalizedstates,
175 : std::array<double, 3> to_z_axis,
176 : std::array<double, 3> to_y_axis);
177 :
178 : Eigen::VectorX<double> getOverlap(const size_t &state_index, std::array<double, 3> to_z_axis,
179 : std::array<double, 3> to_y_axis);
180 :
181 : Eigen::VectorX<double> getOverlap(const std::vector<size_t> &states_indices,
182 : std::array<double, 3> to_z_axis,
183 : std::array<double, 3> to_y_axis);
184 :
185 : Eigen::VectorX<double> getOverlap(const State &generalizedstate, double alpha, double beta,
186 : double gamma);
187 :
188 : Eigen::VectorX<double> getOverlap(const size_t &state_index, double alpha, double beta,
189 : double gamma);
190 :
191 : Eigen::VectorX<double> getOverlap(const std::vector<State> &generalizedstates, double alpha,
192 : double beta, double gamma);
193 :
194 : Eigen::VectorX<double> getOverlap(const std::vector<size_t> &states_indices, double alpha,
195 : double beta, double gamma);
196 :
197 : ////////////////////////////////////////////////////////////////////
198 : /// Methods to get properties of the system ////////////////////////
199 : ////////////////////////////////////////////////////////////////////
200 :
201 : std::vector<State> getStates();
202 :
203 : const typename states_set<State>::type &getStatesMultiIndex();
204 :
205 : Eigen::SparseMatrix<Scalar> &getBasisvectors();
206 :
207 : Eigen::SparseMatrix<Scalar> &getHamiltonian();
208 :
209 : size_t getNumBasisvectors();
210 :
211 : size_t getNumStates();
212 :
213 : std::vector<State> getMainStates();
214 :
215 : std::array<std::vector<size_t>, 2> getConnections(SystemBase<Scalar, State> &system_to,
216 : double threshold);
217 :
218 : ////////////////////////////////////////////////////////////////////
219 : /// Methods to build, transform, and destroy the system ////////////
220 : ////////////////////////////////////////////////////////////////////
221 :
222 : void buildHamiltonian();
223 :
224 : void buildInteraction();
225 :
226 : void buildBasis();
227 :
228 : void diagonalize(double energy_lower_bound, double energy_upper_bound);
229 :
230 : void diagonalize(double energy_lower_bound, double energy_upper_bound, double threshold);
231 :
232 : void diagonalize();
233 :
234 : void diagonalize(double threshold);
235 :
236 : void canonicalize();
237 :
238 : void unitarize();
239 :
240 : void rotate(std::array<double, 3> to_z_axis, std::array<double, 3> to_y_axis);
241 :
242 : void rotate(double alpha, double beta, double gamma);
243 :
244 : void add(SystemBase<Scalar, State> &system);
245 :
246 : void constrainBasisvectors(std::vector<size_t> indices_of_wanted_basisvectors);
247 :
248 : void applySchriefferWolffTransformation(SystemBase<Scalar, State> &system0);
249 :
250 : ////////////////////////////////////////////////////////////////////
251 : /// Methods to manipulate entries of the Hamiltonian ///////////////
252 : ////////////////////////////////////////////////////////////////////
253 :
254 : size_t getStateIndex(const State &searched_state);
255 :
256 : std::vector<size_t> getStateIndex(const std::vector<State> &searched_states);
257 :
258 : size_t getBasisvectorIndex(const State &searched_state);
259 :
260 : std::vector<size_t> getBasisvectorIndex(const std::vector<State> &searched_states);
261 :
262 : void forgetStatemixing();
263 :
264 : Scalar getHamiltonianEntry(const State &state_row, const State &state_col);
265 :
266 : void setHamiltonianEntry(const State &state_row, const State &state_col, Scalar value);
267 :
268 : void addHamiltonianEntry(const State &state_row, const State &state_col, Scalar value);
269 :
270 : protected:
271 : SystemBase();
272 :
273 : SystemBase(MatrixElementCache &cache);
274 :
275 : SystemBase(MatrixElementCache &cache, bool memory_saving);
276 :
277 : virtual void initializeBasis() = 0;
278 : virtual void initializeInteraction() = 0;
279 :
280 : virtual void transformInteraction(const Eigen::SparseMatrix<Scalar> &transformator) = 0;
281 : virtual void addInteraction() = 0;
282 : virtual void deleteInteraction() = 0;
283 :
284 : virtual Eigen::SparseMatrix<Scalar> rotateStates(const std::vector<size_t> &states_indices,
285 : double alpha, double beta, double gamma) = 0;
286 : virtual Eigen::SparseMatrix<Scalar> buildStaterotator(double alpha, double beta,
287 : double gamma) = 0;
288 : virtual void incorporate(SystemBase<Scalar, State> &system) = 0;
289 :
290 0 : virtual void onStatesChange(){};
291 :
292 : MatrixElementCache *m_cache{nullptr};
293 :
294 : double threshold_for_sqnorm{0.05};
295 :
296 : double energy_min, energy_max;
297 : std::set<int> range_n, range_l;
298 : std::set<float> range_j, range_m;
299 : std::set<State> states_to_add;
300 :
301 : bool memory_saving{false};
302 : bool is_interaction_already_contained{false};
303 : bool is_new_hamiltonian_required{false};
304 :
305 : typename states_set<State>::type states;
306 : Eigen::SparseMatrix<Scalar> basisvectors;
307 : Eigen::SparseMatrix<Scalar> hamiltonian;
308 : Eigen::SparseMatrix<Scalar> basisvectors_unperturbed_cache;
309 : Eigen::SparseMatrix<Scalar> hamiltonian_unperturbed_cache;
310 :
311 : ////////////////////////////////////////////////////////////////////
312 : /// Helper method that shoul be called by the derived classes //////
313 : ////////////////////////////////////////////////////////////////////
314 :
315 : void onParameterChange();
316 :
317 : void onSymmetryChange();
318 :
319 : ////////////////////////////////////////////////////////////////////
320 : /// Helper methods to check diagonality and unitarity of a matrix //
321 : ////////////////////////////////////////////////////////////////////
322 :
323 : bool checkIsDiagonal(const Eigen::SparseMatrix<Scalar> &mat);
324 :
325 : bool checkIsUnitary(const Eigen::SparseMatrix<Scalar> &mat);
326 :
327 : ////////////////////////////////////////////////////////////////////
328 : /// Helper methods to check the validity of states /////////////////
329 : ////////////////////////////////////////////////////////////////////
330 :
331 : template <class V>
332 0 : bool checkIsQuantumnumberValid(V q, std::set<V> range_q) {
333 0 : return range_q.empty() || range_q.find(q) != range_q.end();
334 : }
335 :
336 : bool checkIsQuantumstateValid(const State &state);
337 :
338 : bool checkIsQuantumstateValid(const StateOne &state, bool a);
339 :
340 : bool checkIsQuantumstateValid(const StateTwo &state, std::array<bool, 2> a);
341 :
342 : bool checkIsEnergyValid(double e);
343 :
344 : ////////////////////////////////////////////////////////////////////
345 : /// Helper methods to change the set of basis vectors //////////////
346 : ////////////////////////////////////////////////////////////////////
347 :
348 : void applyLeftsideTransformator(std::vector<Eigen::Triplet<Scalar>> &triplets_transformator);
349 :
350 : void applyLeftsideTransformator(Eigen::SparseMatrix<Scalar> &transformator);
351 :
352 : void applyRightsideTransformator(std::vector<Eigen::Triplet<Scalar>> &triplets_transformator);
353 :
354 : void applyRightsideTransformator(Eigen::SparseMatrix<Scalar> &transformator);
355 :
356 : void
357 : removeRestrictedStates(std::function<bool(const enumerated_state<State> &)> checkIsValidEntry);
358 :
359 : ////////////////////////////////////////////////////////////////////
360 : /// Utility methods ////////////////////////////////////////////////
361 : ////////////////////////////////////////////////////////////////////
362 :
363 : template <class V>
364 2428 : void range(std::set<V> &rset, V rmin, V rmax) {
365 2428 : rset.clear();
366 9535 : for (V r = rmin; r <= rmax; ++r) {
367 7107 : rset.insert(r);
368 : }
369 2428 : }
370 :
371 : Eigen::Matrix<double, 3, 3> buildRotator(std::array<double, 3> to_z_axis,
372 : std::array<double, 3> to_y_axis);
373 :
374 : Eigen::Matrix<double, 3, 3> buildRotator(const double &alpha, const double &beta,
375 : const double &gamma);
376 :
377 : Eigen::Matrix<double, 3, 1> getEulerAngles(const std::array<double, 3> &to_z_axis,
378 : const std::array<double, 3> &to_y_axis);
379 :
380 : private:
381 : void forgetRestrictions();
382 :
383 : ////////////////////////////////////////////////////////////////////
384 : /// Method to update the system ////////////////////////////////////
385 : ////////////////////////////////////////////////////////////////////
386 :
387 : void updateEverything();
388 :
389 : ////////////////////////////////////////////////////////////////////
390 : /// Method for serialization ///////////////////////////////////////
391 : ////////////////////////////////////////////////////////////////////
392 :
393 : friend class cereal::access;
394 :
395 : template <class Archive>
396 20 : void serialize(Archive &ar, unsigned int /*version*/) {
397 20 : ar &cereal::make_nvp("cache", *m_cache) & CEREAL_NVP(threshold_for_sqnorm);
398 20 : ar &CEREAL_NVP(energy_min) & CEREAL_NVP(energy_max) & CEREAL_NVP(range_n) &
399 40 : CEREAL_NVP(range_l) & CEREAL_NVP(range_j) & CEREAL_NVP(range_m) &
400 20 : CEREAL_NVP(states_to_add);
401 20 : ar &CEREAL_NVP(memory_saving) & CEREAL_NVP(is_interaction_already_contained) &
402 20 : CEREAL_NVP(is_new_hamiltonian_required);
403 20 : ar &CEREAL_NVP(states) & CEREAL_NVP(basisvectors) & CEREAL_NVP(hamiltonian);
404 20 : ar &CEREAL_NVP(basisvectors_unperturbed_cache) & CEREAL_NVP(hamiltonian_unperturbed_cache);
405 20 : }
406 : };
407 :
408 : extern template class SystemBase<std::complex<double>, StateOne>;
409 : extern template class SystemBase<std::complex<double>, StateTwo>;
410 : extern template class SystemBase<double, StateOne>;
411 : extern template class SystemBase<double, StateTwo>;
412 :
413 : #endif
|
