Now you can combine all studied modules in Tasks 1-10 to organize your DFT program. Follow the following scheme and suggestions.
- Set $r_0=10^{-5}$, $r_f=10$, $N=1000$, and calculate mesh $r_i, i =1,...,N$.
- Calculate analytical solution $\varepsilon^0_{10}$, $P^0_{10}\left(r\right)$ (see Eq. (6) from Task 1), and use it as a first approximation to the solution of the Kohn-Sham (KS) equation $P_1\left(r\right)$, $\varepsilon_1$.
- Calculate nuclear potential $V^{NUC}$ (Eq. (2) in Task 2) and use it as initial approximation to KS potential.
- Solve KS equation using Diff procedure from Task 8 to get the radial wave function $P_2\left(r\right)$.
- Calculate norm of the radial wave function $\lVert P_2\left(r\right) \lVert = \left( 4\pi \int\limits_0^{\infty} \left(P_2\left(r\right)\right)^2dr \right)^{1/2}$ and renormalize the solution $P_2\left(r\right)$.
- Using procedure Vcoul from Task 7 calculate e-e potential ${V^{ee}}$ corresponding to $P_2\left(r\right)$. Using procedure from Task 3 calculate XC potential corresponding to $P_2\left(r\right)$. Calculate KS potential $V_2$ as a sum of ${V^{NUC}}$, ${V^{ee}}\left(P_2\right)$, and ${V^{XC}}\left(P_2\right)$.
- Set the KS potential for next iteration as a mix of potentials $V_2$ and $V_1$ using $\alpha = 0.1$: ${\widetilde{V}}_2={\alpha}V_2+\left(1-\alpha\right)V_1$.
- Using procedure Energy from Task 9 calculate KS energy $\varepsilon_2$ corresponding to $P_2\left(r\right)$.
- Check the difference $\left|{\varepsilon_2-\varepsilon_1}\right|$, if this difference is grater than $\delta=0.0001$, then print energy $\varepsilon_2$, norm of the wave function for current iteration, and set $\varepsilon_2$, $P_2\left(r\right)$, and ${\widetilde{V}}_2$ as initial values for next iteration, go to step 4, and make additional SCF iteration to get new values of $\varepsilon_3$, $P_3\left(r\right)$, and ${\widetilde{V}}_3$, and so on until you reach self-consistent convergence. Set maximum number of iterations to $N^{it}_{max}=300$. Check number of iterations, if number of iterations reached $N^{it}_{max}$, then print error message and stop the program. If the SCF error less than $\delta$, then you reached SCF convergence, go to the next step.
- Calculate the total electronic energy $E$ of the system (see Task 10).
- Print solution $P\left(r\right)$, KS energy $\varepsilon$, total electronic energy $E$, e-e potential, nuclear potential, XC potential, and total KS potential.
- Plot the radial wave function $P\left(r\right)$ and compare with analytical solution $P^0_{10}\left(r\right)$. Plot e-e potential, nuclear potential, XC potential, and total KS potential.
- Change $\alpha $ to 0.2, 0.3,... and repeat the whole procedure from the beginning to see how SCF cycle is working.
Comments.
First, second, and third points are already organized in the code in the Short course. You can use it as a draft for your code.
How to debug your program
Set $V^{ee}=0$ and $V^{KS}=0$. In this case, the KS equation will correspond to the Schrödinger equation for H-like atom. Make several iterations of SCF loop, and check that the solution $P(r)$ corresponds to the solution $P^0_{10}(r)$.