Hartree-fock-slater for single atoms

[sam_bell]

I've been looking at the original Herman-Skillman program that solves for the energy levels of a single atom using Hartree-Fock with Slater exchange orbitals. There is one part I can't understand: the expression of the correction to the energy level for a good guess. What the program *does* is: (1) set $$R^{out}$$(E_guess,r) = 0 for r = 0 and then numerically integrates out to r = r_match = classical turning point; (2) set $$R^{in}$$(E_guess,r) = 0 for r = infinity = end of grid, and then numerically integrates into r = r_match; (3) Computes logarithmic derivatives $$\gamma = R'(r_{match})/R(r_{match})$$ on both sides; and (4) sets dE =

ignore -->(( $$ignore$$ $$ignore$$ ))<---

$$\frac{ \gamma_{out}-\gamma_{in} }{ \int_0^{r_{match}} R_{out}^2/R_{out}(r_{match})^2 + \int_{r_{match}}^{\infty} R_{in}^2/R_{in}(r_{match})^2 }$$

I don't know how I would show that works.

Hoping for some helpful suggestions. Thx.

 

[cgk]

I don't quite get it. What is the meaning of r_out/r_in/r_match? What correction are you talking of? What does it numerically integrate? Is that for making the initial guess of the orbitals? (I take they are represented numerically on a grid and multiplied by some spherical harmonic?)

I have to admit I have not heard of that program before. Is that for some historic research or are you trying to do something specific? We might be able to help if you gave more details of what you are trying to accomplish.

 

[sam_bell]

I believe Herman & Skillman are the first to implement Hartree-Fock with Slater exchange (in 1962). I'm a graduate student and some parts of the code we use bare strong resemblance to that of Herman & Skillman. So I wanted to understand the original program.

R_out(E,r) is the solution of Schrodinger Equation with U(r) = V_centrifigul(r) + V_el-ion(r) + V_el-el(r) + V_exch(r) and with boundary condition R_out(E,0) = 0. It's only defined to an arbitrary overall factor. Similarly, R_in(E,r) is a solution with boundary condition R(E,r) = 0 for r --> +infinity. The "numerical integration" is just a way of finding these functions on a radial grid. It could be done by Runge-Kutta or a predictor-corrector method.

To clarify my question: A priori you don't know what E corresponds to a bound state of the atom. For the Li atom you might initially guess the valance el- to have an energy similar to that of the el- in H, or 13.6 eV (it's closer to half this). So you would take E_guess = 13.6 eV and compute R_in(r) and R_out(r). Then, using the expression in my above post you would get a better energy E_guess' and recompute R_in(E_guess',r) and R_out(E_guess',r). You would continue this until E_guess is very close to E_guess' and take E = E_guess'. It would then also happen that R_in(r) = R_out(r) and you could take this as the radial wave function for the outer el-.

My sticking point is I can't understand how/why the expression for the correction works. It seems to me the crux of the whole procedure. I've tried my hand at it, but I didn't get anywhere. Thx.

 

[cgk]

Okay, thanks for the explanations. I'm not sure, but I guess one can understand this as a fix-point iteration for E. Say, we have E and some function f(E) and we want to achieve E = f(E). In this case E would be the current guess energy, and f(E) would be an expression for the energy in terms of the radial wave function computed from the current guess energy (i.e., take Eguess, use that to calculate R(r), then integrate up R(r) to get the actual energy the wave function would have for the current R(r). The so obtained energy would be f(Eguess)).

For schemes like this the typical update would be to calculate the residual res(E) := E - f(E) (i.e., the amount to which the fixed-point equation is violated), and then update E as E -> E - res(E)/(d/dE f), where the latter is an estimate of the derivative at the current E. This can be understood as a standard zero-point search. I have to admit I'm not quite fit in expressions for the energy and derivatives in terms of radial functions, so I'm not sure if this is what the code is actually doing. But I seem to vaguely remember that this concrete logarithmic derivative is actually an algebraic transformation of some other more complex derivative, which I, however, forgot.

Just as a side remark: is there some reason why you don't simply take some quantum chemistry program working with Gaussian basis sets? Solving either the full Hartree-Fock equations or some DFT analog (e.g., the X-Alpha or the Slater-exchange you have there) is actually not a problem anymore. Most QC programs are also very mature and can use sophisticated techniques to solve the HF equations where simple iterative schemes would fail.

 

[sam_bell]

Thanks for the help cgk. I figured out what was going on by consulting an old book by Hartree himself (The Calculation of Atomic Structures, John Wiley & Sons 1957 Sec 5.3.1). I don't want to reproduce it because it's a page long series of manipulations on the Schroedinger equation. The derivation isn't enlightening, but the result is clever because it avoids needing to estimate any energy derivatives.

We use code similar to the Herman-Skillman program to determine the core energy levels of the atoms in a crystal. (We associate one spherical potential for every atom in the basis of the crystal. For metallic alloys, this approximation works reasonably.) The only major changes needed are to switch the Slater exchange to an LDA exchange-correlation functional and to replace the Schroedinger equation with a scalar relativistic counterpart.