Line data Source code
1 : /*
2 : * Copyright (c) 2017 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 "Wavefunction.hpp"
21 : #include "Constants.hpp"
22 : #include "QuantumDefect.hpp"
23 :
24 : #include "EigenCompat.hpp"
25 : #include <Eigen/Core>
26 :
27 : #include <cctype>
28 : #include <cmath>
29 : #ifdef WITH_GSL
30 : #include <gsl/gsl_sf_hyperg.h>
31 : #endif
32 : #include <string>
33 : #include <vector>
34 :
35 : // --- Numerov's method ---
36 :
37 : namespace model_potential {
38 :
39 174374866 : double V(QuantumDefect const &qd, double x) {
40 : double Z_l =
41 174374866 : 1 + (qd.Z - 1) * std::exp(-qd.a1 * x) - x * (qd.a3 + qd.a4 * x) * std::exp(-qd.a2 * x);
42 174374866 : double V_c = -Z_l / x;
43 174374866 : double V_p = -qd.ac / (2. * x * x * x * x) * (1 - std::exp(-std::pow(x / qd.rc, 6)));
44 174374866 : double V_so = 0.0;
45 174374866 : if (qd.l < 4) {
46 148245211 : double s = 0.5;
47 148245211 : if (std::isdigit(qd.species.back()) != 0) {
48 0 : s = ((qd.species.back() - '0') - 1) / 2.; // TODO think of a better solution
49 : }
50 148245211 : double alpha = 7.2973525664e-3; // ~1/137 fine-structure constant CODATA 2014
51 148245211 : V_so =
52 148245211 : alpha * alpha / (4 * x * x * x) * (qd.j * (qd.j + 1) - qd.l * (qd.l + 1) - s * (s + 1));
53 : }
54 174374866 : return V_c + V_p + V_so;
55 : }
56 :
57 174374864 : double g(QuantumDefect const &qd, double x) {
58 174374864 : return (2. * qd.l + .5) * (2. * qd.l + 1.5) / x + 8 * x * (V(qd, x) - qd.energy / au2GHz);
59 : }
60 :
61 : } // namespace model_potential
62 :
63 6136 : Numerov::Numerov(QuantumDefect const &qd) : qd(qd) {
64 : // augmented classical turning point
65 6136 : double xmin = qd.n * qd.n - qd.n * std::sqrt(qd.n * qd.n - (qd.l - 1) * (qd.l - 1));
66 6136 : if (xmin < 2.08) {
67 4900 : xmin = std::floor(std::sqrt(2.08));
68 : } else {
69 1236 : xmin = std::floor(std::sqrt(xmin));
70 : }
71 :
72 6136 : double const xmax = std::sqrt(2 * qd.n * (qd.n + 15));
73 6136 : double const nsteps = std::ceil((xmax - xmin) / dx);
74 :
75 6136 : xy = Eigen::MatrixX<double>::Zero(nsteps, 2);
76 :
77 58143362 : for (int i = 0; i < nsteps; ++i) {
78 58137226 : xy(i, 0) = (xmin + i * dx);
79 : }
80 6136 : }
81 :
82 6136 : Eigen::MatrixX<double> Numerov::integrate() {
83 : using model_potential::g;
84 :
85 6136 : int const nsteps = xy.rows();
86 :
87 : // Set the initial condition
88 6136 : if ((qd.n - qd.l) % 2 == 0) {
89 2938 : xy(nsteps - 2, 1) = -1e-10;
90 : } else {
91 3198 : xy(nsteps - 2, 1) = 1e-10;
92 : }
93 :
94 : // Perform the integration using Numerov's scheme
95 58131090 : for (int i = nsteps - 3; i >= 0; --i) {
96 58124954 : double A = (2. + 5. / 6. * dx * dx * g(qd, xy(i + 1, 0) * xy(i + 1, 0))) * xy(i + 1, 1);
97 58124954 : double B = (1. - 1. / 12. * dx * dx * g(qd, xy(i + 2, 0) * xy(i + 2, 0))) * xy(i + 2, 1);
98 58124954 : double C = 1. - 1. / 12. * dx * dx * g(qd, xy(i, 0) * xy(i, 0));
99 58124954 : xy(i, 1) = (A - B) / C;
100 : }
101 :
102 : // Normalization
103 6136 : double norm = 0;
104 58143362 : for (int i = 0; i < nsteps; ++i) {
105 58137226 : norm += xy(i, 1) * xy(i, 1) * xy(i, 0) * xy(i, 0) * dx;
106 : }
107 6136 : norm = std::sqrt(2 * norm);
108 :
109 6136 : if (norm > 0.0) {
110 58143362 : for (int i = 0; i < nsteps; ++i) {
111 58137226 : xy(i, 1) /= norm;
112 : }
113 : }
114 :
115 6136 : return xy;
116 : }
117 :
118 : // --- Whittaker method ---
119 :
120 : namespace whittaker_functions {
121 :
122 253288 : double HypergeometricU(double a, double b, double z) {
123 : #ifdef WITH_GSL
124 253288 : if (z == 0) {
125 0 : return NAN;
126 : }
127 253288 : return gsl_sf_hyperg_U(a, b, z);
128 : #else
129 : static_cast<void>(a);
130 : static_cast<void>(b);
131 : static_cast<void>(z);
132 : throw std::runtime_error("Whittaker functions require the GSL library!");
133 : #endif
134 : }
135 :
136 253288 : double WhittakerW(double k, double m, double z) {
137 253288 : return std::exp(-.5 * z) * std::pow(z, m + .5) * HypergeometricU(m - k + .5, 1 + 2 * m, z);
138 : }
139 :
140 253288 : double RadialWFWhittaker(double r, double nu, int l) {
141 253288 : return 1 / std::sqrt(nu * nu * std::tgamma(nu + l + 1) * std::tgamma(nu - l)) *
142 253288 : WhittakerW(nu, l + .5, 2 * r / nu);
143 : }
144 :
145 : } // namespace whittaker_functions
146 :
147 30 : Whittaker::Whittaker(QuantumDefect const &qd) : qd(qd) {
148 30 : double const xmin = 1;
149 30 : double const xmax = std::sqrt(2 * qd.n * (qd.n + 15));
150 30 : double const nsteps = std::ceil((xmax - xmin) / dx);
151 :
152 30 : xy.resize(nsteps, 2);
153 :
154 253318 : for (int i = 0; i < nsteps; ++i) {
155 253288 : xy(i, 0) = (xmin + i * dx);
156 : }
157 30 : }
158 :
159 30 : Eigen::MatrixX<double> Whittaker::integrate() {
160 : using whittaker_functions::RadialWFWhittaker;
161 :
162 : // Set the sign
163 : int sign;
164 30 : if ((qd.n - qd.l) % 2 == 0) {
165 25 : sign = -1;
166 : } else {
167 5 : sign = 1;
168 : }
169 :
170 30 : int const nsteps = xy.rows();
171 :
172 253318 : for (int i = 0; i < nsteps; ++i) {
173 253288 : xy(i, 1) = sign * RadialWFWhittaker(xy(i, 0) * xy(i, 0), qd.nstar, qd.l);
174 : }
175 :
176 30 : return xy;
177 : }
|
