I want to solve the radial Schrodinger equation for a Carbon atom with all six electrons. The equation I want to solve is for $r>0$: $$[\frac{-\hbar^2}{2m}\frac{\partial^2}{\partial r^2}+ \frac{\hbar^2 l(l+1)}{2mr^2}+ V(r)]u(r) = Eu(r),\tag{1}$$

in which:

  • $u(r) = rR(r)$ the wavefunction being $\Phi_l^m(r, \theta, \phi) = R(r)Y_l^m(\theta, \phi)$.
  • $V(r) = V_{proton-electron}(r) + V_{electron-electron}(r) + V_{exchange,correlation}(r)$

The basis function that I choose to decompose $R$ is a simple STO3G basis of the form $\chi_i(r) = r^i e^{-\alpha_i r^2}$ that only depend on $r$. When computing the matrix elements of $V_{proton-electron}$, it is quite the same procedure compared to solving the Schrodinger equation without separating the variables. The matrix element of the Hartree potential doesnt not seem as straightforward: $$V^{Hartree}_{ij} = \sum_{k,l} (ij|kl)P_{kl},\tag{2}$$ with: $$\tag{3}(ij|kl) = \int_{r_1\in \mathcal{R}^+} \int_{r_2\in \mathcal{R}^+}\chi_i(r_1)\chi_j(r_1)\frac{1}{|r_1-r_2|}\chi_k(r_2)\chi_l(r_2),$$ and $P$ is the density matrix.


  • The basis functions in which $R$ is spanned are spherically symmetric, since the angular dependencies have been separated.
  • The two-electron integral as I wrote does not span the entire 3D space and do not take into account cases where the electron density is not spherical (like p orbitals). In addition to that it is wrong.

What are the matrix elements of the electron repulsion potential $V_{electron-electron}$ in when solving the Radial Schrodinger equation?


2 Answers 2


You have not specified sufficiently precisely what is the problem you want to solve. It sounds like you want to solve the radial density functional equations; however, there is still a world of options on how to do that.

  1. Carbon has an open $p$ shell, and due to symmetry breaking, the Kohn-Sham orbitals actually do not obey spherical symmetry but instead you get $s$-$d$ mixing, $p$-$f$ mixing, etc.
  2. The ground state of carbon is also not spin-restricted, so the orbitals are different for spin-up and spin-down electrons.

One can, however, restrict the problem to being spin unrestricted and the orbitals to not break the symmetry; this is achievable by equal occupation of the orbitals (see e.g. doi:10.1103/PhysRevA.101.012516 and doi:10.1021/acs.jctc.3c00183). This will reduce the problem to determining a single Fock matrix for every angular momentum $l$.

The two-electron integrals do account for non-spherical densities; note that your definition still includes the spherical harmonics in the basis functions. To figure out spatial expressions you need to integrate out the angles employing the closure properties of the spherical harmonics, which gives you Gaunt coupling coefficients, see doi:10.1002/qua.25968. You can calculate the nuclear attraction and Coulomb integrals by similar techniques.

The necessary radial matrix elements turn out to be fairly simple in the Gaussian basis, especially if your density is spherically symmetric. Expressions for the radial integrals have been reported by Pitzer in doi:10.1016/j.cpc.2012.02.009.


Solving the radial Schroedinger equation (your Eq. 1) is done with $V(r)$ pre-determined (either it is an analytic expression of a model, such as the Morse/long-range model, or it is known in some other way, such as a point-wise array of $V$ and $r$ values. The most widely used way to solve the radial Schroedinger equation for a diatomic molecule is to use the LEVEL software which originated in the 1960s, but has been maintained by myself since my supervisor (the original author of LEVEL in the 1960s) passed away. LEVEL is also available in OpenMolcas, and you can let me know if you run into any issues when using it.

What you seem to be trying to do is to calculate $V(r)$, and also to solve the radial Schroedinger equation for the calculated $V(r)$. Calculating $V(r)$ can be done with any electronic structure software, such as , or if you would really like to write your own software for this (for many reasons, I do not recommend doing this) you can do so, but basically it has nothing to do with solving the radial Schroedinger equation and everything to do with solving the electronic Schroedinger equation.

  • $\begingroup$ The thing that bugs me is that $V(r)$ depends only on the distance from the nucleus. It is weird when considering orbitals with angular momentum greater than $0$. $\endgroup$
    – L Maxime
    Jan 26 at 14:09
  • $\begingroup$ The radial Schrodinger equation should only be used when $V$ is symmetrically spherical? $\endgroup$
    – L Maxime
    Jan 26 at 14:18

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .