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.


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:

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_] := 
   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
  • $\begingroup$ 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. $\endgroup$
    – joriki
    Commented Mar 4, 2023 at 22:05
  • $\begingroup$ @Roman, thanks for the reply! Could you please comment on my questions in the comment section below your answer ) $\endgroup$
    – Mam Mam
    Commented Apr 18, 2023 at 22:11
  • $\begingroup$ @Roman, could you please look at my questions below $\endgroup$
    – Mam Mam
    Commented Apr 20, 2023 at 9:55

1 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

enter image description here

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}$.

  • 1
    $\begingroup$ 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$. $\endgroup$
    – Roman
    Commented Mar 4, 2023 at 12:18
  • 1
    $\begingroup$ There is no orthogonalization when using a generalized eigenvalue problem. I'll add a derivation. $\endgroup$
    – Roman
    Commented Mar 4, 2023 at 12:18
  • 1
    $\begingroup$ 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. $\endgroup$
    – Roman
    Commented Mar 4, 2023 at 12:37
  • 1
    $\begingroup$ @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. $\endgroup$
    – Roman
    Commented Apr 18, 2023 at 19:55
  • 1
    $\begingroup$ @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. $\endgroup$
    – Roman
    Commented Apr 19, 2023 at 6:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.