From valence to core electrons

[rabbit99]

Dear Exciting-users,
dear Kay,

I have one similar question as ravi-mse
concerning the species-input-files.

I want to change some valence electrons
to core, e.g. the 3s and 3p states in iron.
If I only change the spcore flag from
'F' to 'T' I got strange results for the
energy-volume curve.
Have I to remove some local orbitals?

Thanks in advance,
rabbit

[francesco cricchio]

Hi!
Yes you could remove the local orbitals corresponding to the valence electrons you put in the core.
This will reduce your basis set and will speed up the calculation.

However BE CAREFUL! You can put in the core only electrons that are really core electrons, i.e. deep enough in energy.

So you could follow this procedure

1) Start by including in the core the lowest electrons in the default species and create a new species file, it is convenient to do it in the same dir as your calculation. You can also remove the local orbitals corresponding to such energy levels.

2) Perform 1 iteration (maxscl=1) with new species file.

3) Check EVALCORE.OUT, all the electrons in there must be core electrons, i.e. they must have energy  lower than 1.2 Htr as a rule of thumb.

4) Check the file EIGVAL.OUT or your band structure in BAND.OUT. If they appear ghost bands some of the local orbitals you removed they were actually needed.

5)If not too expensive it is also good to compare your results with the ones obtained with the default species file as a check: the total energy won't be the same but energy differences must be the same (between different volumes as example).

To clarify let's consider Cu example:

The default species file of Cu looks like:

'Cu'                                       : spsymb
'copper'                                   : spname
  -29.0000                                  : spzn
   115837.2717                              : spmass
  0.371391E-06    2.0000   44.7401   500    : sprmin, rmt, sprmax, nrmt
  10                                        : spnst
   1   0   1   2.00000    T                 : spn, spl, spk, spocc, spcore
   2   0   1   2.00000    T
   2   1   1   2.00000    T
   2   1   2   4.00000    T
   3   0   1   2.00000    T
   3   1   1   2.00000    F
   3   1   2   4.00000    F
   3   2   2   4.00000    F
   3   2   3   6.00000    F
   4   0   1   1.00000    F
   1                                        : apword
  0.1500   0  F                             : apwe0, apwdm, apwve
   0                                        : nlx
   5                                        : nlorb
   0   2                                    : lorbl, lorbord
  0.1500   0  F                             : lorbe0, lorbdm, lorbve
  0.1500   1  F
   1   2                                    : lorbl, lorbord
  0.1500   0  F                             : lorbe0, lorbdm, lorbve
  0.1500   1  F
   2   2                                    : lorbl, lorbord
  0.1500   0  F                             : lorbe0, lorbdm, lorbve
  0.1500   1  F
   3   2                                    : lorbl, lorbord
  0.1500   0  F                             : lorbe0, lorbdm, lorbve
  0.1500   1  F
   1   3                                    : lorbl, lorbord
  0.1500   0  F                             : lorbe0, lorbdm, lorbve
  0.1500   1  F
-2.5606   0  T

You can include in the core the 3p electrons since they are quite deep in energy. Then you can also remove the corresponding local orbitals.
You can also remove the f local obitals, since the f levels will be totally unoccupied in Cu. This will reduce  your basis set and speed up the calculation!
The new file Cu_new looks like:

'Cu'                                       : spsymb
'copper'                                   : spname
  -29.0000                                  : spzn
   115837.2717                              : spmass
  0.371391E-06    2.0000   44.7401   500    : sprmin, rmt, sprmax, nrmt
  10                                        : spnst
   1   0   1   2.00000    T                 : spn, spl, spk, spocc, spcore
   2   0   1   2.00000    T
   2   1   1   2.00000    T
   2   1   2   4.00000    T
   3   0   1   2.00000    T
   3   1   1   2.00000    T
   3   1   2   4.00000    T
   3   2   2   4.00000    F
   3   2   3   6.00000    F
   4   0   1   1.00000    F
   1                                        : apword
  0.1500   0  F                             : apwe0, apwdm, apwve
   0                                        : nlx
   3                                        : nlorb
   0   2                                    : lorbl, lorbord
  0.1500   0  F                             : lorbe0, lorbdm, lorbve
  0.1500   1  F
   1   2                                    : lorbl, lorbord
  0.1500   0  F                             : lorbe0, lorbdm, lorbve
  0.1500   1  F
   2   2                                    : lorbl, lorbord
  0.1500   0  F                             : lorbe0, lorbdm, lorbve
  0.1500   1  F
 

