sirius::Radial_solver Class Reference

Finds a solution to radial Schrodinger, Koelling-Harmon or Dirac equation. More...

#include <radial_solver.hpp>

Inherited by sirius::Bound_state, and sirius::Enu_finder.

Public Member Functions
	Radial_solver (int zn__, std::vector< double > const &v__, Radial_grid< double > const &radial_grid__)
std::tuple< int, std::vector< double >, std::vector< double >, std::vector< double >, std::vector< double > >	solve (relativity_t rel__, int dme__, int l__, int k__, double enu__) const
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 solution. More...
int	num_points () const
int	zn () const
double	radial_grid (int i__) const
Radial_grid< double > const &	radial_grid () const

Protected Member Functions
template<relativity_t rel, bool prevent_overflow>
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. More...

Protected Attributes
int	zn_
	Positive charge of the nucleus. More...
Radial_grid< double > const &	radial_grid_
	Radial grid. More...
Spline< double >	ve_
	Electronic part of potential. More...

Detailed Description

Finds a solution to radial Schrodinger, Koelling-Harmon or Dirac equation.

Non-relativistic radial Schrodinger equation:

$$-\frac{1}{2}p''(r) + V_{eff}(r)p(r) = Ep(r)$$

where $V_{eff}(r) = V(r) + \frac{\ell (\ell+1)}{2r^2}$, $p(r) = u(r)r$ and $u(r)$ is the radial wave-function. Energy derivatives of radial solutions obey the slightly different equation:

$$-\frac{1}{2}\dot{p}''(r) + V_{eff}(r)\dot{p}(r) = E\dot{p}(r) + p(r)$$

$$-\frac{1}{2}\ddot{p}''(r) + V_{eff}(r)\ddot{p}(r) = E\ddot{p}(r) + 2\dot{p}(r)$$

and so on. We can generalize the radial Schrodinger equation like this:

$$-\frac{1}{2}p''(r) + \big(V_{eff}(r) - E\big) p(r) = \chi(r)$$

where now $p(r)$ represents m-th energy derivative of the radial solution and $\chi(r) = m \frac{\partial^{m-1} p(r)} {\partial^{m-1}E}$

Let's now decouple second-order differential equation into a system of two first-order euquations. From $p(r) = u(r)r$ we have

$$p'(r) = u'(r)r + u''(r) = 2q(r) + \frac{p(r)}{r}$$

where we have introduced a new variable $q(r) = \frac{u'(r) r}{2}$. Differentiating $p'(r)$ again we arrive to the following equation for $q'(r)$:

$$p''(r) = 2q'(r) + \frac{p'(r)}{r} - \frac{p(r)}{r^2} = 2q'(r) + \frac{2q(r)}{r} + \frac{p(r)}{r^2} - \frac{p(r)}{r^2}$$

$$q'(r) = \frac{1}{2}p''(r) - \frac{q(r)}{r} = \big(V_{eff}(r) - E\big) p(r) - \frac{q(r)}{r} - \chi(r)$$

Final expression for a linear system of differential equations for m-th energy derivative is:

$$\begin{eqnarray*} p'(r) &=& 2q(r) + \frac{p(r)}{r} \\ q'(r) &=& \big(V_{eff}(r) - E\big) p(r) - \frac{q(r)}{r} - \chi(r) \end{eqnarray*}$$

Scalar-relativistic equations look similar. For m = 0 (no energy derivative) we have:

$$\begin{eqnarray*} p'(r) &=& 2Mq(r) + \frac{p(r)}{r} \\ q'(r) &=& \big(V(r) - E + \frac{\ell(\ell+1)}{2Mr^2}\big) p(r) - \frac{q(r)}{r} \end{eqnarray*}$$

where $M = 1 + \frac{\alpha^2}{2}\big(E-V(r)\big)$. Because M depends on energy, the m-th energy derivatives of the scalar-relativistic solution take a slightly different form. For m=1 we have:

$$\begin{eqnarray*} \dot{p}'(r) &=& 2M\dot{q}(r) + \frac{\dot{p}(r)}{r} + \alpha^2 q(r) \\ \dot{q}'(r) &=& \big(V(r) - E + \frac{\ell(\ell+1)}{2Mr^2}\big) \dot{p}(r) - \frac{\dot{q}(r)}{r} - \big(1 + \frac{\ell(\ell+1)\alpha^2}{4M^2r^2}\big)p(r) \end{eqnarray*}$$

For m=2:

$$\begin{eqnarray*} \ddot{p}'(r) &=& 2M\ddot{q}(r) + \frac{\ddot{p}(r)}{r} + 2 \alpha^2 \dot{q}(r) \\ \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 + \frac{\ell(\ell+1)\alpha^2}{4M^2r^2}\big)\dot{p}(r) + \frac{\ell(\ell+1)\alpha^4}{4M^3r^2} p(r) \end{eqnarray*}$$

Derivation of IORA.

First thing to do is to expand the $1/M$ to the first order in $E$:

$$\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) = \frac{1}{M_0} - \frac{\alpha^2}{2} \frac{1}{M_0^2} E$$

From this expansion we can derive the expression for $M$:

$$M = \Big(\frac{1}{M}\Big)^{-1} = \frac{M_0}{1 - \frac{\alpha^2}{2}\frac{E}{M_0}}$$

Now we insert this expressions into the scalar-relativistic equations:

$$\left\{ \begin{array}{ll} \displaystyle p'(r) = 2 \frac{M_0}{1 - \frac{\alpha^2}{2}\frac{E}{M_0}} q(r) + \frac{p(r)}{r} \\ \displaystyle q'(r) = \Big(V(r) - E + \frac{\ell(\ell+1)}{2r^2} \big(\frac{1}{M_0} - \frac{\alpha^2}{2} \frac{E}{M_0^2} \big) \Big) p(r) - \frac{q(r)}{r} = \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) p(r) - \frac{q(r)}{r} \end{array} \right.$$

