[Wien] Problem of HCP Tb (the third time)

Stefaan Cottenier Stefaan.Cottenier at UGent.be
Mon Feb 1 13:20:19 CET 2010


Comments are pasted at the relevant places:

>        After a regular GGA+SO scf, commamd line: *runsp -ec 0.0001 -cc 
> 0.0001 -so -i 1000*
>        then i got some results as follows:
>        :ENE:    -46876.608343 Ry
>        :MMTOT:  11.92563      uB
>        :MMI001: 5.81891       uB
>        :MMINT:  0.28781       uB
>        **
> *     I changed the case.indmc with '1 3' and x lapwdm -c -so -up, and I 
> got the following result:*
> ** 
> case.scfdmup
> *************************************************************
> Spin-polarized + s-o calculation, M||  0.000  0.000  1.000
>   Calculation of <X>, X=c*Xr(r)*Xls(l,s)
>   Xr(r)    =           I
>   Xls(l,s) = L(dzeta)
>   c=  1.00000
>   atom   L        up          dn         total
> :XOP  1  3     0.00332     1.25294     1.25626
> *************************************************************
> 
>       This is the orbital moment with SO only, and it is 1.25294 uB, 
> direction is down.
>  
>       *Here, I am confused: *
> *            <1> The orbital moment is still too small, it should be 
> around 5 uB. I don't know why ?*

This is to be expected. With GGA+SO, the f-electrons are not 
sufficiently localized. All f-orbitals of the minority spin are near 
E_fermi, are more or less equally populated, and this leads to a zero 
(or low) orbital moment.

> *            <2> In case.inst, the initial spin is up, I think the 
> orbital moment should be in the same*
> *                direction as spin. Or **how can the total moment be 
> around 10 uB ??*