Now we run the Cu example with 1 iteration:

tasks
  0

maxscl
1

avec
  0.5  0.5  0.0
  0.5  0.0  0.5
  0.0  0.5  0.5

scale1
  6.83117

scale2
  6.83117

scale3
  6.83117

scale
  0.96

! large cut-off is required for Cu
rgkmax
  8.5

lmaxapw
  10

gmaxvr
  14.0

sppath
  './'

atoms
  1                                   : nspecies
  'Cu_new.in'                             : spfname
  1                                   : natoms
  0.0  0.0  0.0    0.0  0.0  0.0      : atposl, bfcmt

ngridk
  8  8  8

We check EVALCORE.OUT:

Species :    1 (Cu), atom :    1
n = 1, l = 0, k = 1 :   -324.0479415   
n = 2, l = 0, k = 1 :   -38.46759656   
n = 2, l = 1, k = 1 :   -33.68183125   
n = 2, l = 1, k = 2 :   -32.93033142   
n = 3, l = 0, k = 1 :   -3.714155757   
n = 3, l = 1, k = 1 :   -2.232244040   
n = 3, l = 1, k = 2 :   -2.137280134

We are safe since the 3p electrons have energy around -2 htr, deep enough to be treated as core.
We check also the file EIGVAL.OUT for the gamma point:

1   0.000000000       0.000000000       0.000000000     : k-point, vkl
     (state, evalsv, occsv, spnchr below)
     1 -0.7681996212E-01   2.000000000       1.000000000   
     2  0.2517618661       2.000000000       1.000000000   
     3  0.2517618661       2.000000000       1.000000000   
     4  0.2517619241       2.000000000       1.000000000   
     5  0.2952636853       2.000000000       1.000000000   
     6  0.2952636853       2.000000000       1.000000000   
     7   1.198694557       0.000000000       1.000000000   
     8   1.318021212       0.000000000       1.000000000   
     9   1.318021271       0.000000000       1.000000000   
    10   1.318021271       0.000000000       1.000000000   
    11   1.525014128       0.000000000       1.000000000 

It is ok since energy 1 correspond to 4s, the following 2,3,4,5,6 correspond to 3d.
There are no ghost bands. Also there are definetely no f states occupied, that's why we were allowed to remove the f local orbitals in the species file. But we definitely needed the d local orbitals since they are fully occupied.

As final check we could compare the energy difference between Cu fcc with 2 different lattice parameters (scale) obtained with the OLD (default) and the NEW species file, such energy difference must be the same:

OLD Cu scale 0.96 =  -1652.48652634
OLD Cu scale 1.02 =  -1652.48277831  
OLD DELTAE=3.8 mhtr

NEW Cu scale 0.96 = -1652.48692539
NEW Cu scale 1.02 = -1652.48327765
NEW DELTAE=3.7 mhtr

In conclusion same energy difference but computational time reduced about 5%!

OLD species file

Timings (CPU seconds) :
initialisation                        :         5.14
Hamiltonian and overlap matrix set up :         3.39
first-variational secular equation    :         7.02
charge density calculation            :        20.97
potential calculation                 :         5.86
total                                 :        42.39

NEW species file

Timings (CPU seconds) :
initialisation                        :         5.23
Hamiltonian and overlap matrix set up :         3.73
first-variational secular equation    :         6.54
charge density calculation            :        17.21
potential calculation                 :         6.95
total                                 :        39.66

Hope this discussion will help!
Best
Francesco Cricchio

PS Kay please correct me if I'm inaccurate on some point.

[rabbit99]

Dear Exciting-users,
dear Francesco,

first many, many thanks for the fast and detailed answer,
it helps a lot  :).

Maybe you could extend your explanation and
comment on the shortcuts in the species-files i.e.
apword, apwe0, apwdm, apwve
and especially lorbl, lorbord, lorbe0, lorbdm, lorbve.