The first energy derivatives are then:

$$\left\{ \begin{array}{ll} \displaystyle \dot{p}'(r) = 2 \frac{M_0}{1 - \frac{\alpha^2}{2}\frac{E}{M_0}} \dot{q}(r) + \frac{\dot{p}(r)}{r} + \frac{\alpha^2}{(1 - \frac{\alpha^2 E}{2 M_0})^2} q(r) \\ \displaystyle \dot{q}'(r) = \Big(V(r) - E + \frac{\ell(\ell+1)}{2 M_0 r^2} - \frac{\alpha^2}{2} \frac{\ell(\ell+1)}{2 M_0^2 r^2} E \Big) \dot{p}(r) - \frac{\dot{q}(r)}{r} - \Big(1 + \frac{\alpha^2}{2} \frac{\ell(\ell+1)}{2 M_0^2 r^2} \Big) p(r) \end{array} \right.$$

The second energy derivatives are:

$$\left\{ \begin{array}{ll} \displaystyle \ddot{p}'(r) = 2 \frac{M_0}{1 - \frac{\alpha^2}{2}\frac{E}{M_0}} \ddot{q}(r) + \frac{\ddot{p}(r)}{r} + \frac{2 \alpha^2}{(1 - \frac{\alpha^2 E}{2 M_0})^2} \dot{q}(r) + \frac{\alpha^4}{2 M_0 (1 - \frac{\alpha^2 E}{2 M_0})^3} q(r) \\ \displaystyle \ddot{q}'(r) = \Big(V(r) - E + \frac{\ell(\ell+1)}{2 M_0 r^2} - \frac{\alpha^2}{2} \frac{\ell(\ell+1)}{2 M_0^2 r^2} E \Big) \ddot{p}(r) - \frac{\ddot{q}(r)}{r} - 2 \Big(1 + \frac{\alpha^2}{2} \frac{\ell(\ell+1)}{2 M_0^2 r^2} \Big) \dot{p}(r) \end{array} \right.$$

Definition at line 135 of file radial_solver.hpp.

Constructor & Destructor Documentation

Radial_solver()

sirius::Radial_solver::Radial_solver ( int  zn__, std::vector< double > const &  v__, Radial_grid< double > const &  radial_grid__  )

inline

Definition at line 639 of file radial_solver.hpp.

Member Function Documentation

integrate_forward_rk4()

template<relativity_t rel, bool prevent_overflow>

int sirius::Radial_solver::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

inlineprotected

Integrate system of two first-order differential equations forward starting from the origin.

Use Runge-Kutta 4th order method

Definition at line 150 of file radial_solver.hpp.

solve() [1/2]

std::tuple< int, std::vector< double >, std::vector< double >, std::vector< double >, std::vector< double > > sirius::Radial_solver::solve ( relativity_t  rel__, int  dme__, int  l__, int  k__, double  enu__  ) const

inline

Definition at line 652 of file radial_solver.hpp.

solve() [2/2]

int sirius::Radial_solver::solve ( relativity_t  rel__, int  dme__, int  l__, double  enu__, std::vector< double > &  p__, std::vector< double > &  rdudr__, std::array< double, 2 > &  uderiv__  ) const

inline

Integrates the radial equation for a given energy and finds the m-th energy derivative of the radial solution.

Parameters

[in]	rel	Type of relativity
[in]	dme	Order of energy derivative.
[in]	l	Oribtal quantum number.
[in]	enu	Integration energy.
[out]	p	$p(r) = ru(r)$ radial function.
[out]	rdudr	$ru'(r)$.
[out]	uderiv	$u'(R)$ and $u''(R)$.

Returns $p(r) = ru(r)$, $ru'(r)$, $u'(R)$ and $u''(R)$.

Surface derivatives:

From

$$u(r) = \frac{p(r)}{r}$$

we have

$$u'(r) = \frac{p'(r)}{r} - \frac{p(r)}{r^2} = \frac{1}{r}\big(p'(r) - \frac{p(r)}{r}\big)$$

$$u''(r) = \frac{p''(r)}{r} - \frac{p'(r)}{r^2} - \frac{p'(r)}{r^2} + 2 \frac{p(r)}{r^3} = \frac{1}{r}\big(p''(r) - 2 \frac{p'(r)}{r} + 2 \frac{p(r)}{r^2}\big)$$

Definition at line 780 of file radial_solver.hpp.

num_points()

int sirius::Radial_solver::num_points ( ) const

inline

Definition at line 811 of file radial_solver.hpp.

zn()

int sirius::Radial_solver::zn ( ) const

inline

Definition at line 816 of file radial_solver.hpp.

radial_grid() [1/2]

double sirius::Radial_solver::radial_grid ( int  i__ ) const

inline

Definition at line 821 of file radial_solver.hpp.

radial_grid() [2/2]

Radial_grid< double > const & sirius::Radial_solver::radial_grid ( ) const

inline

Definition at line 826 of file radial_solver.hpp.

Member Data Documentation

zn_

int sirius::Radial_solver::zn_

protected

Positive charge of the nucleus.

Definition at line 139 of file radial_solver.hpp.

radial_grid_

Radial_grid<double> const& sirius::Radial_solver::radial_grid_

protected

Radial grid.

Definition at line 142 of file radial_solver.hpp.

ve_

Spline<double> sirius::Radial_solver::ve_

protected

Electronic part of potential.

Definition at line 145 of file radial_solver.hpp.

---

The documentation for this class was generated from the following file:

• radial_solver.hpp