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 : #ifndef WAVEFUNCTION_H
21 : #define WAVEFUNCTION_H
22 :
23 : #include "QuantumDefect.hpp"
24 :
25 : #include "EigenCompat.hpp"
26 : #include <Eigen/Core>
27 :
28 : #include <cstdio>
29 : #include <limits>
30 : #include <string>
31 : #include <vector>
32 :
33 : // --- Numerov's method ---
34 :
35 : namespace model_potential {
36 : /** \brief Calculate the model potentials.
37 : *
38 : * This function calculates the model potentials and returns \f$
39 : * V_{\mathrm{mod}}(r) = V_{\mathrm{C}}(r) + V_{\mathrm{P}}(r) +
40 : * V_{\mathrm{s.o.}}(r) \f$
41 : *
42 : * \f[
43 : * \begin{aligned}
44 : * V_{\mathrm{C}}(r) &= - \frac{e^2}{4\pi\varepsilon_0} \frac{1 +
45 : * (Z-1)\mathrm{e}^{-\alpha_1 r} - r(\alpha_3+\alpha_4 r)\mathrm{e}^{-\alpha_2
46 : * r}}{r} \\
47 : * V_{\mathrm{P}}(r) &= -\frac{e^2}{(4\pi\varepsilon_0)^2}
48 : * \frac{\alpha_d}{2 r^4} \left[ 1 - \mathrm{e}^{-(r/r_c)^6} \right] \\
49 : * V_{\mathrm{s.o.}}(r) &= \frac{1}{2} \left( \frac{e^2}{4\pi\varepsilon_0}
50 : * \right) \left( \frac{g_s}{2 m_e^2 c^2} \right)
51 : * \frac{\boldsymbol{l}\cdot\boldsymbol{s}}{r^3} \end{aligned} \f]
52 : *
53 : * \param[in] qd Quantum defect data (parameters for the potentials)
54 : * \param[in] x Radial position
55 : * \returns Model potentials evaluted at position \p x with parameters \p qd
56 : */
57 : double V(QuantumDefect const &qd, double x);
58 :
59 : /** \brief Helper function for %Numerov's method.s
60 : *
61 : * %Numerov's method solves \f$ y''(x) + g(x) y(x) = 0\f$. This function
62 : * implements \f$ g(x) \f$.
63 : *
64 : * \f[
65 : * g(x) = \frac{(2 l + 1/2)^2}{x} + 8 x (V_{\mathrm{mod}}(x) - E)
66 : * \f]
67 : *
68 : * \param[in] qd Quantum defect data (parameters for the potentials)
69 : * \param[in] x Radial position
70 : * \returns Scaling coefficient evaluted at position \p x with parameters \p qd
71 : */
72 : double g(QuantumDefect const &qd, double x);
73 : } // namespace model_potential
74 :
75 : /** \brief %Numerov's method
76 : *
77 : * This class implements %Numerov's method using a couple of helper functions.
78 : * %Numerov's method solves a differential equation of the form
79 : * \f[
80 : * y''(x) + g(x)y(x) = 0
81 : * \f]
82 : * The equation is solved by iteratively computing
83 : * \f[
84 : * y_{n+1} = \frac{\left( 2 + \frac{5 h^2}{6} g_n \right) y_n - \left( 1 +
85 : * \frac{h^2}{12} g_{n-1} \right) y_{n-1}}{\left( 1 + \frac{h^2}{12} g_n
86 : * \right)} \f]
87 : *
88 : * with \f$ y_n = y(x_n) \f$ and \f$ g_n = g(x_n) \f$. Clearly this equation
89 : * needs two initial conditions \f$ y_0 \f$ and \f$ y_1 \f$. We integrate the
90 : * equation from outside to inside. Since the wavefunction hs to decay to zero
91 : * at infinity we choose \f$ y_0 = 0 \f$. Now depending on the number of knots
92 : * of wavefunction we decide whether to set \f$ y_1 = \pm\varepsilon \f$.
93 : */
94 : class Numerov {
95 : QuantumDefect const &qd;
96 : Eigen::MatrixX<double> xy;
97 :
98 : public:
99 : /** \brief Integration step size */
100 : static constexpr double const dx = 0.01;
101 :
102 : /** \brief Constructor
103 : *
104 : * The constructor initializes the \p qd member and sets the initial
105 : * condition for the integrator. It determines the outer and inner
106 : * integration bounds by semiclassical and empirical arguments.
107 : *
108 : * \param[in] qd Quantum defect data (parameters for the potentials)
109 : */
110 : Numerov(QuantumDefect const &qd);
111 :
112 : /** \brief Perform the integration
113 : *
114 : * This performs the integration using %Numerov's method.
115 : *
116 : * \returns vector with wavefunction amplitude
117 : */
118 : Eigen::MatrixX<double> integrate();
119 :
120 : /** \brief Power kernel for matrix elements
121 : *
122 : * The power kernel accounts for the fact that in %Numerov's method the
123 : * domain is square root scaled. This is important for the calculation of
124 : * matrix elements.
125 : *
126 : * \param[in] power exponent of r
127 : * \returns power kernel
128 : */
129 28766633 : constexpr static inline int power_kernel(int power) { return 2 * power + 2; }
130 : };
131 :
132 : // --- Whittaker method ---
133 :
134 : namespace whittaker_functions {
135 : /** \brief Compute the confluent hypergeometric function
136 : *
137 : * This is merely a wrapper around
138 : * <a
139 : * href="https://www.gnu.org/software/gsl/manual/html_node/Hypergeometric-Functions.html">
140 : * <tt>gsl_sf_hyperg_U</tt>
141 : * </a>
142 : *
143 : * \param[in] a see documentation for <tt>gsl_sf_hyperg_U</tt>
144 : * \param[in] b see documentation for <tt>gsl_sf_hyperg_U</tt>
145 : * \param[in] z see documentation for <tt>gsl_sf_hyperg_U</tt>
146 : * \returns U(a,b,z)
147 : */
148 : double HypergeometricU(double a, double b, double z);
149 :
150 : /** \brief Compute the %Whittaker function
151 : *
152 : * The %Whittaker function is defined in terms of the confluent hypergeometric
153 : * function.
154 : * \f[
155 : * W_{k,m}(z) = \mathrm{e}^{-z/2} z^{m+1/2} U\left(m-k+\frac{1}{2}, 1+2m,
156 : * z\right) \f]
157 : *
158 : * \param[in] k parameter k from %Whittaker's equation
159 : * \param[in] m parameter m from %Whittaker's equation
160 : * \param[in] z radial position
161 : * \returns W(k,m,z)
162 : */
163 : double WhittakerW(double k, double m, double z);
164 :
165 : /** \brief Radial wavefunctions from %Whittaker's function
166 : *
167 : * The radial wavefunction of the Schrödinger equation with the Coulomb
168 : * potential can be expressed in terms of the %Whittaker function.
169 : * \f[
170 : * R_{\nu,l}(r) = \frac{1}{\sqrt{\nu^2 \Gamma(\nu+l+1) \Gamma(\nu-l)}}
171 : * W_{\nu, l+1/2}\left(\frac{2r}{\nu}\right) \f]
172 : *
173 : * \param[in] nu effective principal quantum number \f$ \nu = n -
174 : * \delta_{n,l} \f$ \param[in] l angular quantum number \param[in] r radial
175 : * position \returns R(nu,l,r)
176 : */
177 : double RadialWFWhittaker(double r, double nu, int l);
178 : } // namespace whittaker_functions
179 :
180 : class Whittaker {
181 : QuantumDefect const &qd;
182 : Eigen::MatrixX<double> xy;
183 :
184 : public:
185 : /** \brief Integration step size */
186 : static constexpr double const dx = 0.01;
187 :
188 : /** \brief Constructor
189 : *
190 : * The constructor initializes the \p qd member and sets the initial
191 : * condition for the integrator. It determines the outer integration bound
192 : * by semiclassical and empirical arguments. The inner intergration
193 : * boundary is set to 1.
194 : *
195 : * \param[in] qd Quantum defect data (parameters for the potentials)
196 : */
197 : Whittaker(QuantumDefect const &qd);
198 :
199 : /** \brief Evaluate the radial wavefunction
200 : *
201 : * This does not perform an integration but merely evaluates the %Whittaker
202 : * functions over the domain.
203 : *
204 : * \returns vector with wavefunction amplitude
205 : */
206 : Eigen::MatrixX<double> integrate();
207 :
208 : /** \brief Power kernel for matrix elements
209 : *
210 : * The power kernel accounts for the fact that for the %Whittaker functions
211 : * the domain is linear scaled. This is important for the calculation of
212 : * matrix elements.
213 : *
214 : * \param[in] power exponent of r
215 : * \returns power kernel
216 : */
217 92362 : constexpr static inline double power_kernel(int power) { return 2 * power + 1; }
218 : };
219 :
220 : // --- Matrix element calculation ---
221 :
222 : /** \brief Find and return index
223 : *
224 : * Find a value in an Eigen matrix and return its index in the
225 : * container to the caller. This is intended to work with
226 : * Eigen::Vector or columns of Eigen::Matrix.
227 : *
228 : * \param[in] x An Eigen matrix type
229 : * \param[in] d value to find
230 : * \returns index of d
231 : * \throws std::runtime_error if the value can't be found
232 : */
233 : template <typename T>
234 12308 : int findidx(T const &x, typename T::Scalar const &d,
235 : typename T::Scalar const eps = std::numeric_limits<typename T::Scalar>::epsilon()) {
236 12308 : int L = 0;
237 12308 : int R = x.rows() - 1;
238 147901 : for (;;) {
239 160209 : if (L > R) {
240 : // If we were already restarted crash. It doesn't make sense to
241 : // increase epsilon further.
242 0 : if (eps > std::numeric_limits<typename T::Scalar>::epsilon()) {
243 0 : throw std::runtime_error("Search failed");
244 : }
245 : // Restart search with larger epsilon
246 0 : fprintf(stderr, "WARNING: Restarting search with %.16f\n%s:%d in %s\n", 2 * eps,
247 : __FILE__, __LINE__, __func__);
248 0 : return findidx(x, d, 2 * eps);
249 : }
250 160209 : int m = (L + R) / 2;
251 160209 : if (x(m) < d) {
252 70354 : L = m + 1;
253 89855 : } else if (x(m) > d) {
254 77547 : R = m - 1;
255 : } else {
256 12308 : return m;
257 : }
258 : // x(m) and d might not actually be equal in the floating point sense.
259 : // We do fuzzy compare to account for that.
260 147901 : if (std::abs(x(m) - d) <= std::max(std::abs(x(m)), std::abs(d)) * eps) {
261 0 : return m;
262 : }
263 : }
264 : }
265 :
266 : /** \brief Compute radial matrix elements
267 : *
268 : * The radial matrix element can be calculated from the integral
269 : * \f[
270 : * \langle nlj| \hat{p}^{\mathrm{rad}}_\kappa |n'l'j' \rangle =
271 : * \int \Psi^{\text{rad}}_{nlj}(r) \Psi^{\text{rad}}_{n'l'j'}(r)
272 : * r^{2+\kappa} \; dr \f] The part \f$ r^{2+\kappa} \f$ varies with rescaling
273 : * the domain and is this given by the power_kernel function.
274 : *
275 : * \tparam T Method for calculating the wavefunction
276 : * \param[in] qd1 Quantum defect data for first atom
277 : * \param[in] power Exponent kappa
278 : * \param[in] qd2 Quantum defect data for second atom
279 : * \returns Radial matrix element
280 : */
281 : template <typename T>
282 3077 : double IntegrateRadialElement(QuantumDefect const &qd1, int power, QuantumDefect const &qd2) {
283 6154 : T N1(qd1);
284 6154 : T N2(qd2);
285 :
286 6154 : auto const &xy1 = N1.integrate();
287 3077 : auto const &xy2 = N2.integrate();
288 3077 : auto const dx = T::dx;
289 :
290 3077 : auto const xmin = xy1(0, 0) >= xy2(0, 0) ? xy1(0, 0) : xy2(0, 0);
291 3077 : auto const xmax = xy1(xy1.rows() - 1, 0) <= xy2(xy2.rows() - 1, 0) ? xy1(xy1.rows() - 1, 0)
292 0 : : xy2(xy2.rows() - 1, 0);
293 :
294 3077 : double mu = 0;
295 : // If there is an overlap, calculate the matrix element
296 3077 : if (xmin <= xmax) {
297 3077 : int start1 = findidx(xy1.col(0), xmin);
298 3077 : int end1 = findidx(xy1.col(0), xmax);
299 3077 : int start2 = findidx(xy2.col(0), xmin);
300 3077 : int end2 = findidx(xy2.col(0), xmax);
301 :
302 : int i1, i2;
303 28862072 : for (i1 = start1, i2 = start2; i1 < end1 && i2 < end2; ++i1, ++i2) {
304 28858995 : mu += xy1(i1, 1) * xy2(i2, 1) * std::pow(xy1(i1, 0), T::power_kernel(power)) * dx;
305 : }
306 3077 : mu = 2 * mu;
307 : }
308 :
309 : // The radial matrix element is returned in atomic units
310 6154 : return mu;
311 : }
312 :
313 : #endif // WAVEFUNCTION_H
|
