Why does the minimum energy of motion in a Coulomb field differ from the theoretical value?

Question

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

Actually your code working well for orthogonal polynomials only. For example,

ClearAll["Global`*"]
nmax = 5;
Psi1[r_, n_] := 
  Sqrt[ (n - l - 1)!/(n + l)!] E^(-(r/n)) ((2 r)/n)^l 2/
    n^2 LaguerreL[n - l - 1, 2 l + 1, (2 r)/n] /. l -> 0;

Psi[r_, n_] := Psi1[r, n];
(*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]]

We have out -0.5 as expected. In a case of the Slater type orbitals we need to solve variational problem using representation $\Psi _n=\sum_{m=1}^{nmax}c_{n,m} Psi(r,m)$, for example,

ClearAll["Global`*"]
nmax = 5; var[n_] := Table[c[n, l], {l, nmax}];
Psi[r_, n_] := r^(n - 1) Exp[-r];



(*kinetic energy*)
Kk[r_, n1_, n2_] := 
  FullSimplify[
   Psi[r, n2]*
    Laplacian[Psi[r, n1], {r, \[Theta], \[Phi]}, "Spherical"]];
Kx[n_] := -1/2 Sum[
    c[n, l1] c[n, l2] Integrate[
      Kk[r, l1, l2]*r^2, {r, 0, \[Infinity]}], {l1, nmax}, {l2, nmax}];


(*potential energy*)
VP1[r_] := -1/r;
Px1[n1_, n2_] := 
  Integrate[Psi[r, n2]*VP1[r]*Psi[r, n1]*r^2, {r, 0, \[Infinity]}];
Px[n_] := Sum[c[n, l1] c[n, l2] Px1[l1, l2], {l1, nmax}, {l2, nmax}];



norm[n_] = {Sum[
     c[n, l1] c[n, l2] NIntegrate[
       Psi[r, l1] Psi[r, l2] r^2, {r, 0, \[Infinity]}], {l1, 
      nmax}, {l2, nmax}] == 1};

sol[n_] := NMinimize[{Kx[n] + Px[n], norm[n]}, var[n]];

For sol1=sol[1] we have as expected

{-0.5, {c[1, 1] -> 2., c[1, 2] -> -2.61219*10^-9, 
  c[1, 3] -> 1.96267*10^-9, c[1, 4] -> -5.24162*10^-10, 
  c[1, 5] -> 4.37415*10^-11}} 

We also can check that our state $\Psi_1=\sum c_{1,m}Psi(r,m)$ is normalized

NIntegrate[
 Evaluate[
    Sum[c[1, m] Psi[r, m], {m, 1, nmax}] /. sol1[[2]]]^2 r^2, {r, 0, 
  Infinity}]

Out[]= 1.