Density Functional Theory and the Local Density Approximation

[JMedley]

Im trying to calculate the ground state energy of Helium using a density functional theory approach combined with the local density approximation. So far I have set up universal functionals and I mainly need help with the actual algorithm the evaluation of the Hartree energy functional.

 

[cgk]

You need to be more specific. I understand neither what you are trying to achieve, how you are trying to achieve it, nor what exactly is the problem.

 

[JMedley]

As to what I'm trying to achieve:
"I'm trying to calculate the ground state energy of Helium"
As to how I plan on doing it:
"density functional theory"
And admittedly I wasn't particularly clear on the problems but in particular finding a reasonably accurate exchange-correlation term. I'd read papers where people refer to Monte-Carlo methods for finding a good correlation energy but haven't been able to find any literature on it

 

[cgk]

As to what I'm trying to achieve:
"I'm trying to calculate the ground state energy of Helium"
As to how I plan on doing it:
"density functional theory"

Great. That totally clarifies it. So let me give you a solution to your problem:
(1) Get an account on a computer with molpro installed.
(2) Make a textfile he.com:
memory,50,m;
geometry={He};
basis=vtz
ks,lda
(3) execute "molpro he.com".
...tada, you've just calculated the LDA ground state energy of He. Sure you don't want to go into somewhat more detail about what you are trying to achive and how?

And admittedly I wasn't particularly clear on the problems but in particular finding a reasonably accurate exchange-correlation term.

Your topic says you want to use LDA. LDA implies a specific exchange-correlation term, you don't need to find it. Just get it from the VWN parameterization; or use the one of Perdew et al.

I'd read papers where people refer to Monte-Carlo methods for finding a good correlation energy but haven't been able to find any literature on it

Correlation energy of what? Monte Carlo methods are used to calculate the properties of the homogeneous and non-homegenous electron gas. This data is used in parameterizations of ab initio density functionals (like PBE, for example). If you intend to apply an existing functional, like LDA, you do not need that because it's already parameterized for you.

 

[JMedley]

Ok, I guess I underestimated this one! Thus far I have an exchange-correlation term of:
$$E_{xc}= -\frac{3}{4}\sqrt[3]{\frac{3n(r)}{\pi}}+\frac{0.44}{r_s+7.8}$$
Where $$n(r)$$ is the density of electrons and [ tex ] r_s [/tex] is the Wigner radius. Is this an acceptable expression? I have then tried to use the functional derivative of this to get the exchange-correlation potential as:
$$\frac{\delta E_{xc}[n(r)]}{\delta n(r)} = \frac{d}{dn}(n(r)E_{xc})$$
This is then fed, along with an external Coulombic potential and a hartree potential, back into solving a single-electron Schrodingers equation to get a new density and so on and so on...
Is that the gist of it? Having calculated this differential I know it turns out very messy so I doubt I'm on track! Thanks for your help

 

[cgk]

Ok, I guess I underestimated this one! Thus far I have an exchange-correlation term of:
$$E_{xc}= -\frac{3}{4}\sqrt[3]{\frac{3n(r)}{\pi}}+\frac{0.44}{r_s+7.8}$$
Where $$n(r)$$ is the density of electrons and [ tex ] r_s [/tex] is the Wigner radius. Is this an acceptable expression?

Your exchange correlation term has the right format, but it does look somewhat too simple to me. I don't think it has the correct asymptotics in the RPA limit (shouldn't there be a ln() in there?).

Concrete expressions for the LDA correlation energy were given in
http://dx.doi.org/10.1103/PhysRevB.45.13244 (Perdew-Wang parameterization)
http://dx.doi.org/10.1139/p80-159 (VWN parameterization)

Your expression may be a viable parameterization for your concrete task. In order to check that you could compare it with the expressions given in these^ papers. There was also a nice explanation (including giving all the leading terms) in Schwabl's textbook "Advanced Quantum Mechanics".

I have then tried to use the functional derivative of this to get the exchange-correlation potential as:
$$\frac{\delta E_{xc}[n(r)]}{\delta n(r)} = \frac{d}{dn}(n(r)E_{xc})$$
This is then fed, along with an external Coulombic potential and a hartree potential, back into solving a single-electron Schrodingers equation to get a new density and so on and so on...
Is that the gist of it? Having calculated this differential I know it turns out very messy so I doubt I'm on track! Thanks for your help

Your track seems right. Note that density functional derivatives tend to be very ugly when done correctly. So if you get long and cumbersome expressions I'd rather interpret that as a good sign :). Anyway: In practice density functional derivatives are usually implemented in an automatic fashion. That is, only the energy expression is explicitly programmed, and then code for evaluating it and its derivatives is generated from a computer algebra system or sythesized via automatic differentiation.

If you just want to see if you got the concept right, you could concentrate on just using Slater/Dirac exchange as xc functional. This expression is very simple and allows for a non-messy implementation, including derivatives. And if you get it right, you can just assume that it would work for other functionals too, in principle.

 

[JMedley]

Cool, So I've got it so that it converges on an energy (having used a simple mixing scheme). Unfortunately its the wrong energy! I've solved Possions equation in order to get the Hartree energy of the Helium atom and used the exchange-correlation functional of:

$$V_{xv}[n(r)] = -\left(\frac{3}{2\pi}\right)^{\frac{2}{3}}\frac{1}{r_s}[1+0.0545r_sln(1+11.4/r_s)]$$

Where $$r_s$$ is the Wigner radius. So far as I know I should be getting a value of about 5.9Ry but its converging on about 12.1Ry. The ODE I'm integrating to get the Hartree energy is:

$$\frac{d^2S}{dr^2}=-4\pi rn(r)$$

Where $$S=rV_H$$ and $$V_H$$ is the Hartree potential. The Coulombic potential is still:

$$V_{ext}=-\frac{4}{r}$$

Do these look OK? I can't really figure out what the problem is. Cheers. Jack

 

[cgk]

The He/LDA energy is supposed to be somewhat like -2.834 Hartree (that's the atomic unit; one Hartree is two Ry), with about -0.975 H coming from the functional.

The most common error people make then writing their first SCF program is to sum up the energy wrong. The energy is *not* the sum of the orbital eigenvalues, but rather it is <(1/2)(h + f)> where h is the core Hamiltonian (kinetic energy + nuclear potential) plus the fock matrix f (<f> is the sum of the eigenvalues, <> being the expectation value)

About the external potential: shouldn't this be -2/r for He?

 

[JMedley]

Well because we've been working in Rydberg units I thought the External potential was:
$$-\frac{2Z}{r}$$
In that case I think we may be summing the wrong energy; at the moment we're inputting an initial density of double the Hydrogen 1s density. We then use this to get a hartree potential and an exchange-corellation potential. These are then input to Schrodingers equation with the Coulombic potential to get a new density and so on and so on...Presently we are just looking for convergence of the eigenenergy the comes from Schrodingers equation each time. Is this wrong?