Probably a trivial question, but why is the occupation number of the
4s-state 2.0 in your example and not 1.0?
I would naively expect that this should be 1.0 since there is
just one electron in the 4s-shell.

Thanks a lot :)!
rabbit

[rabbit99]

Hi all,

regarding the second question, I have
overlooked that we are regarding one k-point
of the whole grid.
The average over all kpoints give the correct
occupation, also when just treating only one kpoint
instead of the 8x8x8-grid...

best regards,
rabbit

[francesco cricchio]

Dear Rabbit,
the shortcuts in the species file mean the following:

apwdm is the derivative energy order of the solution of radial schroedinger equation (RDS) included in your basis set.
So in APW basis set there are only included solution of RDS (apwdm=1) with linearization energy "apwe0". In LAPW you have radial solutions (apwdm=1) plus first energy derivative (apwdm=2) to make the basis set energy independent and been able to describe your band structure in only 1 diagonalization. Radial functions and their derivatives can be calculated at different energies apwe0 (number of energies=apword).

Why you need derivatives of radial functions, or alternatively, local orbitals?

In the ordinary APW basis set you have an energy dependent secular equation:
DET(H(E)-E*O(E)=0)
that leads the correct solution only if the energy E was already the one corresponding to your band energy.
You would need to trace the roots of this determinats to find your correct band structure, an horrible job.

The big revolution has been done by making your basis set energy INDEPENDENT in energy (accuracy to 2nd order) by adding first energy derivative (apwdm=2) of the solutions of RDS to ordinary APW, so you get LAPW (L as linearised) basis set.

Or more simply and efficiently, and this is the default approach in EXCITING, you can make your basis set energy INDEPENDENT by adding some local orbitals to the ordinary APW.

These local orbitals are solution of RDS at energy "lorbe0" for l-quantum number "lorbl" inside the muffin tin.
Consequently they have an orbital character.
They are localised, i.e. they are zero outside muffin tin.
They could be also a linear combination of more radial solutions for different energies (number of energies=lorbord), or a linear combination of radial solution (lorbdm 0) plus its energy derivative (lorbdm 1).
Their energy "lorbe0" could be updated during the calculation, this is convenient for semicore states (the flag lorbve=T).

To summarize let's consider the local orbital in Cu.in default species:

1 3 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
-2.5606 0 T

Inside the muffin tin this a linear combination of 3 functions (lorbord=3):

1) the solution of radial schrodinger equation  (lorbdm=0) in nuffin tin for l=1 (lorbl=1) at energy 0.1500 htr

2) the first energy derivative (lorbdm=1) of RDS solutions for l=1 (lorbl=1) at energy 0.1500 htr

3) the solution of radial schrodinger equation for l=1 (lorbl=1)at energy -2.5606 htr (lorbdm=0).
Since there is  lorbve=T, such energy is updated during the calculation to better describe the core state level.
In ordinary calculation you do not actaully need this orbital since it is quite deep in energy. You could put the 3p electrons in the core (spn=3,spl=1) and remove this orbital.
However it is not expensive to keep it since it is a p orbital, then m=-1,0,1, then just few more elements in your basis.
Instead it is good idea for Cu to remove f local orbitals ( lorbl=3) since f states in Cu are far away unoccupied and in this case m=-3,-2,-1,0,1,2,3 ; many more basis states, more time to diagonalize matrices and so on.

FINAL CONSIDERATION: in general you do not need to change the parameters apword, apwe0, apwdm, apwve,lorbe0,lorbl,lorbdm,lorbve of existing local orbitals in the default file.
You need to ADD local orbitals if you move electrons from the core to the valence. Since they are needed to keep your basis set energy independent and give you the correct band structure.

EXAMPLE default Cu.in:

'Cu' : spsymb
'copper' : spname
-29.0000 : spzn
115837.2717 : spmass
0.371391E-06 2.0000 44.7401 500 : sprmin, rmt, sprmax, nrmt
10 : spnst
1 0 1 2.00000 T : spn, spl, spk, spocc, spcore
2 0 1 2.00000 T
2 1 1 2.00000 T
2 1 2 4.00000 T
3 0 1 2.00000 T
3 1 1 2.00000 F
3 1 2 4.00000 F
3 2 2 4.00000 F
3 2 3 6.00000 F
4 0 1 1.00000 F
1 : apword
0.1500 0 F : apwe0, apwdm, apwve
0 : nlx
3 : nlorb
0 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
1 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
2 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
1 3 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
-2.5606 0 T

