SIRIUS 7.5.0
Electronic structure library and applications
radial_solver.hpp
Go to the documentation of this file.
1// Copyright (c) 2013-2016 Anton Kozhevnikov, Thomas Schulthess
2// All rights reserved.
3//
4// Redistribution and use in source and binary forms, with or without modification, are permitted provided that
5// the following conditions are met:
6//
7// 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the
8// following disclaimer.
9// 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions
10// and the following disclaimer in the documentation and/or other materials provided with the distribution.
11//
12// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED
13// WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
14// PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
15// ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
16// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
17// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
18// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
19
20/** \file radial_solver.hpp
21 *
22 * \brief Contains declaration and partial implementation of sirius::Radial_solver class.
23 */
24
25#ifndef __RADIAL_SOLVER_HPP__
26#define __RADIAL_SOLVER_HPP__
27
28#include <tuple>
29#include "spline.hpp"
30#include "core/constants.hpp"
31#include "core/typedefs.hpp"
32
33namespace sirius {
34
35/// Finds a solution to radial Schrodinger, Koelling-Harmon or Dirac equation.
36/** Non-relativistic radial Schrodinger equation:
37 * \f[
38 * -\frac{1}{2}p''(r) + V_{eff}(r)p(r) = Ep(r)
39 * \f]
40 * where \f$ V_{eff}(r) = V(r) + \frac{\ell (\ell+1)}{2r^2} \f$, \f$ p(r) = u(r)r \f$ and \f$ u(r) \f$ is the
41 * radial wave-function. Energy derivatives of radial solutions obey the slightly different equation:
42 * \f[
43 * -\frac{1}{2}\dot{p}''(r) + V_{eff}(r)\dot{p}(r) = E\dot{p}(r) + p(r)
44 * \f]
45 * \f[
46 * -\frac{1}{2}\ddot{p}''(r) + V_{eff}(r)\ddot{p}(r) = E\ddot{p}(r) + 2\dot{p}(r)
47 * \f]
48 * and so on. We can generalize the radial Schrodinger equation like this:
49 * \f[
50 * -\frac{1}{2}p''(r) + \big(V_{eff}(r) - E\big) p(r) = \chi(r)
51 * \f]
52 * where now \f$ p(r) \f$ represents m-th energy derivative of the radial solution and
53 * \f$ \chi(r) = m \frac{\partial^{m-1} p(r)} {\partial^{m-1}E} \f$
54 *
55 * Let's now decouple second-order differential equation into a system of two first-order euquations. From
56 * \f$ p(r) = u(r)r \f$ we have
57 * \f[
58 * p'(r) = u'(r)r + u''(r) = 2q(r) + \frac{p(r)}{r}
59 * \f]
60 * where we have introduced a new variable \f$ q(r) = \frac{u'(r) r}{2} \f$. Differentiating \f$ p'(r) \f$ again
61 * we arrive to the following equation for \f$ q'(r) \f$:
62 * \f[
63 * p''(r) = 2q'(r) + \frac{p'(r)}{r} - \frac{p(r)}{r^2} = 2q'(r) + \frac{2q(r)}{r} + \frac{p(r)}{r^2} -
64 * \frac{p(r)}{r^2} \f] \f[ q'(r) = \frac{1}{2}p''(r) - \frac{q(r)}{r} = \big(V_{eff}(r) - E\big) p(r) - \frac{q(r)}{r}
65 * - \chi(r) \f] Final expression for a linear system of differential equations for m-th energy derivative is:
66 * \f{eqnarray*}{
67 * p'(r) &=& 2q(r) + \frac{p(r)}{r} \\
68 * q'(r) &=& \big(V_{eff}(r) - E\big) p(r) - \frac{q(r)}{r} - \chi(r)
69 * \f}
70 * Scalar-relativistic equations look similar. For m = 0 (no energy derivative) we have:
71 * \f{eqnarray*}{
72 * p'(r) &=& 2Mq(r) + \frac{p(r)}{r} \\
73 * q'(r) &=& \big(V(r) - E + \frac{\ell(\ell+1)}{2Mr^2}\big) p(r) - \frac{q(r)}{r}
74 * \f}
75 * where \f$ M = 1 + \frac{\alpha^2}{2}\big(E-V(r)\big) \f$. Because M depends on energy, the m-th energy
76 * derivatives of the scalar-relativistic solution take a slightly different form. For m=1 we have:
77 * \f{eqnarray*}{
78 * \dot{p}'(r) &=& 2M\dot{q}(r) + \frac{\dot{p}(r)}{r} + \alpha^2 q(r) \\
79 * \dot{q}'(r) &=& \big(V(r) - E + \frac{\ell(\ell+1)}{2Mr^2}\big) \dot{p}(r) - \frac{\dot{q}(r)}{r} -
80 * \big(1 + \frac{\ell(\ell+1)\alpha^2}{4M^2r^2}\big)p(r)
81 * \f}
82 * For m=2:
83 * \f{eqnarray*}{
84 * \ddot{p}'(r) &=& 2M\ddot{q}(r) + \frac{\ddot{p}(r)}{r} + 2 \alpha^2 \dot{q}(r) \\
85 * \ddot{q}'(r) &=& \big(V(r) - E + \frac{\ell(\ell+1)}{2Mr^2}\big) \ddot{p}(r) - \frac{\ddot{q}(r)}{r} - 2 \big(1 +
86 * \frac{\ell(\ell+1)\alpha^2}{4M^2r^2}\big)\dot{p}(r)
87 * + \frac{\ell(\ell+1)\alpha^4}{4M^3r^2} p(r)
88 * \f}
89 *
90 * Derivation of IORA.
91 *
92 * First thing to do is to expand the \f$ 1/M \f$ to the first order in \f$ E \f$:
93 * \f[
94 * \frac{1}{M} = \frac{1}{1-\frac{\alpha^2}{2} V(r)} - \frac{\alpha^2 E}{2 (1-\frac{\alpha^2}{2} V(r))^2} + O(E^2) =
95 * \frac{1}{M_0} - \frac{\alpha^2}{2} \frac{1}{M_0^2} E
96 * \f]
97 * From this expansion we can derive the expression for \f$ M \f$:
98 * \f[
99 * M = \Big(\frac{1}{M}\Big)^{-1} = \frac{M_0}{1 - \frac{\alpha^2}{2}\frac{E}{M_0}}
100 * \f]
101 * Now we insert this expressions into the scalar-relativistic equations:
102 * \f[
103 * \left\{
104 * \begin{array}{ll}
105 * \displaystyle p'(r) = 2 \frac{M_0}{1 - \frac{\alpha^2}{2}\frac{E}{M_0}} q(r) + \frac{p(r)}{r} \\
106 * \displaystyle q'(r) = \Big(V(r) - E + \frac{\ell(\ell+1)}{2r^2} \big(\frac{1}{M_0} - \frac{\alpha^2}{2}
107 * \frac{E}{M_0^2} \big) \Big) p(r) - \frac{q(r)}{r}
108 * = \Big(V(r) - E + \frac{\ell(\ell+1)}{2 M_0 r^2} - \frac{\alpha^2}{2} \frac{\ell(\ell+1) E}{2 M_0^2 r^2} \Big)
109 * p(r) - \frac{q(r)}{r} \end{array} \right. \f] The first energy derivatives are then: \f[ \left\{ \begin{array}{ll}
110 * \displaystyle \dot{p}'(r) = 2 \frac{M_0}{1 - \frac{\alpha^2}{2}\frac{E}{M_0}} \dot{q}(r) +
111 * \frac{\dot{p}(r)}{r} + \frac{\alpha^2}{(1 - \frac{\alpha^2 E}{2 M_0})^2} q(r) \\
112 * \displaystyle \dot{q}'(r) = \Big(V(r) - E + \frac{\ell(\ell+1)}{2 M_0 r^2} - \frac{\alpha^2}{2}
113 * \frac{\ell(\ell+1)}{2 M_0^2 r^2} E \Big) \dot{p}(r) -
114 * \frac{\dot{q}(r)}{r} -
115 * \Big(1 + \frac{\alpha^2}{2} \frac{\ell(\ell+1)}{2 M_0^2 r^2} \Big) p(r)
116 * \end{array}
117 * \right.
118 * \f]
119 * The second energy derivatives are:
120 * \f[
121 * \left\{
122 * \begin{array}{ll}
123 * \displaystyle \ddot{p}'(r) = 2 \frac{M_0}{1 - \frac{\alpha^2}{2}\frac{E}{M_0}} \ddot{q}(r) +
124 * \frac{\ddot{p}(r)}{r} +
125 * \frac{2 \alpha^2}{(1 - \frac{\alpha^2 E}{2 M_0})^2} \dot{q}(r) +
126 * \frac{\alpha^4}{2 M_0 (1 - \frac{\alpha^2 E}{2 M_0})^3} q(r) \\
127 * \displaystyle \ddot{q}'(r) = \Big(V(r) - E + \frac{\ell(\ell+1)}{2 M_0 r^2} - \frac{\alpha^2}{2}
128 * \frac{\ell(\ell+1)}{2 M_0^2 r^2} E \Big) \ddot{p}(r) -
129 * \frac{\ddot{q}(r)}{r} -
130 * 2 \Big(1 + \frac{\alpha^2}{2} \frac{\ell(\ell+1)}{2 M_0^2 r^2} \Big) \dot{p}(r)
131 * \end{array}
132 * \right.
133 * \f]
134 */
135class Radial_solver
136{
137 protected:
138 /// Positive charge of the nucleus.
139 int zn_;
140
141 /// Radial grid.
142 Radial_grid<double> const& radial_grid_;
143
144 /// Electronic part of potential.
145 Spline<double> ve_;
146
147 /// Integrate system of two first-order differential equations forward starting from the origin.
148 /** Use Runge-Kutta 4th order method */
149 template <relativity_t rel, bool prevent_overflow>
150 int integrate_forward_rk4(double enu__, int l__, int k__, Spline<double> const& chi_p__,
151 Spline<double> const& chi_q__, std::vector<double>& p__, std::vector<double>& dpdr__,
152 std::vector<double>& q__, std::vector<double>& dqdr__) const
153 {
154 /* number of mesh points */
155 int nr = num_points();
156
157 double rest_energy = std::pow(speed_of_light, 2);
158
159 double alpha = 1.0 / speed_of_light;
160
161 double sq_alpha_half = 0.5 / rest_energy;
162
163 if (rel == relativity_t::none) {
164 sq_alpha_half = 0;
165 }
166
167 double ll_half = l__ * (l__ + 1) / 2.0;
168
169 double kappa{0};
170 if (rel == relativity_t::dirac) {
171 if (k__ == l__) {
172 kappa = k__;
173 } else if (k__ == l__ + 1) {
174 kappa = -k__;
175 } else {
176 throw std::runtime_error("wrong k");
177 }
178 }
179
180 auto rel_mass = [sq_alpha_half](double enu__, double v__) -> double {
181 switch (rel) {
182 case relativity_t::none: {
183 return 1.0;
184 }
185 case relativity_t::koelling_harmon: {
186 return 1.0 + sq_alpha_half * (enu__ - v__);
187 }
188 case relativity_t::zora: {
189 return 1.0 - sq_alpha_half * v__;
190 }
191 case relativity_t::iora: {
192 double m0 = 1.0 - sq_alpha_half * v__;
193 return m0 / (1 - sq_alpha_half * enu__ / m0);
194 }
195 default: {
196 return 1.0;
197 }
198 }
199 };
200
201 /* try to find classical turning point */
202 int idx_ctp{-1};
203 for (int ir = 0; ir < nr; ir++) {
204 if (ve_(ir) - zn_ * radial_grid_.x_inv(ir) > enu__) {
205 idx_ctp = ir;
206 break;
207 }
208 }
209 /* if we didn't fint the classical turning point, take half of the grid */
210 // int idx_tail{-1};
211 if (idx_ctp == -1) {
212 double rmid = radial_grid_.last() / 2;
213 for (int ir = 0; ir < nr; ir++) {
214 if (radial_grid_[ir] > rmid) {
215 idx_ctp = ir;
216 break;
217 }
218 }
219 }
220
221 /* here and below var0 means var(x), var2 means var(x+h) and var1 means var(x+h/2) */
222 double x2 = radial_grid_[0];
223 double xinv2 = radial_grid_.x_inv(0);
224 double v2 = ve_(0) - zn_ / x2;
225 double M2 = rel_mass(enu__, v2);
226 double chi_p2 = chi_p__(0);
227 double chi_q2 = chi_q__(0);
228
229 /* r->0 asymptotics */
230 if (rel != relativity_t::dirac) {
231 if (l__ == 0) {
232 p__[0] = 2 * zn_ * x2;
233 q__[0] = -std::pow(zn_, 2) * x2;
234 } else {
235 p__[0] = std::pow(x2, l__ + 1);
236 q__[0] = std::pow(x2, l__) * l__ / 2;
237 }
238 } else {
239 double b = std::sqrt(std::pow(kappa, 2) - std::pow(zn_ / speed_of_light, 2));
240 p__[0] = std::pow(x2, b);
241 q__[0] = std::pow(x2, b) * speed_of_light * (b + kappa) / zn_;
242 }
243
244 // p__[0] = std::pow(radial_grid_[0], l_ + 1);
245 // if (l_ == 0)
246 //{
247 // q__[0] = -zn_ * radial_grid_[0] / M2 / 2;
248 //}
249 // else
250 //{
251 // q__[0] = std::pow(radial_grid_[0], l_) * l_ / M2 / 2;
252 //}
253
254 double p2 = p__[0];
255 double q2 = q__[0];
256
257 double pk[4];
258 double qk[4];
259
260 int last{0};
261
262 for (int i = 0; i < nr - 1; i++) {
263 /* copy previous values */
264 double x0 = x2;
265 double xinv0 = xinv2;
266 double M0 = M2;
267 double chi_p0 = chi_p2;
268 double chi_q0 = chi_q2;
269 double v0 = v2;
270 double p0 = p2;
271 double q0 = q2;
272 /* radial grid step */
273 double h = radial_grid_.dx(i);
274 /* mid-point */
275 double h_half = h / 2;
276 double x1 = x0 + h_half;
277 double xinv1 = 1.0 / x1;
278 double chi_p1 = chi_p__(i, h_half);
279 double chi_q1 = chi_q__(i, h_half);
280 double v1 = ve_(i, h_half) - zn_ * xinv1;
281 double M1 = rel_mass(enu__, v1);
282
283 /* next point */
284 x2 = radial_grid_[i + 1];
285 xinv2 = radial_grid_.x_inv(i + 1);
286 chi_p2 = chi_p__(i + 1);
287 chi_q2 = chi_q__(i + 1);
288 v2 = ve_(i + 1) - zn_ * xinv2;
289 M2 = rel_mass(enu__, v2);
290
291 if (rel == relativity_t::none || rel == relativity_t::koelling_harmon || rel == relativity_t::zora) {
292 /* k0 = F(Y(x), x) */
293 pk[0] = 2 * M0 * q0 + p0 * xinv0;
294 qk[0] = (v0 - enu__ + ll_half / M0 / std::pow(x0, 2)) * p0 - q0 * xinv0 + chi_q0;
295
296 /* k1 = F(Y(x) + k0 * h/2, x + h/2) */
297 pk[1] = 2 * M1 * (q0 + qk[0] * h_half) + (p0 + pk[0] * h_half) * xinv1;
298 qk[1] = (v1 - enu__ + ll_half / M1 / std::pow(x1, 2)) * (p0 + pk[0] * h_half) -
299 (q0 + qk[0] * h_half) * xinv1 + chi_q1;
300
301 /* k2 = F(Y(x) + k1 * h/2, x + h/2) */
302 pk[2] = 2 * M1 * (q0 + qk[1] * h_half) + (p0 + pk[1] * h_half) * xinv1;
303 qk[2] = (v1 - enu__ + ll_half / M1 / std::pow(x1, 2)) * (p0 + pk[1] * h_half) -
304 (q0 + qk[1] * h_half) * xinv1 + chi_q1;
305
306 /* k3 = F(Y(x) + k2 * h, x + h) */
307 pk[3] = 2 * M2 * (q0 + qk[2] * h) + (p0 + pk[2] * h) * xinv2;
308 qk[3] = (v2 - enu__ + ll_half / M2 / std::pow(x2, 2)) * (p0 + pk[2] * h) - (q0 + qk[2] * h) * xinv2 +
309 chi_q2;
310 }
311 if (rel == relativity_t::iora) {
312 double m0 = 1 - sq_alpha_half * v0;
313 double m1 = 1 - sq_alpha_half * v1;
314 double m2 = 1 - sq_alpha_half * v2;
315
316 double a0 = ll_half / m0 / std::pow(x0, 2);
317 double a1 = ll_half / m1 / std::pow(x1, 2);
318 double a2 = ll_half / m2 / std::pow(x2, 2);
319
320 /* k0 = F(Y(x), x) */
321 pk[0] = 2 * M0 * q0 + p0 * xinv0 + chi_p0;
322 qk[0] = (v0 - enu__ + a0 - sq_alpha_half * a0 * enu__ / m0) * p0 - q0 * xinv0 + chi_q0;
323
324 /* k1 = F(Y(x) + k0 * h/2, x + h/2) */
325 pk[1] = 2 * M1 * (q0 + qk[0] * h_half) + (p0 + pk[0] * h_half) * xinv1 + chi_p1;
326 qk[1] = (v1 - enu__ + a1 - sq_alpha_half * a1 * enu__ / m1) * (p0 + pk[0] * h_half) -
327 (q0 + qk[0] * h_half) * xinv1 + chi_q1;
328
329 /* k2 = F(Y(x) + k1 * h/2, x + h/2) */
330 pk[2] = 2 * M1 * (q0 + qk[1] * h_half) + (p0 + pk[1] * h_half) * xinv1 + chi_p1;
331 qk[2] = (v1 - enu__ + a1 - sq_alpha_half * a1 * enu__ / m1) * (p0 + pk[1] * h_half) -
332 (q0 + qk[1] * h_half) * xinv1 + chi_q1;
333
334 /* k3 = F(Y(x) + k2 * h, x + h) */
335 pk[3] = 2 * M2 * (q0 + qk[2] * h) + (p0 + pk[2] * h) * xinv2 + chi_p2;
336 qk[3] = (v2 - enu__ + a2 - sq_alpha_half * a2 * enu__ / m2) * (p0 + pk[2] * h) -
337 (q0 + qk[2] * h) * xinv2 + chi_q2;
338 }
339 if (rel == relativity_t::dirac) {
340 /* k0 = F(Y(x), x) */
341 pk[0] = alpha * (2 * rest_energy + enu__ - v0) * q0 - kappa * p0 * xinv0;
342 qk[0] = alpha * (v0 - enu__) * p0 + kappa * q0 * xinv0;
343
344 /* k1 = F(Y(x) + k0 * h/2, x + h/2) */
345 pk[1] = alpha * (2 * rest_energy + enu__ - v1) * (q0 + qk[0] * h_half) -
346 kappa * (p0 + pk[0] * h_half) * xinv1;
347 qk[1] = alpha * (v1 - enu__) * (p0 + pk[0] * h_half) + kappa * (q0 + qk[0] * h_half) * xinv1;
348
349 /* k2 = F(Y(x) + k1 * h/2, x + h/2) */
350 pk[2] = alpha * (2 * rest_energy + enu__ - v1) * (q0 + qk[1] * h_half) -
351 kappa * (p0 + pk[1] * h_half) * xinv1;
352 qk[2] = alpha * (v1 - enu__) * (p0 + pk[1] * h_half) + kappa * (q0 + qk[1] * h_half) * xinv1;
353
354 /* k3 = F(Y(x) + k2 * h, x + h) */
355 pk[3] = alpha * (2 * rest_energy + enu__ - v2) * (q0 + qk[2] * h) - kappa * (p0 + pk[2] * h) * xinv2;
356 qk[3] = alpha * (v2 - enu__) * (p0 + pk[2] * h) + kappa * (q0 + qk[2] * h) * xinv2;
357 }
358 /* Y(x + h) = Y(x) + h * (k0 + 2 * k1 + 2 * k2 + k3) / 6 */
359 p2 = p0 + (pk[0] + 2 * (pk[1] + pk[2]) + pk[3]) * h / 6.0;
360 q2 = q0 + (qk[0] + 2 * (qk[1] + qk[2]) + qk[3]) * h / 6.0;
361
362 /* check overflow */
363 // if (std::abs(p2) > 1e10 || std::abs(q2) > 1e10) {
364 if (std::abs(p2) > 1e4) {
365 /* if we didn't expect the overflow and it happened, or it happened before the
366 * classical turning point, it's bad */
367 if (!prevent_overflow || (prevent_overflow && i < idx_ctp)) {
368 std::stringstream s;
369 if (!prevent_overflow) {
370 s << "unexpected overflow ";
371 } else {
372 s << "overflow before the classical turning point ";
373 }
374 s << "for atom type with zn = " << zn_ << ", l = " << l__ << ", enu = " << enu__ << ", ir = " << i
375 << ", idx_ctp: " << idx_ctp;
376 for (int j = 0; j <= i; j++) {
377 p__[j] /= 1e4;
378 q__[j] /= 1e4;
379 }
380 p2 /= 1e4;
381 q2 /= 1e4;
382
383 // WARNING(s);
384 } else { /* if overflow happened after the classical turning point, it's ok */
385 last = i;
386 break;
387 }
388 }
389
390 p__[i + 1] = p2;
391 q__[i + 1] = q2;
392 }
393
394 if (prevent_overflow && last) {
395 /* find the minimum value of the "tail" */
396 double pmax = std::abs(p__[last]);
397 /* go towards the origin */
398 for (int j = last; j >= 0; j--) {
399 if (std::abs(p__[j]) < pmax) {
400 pmax = std::abs(p__[j]);
401 } else {
402 /* we may go through zero here and miss one node,
403 * so stay on the safe side with one extra point */
404 last = j + 1;
405 break;
406 }
407 }
408 /* zero the tail */
409 for (int j = last; j < nr; j++) {
410 p__[j] = 0;
411 q__[j] = 0;
412 }
413 }
414
415 /* get number of nodes */
416 int nn{0};
417 for (int i = 0; i < nr - 1; i++) {
418 if (p__[i] * p__[i + 1] < 0.0) {
419 nn++;
420 }
421 }
422
423 for (int i = 0; i < nr; i++) {
424 if (rel == relativity_t::none || rel == relativity_t::koelling_harmon || rel == relativity_t::zora) {
425 double V = ve_(i) - zn_ * radial_grid_.x_inv(i);
426 double M = rel_mass(enu__, V);
427 double v1 = ll_half / M / std::pow(radial_grid_[i], 2);
428
429 /* P' = 2MQ + \frac{P}{r} */
430 dpdr__[i] = 2 * M * q__[i] + p__[i] * radial_grid_.x_inv(i);
431
432 /* Q' = (V - E + \frac{\ell(\ell + 1)}{2 M r^2}) P - \frac{Q}{r} */
433 dqdr__[i] = (V - enu__ + v1) * p__[i] - q__[i] * radial_grid_.x_inv(i) + chi_q__(i);
434 }
435 if (rel == relativity_t::iora) {
436 double V = ve_(i) - zn_ * radial_grid_.x_inv(i);
437 double M = rel_mass(enu__, V);
438 double M0 = 1 - sq_alpha_half * V;
439 double v1 = ll_half / M0 / std::pow(radial_grid_[i], 2);
440 double v2 = sq_alpha_half * ll_half * enu__ / std::pow(radial_grid_[i] * M0, 2);
441
442 /* P' = 2MQ + \frac{P}{r} */
443 dpdr__[i] = 2 * M * q__[i] + p__[i] * radial_grid_.x_inv(i) + chi_p__(i);
444
445 /* Q' = (V - E + \frac{\ell(\ell + 1)}{2 M r^2}) P - \frac{Q}{r} */
446 dqdr__[i] = (V - enu__ + v1 - v2) * p__[i] - q__[i] * radial_grid_.x_inv(i) + chi_q__(i);
447 }
448 if (rel == relativity_t::dirac) {
449 double V = ve_(i) - zn_ * radial_grid_.x_inv(i);
450 /* P' = ... */
451 dpdr__[i] = alpha * (2 * rest_energy + enu__ - V) * q__[i] - kappa * p__[i] * radial_grid_.x_inv(i);
452 /* Q' = ... */
453 dqdr__[i] = alpha * (V - enu__) * p__[i] + kappa * q__[i] * radial_grid_.x_inv(i);
454 }
455 }
456
457 return nn;
458 }
459
460 //== inline double extrapolate_to_zero(int istep, double y, double* x, double* work) const
461 //== {
462 //== double dy = y;
463 //== double result = y;
464 //==
465 //== if (istep == 0)
466 //== {
467 //== work[0] = y;
468 //== }
469 //== else
470 //== {
471 //== double c = y;
472 //== for (int k = 1; k < istep; k++)
473 //== {
474 //== double delta = 1.0 / (x[istep - k] - x[istep]);
475 //== double f1 = x[istep] * delta;
476 //== double f2 = x[istep - k] * delta;
477 //==
478 //== double q = work[k];
479 //== work[k] = dy;
480 //== delta = c - q;
481 //== dy = f1 * delta;
482 //== c = f2 * delta;
483 //== result += dy;
484 //== }
485 //== }
486 //== work[istep] = dy;
487 //== return result;
488 //== }
489
490 //== /// Integrate system of two first-order differential equations forward starting from the origin.
491 //== /** Use Bulirsch-Stoer technique with Gragg's modified midpoint method */
492 //== template <bool check_overflow>
493 //== int integrate_forward_gbs(double enu__,
494 //== Spline<double> const& ve__,
495 //== Spline<double> const& mp__,
496 //== std::vector<double>& p__,
497 //== std::vector<double>& dpdr__,
498 //== std::vector<double>& q__,
499 //== std::vector<double>& dqdr__) const
500 //== {
501 //== /* number of mesh points */
502 //== int nr = num_points();
503 //==
504 //== double alpha2 = 0.5 * std::pow((1 / speed_of_light), 2);
505 //== if (!relativistic_) alpha2 = 0.0;
506
507 //== double enu0 = 0.0;
508 //== if (relativistic_) enu0 = enu__;
509
510 //== double ll2 = 0.5 * l_ * (l_ + 1);
511
512 //== double x2 = radial_grid_[0];
513 //== double v2 = ve__[0] - zn_ / x2;
514 //== double M2 = 1 - (v2 - enu0) * alpha2;
515
516 //== p__[0] = std::pow(radial_grid_[0], l_ + 1);
517 //== if (l_ == 0)
518 //== {
519 //== q__[0] = -zn_ * radial_grid_[0] / M2 / 2;
520 //== }
521 //== else
522 //== {
523 //== q__[0] = std::pow(radial_grid_[0], l_) * l_ / M2 / 2;
524 //== }
525 //==
526 //== double step_size2[20];
527 //== double work_p[20];
528 //== double work_q[20];
529
530 //== int last = 0;
531
532 //== for (int ir = 0; ir < nr - 1; ir++)
533 //== {
534 //== double H = radial_grid_.dx(ir);
535 //== double p_est, q_est, p_old, q_old;
536 //== for (int j = 0; j < 12; j++)
537 //== {
538 //== int num_steps = 2 * (j + 1);
539 //== double h = H / num_steps;
540 //== double h2 = h + h;
541 //== step_size2[j] = std::pow(h, 2);
542 //== double x0 = radial_grid_[ir];
543 //== double x0inv = radial_grid_.x_inv(ir);
544 //== double x1inv = radial_grid_.x_inv(ir + 1);
545
546 //== p_old = p__[ir + 1];
547 //== q_old = q__[ir + 1];
548
549 //== double p0 = p__[ir];
550 //== double q0 = q__[ir];
551 //== double p1 = p0 + h * (2 * q0 + p0 * x0inv);
552 //== double q1 = q0 + h * ((ve__[ir] + x0inv * (ll2 * x0inv - zn_) - enu__) * p0 - q0 * x0inv -
553 // mp__[ir]);
554
555 //== for (int step = 1; step < num_steps; step++)
556 //== {
557 //== double x = x0 + h * step;
558 //== double xinv = 1.0 / x;
559 //== double p2 = p0 + h2 * (2 * q1 + p1 * xinv);
560 //== double q2 = q0 + h2 * ((ve__(ir, h * step) + xinv * (ll2 * xinv - zn_) - enu__) * p1 -
561 //== q1 * xinv - mp__(ir, h * step));
562 //== p0 = p1;
563 //== p1 = p2;
564 //== q0 = q1;
565 //== q1 = q2;
566 //== }
567 //== p_est = 0.5 * (p0 + p1 + h * (2 * q1 + p1 * x1inv));
568 //== q_est = 0.5 * (q0 + q1 + h * ((ve__[ir + 1] + x1inv * (ll2 * x1inv - zn_) - enu__) * p1 -
569 //== q1 * x1inv - mp__[ir + 1]));
570
571 //== p__[ir + 1] = extrapolate_to_zero(j, p_est, step_size2, work_p);
572 //== q__[ir + 1] = extrapolate_to_zero(j, q_est, step_size2, work_q);
573 //==
574 //== if (j > 1)
575 //== {
576 //== if (std::abs(p__[ir + 1] - p_old) < 1e-12 && std::abs(q__[ir + 1] - q_old) < 1e-12) break;
577 //== }
578 //== }
579
580 //== /* don't allow overflow */
581 //== if (check_overflow && std::abs(p__[ir + 1]) > 1e10)
582 //== {
583 //== last = ir;
584 //== break;
585 //== }
586
587 //== if (!check_overflow && std::abs(p__[ir + 1]) > 1e10)
588 //== {
589 //== TERMINATE("overflow");
590 //== }
591 //== }
592
593 //== if (check_overflow && last)
594 //== {
595 //== /* find the minimum value of the "tail" */
596 //== double pmax = std::abs(p__[last]);
597 //== for (int j = last - 1; j >= 0; j++)
598 //== {
599 //== if (std::abs(p__[j]) < pmax)
600 //== {
601 //== pmax = std::abs(p__[j]);
602 //== }
603 //== else
604 //== {
605 //== /* we may go through zero here and miss one node,
606 //== * so stay on the safe side with one extra point */
607 //== last = j + 1;
608 //== break;
609 //== }
610 //== }
611 //== for (int j = last; j < nr; j++)
612 //== {
613 //== p__[j] = 0;
614 //== q__[j] = 0;
615 //== }
616 //== }
617 //==
618 //== /* get number of nodes */
619 //== int nn = 0;
620 //== for (int i = 0; i < nr - 1; i++) if (p__[i] * p__[i + 1] < 0.0) nn++;
621
622 //== for (int i = 0; i < nr; i++)
623 //== {
624 //== double V = ve__[i] - zn_ * radial_grid_.x_inv(i);
625 //== double M = 1.0 - (V - enu0) * alpha2;
626
627 //== /* P' = 2MQ + \frac{P}{r} */
628 //== dpdr__[i] = 2 * M * q__[i] + p__[i] * radial_grid_.x_inv(i);
629
630 //== /* Q' = (V - E + \frac{\ell(\ell + 1)}{2 M r^2}) P - \frac{Q}{r} */
631 //== dqdr__[i] = (V - enu__ + double(l_ * (l_ + 1)) / (2 * M * std::pow(radial_grid_[i], 2))) * p__[i] -
632 //== q__[i] * radial_grid_.x_inv(i) - mp__[i];
633 //== }
634
635 //== return nn;
636 //== }
637
638 public:
639 Radial_solver(int zn__, std::vector<double> const& v__, Radial_grid<double> const& radial_grid__)
640 : zn_(zn__)
641 , radial_grid_(radial_grid__)
642 {
643 ve_ = Spline<double>(radial_grid__);
644
645 for (int i = 0; i < num_points(); i++) {
646 ve_(i) = v__[i] + zn_ * radial_grid_.x_inv(i);
647 }
648 ve_.interpolate();
649 }
650
651 std::tuple<int, std::vector<double>, std::vector<double>, std::vector<double>, std::vector<double>>
652 solve(relativity_t rel__, int dme__, int l__, int k__, double enu__) const
653 {
654 int nr = num_points();
655 std::vector<std::vector<double>> p;
656 std::vector<std::vector<double>> q;
657 std::vector<std::vector<double>> dpdr;
658 std::vector<std::vector<double>> dqdr;
659
660 Spline<double> chi_p(radial_grid_);
661 Spline<double> chi_q(radial_grid_);
662
663 int nn{0};
664
665 for (int j = 0; j <= dme__; j++) {
666 p.push_back(std::vector<double>(nr));
667 q.push_back(std::vector<double>(nr));
668 dpdr.push_back(std::vector<double>(nr));
669 dqdr.push_back(std::vector<double>(nr));
670
671 if (j) {
672 if (rel__ == relativity_t::none || rel__ == relativity_t::zora) {
673 for (int i = 0; i < nr; i++) {
674 chi_q(i) = -j * p[j - 1][i];
675 }
676 chi_q.interpolate();
677 } else if (rel__ == relativity_t::koelling_harmon) {
678 double sq_alpha = std::pow(speed_of_light, -2);
679 double ll_half = l__ * (l__ + 1) / 2.0;
680 if (j == 1) {
681 for (int i = 0; i < nr; i++) {
682 double x = radial_grid_[i];
683 double V = ve_(i) - zn_ * radial_grid_.x_inv(i);
684 double M = 1 + 0.5 * sq_alpha * (enu__ - V);
685 chi_p(i) = sq_alpha * q[j - 1][i];
686 chi_q(i) = -p[j - 1][i] * (1 + 0.5 * sq_alpha * ll_half / std::pow(M * x, 2));
687 }
688 } else {
689 throw std::runtime_error("not implemented");
690 }
691 } else if (rel__ == relativity_t::iora) {
692 double sq_alpha = std::pow(speed_of_light, -2);
693 double ll_half = l__ * (l__ + 1) / 2.0;
694
695 if (j == 1) {
696 for (int i = 0; i < nr; i++) {
697 double V = ve_(i) - zn_ * radial_grid_.x_inv(i);
698 double M = 1 - 0.5 * sq_alpha * V;
699 double x = radial_grid_[i];
700 chi_p(i) = q[j - 1][i] * sq_alpha / std::pow(1 - sq_alpha * enu__ / 2 / M, 2);
701 chi_q(i) = -p[j - 1][i] * (1 + 0.5 * sq_alpha * ll_half / std::pow(M * x, 2));
702 }
703 } else if (j == 2) {
704 for (int i = 0; i < nr; i++) {
705 double V = ve_(i) - zn_ * radial_grid_.x_inv(i);
706 double M = 1 - 0.5 * sq_alpha * V;
707 double x = radial_grid_[i];
708 chi_p(i) =
709 q[j - 1][i] * 2 * sq_alpha / std::pow(1 - sq_alpha * enu__ / 2 / M, 2) +
710 q[j - 2][i] * std::pow(sq_alpha, 2) / 2 / M / std::pow(1 - sq_alpha * enu__ / 2 / M, 3);
711 chi_q(i) = -p[j - 1][i] * 2 * (1 + 0.5 * sq_alpha * ll_half / std::pow(M * x, 2));
712 }
713 } else {
714 throw std::runtime_error("not implemented");
715 }
716 chi_p.interpolate();
717 chi_q.interpolate();
718 } else {
719 throw std::runtime_error("not implemented");
720 }
721 }
722
723 switch (rel__) {
724 case relativity_t::none: {
725 nn = integrate_forward_rk4<relativity_t::none, false>(enu__, l__, 0, chi_p, chi_q, p[j], dpdr[j],
726 q[j], dqdr[j]);
727 break;
728 }
729 case relativity_t::koelling_harmon: {
730 nn = integrate_forward_rk4<relativity_t::koelling_harmon, false>(enu__, l__, 0, chi_p, chi_q, p[j],
731 dpdr[j], q[j], dqdr[j]);
732 break;
733 }
734 case relativity_t::zora: {
735 nn = integrate_forward_rk4<relativity_t::zora, false>(enu__, l__, 0, chi_p, chi_q, p[j], dpdr[j],
736 q[j], dqdr[j]);
737 break;
738 }
739 case relativity_t::iora: {
740 nn = integrate_forward_rk4<relativity_t::iora, false>(enu__, l__, 0, chi_p, chi_q, p[j], dpdr[j],
741 q[j], dqdr[j]);
742 break;
743 }
744 default: {
745 throw std::runtime_error("not implemented");
746 }
747 }
748 }
749
750 return std::make_tuple(nn, p.back(), dpdr.back(), q.back(), dqdr.back());
751 }
752
753 /// Integrates the radial equation for a given energy and finds the m-th energy derivative of the radial solution.
754 /** \param [in] rel Type of relativity
755 * \param [in] dme Order of energy derivative.
756 * \param [in] l Oribtal quantum number.
757 * \param [in] enu Integration energy.
758 * \param [out] p \f$ p(r) = ru(r) \f$ radial function.
759 * \param [out] rdudr \f$ ru'(r) \f$.
760 * \param [out] uderiv \f$ u'(R) \f$ and \f$ u''(R) \f$.
761 *
762 * Returns \f$ p(r) = ru(r) \f$, \f$ ru'(r)\f$, \f$ u'(R) \f$ and \f$ u''(R) \f$.
763 *
764 * Surface derivatives:
765 *
766 * From
767 * \f[
768 * u(r) = \frac{p(r)}{r}
769 * \f]
770 * we have
771 * \f[
772 * u'(r) = \frac{p'(r)}{r} - \frac{p(r)}{r^2} = \frac{1}{r}\big(p'(r) - \frac{p(r)}{r}\big)
773 * \f]
774 *
775 * \f[
776 * u''(r) = \frac{p''(r)}{r} - \frac{p'(r)}{r^2} - \frac{p'(r)}{r^2} + 2 \frac{p(r)}{r^3} =
777 * \frac{1}{r}\big(p''(r) - 2 \frac{p'(r)}{r} + 2 \frac{p(r)}{r^2}\big)
778 * \f]
779 */
780 int solve(relativity_t rel__, int dme__, int l__, double enu__, std::vector<double>& p__,
781 std::vector<double>& rdudr__, std::array<double, 2>& uderiv__) const
782 {
783 auto result = solve(rel__, dme__, l__, 0, enu__);
784 int nr = num_points();
785
786 auto p0 = std::get<1>(result);
787 auto p1 = std::get<2>(result);
788 auto q0 = std::get<3>(result);
789 auto q1 = std::get<4>(result);
790
791 /* save the results */
792 p__.resize(nr);
793 rdudr__.resize(nr);
794 for (int i = 0; i < radial_grid_.num_points(); i++) {
795 p__[i] = p0[i];
796 rdudr__[i] = p1[i] - p0[i] * radial_grid_.x_inv(i);
797 }
798
799 double R = radial_grid_.last();
800
801 /* 1st radial derivative of u(r) */
802 uderiv__[0] = (p1.back() - p0.back() / R) / R;
803 /* 2nd radial derivative of u(r) */
804 Spline<double> sdpdr(radial_grid_, p1);
805 sdpdr.interpolate();
806 uderiv__[1] = (sdpdr.deriv(1, nr - 1) - 2 * p1.back() / R + 2 * p0.back() / std::pow(R, 2)) / R;
807
808 return std::get<0>(result);
809 }
810
811 inline int num_points() const
812 {
813 return radial_grid_.num_points();
814 }
815
816 inline int zn() const
817 {
818 return zn_;
819 }
820
821 inline double radial_grid(int i__) const
822 {
823 return radial_grid_[i__];
824 }
825
826 inline Radial_grid<double> const& radial_grid() const
827 {
828 return radial_grid_;
829 }
830};
831
832class Bound_state : public Radial_solver
833{
834 private:
835 int n_;
836
837 int l_;
838
839 int k_;
840
841 /// Tolerance of bound state energy.
842 double enu_tolerance_;
843
844 double enu_;
845
846 Spline<double> p_;
847
848 Spline<double> q_;
849
850 Spline<double> u_;
851
852 Spline<double> rdudr_;
853
854 Spline<double> rho_;
855
856 std::vector<double> dpdr_;
857
858 void solve(relativity_t rel__, double enu_start__)
859 {
860 int np = num_points();
861
862 Spline<double> chi_p(radial_grid());
863 Spline<double> chi_q(radial_grid());
864
865 std::vector<double> p(np);
866 std::vector<double> q(np);
867 std::vector<double> dqdr(np);
868 std::vector<double> rdudr(np);
869 dpdr_ = std::vector<double>(np);
870
871 int s{1};
872 int sp;
873 enu_ = enu_start__;
874 double denu{0.1};
875
876 /* search for the bound state */
877 for (int iter = 0; iter < 1000; iter++) {
878 int nn{0};
879
880 switch (rel__) {
881 case relativity_t::none: {
882 nn = integrate_forward_rk4<relativity_t::none, true>(enu_, l_, k_, chi_p, chi_q, p, dpdr_, q, dqdr);
883 break;
884 }
885 case relativity_t::koelling_harmon: {
886 nn = integrate_forward_rk4<relativity_t::koelling_harmon, true>(enu_, l_, k_, chi_p, chi_q, p,
887 dpdr_, q, dqdr);
888 break;
889 }
890 case relativity_t::zora: {
891 nn = integrate_forward_rk4<relativity_t::zora, true>(enu_, l_, k_, chi_p, chi_q, p, dpdr_, q, dqdr);
892 break;
893 }
894 case relativity_t::dirac: {
895 nn = integrate_forward_rk4<relativity_t::dirac, true>(enu_, l_, k_, chi_p, chi_q, p, dpdr_, q,
896 dqdr);
897 break;
898 }
899 default: {
900 throw std::runtime_error("not implemented");
901 }
902 }
903
904 sp = s;
905 s = (nn > (n_ - l_ - 1)) ? -1 : 1;
906 denu = s * std::abs(denu);
907 denu = (s != sp) ? denu * 0.5 : denu * 1.25;
908 enu_ += denu;
909
910 if (std::abs(denu) < enu_tolerance_ && iter > 4) {
911 break;
912 }
913 }
914
915 if (std::abs(denu) >= enu_tolerance_) {
916 std::stringstream s;
917 s << "enu is not converged for n = " << n_ << " and l = " << l_ << std::endl
918 << "enu = " << enu_ << ", denu = " << denu;
919 throw std::runtime_error(s.str());
920 }
921
922 /* compute r * u'(r) */
923 for (int i = 0; i < num_points(); i++) {
924 rdudr[i] = dpdr_[i] - p[i] / radial_grid(i);
925 }
926
927 /* search for the turning point */
928 int idxtp = np - 1;
929 for (int i = 0; i < np; i++) {
930 if (ve_(i) - zn_ * radial_grid_.x_inv(i) > enu_) {
931 idxtp = i;
932 break;
933 }
934 }
935
936 /* zero the tail of the wave-function */
937 double t1 = 1e100;
938 for (int i = idxtp; i < np; i++) {
939 if (std::abs(p[i]) < t1 && p[i - 1] * p[i] > 0) {
940 t1 = std::abs(p[i]);
941 } else {
942 t1 = 0.0;
943 p[i] = 0.0;
944 q[i] = 0.0;
945 rdudr[i] = 0.0;
946 }
947 }
948
949 for (int i = 0; i < np; i++) {
950 p_(i) = p[i];
951 q_(i) = q[i];
952 u_(i) = p[i] * radial_grid_.x_inv(i);
953 rdudr_(i) = rdudr[i];
954 }
955 p_.interpolate();
956 q_.interpolate();
957 u_.interpolate();
958 rdudr_.interpolate();
959
960 /* p is not divided by r, so we integrate with r^0 prefactor */
961 double norm = inner(p_, p_, 0);
962 if (rel__ == relativity_t::dirac) {
963 norm += inner(q_, q_, 0);
964 }
965
966 norm = 1.0 / std::sqrt(norm);
967 p_.scale(norm);
968 q_.scale(norm);
969 u_.scale(norm);
970 rdudr_.scale(norm);
971
972 /* count number of nodes of the function */
973 int nn{0};
974 for (int i = 0; i < np - 1; i++) {
975 if (p_(i) * p_(i + 1) < 0.0) {
976 nn++;
977 }
978 }
979
980 if (nn != (n_ - l_ - 1)) {
981 FILE* fout = fopen("p.dat", "w");
982 for (int ir = 0; ir < np; ir++) {
983 double x = radial_grid(ir);
984 fprintf(fout, "%12.6f %16.8f\n", x, p_(ir));
985 }
986 fclose(fout);
987
988 std::stringstream s;
989 s << "n = " << n_ << std::endl
990 << "l = " << l_ << std::endl
991 << "enu = " << enu_ << std::endl
992 << "wrong number of nodes : " << nn << " instead of " << (n_ - l_ - 1);
993 throw std::runtime_error(s.str());
994 }
995
996 for (int i = 0; i < np - 1; i++) {
997 rho_(i) += std::pow(u_(i), 2);
998 if (rel__ == relativity_t::dirac) {
999 rho_(i) += std::pow(q_(i) * radial_grid_.x_inv(i), 2);
1000 }
1001 }
1002 }
1003
1004 public:
1005 Bound_state(relativity_t rel__, int zn__, int n__, int l__, int k__, Radial_grid<double> const& radial_grid__,
1006 std::vector<double> const& v__, double enu_start__)
1007 : Radial_solver(zn__, v__, radial_grid__)
1008 , n_(n__)
1009 , l_(l__)
1010 , k_(k__)
1011 , enu_tolerance_(1e-12)
1012 , p_(radial_grid__)
1013 , q_(radial_grid__)
1014 , u_(radial_grid__)
1015 , rdudr_(radial_grid__)
1016 , rho_(radial_grid__)
1017 {
1018 solve(rel__, enu_start__);
1019 }
1020
1021 inline double enu() const
1022 {
1023 return enu_;
1024 }
1025
1026 /// Return charge density, corresponding to a radial solution.
1027 Spline<double> const& rho() const
1028 {
1029 return rho_;
1030 }
1031
1032 /// Return radial function.
1033 Spline<double> const& u() const
1034 {
1035 return u_;
1036 }
1037
1038 /// Return radial function multiplied by x.
1039 Spline<double> const& p() const
1040 {
1041 return p_;
1042 }
1043
1044 std::vector<double> const& dpdr() const
1045 {
1046 return dpdr_;
1047 }
1048};
1049
1050class Enu_finder : public Radial_solver
1051{
1052 private:
1053 int n_;
1054
1055 int l_;
1056
1057 double enu_;
1058
1059 double etop_;
1060 double ebot_;
1061
1062 void find_enu(relativity_t rel__, double enu_start__)
1063 {
1064 int np = num_points();
1065
1066 Spline<double> chi_p(radial_grid());
1067 Spline<double> chi_q(radial_grid());
1068
1069 std::vector<double> p(np);
1070 std::vector<double> q(np);
1071 std::vector<double> dpdr(np);
1072 std::vector<double> dqdr(np);
1073
1074 double enu = enu_start__;
1075 double de = 0.001;
1076 bool found = false;
1077 int nndp = 0;
1078
1079 /* We want to find enu such that the wave-function at the muffin-tin boundary is zero
1080 * and the number of nodes inside muffin-tin is equal to n-l-1. This will be the top
1081 * of the band. */
1082 for (int i = 0; i < 1000; i++) {
1083 int nnd{0};
1084
1085 switch (rel__) {
1086 case relativity_t::none: {
1087 nnd = integrate_forward_rk4<relativity_t::none, false>(enu, l_, 0, chi_p, chi_q, p, dpdr, q, dqdr);
1088 break;
1089 }
1090 case relativity_t::koelling_harmon: {
1091 nnd = integrate_forward_rk4<relativity_t::koelling_harmon, false>(enu, l_, 0, chi_p, chi_q, p, dpdr,
1092 q, dqdr);
1093 break;
1094 }
1095 case relativity_t::zora: {
1096 nnd = integrate_forward_rk4<relativity_t::zora, false>(enu, l_, 0, chi_p, chi_q, p, dpdr, q, dqdr);
1097 break;
1098 }
1099 case relativity_t::iora: {
1100 nnd = integrate_forward_rk4<relativity_t::iora, false>(enu, l_, 0, chi_p, chi_q, p, dpdr, q, dqdr);
1101 break;
1102 }
1103 default: {
1104 throw std::runtime_error("not implemented");
1105 }
1106 }
1107 nnd -= (n_ - l_ - 1);
1108
1109 enu = (nnd > 0) ? enu - de : enu + de;
1110
1111 if (i) {
1112 de = (nnd != nndp) ? de * 0.5 : de * 1.25;
1113 }
1114 if (std::abs(de) < 1e-10) {
1115 found = true;
1116 break;
1117 }
1118 nndp = nnd;
1119 }
1120 etop_ = (!found) ? enu_start__ : enu;
1121
1122 auto surface_deriv = [this, &dpdr, &p]() {
1123 if (true) {
1124 /* return p'(R) */
1125 return dpdr.back();
1126 } else {
1127 /* return R*u'(R) */
1128 return dpdr.back() - p.back() / radial_grid_.last();
1129 }
1130 };
1131
1132 double sd = surface_deriv();
1133
1134 /* Now we go down in energy and search for enu such that the wave-function derivative is zero
1135 * at the muffin-tin boundary. This will be the bottom of the band. */
1136 de = 1e-4;
1137 for (int i = 0; i < 100; i++) {
1138 de *= 1.1;
1139 enu -= de;
1140 switch (rel__) {
1141 case relativity_t::none: {
1142 integrate_forward_rk4<relativity_t::none, false>(enu, l_, 0, chi_p, chi_q, p, dpdr, q, dqdr);
1143 break;
1144 }
1145 case relativity_t::koelling_harmon: {
1146 integrate_forward_rk4<relativity_t::koelling_harmon, false>(enu, l_, 0, chi_p, chi_q, p, dpdr, q,
1147 dqdr);
1148 break;
1149 }
1150 case relativity_t::zora: {
1151 integrate_forward_rk4<relativity_t::zora, false>(enu, l_, 0, chi_p, chi_q, p, dpdr, q, dqdr);
1152 break;
1153 }
1154 case relativity_t::iora: {
1155 integrate_forward_rk4<relativity_t::iora, false>(enu, l_, 0, chi_p, chi_q, p, dpdr, q, dqdr);
1156 break;
1157 }
1158 default: {
1159 throw std::runtime_error("not implemented");
1160 }
1161 }
1162 if (surface_deriv() * sd <= 0) {
1163 break;
1164 }
1165 }
1166
1167 /* refine bottom energy */
1168 double e1 = enu;
1169 double e0 = enu + de;
1170
1171 for (int i = 0; i < 100; i++) {
1172 enu = (e1 + e0) / 2.0;
1173 switch (rel__) {
1174 case relativity_t::none: {
1175 integrate_forward_rk4<relativity_t::none, false>(enu, l_, 0, chi_p, chi_q, p, dpdr, q, dqdr);
1176 break;
1177 }
1178 case relativity_t::koelling_harmon: {
1179 integrate_forward_rk4<relativity_t::koelling_harmon, false>(enu, l_, 0, chi_p, chi_q, p, dpdr, q,
1180 dqdr);
1181 break;
1182 }
1183 case relativity_t::zora: {
1184 integrate_forward_rk4<relativity_t::zora, false>(enu, l_, 0, chi_p, chi_q, p, dpdr, q, dqdr);
1185 break;
1186 }
1187 case relativity_t::iora: {
1188 integrate_forward_rk4<relativity_t::iora, false>(enu, l_, 0, chi_p, chi_q, p, dpdr, q, dqdr);
1189 break;
1190 }
1191 default: {
1192 throw std::runtime_error("not implemented");
1193 }
1194 }
1195 /* derivative at the boundary */
1196 if (std::abs(surface_deriv()) < 1e-10) {
1197 break;
1198 }
1199
1200 if (surface_deriv() * sd > 0) {
1201 e0 = enu;
1202 } else {
1203 e1 = enu;
1204 }
1205 }
1206
1207 ebot_ = enu;
1208 /* last check */
1209 int nn{0};
1210 switch (rel__) {
1211 case relativity_t::none: {
1212 nn = integrate_forward_rk4<relativity_t::none, false>(enu, l_, 0, chi_p, chi_q, p, dpdr, q, dqdr);
1213 break;
1214 }
1215 case relativity_t::koelling_harmon: {
1216 nn = integrate_forward_rk4<relativity_t::koelling_harmon, false>(enu, l_, 0, chi_p, chi_q, p, dpdr, q,
1217 dqdr);
1218 break;
1219 }
1220 case relativity_t::zora: {
1221 nn = integrate_forward_rk4<relativity_t::zora, false>(enu, l_, 0, chi_p, chi_q, p, dpdr, q, dqdr);
1222 break;
1223 }
1224 case relativity_t::iora: {
1225 nn = integrate_forward_rk4<relativity_t::iora, false>(enu, l_, 0, chi_p, chi_q, p, dpdr, q, dqdr);
1226 break;
1227 }
1228 default: {
1229 throw std::runtime_error("not implemented");
1230 }
1231 }
1232
1233 if (nn != n_ - l_ - 1) {
1234 FILE* fout = fopen("p.dat", "w");
1235 for (int ir = 0; ir < np; ir++) {
1236 double x = radial_grid(ir);
1237 fprintf(fout, "%16.8f %16.8f %16.8f\n", x, p[ir], q[ir]);
1238 }
1239 fclose(fout);
1240
1241 // printf("n: %i, l: %i, nn: %i", n_, l_, nn);
1242 std::stringstream s;
1243 s << "wrong number of nodes: " << nn << " instead of " << n_ - l_ - 1 << std::endl
1244 << "n: " << n_ << ", l: " << l_ << std::endl
1245 << "etop: " << etop_ << " ebot: " << ebot_ << std::endl
1246 << "initial surface derivative: " << sd;
1247
1248 throw std::runtime_error(s.str());
1249 }
1250
1251 enu_ = (ebot_ + etop_) / 2.0;
1252 }
1253
1254 public:
1255 /// Constructor
1256 Enu_finder(relativity_t rel__, int zn__, int n__, int l__, Radial_grid<double> const& radial_grid__,
1257 std::vector<double> const& v__, double enu_start__)
1258 : Radial_solver(zn__, v__, radial_grid__)
1259 , n_(n__)
1260 , l_(l__)
1261 {
1262 assert(l_ < n_);
1263 find_enu(rel__, enu_start__);
1264 }
1265
1266 inline double enu() const
1267 {
1268 return enu_;
1269 }
1270
1271 inline double ebot() const
1272 {
1273 return ebot_;
1274 }
1275
1276 inline double etop() const
1277 {
1278 return etop_;
1279 }
1280};
1281
1282}; // namespace sirius
1283
1284#endif // __RADIAL_SOLVER_HPP__
sirius::Bound_state
Definition: radial_solver.hpp:833
sirius::Bound_state::u
Spline< double > const & u() const
Return radial function.
Definition: radial_solver.hpp:1033
sirius::Bound_state::p
Spline< double > const & p() const
Return radial function multiplied by x.
Definition: radial_solver.hpp:1039
sirius::Bound_state::rho
Spline< double > const & rho() const
Return charge density, corresponding to a radial solution.
Definition: radial_solver.hpp:1027
sirius::Bound_state::enu_tolerance_
double enu_tolerance_
Tolerance of bound state energy.
Definition: radial_solver.hpp:842
sirius::Enu_finder
Definition: radial_solver.hpp:1051
sirius::Enu_finder::Enu_finder
Enu_finder(relativity_t rel__, int zn__, int n__, int l__, Radial_grid< double > const &radial_grid__, std::vector< double > const &v__, double enu_start__)
Constructor.
Definition: radial_solver.hpp:1256
sirius::Radial_grid< double >
sirius::Radial_grid::last
T last() const
Last point of the grid.
Definition: radial_grid.hpp:133
sirius::Radial_grid::dx
T dx(const int i) const
Return .
Definition: radial_grid.hpp:152
sirius::Radial_grid::x_inv
T x_inv(const int i) const
Return .
Definition: radial_grid.hpp:158
sirius::Radial_grid::num_points
int num_points() const
Number of grid points.
Definition: radial_grid.hpp:121
sirius::Radial_solver
Finds a solution to radial Schrodinger, Koelling-Harmon or Dirac equation.
Definition: radial_solver.hpp:136
sirius::Radial_solver::radial_grid_
Radial_grid< double > const & radial_grid_
Radial grid.
Definition: radial_solver.hpp:142
sirius::Radial_solver::zn_
int zn_
Positive charge of the nucleus.
Definition: radial_solver.hpp:139
sirius::Radial_solver::solve
int solve(relativity_t rel__, int dme__, int l__, double enu__, std::vector< double > &p__, std::vector< double > &rdudr__, std::array< double, 2 > &uderiv__) const
Integrates the radial equation for a given energy and finds the m-th energy derivative of the radial ...
Definition: radial_solver.hpp:780
sirius::Radial_solver::ve_
Spline< double > ve_
Electronic part of potential.
Definition: radial_solver.hpp:145
sirius::Radial_solver::integrate_forward_rk4
int integrate_forward_rk4(double enu__, int l__, int k__, Spline< double > const &chi_p__, Spline< double > const &chi_q__, std::vector< double > &p__, std::vector< double > &dpdr__, std::vector< double > &q__, std::vector< double > &dqdr__) const
Integrate system of two first-order differential equations forward starting from the origin.
Definition: radial_solver.hpp:150
sirius::Spline::interpolate
Spline< T, U > & interpolate()
Definition: spline.hpp:275
constants.hpp
Various constants.
sirius::wf::inner
std::enable_if_t< std::is_same< T, real_type< F > >::value, void > inner(::spla::Context &spla_ctx__, sddk::memory_t mem__, spin_range spins__, W const &wf_i__, band_range br_i__, Wave_functions< T > const &wf_j__, band_range br_j__, la::dmatrix< F > &result__, int irow0__, int jcol0__)
Compute inner product between the two sets of wave-functions.
Definition: wave_functions.hpp:1622
sirius
Namespace of the SIRIUS library.
Definition: sirius.f90:5
sirius::relativity_t
relativity_t
Type of relativity treatment in the case of LAPW.
Definition: typedefs.hpp:95
sirius::speed_of_light
const double speed_of_light
NIST value for the inverse fine structure (http://physics.nist.gov/cuu/Constants/index....
Definition: constants.hpp:33
spline.hpp
Contains definition and partial implementation of sirius::Spline class.
typedefs.hpp
Contains typedefs, enums and simple descriptors.
Generated on Wed Nov 22 2023 17:00:16 for SIRIUS by doxygen 1.9.3