LCOV - code coverage report
Current view: top level - pairinteraction - MatrixElementCache.hpp (source / functions) Hit Total Coverage
Test: coverage.info Lines: 30 30 100.0 %
Date: 2024-04-29 00:41:50 Functions: 5 20 25.0 %

          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 MATRIXELEMENTCACHE_H
      21             : #define MATRIXELEMENTCACHE_H
      22             : 
      23             : #include "Basisnames.hpp"
      24             : #include "State.hpp"
      25             : #include "Wavefunction.hpp"
      26             : #include "utils.hpp"
      27             : 
      28             : #include <cereal/types/unordered_map.hpp>
      29             : #include <cereal/types/unordered_set.hpp>
      30             : #include <wignerSymbols/wignerSymbols-cpp.h>
      31             : 
      32             : #include <memory>
      33             : #include <sstream>
      34             : #include <string>
      35             : #include <tuple>
      36             : #include <unordered_map>
      37             : 
      38             : enum method_t {
      39             :     NUMEROV = 0,
      40             :     WHITTAKER = 1,
      41             : };
      42             : 
      43             : class StateOne;
      44             : class StateTwo;
      45             : 
      46             : bool selectionRulesMomentumNew(
      47             :     StateOne const &state1, StateOne const &state2,
      48             :     int q); // TODO rename, integrate into the MatrixElementCache namespace
      49             : bool selectionRulesMomentumNew(StateOne const &state1, StateOne const &state2);
      50             : bool selectionRulesMultipoleNew(StateOne const &state1, StateOne const &state2, int kappa, int q);
      51             : bool selectionRulesMultipoleNew(StateOne const &state1, StateOne const &state2, int kappa);
      52             : 
      53             : class MatrixElementCache {
      54             : public:
      55             :     MatrixElementCache();
      56             :     MatrixElementCache(std::string const &cachedir);
      57             : 
      58             :     double getElectricDipole(StateOne const &state_row,
      59             :                              StateOne const &state_col); // return value in GHz/(V/cm)
      60             :     double getElectricMultipole(StateOne const &state_row, StateOne const &state_col,
      61             :                                 int k); // return value in GHz/(V/cm)*um^(kappa-1)
      62             :     double getDiamagnetism(StateOne const &state_row, StateOne const &state_col,
      63             :                            int k); // return value in GHz/(V/cm)*um
      64             :     double getMagneticDipole(StateOne const &state_row,
      65             :                              StateOne const &state_col); // return value in GHz/G
      66             :     double
      67             :     getElectricMultipole(StateOne const &state_row, StateOne const &state_col, int kappa_radial,
      68             :                          int kappa_angular); // return value in GHz/(V/cm)*um^(kappa_radial-1)
      69             :     double getRadial(StateOne const &state_row, StateOne const &state_col,
      70             :                      int kappa); // return value in um^kappa
      71             : 
      72             :     void precalculateElectricMomentum(const std::vector<StateOne> &basis_one, int q);
      73             :     void precalculateMagneticMomentum(const std::vector<StateOne> &basis_one, int q);
      74             :     void precalculateDiamagnetism(const std::vector<StateOne> &basis_one, int k, int q);
      75             :     void precalculateMultipole(const std::vector<StateOne> &basis_one, int k);
      76             :     void precalculateRadial(const std::vector<StateOne> &basis_one, int k);
      77             : 
      78             :     void setDefectDB(std::string const &path);
      79             :     const std::string &getDefectDB() const;
      80             :     void setMethod(method_t const &m);
      81             :     void loadElectricDipoleDB(std::string const &path, std::string const &species);
      82             : 
      83             :     size_t size();
      84             : 
      85             : private:
      86             :     int update();
      87             :     void precalculate(std::shared_ptr<const BasisnamesOne> basis_one, int kappa, int q, int kappar,
      88             :                       bool calcMultipole, bool calcMomentum, bool calcRadial);
      89             :     double calcRadialElement(const QuantumDefect &qd1, int power, const QuantumDefect &qd2);
      90             :     void precalculate(const std::vector<StateOne> &basis_one, int kappa_angular, int q,
      91             :                       int kappa_radial, bool calcElectricMultipole, bool calcMagneticMomentum,
      92             :                       bool calcRadial);
      93             : 
      94             :     ////////////////////////////////////////////////////////////////////
      95             :     /// Keys ///////////////////////////////////////////////////////////
      96             :     ////////////////////////////////////////////////////////////////////
      97             : 
      98             :     struct CacheKey_cache_radial {
      99             :         CacheKey_cache_radial(method_t method, std::string species, int kappa, int n1, int n2,
     100             :                               int l1, int l2, float j1, float j2);
     101             :         CacheKey_cache_radial();
     102             :         bool operator==(const CacheKey_cache_radial &rhs) const;
     103             :         std::string species;
     104             :         method_t method;
     105             :         int kappa;
     106             :         std::array<int, 2> n, l;
     107             :         std::array<float, 2> j;
     108             : 
     109             :     private:
     110             :         friend class cereal::access;
     111             :         template <class Archive>
     112         570 :         void serialize(Archive &ar, unsigned int /* version */) {
     113         570 :             ar &CEREAL_NVP(method) & CEREAL_NVP(species) & CEREAL_NVP(kappa) & CEREAL_NVP(n) &
     114        1140 :                 CEREAL_NVP(l) & CEREAL_NVP(j);
     115         570 :         }
     116             :     };
     117             : 
     118             :     struct CacheKey_cache_angular {
     119             :         CacheKey_cache_angular(int kappa, float j1, float j2, float m1, float m2);
     120             :         CacheKey_cache_angular();
     121             :         bool operator==(const CacheKey_cache_angular &rhs) const;
     122             :         int kappa;
     123             :         std::array<float, 2> j, m;
     124             :         int sgn;
     125             : 
     126             :     private:
     127             :         friend class cereal::access;
     128             :         template <class Archive>
     129         920 :         void serialize(Archive &ar, unsigned int /* version */) {
     130         920 :             ar &CEREAL_NVP(kappa) & CEREAL_NVP(j) & CEREAL_NVP(m) & CEREAL_NVP(sgn);
     131         920 :         }
     132             :     };
     133             : 
     134             :     struct CacheKey_cache_reduced_commutes {
     135             :         CacheKey_cache_reduced_commutes(float s, int kappa, int l1, int l2, float j1, float j2);
     136             :         CacheKey_cache_reduced_commutes();
     137             :         bool operator==(const CacheKey_cache_reduced_commutes &rhs) const;
     138             :         float s;
     139             :         int kappa;
     140             :         std::array<int, 2> l;
     141             :         std::array<float, 2> j;
     142             :         int sgn;
     143             : 
     144             :     private:
     145             :         friend class cereal::access;
     146             :         template <class Archive>
     147         120 :         void serialize(Archive &ar, unsigned int /* version */) {
     148         120 :             ar &CEREAL_NVP(s) & CEREAL_NVP(kappa) & CEREAL_NVP(l) & CEREAL_NVP(j) & CEREAL_NVP(sgn);
     149         120 :         }
     150             :     };
     151             : 
     152             :     struct CacheKey_cache_reduced_multipole {
     153             :         CacheKey_cache_reduced_multipole(int kappa, int l1, int l2);
     154             :         CacheKey_cache_reduced_multipole();
     155             :         bool operator==(const CacheKey_cache_reduced_multipole &rhs) const;
     156             :         int kappa;
     157             :         std::array<int, 2> l;
     158             :         int sgn;
     159             : 
     160             :     private:
     161             :         friend class cereal::access;
     162             :         template <class Archive>
     163          40 :         void serialize(Archive &ar, unsigned int /* version */) {
     164          40 :             ar &CEREAL_NVP(kappa) & CEREAL_NVP(l) & CEREAL_NVP(sgn);
     165          40 :         }
     166             :     };
     167             : 
     168             :     struct CacheKeyHasher_cache_radial {
     169             :         std::size_t operator()(const CacheKey_cache_radial &c) const;
     170             :     };
     171             :     struct CacheKeyHasher_cache_angular {
     172             :         std::size_t operator()(const CacheKey_cache_angular &c) const;
     173             :     };
     174             :     struct CacheKeyHasher_cache_reduced_commutes {
     175             :         std::size_t operator()(const CacheKey_cache_reduced_commutes &c) const;
     176             :     };
     177             :     struct CacheKeyHasher_cache_reduced_multipole {
     178             :         std::size_t operator()(const CacheKey_cache_reduced_multipole &c) const;
     179             :     };
     180             : 
     181             :     std::unordered_map<CacheKey_cache_radial, double, CacheKeyHasher_cache_radial> cache_radial;
     182             :     std::unordered_map<CacheKey_cache_angular, double, CacheKeyHasher_cache_angular> cache_angular;
     183             :     std::unordered_map<CacheKey_cache_reduced_commutes, double,
     184             :                        CacheKeyHasher_cache_reduced_commutes>
     185             :         cache_reduced_commutes_s;
     186             :     std::unordered_map<CacheKey_cache_reduced_commutes, double,
     187             :                        CacheKeyHasher_cache_reduced_commutes>
     188             :         cache_reduced_commutes_l;
     189             :     std::unordered_map<CacheKey_cache_reduced_multipole, double,
     190             :                        CacheKeyHasher_cache_reduced_multipole>
     191             :         cache_reduced_multipole;
     192             : 
     193             :     std::unordered_set<CacheKey_cache_radial, CacheKeyHasher_cache_radial> cache_radial_missing;
     194             :     std::unordered_set<CacheKey_cache_angular, CacheKeyHasher_cache_angular> cache_angular_missing;
     195             :     std::unordered_set<CacheKey_cache_reduced_commutes, CacheKeyHasher_cache_reduced_commutes>
     196             :         cache_reduced_commutes_s_missing;
     197             :     std::unordered_set<CacheKey_cache_reduced_commutes, CacheKeyHasher_cache_reduced_commutes>
     198             :         cache_reduced_commutes_l_missing;
     199             :     std::unordered_set<CacheKey_cache_reduced_multipole, CacheKeyHasher_cache_reduced_multipole>
     200             :         cache_reduced_multipole_missing;
     201             : 
     202             :     method_t method{NUMEROV};
     203             :     std::string defectdbname;
     204             :     std::string dbname;
     205             :     std::unique_ptr<sqlite::handle> db;
     206             :     std::unique_ptr<sqlite::statement> stmt;
     207             :     long pid_which_created_db;
     208             : 
     209             :     ////////////////////////////////////////////////////////////////////
     210             :     /// Method for serialization ///////////////////////////////////////
     211             :     ////////////////////////////////////////////////////////////////////
     212             : 
     213             :     friend class cereal::access;
     214             :     template <class Archive>
     215          20 :     void serialize(Archive &ar, unsigned int /* version */) {
     216          20 :         ar &CEREAL_NVP(method);
     217          20 :         ar &CEREAL_NVP(dbname);
     218          20 :         ar &CEREAL_NVP(cache_radial) & CEREAL_NVP(cache_angular) &
     219          40 :             CEREAL_NVP(cache_reduced_commutes_s) & CEREAL_NVP(cache_reduced_commutes_l) &
     220          20 :             CEREAL_NVP(cache_reduced_multipole);
     221          20 :         ar &CEREAL_NVP(cache_radial_missing) & CEREAL_NVP(cache_angular_missing) &
     222          40 :             CEREAL_NVP(cache_reduced_commutes_s_missing) &
     223          40 :             CEREAL_NVP(cache_reduced_commutes_l_missing) &
     224          20 :             CEREAL_NVP(cache_reduced_multipole_missing);
     225             : 
     226          20 :         if (Archive::is_loading::value && !dbname.empty()) {
     227             :             // Open database
     228          20 :             db = std::make_unique<sqlite::handle>(dbname);
     229          20 :             stmt = std::make_unique<sqlite::statement>(*db);
     230          20 :             pid_which_created_db = utils::get_pid();
     231             : 
     232             :             // Speed up database access
     233          20 :             stmt->exec(
     234             :                 "PRAGMA synchronous = OFF"); // do not wait on write, hand off to OS and carry on
     235          20 :             stmt->exec("PRAGMA journal_mode = MEMORY"); // keep rollback journal in memory during
     236             :                                                         // transaction
     237             :         }
     238          20 :     }
     239             : };
     240             : 
     241             : #endif

Generated by: LCOV version 1.14