I am trying to find a solution to $$\displaystyle \left[-\frac{1}{2}\nabla^2 - \frac{2}{r} + C(r)\right]\phi(r) = E\phi(r)$$ where $C(r)$ is a known function of r. I am just looking for some help on the strategy to solve explicitly for $\phi(r)$. I have read through the solution for hydrogen using separation of variables but

1) I don't really understand it, they make a lot of handwavy arguments about terms disappearing :)

2) How does one handle the extra $C(r)$?

FYI: this is for the Hartree Fock method on Helium so the actual equation is the one posted in the answer here: Hartree Fock equations

P.S. I posted this in math.stackechange first should I close that one?

  • $\begingroup$ Do you have mathematica? $\endgroup$ – Sponge Bob Nov 28 '15 at 3:02
  • $\begingroup$ Haha no I don't :( but can Mathematica handle it? Also I am trying to implement my own HF algorithm in Python so if I could come up with an analytical solution that would be awesome. $\endgroup$ – user2879934 Nov 28 '15 at 3:04
  • 3
    $\begingroup$ There's no general algorithm that can provide solutions for a generic $C(r)$. If only it were that simple! Wolfram Alpha's 'DSolve' can provide solace for some simple $C(r)$ but will quickly run out of computing steam for more complicated $C(r)$. $\endgroup$ – Gert Nov 28 '15 at 3:07
  • 2
    $\begingroup$ @user2879934: HF-based methods provide numerical solutions, not analytical ones. $\endgroup$ – Gert Nov 28 '15 at 3:10
  • $\begingroup$ @Gert well kind of this is just one step in HF for helium...we have already turned it into a one body problem that should be solvable by integrating out the effect of the second electron and turned it into an effective field. The step in HF that I am stuck on is I guessed my orbitals as Guassian and I calculated C(r) but how the heck to I generate a new psi using the Hamiltonian? Sorry for typos on phone $\endgroup$ – user2879934 Nov 28 '15 at 3:14

To do Hartree-Fock to Helium atom, you just need to calculate one orbital, which for helium is spherically symmetric. The Hartree-Fock 'integro-differential' equation for spherically symmetric atom with one eigenstate can be written as

$$\left(-\frac{1}{2}\nabla^2 - \frac{2}{r} + V_{Hx}(r) \right) u(r) = \epsilon u(r),$$

where $u(r) = \psi(r) r$ and

$$V_{Hx}(r) = \frac{1}{2} \int d{\bf r}' \frac{n({\bf r'})}{|{\bf r-r}'|}$$ The factor $\frac{1}{2}$ comes from exchange term cancelling half of the Hartree-potential.

I would recommend shooting method. (Alternatives: gaussian basis set, finite element, finite difference methods).

Shooting method is simple. Guess an eigenvalue and integrate the solution of radial Schrödinger equation from r=0 to r=r_max. Require that boundary condition is fulfilled at r_max (zero). If not, change the eigenvalue accordingly.


Solving the radial Hartree-Fock potential equation is also super simple.

$$ \int dr' \frac{n(r')}{|r-r'|} = \frac{1}{r} \int_0^r dr' 4\pi n(r) r^2 + \int_r^\infty 4\pi n(r) r $$

(Eq. 7.6 here) https://wiki.fysik.dtu.dk/gpaw/_static/rostgaard_master.pdf

The latter part of my answer here explains how this equation works. Dipole in a spherical cavity in an infinite dielectric

Beware, that the Fock-operator is only multiplicative for the Helium 1s orbitals. Rest of the eigenvalue spectrum of the given equation does not match to Fock-operator spectrum.

edit: I just read a comment of yours, which stated that you are using gaussian orbitals. You can solve it easily. You just need to calculate some matrix elements. Given a set of gaussian orbitals $\phi_i$, you calculate two matrices

$$H_{ij} = < \phi_i | H | \phi_j>$$


$$S_{ij} = < \phi_i | \phi_j >$$.

Then you solve this generalized eigenvalue equation

$$Hc=\epsilon Sc$$, and you get the new coefficients for the wave function as the eigenvector. Your wave function is now

$$ \psi(r) = \sum_n c_n \phi_n(r) $$

The method is called subspace diagonalization.

To solve the generalized eigenvalue equation with octave or Matlab,

octave:1> H = [ -2 0.1 ; 0.1 -3 ];
octave:2> S = [ 1 0.1 ; 0.1 1 ];
octave:5> [psi, lambda] = eig(H,S)
psi =

  -0.35088  -0.94180
   0.97217  -0.25494

lambda =

Diagonal Matrix

  -3.1498        0
        0  -1.9209

Here you obtain the eigenvalues (-3.14 is the lowest), and the corresponding eigenvector [-0.35, 0.97]. Now, all that there is left, is to fill matrices H and S with the real matrix elements.

To be more spesific, following integrals must be evaluated:

$$H_{ij} = \int dr \phi_i(r) (-1/2 \nabla^2 + V(r) ) \phi_j(r) \\ S_{ij} = \int dr \phi_i(r) \phi_j(r) $$

  • $\begingroup$ In $$Hc=\epsilon Sc$$ are you solving for $\epsilon$ to update the coefficient? Also how does one solve that? Thanks so much Mikael! $\endgroup$ – user2879934 Nov 29 '15 at 2:02
  • $\begingroup$ I updated the answer accordingly. One doesn't update epsilon, but obtains directly the correct value. This is called direct diagonalization. For large problems, one uses iterative diagonalization where $\epsilon$ is estimated and converged iteratively. $\endgroup$ – Mikael Kuisma Nov 29 '15 at 2:09
  • $\begingroup$ oh man! I didn't realize math packages had the ability to solve eigien-value equations. That is awesome thank you! $\endgroup$ – user2879934 Nov 29 '15 at 2:15
  • $\begingroup$ Actually one last thing...what "does fill matrices H and S with the real matrix elements" mean? Sorry I am just a computer programmer :) $\endgroup$ – user2879934 Nov 29 '15 at 2:52
  • $\begingroup$ Haven't we already filled them by evaluating the two given integrals? $\endgroup$ – user2879934 Nov 29 '15 at 14:32

I feel like this question is a little too open ended and I'm sorry for that...but I did find a WONDERFUL HF walk through http://www.phys.sinica.edu.tw/TIGP-NANO/Course/2011_Spring/classnotes/CMS_20110511.pdf that clears everything up.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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