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
    Mar 4 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
    Apr 18 at 22:11
  • $\begingroup$ @Roman, could you please look at my questions below $\endgroup$
    – Mam Mam
    Apr 20 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
    Mar 4 at 12:18
  • 1
    $\begingroup$ There is no orthogonalization when using a generalized eigenvalue problem. I'll add a derivation. $\endgroup$
    – Roman
    Mar 4 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
    Mar 4 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
    Apr 18 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
    Apr 19 at 6:56

Your Answer

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

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