[Wien] Electron density at the nucleus

[Stefaan Cottenier]

> I have run BeO lattice case (space group P63mc) using WIEN2K code and 
> found that the Be 1s state energy is = -6.204169219 Ry and the electron 
> density at Be nucleus (RTO001) due to 1s core state is = 34.428627. 

(just to fill out a small detail: the fact that you are able in this 
case to separate the contribution by 1s and 2s electrons, is because the 
1s is the only core state, and the 2s the only valence state)

> Then 
> the code was run again by reducing the BeO lattice parameters by 5%. As 
> a result, Be 1s state energy increases to -5.971227606 Ry as 
> qualitatively expected from uncertainty principle considerations. The 
> core force also increases as expected and the total charge in sphere 1 
> also increases by about 0.3%. However as a result of the compression, 
> the electron density at the nucleus due to the core 1s state decreases 
> to 34.331679 from the earlier value of 34.428627 by about 0.3%. This 
> result looks puzzling. As a result of compression, the kinetic energy of 
> Be 1s state should increase and the 1s state electrons should be 
> confined to a smaller volume thus increasing the electron density at the 
> nucleus. The electron density of 2s states of Be at the nucleus 
> increases due to the compression as expected. So the puzzle is why the 
> electron density of the core 1s state of Be is decreasing due to the 
> compression.

I'm not sure on this, just a try: your argument probably holds true for 
a hydrogen atom. But in Be, there is 1s as well as 2s at the nucleus. 
Doesn't the interaction between 1s (no nodes) and 2s (1 node) prevent 
you from making the same conclusions as for a single hydrogen orbital?

> I am interested to know how the calculation regarding the electron 
> density at the nucleus is being done in WIEN2K code and what are the 
> relevant subroutines to look at. If there is any published literature 
> then please also refer me to those papers.

The :RTO is found by taking the density at the nearest radial grid point 
(R0), and considering it to be constant over a sphere with radius R0. 
This is the R0 you find at the corresponding atom in case.struct.

It might be worthwhile to repeat the calculation with a R0 that is 
comparable to a nuclear radius. Does that lead to similar behaviour?

What happens if you take both 1s and 2s as valence states? (-8 as input 
for lstart)

There are issues about the divergence of the density at a point nucleus, 
and an incorrect behaviour near the nucleus that is due to the 
scalar-relativistic approximation. That makes the :RTO numbers 
suspicious to some extent (they do correlate reasonably with Mossbauer 
isomer shifts, but perhaps now with the decay rates you are looking at). 
These problems are discussed a.o. in the following two papers:

M. Filatov, Coordination Chemistry Reviews
http://dx.doi.org/10.1016/j.ccr.2008.05.002
(on an alternative method to determine electron density at the nucleus)

K. Koch et al., Phys. Rev. A
http://dx.doi.org/10.1103/PhysRevA.81.032507
(on p1/2 density at the nucleus by finite nucleus calculations)

Stefaan