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

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

• Cross-posted to math.SE. Please don't do that, especially not without linking the questions to each other and explaining your reason for cross-posting. It wastes everyone's time. Commented Mar 4, 2023 at 22:05
• @Roman, thanks for the reply! Could you please comment on my questions in the comment section below your answer ) Commented Apr 18, 2023 at 22:11
• @Roman, could you please look at my questions below Commented Apr 20, 2023 at 9:55

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 gives

$$\left=E\left$$

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=\sum_{j'}c_{j'}E\left$$

The "usual" assumption of $$\left=\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$$, and $$B$$ the matrix with elements $$B_{ij}=\left$$. 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}$$.

• With 20 orbitals there seems to be a bit of a numerical problem. I've modified the code a bit to work with the normalized (but not orthogonalized) orbitals: it makes the problem a bit easier but there are still complex numbers appearing for $n_{\text{max}}=20$. Commented Mar 4, 2023 at 12:18
• There is no orthogonalization when using a generalized eigenvalue problem. I'll add a derivation. Commented Mar 4, 2023 at 12:18
• For $n_{\text{max}}=20$ we seem indeed to run into numerical difficulties. The matrix condition number of $B$ is larger than $10^{22}$, exceeding the numerical precision of machine-double numbers, and so the condition that $B$ is positive-definite is no longer strictly true at the numerical level. Increasing the working precision (beyond machine-double numbers) helps: Min@Eigenvalues[N[{KK + PP, B}, 100]] still gives a good answer. Commented Mar 4, 2023 at 12:37
• @MamMam please check the given derivation of the general eigenvalue problem: there are no assumptions on the basis functions. They don't need to be normalized, and they don't need to be orthogonal. Commented Apr 18, 2023 at 19:55
• @MamMam Once you have {eval, evec} = Eigensystem[{KK + PP, B} // N] you can normalize the eigenvectors with nevec = Normalize[#, Sqrt[# . B . #] &] & /@ evec. Then you have, for example, f1[r_] = v[[1]] . Table[Psi[r, n], {n, nmax}] as the first eigen-function. Check that it is normalized properly: Integrate[f1[r]^2*r^2, {r, 0, ∞}]` gives 1. Commented Apr 19, 2023 at 6:56