If you move the 3s electron (spn=3,spl=0) from valence to core you need an additional local orbital with (changing) energy approximately corresponding to his energy level.
So you can estimate such energy by looking at the file EVALCORE.OUT in the calculation with the default species (1 iteration,maxscl=1, is enough). We get

Species :    1 (Cu), atom :    1
n = 1, l = 0, k = 1 :   -324.0479415
n = 2, l = 0, k = 1 :   -38.46759656
n = 2, l = 1, k = 1 :   -33.68183125
n = 2, l = 1, k = 2 :   -32.93033142
n = 3, l = 0, k = 1 :   -3.714155757

Ok, so let's estimate a lower bound as -3.75 htr, the code will update it upward.
So we now add this local orbital in Cu_new.in:

'Cu' : spsymb
'copper' : spname
-29.0000 : spzn
115837.2717 : spmass
0.371391E-06 2.0000 44.7401 500 : sprmin, rmt, sprmax, nrmt
10 : spnst
1 0 1 2.00000 T : spn, spl, spk, spocc, spcore
2 0 1 2.00000 T
2 1 1 2.00000 T
2 1 2 4.00000 T
3 0 1 2.00000 F
3 1 1 2.00000 F
3 1 2 4.00000 F
3 2 2 4.00000 F
3 2 3 6.00000 F
4 0 1 1.00000 F
1 : apword
0.1500 0 F : apwe0, apwdm, apwve
0 : nlx
5 : nlorb !!!!CHANGED TO 5
0 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
1 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
2 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
1 3 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
-2.5606 0 T
0 3 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
-3.7500 0 T

In this file I also removed the f local orbitals (lorbl=3) since they are useless, and I updated the number of local orbitals nlorb=5 (IMPORTANT!!!!, if not the last local orbital will not be read by the code).

Let's now look at the eigenvalues in EIGVAL.OUT for Gamma point
1   0.000000000       0.000000000       0.000000000     : k-point, vkl
     (state, evalsv, occsv, spnchr below)
     1  -3.777120563       2.000000000       1.000000000
     2  -2.223187476       2.000000000       1.000000000
     3  -2.223187476       2.000000000       1.000000000
     4  -2.223187476       2.000000000       1.000000000
     5 -0.3782050929E-01   2.000000000       1.000000000
     6  0.2158364058       2.000000000       1.000000000
     7  0.2158364058       2.000000000       1.000000000
     8  0.2158364058       2.000000000       1.000000000
     9  0.2531067813       2.000000000       1.000000000
    10  0.2531067813       2.000000000       1.000000000
    11   1.247389266       0.000000000       1.000000000
    12   1.360060769       0.000000000       1.000000000
    13   1.360060769       0.000000000       1.000000000
    14   1.360060769       0.000000000       1.000000000
    15   1.541151569       0.000000000       1.000000000

as you see we get the additional -3.77120563 corresponding to the 3s state we added in the valence.
If we didn't have added the additional local orbitals our basis set. we would not have been able to recover such deep energy state.
You would not have got any 3s state ,or worse, you would have got a wrong one shifted up in energy.
We can try it:

1   0.000000000       0.000000000       0.000000000     : k-point, vkl
     (state, evalsv, occsv, spnchr below)
     1  -2.526193107       2.000000000       1.000000000
     2  -2.526193107       2.000000000       1.000000000
     3  -2.526193107       2.000000000       1.000000000
     4 -0.3672795995E-01   2.000000000       1.000000000
     5 -0.2497265583E-01   2.000000000       1.000000000
     6 -0.2497265582E-01   2.000000000       1.000000000
     7 -0.2497265582E-01   2.000000000       1.000000000
     8 -0.2978748066E-02   2.000000000       1.000000000
     9 -0.2978748066E-02   2.000000000       1.000000000
    10   1.276216776       0.000000000       1.000000000
    11   1.346818776       0.000000000       1.000000000
  