For the second half of the lanthanide series, the orbital and spin 
moment have the same direction, indeed (Hund's third rule).

>      *(3)*based on the GGA+SO scf, then i added the U parameter.
>         case.inorb and case.indm were prepared as follows:
>  
> *case.indm*

Are you sure you need case.indm? Do "grep lapw2 :log" : if lapw2 has the 
-c option, then lapwdm needs it too (and then you need case.indmc). If 
an old case.indmc is around with different content, that might spoil 
your calculation.

> *************************************************************
> -9.                      Emin cutoff energy
>  1                       number of atoms for which density matrix is 
> calculated
>  1  1  3      index of 1st atom, number of L's, L1
>  0 0           r-index, (l,s)index
> *************************************************************
> *case.inorb  *
> *************************************************************
>   1  1  0                     nmod, natorb, ipr
> PRATT  1.0                    BROYD/PRATT, mixing
>   1 1 3                          iatom nlorb, lorb
>   1                              nsic 0..AFM, 1..SIC, 2..HFM
>   0.50 0.00        U J (Ry)   Note: we recommend to use U_eff = U-J and J=0
> *************************************************************
>  
>        run a regular GGA+SO scf, commamd line: *runsp -ec 0.0001 -cc 
> 0.0001 -so -orb -i 1000*
>        then i got some results as follows:
>        :ENE:    -46876.47444 Ry
>        :MMTOT:  13.45593      uB
>        :MMI001: 6.25394       uB
>        :MMINT:  0.94806       uB
>        *:ORB001: 0.00439       uB *
> ** 
> *      After the three steps above, I found orbital moment in case.scf 
> which is 0.00439 uB. comparing with the " 1.25294 uB" generated last 
> step, why it's so small after plus U( Ueff = 0.5 Ry) ??*

Impossible to tell with this information (but see hereafter, the dmat 
files). Apparently adding U leads to a redistribution of the electrons 
over the f-orbitals, with a zero orbital moment as a result.

>       Then, I changed the last line of case.indmc from '0 0' to '1 3' 
> after convergence,
> *case.indmc*
> *************************************************************
> -9.                      Emin cutoff energy
>  1                       number of atoms for which density matrix is 
> calculated
>  1  1  3      index of 1st atom, number of L's, L1
>  1 3           r-index, (l,s)index
> *************************************************************
>  
>     then I run *'x lapwdm -c -so -up'*, and I find the orbital moment in 
> case.scfdmup is as follows:
>  
> *case.scfdmup*
> *************************************************************
>  Spin-polarized + s-o calculation, M||  0.000  0.000  1.000
>   Calculation of <X>, X=c*Xr(r)*Xls(l,s)
>   Xr(r)    =           I
>   Xls(l,s) = L(dzeta)
>   c=  1.00000
>   atom   L        up          dn         total
> :XOP  1  3    -0.00187     0.00626     0.00439
> *************************************************************
>  
>      * Orbital moment is 0.00439 uB, I was completely confused ??????*

This is just another way of calculating the number that is reported as 
:ORB in case.scf, hence you should find the same number, indeed.

> *   *The following files are the dmatup/dmatdn files: ( I didn't get it)
>  
> case.dmatup
> *************************************************************
>     1 atom density matrix
>     3  0.000000  0.000000 -0.001873 L, Lx,Ly,Lz in global orthogonal system
>   0.99122776E+00 -0.11754944E-36    0.15671123E-12 -0.12705494E-20
>  -0.33034285E-19  0.84703295E-21    0.12447650E-12  0.13764285E-20
>  -0.29646153E-19  0.00000000E+00   -0.58634492E-12  0.84703295E-21
>  -0.51769973E-03  0.19058241E-20
>   0.15671123E-12  0.12705494E-20    0.98989657E+00  0.42317797E-36
>  -0.16740146E-12 -0.16940659E-20   -0.72844833E-19  0.00000000E+00
>   0.78433837E-12  0.42351647E-21   -0.12705494E-19  0.00000000E+00
>   0.51044993E-12 -0.25410988E-20
>  -0.33034285E-19 -0.84703295E-21   -0.16740146E-12  0.16940659E-20
>   0.98949500E+00  0.56423729E-36   -0.11551635E-13  0.21175824E-21
>   0.51669010E-19  0.33881318E-20   -0.97451913E-12 -0.21175824E-20
>  -0.38963516E-19 -0.33881318E-20
>   0.12447650E-12 -0.13764285E-20   -0.72844833E-19  0.00000000E+00
>  -0.11551635E-13 -0.21175824E-21    0.98721653E+00 -0.29779190E-36
>   0.68056484E-14  0.21175824E-21   -0.78668185E-19  0.33881318E-20
>  -0.10434050E-12  0.21175824E-21
>  -0.29646153E-19  0.00000000E+00    0.78433837E-12 -0.42351647E-21
>   0.51669010E-19 -0.33881318E-20    0.68056484E-14 -0.21175824E-21
>   0.98801066E+00  0.31346516E-36    0.18533121E-12 -0.42351647E-21
>  -0.39069395E-19  0.16940659E-20
>  -0.58634492E-12 -0.84703295E-21   -0.12705494E-19  0.00000000E+00
>  -0.97451913E-12  0.21175824E-20   -0.78668185E-19 -0.33881318E-20
>   0.18533121E-12  0.42351647E-21    0.98954264E+00 -0.20375235E-36
>  -0.15970123E-12 -0.21175824E-21
>  -0.51769973E-03 -0.19058241E-20    0.51044993E-12  0.25410988E-20
>  -0.38963516E-19  0.33881318E-20   -0.10434050E-12 -0.21175824E-21
>  -0.39069395E-19 -0.16940659E-20   -0.15970123E-12  0.21175824E-21
>   0.99133429E+00  0.10971281E-36
> *************************************************************

This is a 7x7 matrix of complex numbers. The horizontal index runs over 
m (-3, -2, -1, 0, 1, 2, 3 -- or the other way around, I always forget), 
and the vertical index does the same. The diagonal elements give the 
occupation of the corresponding m-orbital. Hence, from this matrix you 
read the following occupations for the m-orbitals (f-up) :

m=-3 : 0.99122776E+00
m=-2 : 0.98989657E+00
(...)
m=+3 : 0.99133429E+00

Hence, all orbital occupied by exactly 1 electron -- as it should be, 
the up-shell is completely filled.

The dn-orbitals will be more interesting:

> case.dmatdn
> *************************************************************
>   1 atom density matrix
>     3  0.000000  0.000000  0.006263 L, Lx,Ly,Lz in global orthogonal system
>   0.76999012E-03  0.14497764E-36   -0.82982858E-13  0.44998625E-21
>  -0.72844833E-19  0.35998900E-20   -0.86064816E-13 -0.16940659E-20
>  -0.30969642E-19 -0.50821977E-20    0.19253073E-11  0.19905274E-19
>  -0.23424275E-03 -0.13552527E-19
>  -0.82982858E-13 -0.44998625E-21    0.16025310E-02  0.38399482E-36
>   0.21974289E-12  0.84703295E-21    0.18973538E-18  0.00000000E+00
>   0.48535084E-12 -0.16940659E-20   -0.27952087E-19 -0.25410988E-20
>  -0.18984536E-11 -0.15246593E-19
>  -0.72844833E-19 -0.35998900E-20    0.21974289E-12 -0.84703295E-21
>   0.13365436E-02 -0.13322269E-36    0.21297389E-12  0.14823077E-20
>  -0.75385932E-19  0.33881318E-20   -0.17868594E-12  0.78085850E-21
>  -0.32187252E-19 -0.16940659E-20
>  -0.86064816E-13  0.16940659E-20    0.18973538E-18  0.00000000E+00
>   0.21297389E-12 -0.14823077E-20    0.98011221E+00  0.51721751E-36
>  -0.20662598E-12  0.00000000E+00    0.19989978E-18  0.13552527E-19
>   0.23065083E-13 -0.16940659E-20
>  -0.30969642E-19  0.50821977E-20    0.48535084E-12  0.16940659E-20
>  -0.75385932E-19 -0.33881318E-20   -0.20662598E-12  0.00000000E+00
>   0.20209675E-02 -0.76798964E-36   -0.22783686E-12 -0.15881868E-20
>  -0.92008954E-19 -0.16940659E-20
>   0.19253073E-11 -0.19905274E-19   -0.27952087E-19  0.25410988E-20
>  -0.17868594E-12 -0.78085850E-21    0.19989978E-18 -0.13552527E-19
>  -0.22783686E-12  0.15881868E-20    0.26211585E-02 -0.65827684E-36
>  -0.67107735E-13  0.21175824E-20
>  -0.23424275E-03  0.13552527E-19   -0.18984536E-11  0.15246593E-19
>  -0.32187252E-19  0.16940659E-20    0.23065083E-13  0.16940659E-20
>  -0.92008954E-19  0.16940659E-20   -0.67107735E-13 -0.21175824E-20
>   0.19503738E-02 -0.34481168E-36
> *************************************************************

m=-3 : 0.76999012E-03  == 0
m=-2 : 0.16025310E-02  == 0
m=-1 : 0.13365436E-02  == 0
m= 0 : 0.98011221E+00  == 1
m=+1 : 0.20209675E-02  == 0
m=+2 : 0.26211585E-02  == 0
m=+3 : 0.19503738E-02  == 0

This shows you have 1 f-dn electron, and it is the m=0 orbital. This 
orbital contributes 0 to the orbital moment, consistent with what you 
have found.

You should try now to move this single electron to the m=3 orbital. To 
do this, interchange 0.98011221E+00 and 0.19503738E-02 in case.dmatdn, 
run orb (x orb, x orb -up, x orb -dn) to create orbital potentials 
consistent with these newly made dmat files, and then do an scf cycle 
with frozen orbital potentials (= -orbc). When this is roughly 
converged, do a normal scf cycle with very small mixing value (at least, 
that was how it worked a few years ago, no idea how the improved mixer 
will interfere with this). Keep your fingers crossed and hope that the 
electron will stay in m=+3 (it might as well go back to m=0, to another 
orbital, or get split over several orbitals).

The five different cases for Tm in Table I of 
http://dx.doi.org/10.1103/PhysRevB.74.014409 are produced this way.

The total energy might tell you which occupation is the best one, 
although that might not always be what you find in reality (see 
discussion in that paper).

If this works, you'll get (7-1)+(+3)=9 mu_B for the f-electron moment 
inside the spheres (6 spin moment, 3 orbital moment). And you had 0.94 
mu_B in the interstitial region (d-electrons?). That gives you the 10 
mu_B you're looking for (note that the configuration which you find here 
is different from the f^9s^2 for the free Tb atom: you find f^8d^1s^2. 
The free atom has (7-2)+(+3+2)=10 mu_B for f-electrons only, and an 
extra 1 mu_B by the d-electron (ferromagnetic intra-atomic 5d-4f coupling)).

Stefaan



More information about the Wien mailing list