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 "QuantumDefect.hpp"
      21             : #include "Cache.hpp"
      22             : #include "Constants.hpp"
      23             : #include "EmbeddedDatabase.hpp"
      24             : #include "SQLite.hpp"
      25             : 
      26             : #include <algorithm>
      27             : #include <cmath>
      28             : #include <exception>
      29             : #include <fstream>
      30             : #include <iostream>
      31             : #include <sstream>
      32             : #include <string>
      33             : #include <utility>
      34             : 
      35             : struct no_defect : public std::exception {
      36             : private:
      37             :     std::string msg;
      38             : 
      39             : public:
      40           2 :     no_defect(QuantumDefect const &qd)
      41           4 :         : msg{"There is no defect available for " + qd.species + ", l = " + std::to_string(qd.l) +
      42           6 :               ", j = " + std::to_string(qd.j)} {}
      43             : 
      44           2 :     const char *what() const noexcept override { return msg.c_str(); }
      45             : };
      46             : 
      47             : struct no_potential : public std::exception {
      48             : private:
      49             :     std::string msg;
      50             : 
      51             : public:
      52           2 :     no_potential(QuantumDefect const &qd)
      53           4 :         : msg{"There is no model potential available for " + qd.species +
      54           6 :               ", l = " + std::to_string(qd.l)} {}
      55             : 
      56           2 :     const char *what() const noexcept override { return msg.c_str(); }
      57             : };
      58             : 
      59      166114 : QuantumDefect::QuantumDefect(std::string _species, int _n, int _l, double _j, std::nullptr_t)
      60      332228 :     : e(), species(std::move(_species)), n(_n), l(_l), j(_j), ac(e.ac), Z(e.Z), a1(e.a1), a2(e.a2),
      61      166114 :       a3(e.a3), a4(e.a4), rc(e.rc), nstar(e.nstar), energy(e.energy) {}
      62             : 
      63        1237 : QuantumDefect::QuantumDefect(std::string const &species, int n, int l, double j)
      64        1241 :     : QuantumDefect(species, n, l, j, nullptr) {
      65        1237 :     static thread_local EmbeddedDatabase embedded_database{};
      66        1241 :     setup(embedded_database, "");
      67        1233 : }
      68             : 
      69      164877 : QuantumDefect::QuantumDefect(std::string const &species, int n, int l, double j,
      70      164877 :                              std::string const &database)
      71      164877 :     : QuantumDefect(species, n, l, j, nullptr) {
      72      164877 :     if (database.empty()) {
      73      164877 :         static thread_local EmbeddedDatabase embedded_database{};
      74      164877 :         setup(embedded_database, database);
      75             :     } else {
      76           0 :         static thread_local sqlite::handle db(":memory:");
      77           0 :         sqlite::statement stmt(db);
      78           0 :         std::ifstream ifs(database);
      79           0 :         std::stringstream buffer;
      80           0 :         buffer << ifs.rdbuf();
      81           0 :         stmt.exec(buffer.str());
      82           0 :         setup(db, database);
      83             :     }
      84      164877 : }
      85             : 
      86      166114 : void QuantumDefect::setup(sqlite3 *db, std::string const &db_name) {
      87      166114 :     static Cache<Key, Element, Hash> cache;
      88      166114 :     static std::string used_db_name;
      89             : 
      90      166114 :     if (used_db_name != db_name) {
      91           0 :         used_db_name = db_name;
      92           0 :         cache.clear(); // Clear cache
      93             :     }
      94             : 
      95      166118 :     Key const key{species, n, l, j};
      96      166114 :     if (auto oe = cache.restore(key)) {
      97      164674 :         e = oe.value(); // Restore cache
      98      164674 :         return;         // Return early
      99             :     }
     100             : 
     101        2880 :     std::stringstream ss;
     102        2880 :     sqlite::statement stmt(db);
     103             :     int pot_max_l, ryd_max_l;
     104             :     int pot_l, ryd_l;
     105             :     double ryd_max_j;
     106             :     double ryd_j;
     107             : 
     108             :     // Determine maximal L for model potentials
     109        1440 :     stmt.set("select MAX(L) from model_potential where (element = ?);");
     110        1440 :     stmt.prepare();
     111        1440 :     stmt.bind(1, species);
     112        1440 :     if (stmt.step()) {
     113        1440 :         pot_max_l = stmt.get<int>(0);
     114             :     } else {
     115           0 :         throw no_potential(*this);
     116             :     }
     117             : 
     118             :     // The l to be used is the minimum of the two below
     119        1440 :     pot_l = std::min(l, pot_max_l);
     120             : 
     121             :     // Determine maximal L for Rydberg-Ritz coefficients
     122        1440 :     stmt.reset();
     123        1440 :     stmt.set("select MAX(L) from rydberg_ritz where (element = ?);");
     124        1440 :     stmt.prepare();
     125        1440 :     stmt.bind(1, species);
     126        1440 :     if (stmt.step()) {
     127        1440 :         ryd_max_l = stmt.get<int>(0);
     128             :     } else {
     129           0 :         throw no_defect(*this);
     130             :     }
     131             : 
     132             :     // The l to be used is the minimum of the two below
     133        1440 :     ryd_l = std::min(l, ryd_max_l);
     134             : 
     135             :     // Determine maximal J for Rydberg-Ritz coefficients
     136        1440 :     stmt.reset();
     137        1440 :     stmt.set("select MAX(J) from rydberg_ritz where  ((element = ?1) and (L = "
     138        2880 :              "?2));");
     139        1440 :     stmt.prepare();
     140        1440 :     stmt.bind(1, species);
     141        1440 :     stmt.bind(2, ryd_l);
     142        1440 :     if (stmt.step()) {
     143        1440 :         ryd_max_j = stmt.get<double>(0);
     144             :     } else {
     145           0 :         throw no_defect(*this);
     146             :     }
     147             : 
     148             :     // The j to be used is the minimum of the two below
     149        1440 :     ryd_j = std::min(j, ryd_max_j);
     150             : 
     151             :     // Load model potentials from database
     152        1440 :     stmt.reset();
     153        1440 :     stmt.set("select ac,Z,a1,a2,a3,a4,rc from model_potential where ((element "
     154        2880 :              "= ?1) and (L = ?2));");
     155        1440 :     stmt.prepare();
     156        1440 :     stmt.bind(1, species);
     157        1440 :     stmt.bind(2, pot_l);
     158        1440 :     if (stmt.step()) {
     159        1438 :         e.ac = stmt.get<double>(0);
     160        1438 :         e.Z = stmt.get<int>(1);
     161        1438 :         e.a1 = stmt.get<double>(2);
     162        1438 :         e.a2 = stmt.get<double>(3);
     163        1438 :         e.a3 = stmt.get<double>(4);
     164        1438 :         e.a4 = stmt.get<double>(5);
     165        1438 :         e.rc = stmt.get<double>(6);
     166             :     } else {
     167           2 :         throw no_potential(*this);
     168             :     }
     169             : 
     170             :     // Load Rydberg-Ritz coefficients from database
     171        1438 :     stmt.reset();
     172        1438 :     stmt.set("select d0,d2,d4,d6,d8,Ry from rydberg_ritz where ((element = ?1) "
     173        2876 :              "and (L = ?2) and (J = ?3));");
     174        1438 :     stmt.prepare();
     175        1438 :     stmt.bind(1, species);
     176        1438 :     stmt.bind(2, ryd_l);
     177        1438 :     stmt.bind(3, ryd_j);
     178        1438 :     e.nstar = n;
     179        1438 :     double Ry_inf = 109737.31568525;
     180             :     double Ry;
     181        1438 :     if (stmt.step()) {
     182        1436 :         double d0 = stmt.get<double>(0);
     183        1436 :         double d2 = stmt.get<double>(1);
     184        1436 :         double d4 = stmt.get<double>(2);
     185        1436 :         double d6 = stmt.get<double>(3);
     186        1436 :         double d8 = stmt.get<double>(4);
     187        1436 :         Ry = stmt.get<double>(5);
     188        1436 :         e.nstar -= d0 + d2 / std::pow(n - d0, 2) + d4 / std::pow(n - d0, 4) +
     189        1436 :             d6 / std::pow(n - d0, 6) + d8 / std::pow(n - d0, 8);
     190             :     } else {
     191           2 :         throw no_defect(*this);
     192             :     }
     193             : 
     194        1436 :     e.energy = -.5 * (Ry / Ry_inf) / (e.nstar * e.nstar) * au2GHz;
     195             : 
     196        1436 :     cache.save(key, e);
     197             : }
     198             : 
     199      144988 : double energy_level(std::string const &species, int n, int l, double j,
     200             :                     std::string const &database) {
     201      144988 :     QuantumDefect qd(species, n, l, j, database);
     202      289976 :     return qd.energy;
     203             : }
     204             : 
     205       14961 : double nstar(std::string const &species, int n, int l, double j, std::string const &database) {
     206       14961 :     QuantumDefect qd(species, n, l, j, database);
     207       29922 :     return qd.nstar;
     208             : }