How to use Slater Type Orbitals as a basis functions in matrix method correctly?

Question

This question is a continuation of my previous series of questions about basis functions.

I would like to find the minimum energy of Coulomb potential motion using matrix method.

$H=-\frac{1}{2}\Delta-\frac{1}{r}$

I have chosen Slater Type Orbitals as a basis functions $R(r)=Nr^{n-1}e^{-r}$ ,$ \quad n = 1, 2, 3, 4, 5...$ (formula 11.2.2 from here https://chem.libretexts.org/Courses/Pacific_Union_College/Quantum_Chemistry/11%3A_Computational_Quantum_Chemistry/11.02%3A_Gaussian_Basis_Sets )

The answer is known from the general theory ($E_{min} = -0.5$)

My Wolfram Mathematica code give the next values:

$E_{min} = -0.874888131489401$ - 5 basis functions
$E_{min} = -1.7871262624523565$ - 10 basis functions
$E_{min} = -2.407853444926228$ - 15 basis functions
$E_{min} = -2.858626490878647$ - 20 basis functions

The increase in basis functions number gives increase in the absolute value of the energy and it is so far from -0.5. But with increase of the basis functions number, the energy value should get closer to the correct value. What have I done wrong?

The code:

ClearAll["Global`*"]
nmax = 5;
Psi1[r_, n_] := r^(n - 1) Exp[-1*r^2];
NN[n_] := 
  1/Sqrt[Integrate[Psi1[r, n]*Psi1[r, n]*r^2, {r, 0, Infinity}]];
Psi[r_, n_] := NN[n]*r^(n - 1) Exp[-1*r^2];
(*kinetic energy*)
Kk[r_, n1_, n2_] := 
  FullSimplify[
   Psi[r, n2]*
    Laplacian[Psi[r, n1], {r, \[Theta], \[Phi]}, "Spherical"]];
Kx[n1_, n2_] := -0.5*Integrate[Kk[r, n1, n2]*r^2, {r, 0, \[Infinity]}];
KK = Table[Kx[n1, n2], {n1, 1, nmax}, {n2, 1, nmax}];

(*potential energy*)
VP1[r_] := -1/r;
Px1[n1_, n2_] := 
  Integrate[Psi[r, n2]*VP1[r]*Psi[r, n1]*r^2, {r, 0, \[Infinity]}];
PP = Table[Px1[n1, n2], {n1, 1, nmax}, {n2, 1, nmax}];
Coulomb = Min[Eigenvalues[KK + PP]]

Out[2009]= -0.874888

Answer 1

Your orbitals are not orthogonal to each other: here I've normalized them with an explicit formula (using FindSequenceFunction to discover it):

Psi[r_, n_] = 2^(n + 3/4)/Sqrt[(2 n - 1)!! Sqrt[π]] * r^(n - 1) E^(-r^2);

Table[Integrate[Psi[r, n1] Psi[r, n2] r^2, {r, 0, ∞}], {n1, 0, 5}, {n2, 0, 5}] // MatrixForm

Therefore, you need to solve a generalized eigenvalue problem. Alternatively, you can ortho-normalize your basis vectors with the Gram–Schmidt process.

Corrected code (generalized eigenvalue problem):

nmax = 10;
Psi[r_, n_] = 2^(n + 3/4)/Sqrt[(2 n - 1)!! Sqrt[π]] * r^(n - 1) E^(-r^2);

(*kinetic energy*)
Kk[r_, n1_, n2_] := FullSimplify[Psi[r, n2]*
  Laplacian[Psi[r, n1], {r, \[Theta], \[Phi]}, "Spherical"]];
Kx[n1_, n2_] := (-1/2) Integrate[Kk[r, n1, n2]*r^2, {r, 0, ∞}];
KK = Table[Kx[n1, n2], {n1, 1, nmax}, {n2, 1, nmax}];

(*potential energy*)
VP1[r_] = -1/r;
Px1[n1_, n2_] := Integrate[Psi[r, n2]*VP1[r]*Psi[r, n1]*r^2, {r, 0, ∞}];
PP = Table[Px1[n1, n2], {n1, 1, nmax}, {n2, 1, nmax}];

(*basis function overlaps*)
B = Table[Integrate[Psi[r, n2]*Psi[r, n1]*r^2, {r, 0, ∞}],
          {n1, 1, nmax}, {n2, 1, nmax}];

Min@Eigenvalues[{KK + PP, B} // N]
(*    -0.480814    *)

derivation of the generalized eigenvalue problem

The time-independent Schrödinger equation in operator form is

$$ \hat{H}\left|\psi\right>=E\left|\psi\right> $$

Multiplying from the left with an orbital $\left<i\right|$ gives

$$ \left<i\left|\hat{H}\right|\psi\right>=E\left<i|\psi\right> $$

Now we write $\left|\psi\right>=\sum_jc_j\left|j\right>$, which is possible even if the orbitals $\left|j\right>$ are not ortho-normalized. Inserting this expression we get

$$ \sum_j c_j \left<i\left|\hat{H}\right|j\right>=\sum_{j'}c_{j'}E\left<i|j'\right> $$

The "usual" assumption of $\left<i|j\right>=\delta_{ij}$ is not satisfied here (we don't have ortho-normalized basis functions).

Let's call $H$ the matrix with elements $H_{ij} = \left<i\left|\hat{H}\right|j\right>$, and $B$ the matrix with elements $B_{ij}=\left<i|j\right>$. The Schrödinger equation becomes

$$ \sum_j H_{ij} c_j = E \sum_{j'} B_{i j'} c_{j'} $$

or, in matrix writing:

$$ H\cdot\vec{c} = E B\cdot\vec{c} $$

This is a generalized eigenvalue problem. For ortho-normalized basis functions, $B$ is the unit matrix and the problem becomes a regular eigenvalue problem $H\cdot\vec{c} = E \vec{c}$.