Non-Collinear Spin Density Functional: Ver. 1.0
Taisuke Ozaki, RCIS, JAIST
A two component spinor wave function is defined by
where
with a spatial function
and a spin function
.
In the notes we consider non-Bloch functions, but the generalization
of the description to the Bloch function is straightforward.
Then, a density operator is given by
where
should be a step function, but it is replaced
by the Fermi function in the implementation of OpenMX.
With the definition of density operator
, a non-collinear
electron density in real space is given by
where
, and
is
a position eigenvector. The up- and down-spin densities
,
at each point are defined
by diagonalizing a matrix consisting of a non-collinear
electron densities as follows:
Based on the spinor wave function Eq. (1), the non-collinear
electron density Eq. (3), and the up- and down-spin densities,
the total energy non-collinear functional [1,2]
could be written by
where the first term is the kinetic energy, the second the electron-core
Coulomb energy, the third term the electron-electron Coulomb energy, and
the fourth term the exchange-correlation energy, respectively.
Also the total electron density
at each point is the sum of
up- and down-spin densities
,
.
Alternatively, the total energy
can be expressed in terms of
the Kohn-Sham eigenenergies
as follows:
where
is a non-collinear exchange-correlation potential
which will be discussed later on.
Considering an orthogonality relation among spinor wave functions,
let us introduce a functional
:
The variation of
with respect to the spatial wave function
is found as:
with
By setting the variation of
with respect to the spatial wave
function
to zero, and considering a unitary transformation
of
so that
can be diagonalized,
we can obtain the non-collinear Kohn-Sham equation as follows:
We see that the off-diagonal potentials produce explicitly a direct interaction
between
and
spin components in this
-
coupled
equation. The off-diagonal potentials consist of the exchange-correlation
potential
and the other contributions
such as
spin-orbit interactions.
The
-matrix in Eq. (4) which relates the non-collinear electron densities
to the up- and down-spin densities is expressed by a rotation
operator
[4]:
with Pauli matrices
where
is a unit vector along certain direction, and
a rotational angle around
.
Then, consider the following two-step rotation of a unit vector
(1,0) along the z-axis:
- First, rotate
on the y-axis
- Second, rotate
on the z-axis
The unit vector (1,0) along the z-axis is then transformed as follows:
where
Thus, if the direction of the spin is specified by the Euler
angle (
), the
-matrix in Eq. (4) is given by
the conjugate transposed matrix of Eq. (15).
The meaning of Eq. (4) becomes more clear when it is written
in a matrix form as follows:
We see that the U-matrix diagonalizes the total (average) non-collinear
spin matrix rather than the non-collinear spin matrix of each state
.
Since the exchange-correlation term is approximated by the LDA or GGA,
once the non-collinear spin matrix
is diagonalized, the diagonal up- and
down-densities are used to evaluate the exchange-correlation potentials
within LDA or GGA:
Then, the potential
is transformed to the non-collinear
exchange-correlation potential
as follows:
where
The Euler angle (
) and the up- and down-spin densities
(
,
) are determined
from the non-collinear electron densities so that the following
relation can be satisfied:
After some algebra, they are given by
Then, it is noted that the effective potential
in Eq. (11)
can be written in Pauli matrices as follows:
where
As well, the non-collinear spin density can be also written in Pauli matrices
as follows:
with
where
is a
unit matrix.
In OpenMX, the spin-orbit coupling is incorporated through
j-dependent pseudo potentials [3].
Under a spherical potential,
a couple of Dirac equations for the radial part is given by
|
|
![$\displaystyle \frac{dG_{nlj}}{dr}
+ \frac{\kappa}{r} G_{nlj}
- a
\left[
\frac{2}{a^2}
+
\varepsilon_{nlj}
-
V(r)
\right]
F_{nlj}
= 0,$](img117.png) |
(36) |
|
|
![$\displaystyle \frac{dF_{nlj}}{dr}
- \frac{\kappa}{r} F_{nlj}
+ a
\left[
\varepsilon_{nlj}
-
V(r)
\right]
G_{nlj}
= 0,$](img118.png) |
(37) |
where
and
are the majority and minority components of
the radial wave function.
(1/137.036 in a.u.).
and
for
and
, respectively.
Combining both Eqs. and eliminating
, we have
the following equation for
:
|
|
![$\displaystyle \left[
\frac{1}{2M(r)}
\left(
\frac{d^2}{dr^2}
+ \frac{a^2}{2M(r)...
...frac{\kappa(\kappa+1)}{r^2}
\right)
+ \varepsilon_{nlj}
- V
\right]
G_{nlj}
= 0$](img126.png) |
(38) |
with
By solving numerically Eq. (38) and generating j-dependent
pseudo potential
by the Troullier and Martine
(TM) scheme, we can define a general pseudopotential by
where
for
and
and for
and
The
and
are constituents of
the eigenfunction of Dirac equation.
Since
and
, the degeneracies
of
and
are
and
, respectively.
In the use of the pseudopotential defined by Eq. (40),
it is transformed to a separable form.
By introducing a local potential
which approaches
as
increases,
the j-dependent pseudo potential is divided into two contributions:
The non-local potentials
and
are non-zero within a certain radius.
Then, the pseudopotential defined by Eq. (40) is written by
The non-local part is transformed by the Blochl projector
into a separable form:
with
where
and
is an orthonormal set defined by
a norm
,
and is calculated by the following Gram-Schmidt orthogonalization:
Similarly,
with
Moreover, by unitary transforming the complex spherical harmonics functions
into the real spherical harmonics function
,
we obtain the following expressions:
 |
 |
 |
(59) |
 |
 |
 |
(60) |
 |
 |
 |
(61) |
 |
 |
 |
(62) |
 |
 |
 |
(63) |
 |
 |
 |
(64) |
 |
 |
 |
(65) |
 |
 |
 |
(66) |
with
where the real spherical harmonics functions
are
denoted by
,
, and
for
-,
-, and
-orbitals, respectively.
In conjunction with the on-site exchange term of the unrestricted
Hartree-Fock theory, the total energy of a non-collinear LDA+U method
could be defined by
with
where
is a site index,
an angular momentum quantum number,
a multiplicity number of radial basis functions, and
an
organized index of
.
is an diagonalized
occupation matrix with the size of
.
The
is the effective Coulomb electron-electron interaction energy.
Also,
is given by Eq. (5).
It should be noted that the occupation matrix is twice as size as
the collinear case. In this definition it is assumed that the exchange
interaction arises when an electron is occupied with a certain spin
direction in each localized orbital.
Considering the rotational invariance of total energy with respect
to each sub-shell
, Eq. (100) can be transformed as follows:
where
.
In this Eq. (101), although off-diagonal occupation terms in each sub-shell
are taken into account, however, those between sub-shells are neglected.
This treatment is consistent with their rotational invariant functional
by Dudarev et al.[5], and is a simple extension
of the rotational invariant functional for the case that a different
U-value is given for each basis orbital indexed with
.
In addition, the functional is rotationally invariant in the spin-space.
In this simple extension, we can not only include multiple d-orbitals
as basis set, but also can easily derive the force on atoms in a simple
form as discussed later on.
The total energy
can be expressed in terms of the
Kohn-Sham eigenenergies
as follows:
where
and
are the double counting
corrections of LDA- and U-energies, respectively.
The occupation number
(which is written by an italic font, while
the electron density, appears in Sec. 1, is denoted by a roman font)
is defined by
where, to count the occupation number
, we define three occupation number
operators given by
on-site
full
dual
where
is the dual orbital of a original
non-orthogonal basis orbital
, and is defined by
with the overlap matrix
between non-orthogonal basis orbitals.
Then, the following bi-orthogonal relation is verified:
The on-site and full occupation number operators have been
proposed by Eschrig et al. [6]
and Pickett et al. [7], respectively.
It is noted that these definitions do not satisfy a sum rule that the
trace of the occupation number matrix is equivalent to
the total number of electrons, while only the dual occupation number
operator fulfills the sum rule as follows:
where
is the density matrix defined by
with a density operator:
For three definition of occupation number operators,
on-site, full, and dual,
the occupation numbers are given by
on-site
full
dual
The derivative of the total energy Eq. (99) with respect to LCAO coefficient
is given by
with
on-site
full
dual
Substituting Eqs. (116)-(118) for the second term of Eq. (115),
we see
on-site
full
dual
Therefore, the effective projector potentials
can be expressed by
on-site
full
dual
It is clear that the effective potentials of on-site and full
are Hermitian. Also, it is verified that the effective potential
of dual is Hermitian as follows:
It should be noted that in the full and dual the
of the
site
can affect the different sites by the projector potentials
Eqs. (123) and (124) because of the overlap.
The force on atom is evaluated by
The first term can be calculated in the same way as in the collinear case.
The second term is evaluated as follows:
Considering
and
,
the first and second terms in Eq. (127) can be transformed into
derivatives of the overlap matrix. The third term in Eq. (127) means
that only the differentiation for the overlap matrix is considered,
And it is analytically differentiated, since it contains just two-center
integrals.
The LDA+U functional can possess multiple stationary points
due to the degree of freedom in the configuration space
of occupation ratio for degenerate orbitals. If electrons are
occupied with a nearly same occupancy ratio in degenerate
orbitals at the first stage of SCF steps,
the final electronic state often converges a stationary
minimum with non-orbital polarization after the SCF iteration.
Also, it is often likely that electrons are disproportionately
occupied in some of degenerate orbitals due to the exchange
interaction, which is so-called 'orbital polarization'.
As an example of the multiple minima,
we can point out a cobalt oxide (CoO) bulk in which d-orbitals
of the cobalt atom are split to
and
states, and
the five of seven d-electrons are occupied in
and
states of the majority spin, and remaining two d-electrons are
occupied in the
state of the minority spin.
Then, it depends on the initial occupancy ratios for the
states of the minority spin how the remaining two d-electrons
are occupied in three
states.
If the initial occupancy ratios are uniform, we may arrive at
the non-orbital polarized state. In fact, unless any special
treatment is considered for the initial occupancy ratios,
we see the non-orbital polarized state of the CoO bulk.
In order to explore the degree of freedom for the orbital
occupation, therefore, it is needed to develop a general method
which explicitly induces the orbital polarization.
To induce the orbital polarization, a polarized redistribution scheme
is proposed as follows:
 |
|
 |
(128) |
 |
|
 |
(129) |
 |
|
 |
|
|
|
 |
|
|
|
 |
|
|
|
 |
|
|
|
 |
(130) |
|
|
 |
(131) |
 |
|
 |
(132) |
After diagonalizing each sub-shell matrix consisting of occupation numbers,
we introduce a polarized redistribution scheme given by Eq. (130) while keeping
Eq. (129). Then, by a back transformation Eq. (130), we can obtain
a polarized occupation matrix for each sub-shell. This polarized redistribution
scheme is applied during the first few SCF steps, and then no modification
is made during subsequent SCF steps.
This proposed scheme maybe applicable to a general case:
any crystal field, any number of electrons
in the sub-shell, and any orbitals: p,d,f,...
Define
Then, the density of states,
, is given by
Also, the Mulliken populations,
, are given by
The local spin direction is determined by
The contribution to the total energy arising the Zeeman term is given by
where
The vector components of the spin magnetic moment are given by
with
where
is given by
After some alegebra we have
The vector components of the orbital magnetic moment are given by
where
.
Noting that
 |
 |
 |
(152) |
 |
 |
 |
(153) |
 |
 |
 |
(154) |
 |
 |
 |
(155) |
 |
 |
 |
(156) |
and considering a unitary transformation of the spherical harmonic functions
into a set of real harmonic functions defined by
it can be shown that
It is noted that the expectation values of
in terms of the real
harmonic functions are purely imaginary numbers.
The unitary transformation for the other
-components can be found
in a subroutine 'Set_Comp2Real()' in 'SetPara_DFT.c'. Thus, one can obtain
the matrix representation for
in terms of the real harmonic functions.
-
- 1
-
U. Von Barth and L. Hedin, J. Phys. C: Solid State Phys. 5, 1629 (1972).
- 2
-
J. Kubler, K-H. Hock, J. Sticht, and A. R. Williams,
J. Phys. F: Met. Phys. 18, 469 (1988).
- 3
-
G. Theurich and N. A. Hill, Phys. Rev. B 64, 073106 (2001).
- 4
-
J. J. Sakurai, Modern Quantum Mechanics.
- 5
-
S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and
A. P. Sutton, Phys. Rev. B 57, 1505 (1998).
- 6
-
H. Eschrig, K. Koepernik, and I. Chaplygin, J. Solid State Chem. 176, 482 (2003).
- 7
-
W. E. Pickett, SC. Erwin, E. C. Ethridge, Phy. Rev. B 58, 1201 (1998).
2007-08-14