[Wien] Help on Core hole calculations

Fri May 17 09:56:56 CEST 2019

> I have read some of the discussions in the mailing list (https://www.mail-archive.com/wien@zeus.theochem.tuwien.ac.at/msg11037.html) and somewhere else but I couldn’t find some useful information to solve my doubts; I actually got more confused. I am doing core hole calculations for alpha-Fe2O3 (hematite) and I want to compute just the DOS for this system for holes from Fe 1s, 2p, and 3s (and possibly 3p). I have some doubts regarding the correct procedure. I would appreciate any feedback on the following points.
> 
> 1) I am not sure if the crystal structure I made is correct. As a starting point I used the .cif file downloaded from AMCSD web site with space group R-3c (167), then converted it to .struct file with RMT Fe=1.91, RMT O=1.64. This structure is corundum and all Fe atoms are equivalent. However, to reproduce experimental observations, spin-polarized calculations with LDA+U are required. So I need to make Fe atoms nonequivalent and assign the spin orientation along the [0001] direction, with Fe atoms alternating spin in planes {0001}. For this I created a 1X1X1 supercell (just to test) from the original file and later execute x nn, x sgroup and x symmetry. Sgroup suggested a new space group P3c1 (158) with 6 nonquivalent Fe atoms and 3 O equivalent atoms which I accepted and then labelled the Fe atoms as Fe1,... Fe6. Then I called instgen and assigned the spin orientation to each Fe and O atom (oxygen is non-magnetic), and run the rest of the initialization programs with 100 k points, ecut= -6 Ry, RKmax=7, and for the functional PBE. Thus, since now I have six nonequivalent Fe atoms, I would have, after the scf cycle, 6 partial DOS (actually 12, 6 up and 6 down), one for each Fe atom. Is this correct?

No. You should not run a supercell for a R structure, since this will 
even with 1x1x1 create a 3 times (conventional hexagonal) cell.
Instead,
i) remove the spacegroup and change to "R"-lattice.
ii) split the 4 Fe positions into 2 and 2 (Fe1 and Fe2). This splitting 
can be done in different ways (up-up-dn-dn along z; or in other ways. 
Please check literature which one is correct). After nn (or sgroup) you 
should have only 2 Fe positions.

> 2) For core hole calculations there exist at least two options to send the core electron. After removing a core hole (from a core level: 1s, 2s, 2p, or 3s) in the file case.inc we can send this electron to:
> 
> a) the valence electrons in case.in2c; or,
> 
> b) to the background charge in case.inm
> 
> I’m not sure where to add the electron, what physical criterion should I apply? In any case, I decided to try both cases. For a), I removed the electron from the file case.inc for the levels 1s, 2s and 3s (of course, in three different computing sessions). For b), I removed one electron from the same three levels in case.inc. The core electron (for each calculation) was removed from only one Fe atom, say Fe1 (the remaining atoms kept the same configuration), as follows:

Much more critical than where to add the extra electron is, that you now 
create a big supercell (as big as you can handle), break symmetry and 
only one ONE Fe atom should have the core hole. This ways we "simulate" 
experiment, where never all atoms of the crystal are excited at the same 
time, but only a tiny fraction, while the rest remains in the 
groundstate and can contribute to screening.

> 3) Here I have a question, where exactly is the electron added? The 3d band, the 4s band? In the user’s guide, it reads that the -12 stands for the lower energy cut-off for defining the range of occupied states and NE is the number of electrons (per unit cell) in that energy range. Does it mean that those 277 electrons are in states with energies higher than -12 Ry. In other words, is it possible to decide where exactly in the valence to put the excited electron?

No. It will add the electron in the first few conduction band states - 
whatever these states are (in your example of course mainly into Fe-3d 
states of the atom with the core hole.

> 4) Here I have another question. The user’s guide reads that the second number in the line, bgch, means “background charge to apply to the cell”. Normally, this is set to 0.0, meaning that the cell is not charged??? I do not have clear what is called background charge. Could you please explain what is the nature of this background charge. As I understand, the cell/material is electrically neutral, and we only have charges due to nuclei and electrons, where does this background charge come from? This question is connected with the following.
> 
> Since I want LDA+U, I also included the files case.orb and case.indmc. Then I executed:
> 
> runsp_lapw -orb -ec 0.001
> 
> After the scf cycle, Prof. Blaha suggested to check charge neutrality. I did in all cases using grep :NEC01 in case.scf and obtained:
> 
> NEC01: NUCLEAR AND ELECTRONIC CHARGE 456.00000 454.90177
> 
> so, I have no charge neutrality. I added the electron to either the valence or background before running the scf cycle, as suggested everywhere, why there is no charge neutrality?

Cells must always be neutral. Either you have a huge core leakage (NEVER 
change cutt-off to -4Ry, it means huge core leakage), or you have done 
something wrong when putting the extra electron or the background 
charge. It could also be that your core-hole atom has MULT=2, then 
obviously you would have to add 2 electrons - but you should never have 
core-hole atoms with MULT=2 (see above).

> 5) After the scf-cycle, for the case in which the electron was added to the valence electrons, I removed the electron from the file case.in2c and computed the DOS. But, what should I do for the case in which the electron was added to the background? Should I reset to 0.0 in case.inm before calculating the DOS?

It does not matter. A DOS calculation does not involve mixer !!??

> 6) My final question. When we run lstart we define the cut-off energy, say – 6 Ry, that separates core from valence states. In my case, with this energy, I have, in case.outputst, as core states 1s, 2s, 2p, and 3s; and as valence states, 3p, 3d, 4s. However, in the file case.inst core states are considered those of Ar and 3d and 4s as valence states. So, during the computation of the scf, how is the 3p level taken, as a core or a valence state? If I changed cut-off energy=-4 Ry to consider the 3p as valence state (and so do core hole calculations for this level) then I get a charge leakage warning! I am confused about this apparent contradiction in these two files.

The case.inst file does NOT decide which states are core/valence. It is 
only for convenience that one can put Ne,Ar,.. as a shortcut for the 
specification of the total electronic configuration.

> 
> 
> I would be happy if you could help me.
> 
> Best Regards
> 
> Israel Perez
> 
> 
> Institute of Engineering and Technology
> 
> Department of Physics and Mathematics,
> 
> Universidad Autónoma de Ciudad Juárez
> 
> Av. del Charro 450 Nte., Col. Partido Romero,
> 
> Ciudad Juárez, Juarez Chihuahua. Mexico C. P. 32310
> 
> Tel: +52 (656) 688 4887
> 
> 
> National Council of Science and Technology
> 
> Insurgentes Sur No. 1582,
> 
> Col. Crédito Constructor,  C.P. 03940
> 
> Del. Benito Juárez, México D. F.
> 
> 
> _______________________________________________
> Wien mailing list
> Wien at zeus.theochem.tuwien.ac.at
> http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien
> SEARCH the MAILING-LIST at:  http://www.mail-archive.com/wien@zeus.theochem.tuwien.ac.at/index.html
> 

-- 

                                       P.Blaha
--------------------------------------------------------------------------
Peter BLAHA, Inst.f. Materials Chemistry, TU Vienna, A-1060 Vienna
Phone: +43-1-58801-165300             FAX: +43-1-58801-165982
Email: blaha at theochem.tuwien.ac.at    WIEN2k: http://www.wien2k.at
WWW:   http://www.imc.tuwien.ac.at/TC_Blaha
--------------------------------------------------------------------------