Hf Atomic Calc: electronic structure problems

[grassman]

Hello all. We're doing some adsorption studies of Hf on semiconductors in my group, but we've run into a curiosity during the Hf atomic calcs, a problem we've never seen before. We're getting various electronic configurations of the free Hf atom depending on a few different input parameters, which we would have thought to be pretty flexible (a general INCAR is given below, with the parameters that had an unpredicted effect indicated; NOTE: Since we are doing these calcs for adsorptions on semiconductor slabs, we're trying to keep all of the settings exactly the same). This has been tested both with Gamma k-pt vs. k-pt mesh, cubic supercell (20x20x20A) vs. orthorhombic (~9x18x31A), but these were found to basically have no effect. Also, we are using the PAW-PBE potentials.

Running the spin-polarized (ISPIN=2) calc we get the correct overall electron configuration (6s2,5d2; triplet state), but we're seeing variations in the eigenvalues for the 5d orbitals depending on SIGMA and NELMDL, as well as the starting CHGCAR/WAVECAR (i.e. new calc vs. continuation from spin-restricted calc), which of course also gives different energies (all energies given are without entropy).

We're getting three basic variations of the Hf atomic electronic structure: 1) all 5d orbitals degenerate with equal charge, 2) 5d orbitals split 2-3 (deltaE ~ 0.6 eV), with electrons in lower energy 2-level set, and 3) 5d orbitals split 1-2-2 (deltaE1 ~ 0.27 eV, deltaE2 ~ 0.35 eV), with lowest energy 1-level filled, and middle-energy 2-level set half-filled. The lowest total energy configuration is the 2-3 split, followed by the 1-2-2 split, and finally the fully degenerate 5d system (with a total energy spread of almost 0.4 eV).

One of the main problems is that we can't seem to really find any data on what the exact configuration of atomic Hf should be. For an free atom with no external fields, it seems odd that the level would split like that, except maybe because of strong spin-orbit coupling? Does anyone know? The absolute lowest energy calc I got was with using SIGMA=0.01 and no NELMDL, which makes me think it's the most solid of the bunch, but the fact that it's a 2-3 split seems odd to me. However, I have very little experience working with transition metals, and especially with heavy atoms such as Hf, so it's entirely possible that I'm missing something. So, if anyone has any info or suggestions, I'd love to hear them. Thanks. (Also, I can provide more results/details for anyone interested, but since this message is already quite large, I decided to leave them out.)

Regards,
Tyler Grassman

----------------------
INCAR:
System = Hf

Startparameter for this run:
ISTART = 0
ICHARG = 2
INIWAV = 1

Electronic Relaxation
ENCUT = 500 eV
ENAUG = 2000 eV
ADDGRID = .TRUE.
NELM = 60
NELMIN = 4
NELMDL = -12 >> -12 vs. -3 vs. commented out
EDIFF = 1E-06
LREAL = Auto
PREC = Normal
ISMEAR = 0
SIGMA = 0.15 >> 0.15 vs. 0.01
ALGO = Normal
ISYM = 0
ISPIN = 2
WEIMIN = 0
LORBIT = 10
MAXMIX = 50
LMAXMIX = 4

Ionic relaxation
IBRION = -1

Parallelization Inputs
LPLANE = .TRUE.
LSCALU = .FALSE.
NSIM = 8
NPAR = 8
IMAGES = 0

[admin]

your problems are most probably due to the fact that the pseudopotentials are generated for the bulk electronic ground state configuration (6s1 5d3), not for the free atom ground state (6s2 5d2).
If the calculation does not converge to the correct Hund's rule result,
you have to enforce that result by explicitely occupying the levels
please have a look into the Examples:Atoms chapter of the online handbook:
use
KPOINTS Gamma only
INCAR :
SIGMA=0.0001 # sharp discrete atomic levels
ISMEAR= -2
LDIAG=.False.
ISPIN=2; NUPDOWN=2
FERWE= 3*1.0 (+ fill up to NBANDS with 0.0)
FERDO= 1*1.0 (+fill up to NBANDS with 0.0)
please check the levels occupancies in OUTCAR to verify whether the spin state and the occupancies of the levels are correct.

[grassman]

Thank you for the response. Yes, I didn't realize the POTCAR was generated with respect to the bulk ground state (s1d3), since the PAW-PBE POTCAR file only states 6s5d (but I just checked the GGA US-PP POTCAR and it does actually state s1d3). So, it makes sense that I'm getting odds results. It seems like I may not necessarily need to force a particular electronic configuration, however, as the calculation is basically getting it right (except for the splitting of the d orbital). I'm still guessing that the d orbital should be fully degenerate, though, so I guess I can force that configuration, unless anyone else knows this to be untrue? Thank you.

[admin]

you are right, for a really 'free atom' (with the corresponding spheric symmetry) the d-levels should be degenerate. However, as vasp is based on periodic boundary conditions, the highest possible symmetry is that of a cubic box (Oh), because space has to be filled by unit cells without leaving void space, which automaticaly induces the eg, t2g splitting oif the levels. (2-3). If the symmetry is reduced (ISYM=0 for SOC calculations), the splitting may even be increased to 1-2-2 .

[grassman]

Ah yes, that makes perfect sense. Okay, I will work on setting the occupancies manually. Thank you very much for your help and the explanation.