Our task is to find the electron-electron repulsion potential ${V^{ee}}\left(r\right)$, if we know the radial wave function $P\left( r \right)$. For that we use Eq. (3) from e-e potential, discretize it using finite differences to get the matrix equation (see Numerov method, Eq. (7)), and then solve the matrix equation by Thomas method.
Let us introduce new function
- $$Y\left(r\right) = r{V^{ee}}\left(r\right)$$
Our task is to solve numerically the following equation (see Eq. (3) from e-e potential)
- $${{{d^2}} \over {d{r^2}}}\left( {Y\left( r \right)} \right) = - {{4\pi } \over r}{P^2}\left( r \right)$$
where $P\left( r \right)$ is known radial wave function.
We will solve equation (2) for $r_i, i=1,...,N$, where points $r_i$ equally separate the region $\left[ {{r_0},{r_f}} \right]$, where r0 is a very small number and rf is sufficiently big number.
The boundary condition is (see Eq. (6) from e-e potential)
- $$\eqalign{
& {\left. {Y\left( r \right)} \right|_{r \to 0}} = 0 \cr
& {\left. {Y\left( r \right)} \right|_{r \to + \infty }} = 1 \cr} $$
Since we use the finite region $\left[ {{r_0},{r_f}} \right]$, we can write the boundary condition in the form
- $$\eqalign{
& Y\left( {{r_0}} \right) = Z{r_0} \cr
& Y\left( {{r_f}} \right) = 1 \cr} $$
For $Y\left( {{r_0}} \right)$ see Eqs. (5) and (6) from e-e potential.
Then using Thomas method we will look for solution in the form ${y_i} = {\alpha _{i + 1}}{y_{i + 1}} + {\beta _{i + 1}}$. For that, as a first step, we find coefficients $\alpha_i$ and $\beta_i$ (see formula (4) in Thomas method), and then, for second step, calculate the solution $Y_i$ using formula (3) from Thomas method.
For that using Eq. (4) for $Y\left( {{r_0}} \right)$ we can define
- $$\eqalign{
& \alpha_2 = 0 \cr
& \beta_2 = Z{r_0} \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 $Y_i$. For that, using Eq. (4) we set $Y_N = 1$ and then using recurrent formula (3) from Thomas method, we can calculate
- $${Y_i} = {\alpha _{i + 1}}{Y_{i + 1}} + {\beta _{i + 1}}$$
for $i=N-1,...,2$ (backward sweep), whereas $Y_1$ is defined from Eq. (4).
The electron-electron repulsion potential then can be calculated using Eq. (1)
- $${V^{ee}_i} = Y_i / r_i $$
Download the current page as Word document, and the developed procedure in Fortran, C++, or Python modules.