[Wien] Problems of HCP Tb!!! (the second time)

[Stefaan Cottenier]

>     The reason I use *kpoints = 1000* is because I refered to S. 
> Cottenier's book ---Density Functional Theory and the Family of 
> (L)APW-methods: a step-by-step introduction, it uses kpoints = 912 for 
> hcp-Cd 

That is in the irreducible part of the Brillouin zone (= "number of
lines in case.klist"). That corresponds to (roughly) 20000 as input for
kgen.

> reasonable. I also use kpoints = 4000 for HCP Tb, the total energy only 
> changed *0.4 mRy*, so kpoints = 1000 is enough.

That might be fine -- depends on the properties you are finally
interested in.

>     According to Stefaan's suggestions, I recalculate HCP Tb. ( 
> Userguide says first use +SO, then use +U in page 85 of 2009's userguide.)
>  
> *order: GGA==>>GGA+SO==>>GGA+SO+U*
> ** 
> *    **(1)*run a regular GGA scf, commamd line: *runsp -ec 0.0001 -cc 
> 0.0001 -i 1000*
>        then i got some results as follows:
>        :ENE:    -46876.362493 Ry
>        :MMTOT:  12.01107      uB
>        :MMI001: 5.86470       uB
>        :MMINT:  0.28167       uB
>  
> *    **(2)*based on the GGA scf, i added the spin-orbit coupling.
>        case.indmc, case.inso and case.inorb were prepared as follows:
>  
> *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
>  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
> *************************************************************
> *case.inso*
> *************************************************************
> WFFIL
>  4  1  0                      llmax,ipr,kpot
>  -9.0000   2.0000           emin,emax (output energy window)
>    0.  0.  1.                 direction of magnetization (lattice vectors)
>  1                           number of atoms for which RLO is added
>  1   -1.7      0.01       atom number,e-lo,de (case.in1), repeat NX times
>  0 0 0 0 0                    number of atoms for which SO is switch 
> off; atoms
> *************************************************************
>  
>        run 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 added spin orbit coupling, but there is no :ORB information, i 
> don't why ?*
> ** 

That is just an output feature. :ORB is written only when LDA+U is used.
If you want to know the orbital moment with SO only, you have to run
lapwdm once after the scf cycle with '1 3'.

>      *(3)*based on the GGA+SO scf, then i added the U parameter.
>         case.inorb and case.indm were prepared as follows:
>  
> *case.indm*
> *************************************************************
> -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 ( I don't know why it is so small ?)*
> ** 

That can be due to the m-orbitals occupied by the f-electrons of the
minority spin. If one is in m=-2 and the other in m=+2, the orbital
moment is +2-2=0. If they are equally distributed over all m-values
(fractional occupation), the sum is zero too.

You can inspect this by looking at the diagonal elements of the f
density matrix -- you find the full matrix (complex numbers) in
case.dmatup and case.dmatdn.

If you do have indeed a case where the occupation is such that you have
zero orbital moment, you might try to write an occupation in the dmat
files that gives the occupation you want, run orb, then do the scf cycle
with -orbc (to freeze the orbital potential corresponding to that
density matrix), and finally replace -orbc by -orb again -- and hope
that now it will stay in that type of solution (some examples of this
are given in the papers I refered to before).

Stefaan