No 4s state!

If you would like to move back the electron from valence to core, you may remove (not compulsory) this additional local orbital to speed up calculation time.

In general, as I pointed out in my previous answer, you only need in the valence electrons with energy higher than 1.2 htr as rule of thumb.

To understand better the great power and accuracy of this method I suggest you to read:

APW+LO
E. Sjostedt and L. Nordstrom J. Phys Cond Matt 14 (2002) 12485
E. Sjostedt and L. Nordstrom Solid State Comm.114 (2000) 15-20
Madsen et al Phys Rev B 64, 195134

LAPW
O. K. Andersen , Phys Rev B 12 (1975) 3060

Original APW
J.C. Slater. Phys Rev, 51 (1937) 846

Hope this will help!
Best
Francesco

[rabbit99]

Dear Francesco,

again many, many thanks for this nice and detailed answer and
the references :)!

best regards,
rabbit

[francesco cricchio]

Dear Rabbit,
the computational cost of adding local orbitals (LO) does not come from matrix diagonalization as I wrongly stated in my prevoius answer. In fact the number of local orbitals does not affect the number of plane waves (the size of your hamiltonian). Such computational cost simply comes from the fact that with more LO you calculate more radial functions in muffin tin and you need to match them with plane waves at the boundary.
Best
Francesco

PS Kay please correct me if I'm wrong!

[francesco cricchio]

Dear Rabbit,

I forgot to answer about your other question regarding occupancies.
In case of Cu example, that is spin unpolarised calculation since Cu is not magnetic, you have 2 electrons per energy (spin up/spin down)

In other kpoints different than Gamma you do not get integer occupancies, but fractional, since, of course we are in a solid.

If we force the calculation to be spin-polarised (spinpol=.true. and a small bfieldc to break symmetry) we do get one electron per energy, but then you get TWICE the same energy for Cu that is non magnetic.

You get something like:
1   0.000000000       0.000000000       0.000000000     : k-point, vkl
     (state, evalsv, occsv, spnchr below)
     1  -2.165782961       1.000000000       1.000000000       0.000000000
     2  -2.165782961       1.000000000       1.000000000       0.000000000
     3  -2.165782772       1.000000000       1.000000000       0.000000000
     4 -0.7685649052E-01   1.000000000       1.000000000       0.000000000
     5  0.2517253328       1.000000000       1.000000000       0.000000000
     6  0.2517253328       1.000000000       1.000000000       0.000000000
     7  0.2517253908       1.000000000       1.000000000       0.000000000
     8  0.2952271516       1.000000000       1.000000000       0.000000000
     9  0.2952271516       1.000000000       1.000000000       0.000000000
    10   1.198657937       0.000000000       1.000000000       0.000000000
    11   1.322047742       0.000000000       1.000000000       0.000000000
    12   1.322047808       0.000000000       1.000000000       0.000000000
    13   1.322047808       0.000000000       1.000000000       0.000000000
    14   1.524977483       0.000000000       1.000000000       0.000000000
    15  -2.165709898       1.000000000       0.000000000       1.000000000
    16  -2.165709898       1.000000000       0.000000000       1.000000000
    17  -2.165709709       1.000000000       0.000000000       1.000000000
    18 -0.7678343371E-01   1.000000000       0.000000000       1.000000000
    19  0.2517983994       1.000000000       0.000000000       1.000000000
    20  0.2517983994       1.000000000       0.000000000       1.000000000
    21  0.2517984574       1.000000000       0.000000000       1.000000000
    22  0.2953002190       1.000000000       0.000000000       1.000000000
    23  0.2953002190       1.000000000       0.000000000       1.000000000
    24   1.198731176       0.000000000       0.000000000       1.000000000
    25   1.322120826       0.000000000       0.000000000       1.000000000
    26   1.322120891       0.000000000       0.000000000       1.000000000
    27   1.322120891       0.000000000       0.000000000       1.000000000
    28   1.525050773       0.000000000       0.000000000       1.000000000

More info on magnetism on
Phys Rev B 66, 014447 (2002)
for example.

Best
Francesco