Kohn-Sham equation solver: eigenvector

Here we consider how to find eigenvector of Kohn-Sham equation, if we know the eigenvalue.

Let us write the radial Kohn-Sham equation for $1s$ state (see Eq. (4) from Task 2) in the form:

  1. $$ - {1 \over 2} {d^2 \over{d{r^2}}} P\left( r \right) - \left( \varepsilon - V^{KS}\left(r\right)\right) P\left(r\right) = 0$$

Eq. (1) can be solved by Numerov method.

We will solve equation (1) on the region $\left[ {{r_0},{r_f}} \right]$, where r0 is a very small number and rf is sufficiently big number. We separate this region evenly by points $r_i, i=1,...,N$. 

The boundary condition is

  1. $$\eqalign{
      &  P\left( r_0 \right)  = P_{10}^0\left( r_0 \right)  \cr 
      &  P\left( r_f \right)  = P_{10}^0\left( r_f \right) \cr} $$

where $P_{10}^0\left( r \right)$ is the analytical solution considered in Task 1 (see Eq. (6)).

Using Thomas method we look for solution in the form ${P_i} = {\alpha _{i + 1}}{P_{i + 1}} + {\beta _{i + 1}}$ (see Thomas method, formula (3)). As a first step, we find coefficients $\alpha_i$ and $\beta_i$ (see Thomas method, formula (4)), and for second step, calculate the solution $P_i$. 

Using Eq. (2) for $P\left( {{r_0}} \right)$ we can define

  1. $$\eqalign{
      & \alpha_2 = 0  \cr 
      & \beta_2 =   P_{10}^0\left( r_0 \right)\cr} $$

Using Eqs. (4) from Thomas method we can calculate $\alpha_{i+1}$, $\beta_{i+1}$ for $i=2,...,N-1$.

For step two of Thomas method we calculate the functions $P_i$. For that, using Eq. (2) we set $P_N = P_{10}^0\left( r_f \right)$ and then using recurrent formula (3) from Thomas method, we can calculate 

  1. $${P_i} = {\alpha _{i + 1}}{P_{i + 1}} + {\beta _{i + 1}}$$

for $i=N-1,...,2$, whereas $P_1$ is defined from Eq. (2).

Download the current page as Word document, and the developed procedure in Fortran, C++, or Python module.