# Non-Collinear Spin Density Functional: Ver. 1.0

Taisuke Ozaki, RCIS, JAIST

# Non-collinear spin density functional

A two component spinor wave function is defined by

 (1)

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
 (2)

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
 (3)

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:
 (4)

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

 (5)

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:
 (6)

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 :
 (7)

The variation of with respect to the spatial wave function is found as:
 (8)

with
 (9) (10)

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:
 (11)

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]:

 (12)

with Pauli matrices
 (13)

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:
 (14)

where
 (15)

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).
 (16)

The meaning of Eq. (4) becomes more clear when it is written in a matrix form as follows:
 (17)

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:
 (18)

Then, the potential is transformed to the non-collinear exchange-correlation potential as follows:
 (19)

where
 (20)

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:
 (21)

After some algebra, they are given by
 (22) (23) (24) (25)

Then, it is noted that the effective potential in Eq. (11) can be written in Pauli matrices as follows:

 (26)

where
 (27)

 (28)

 (29) (30) (31)

 (32)

As well, the non-collinear spin density can be also written in Pauli matrices as follows:
 (33)

with
 (34) (35)

where is a unit matrix.

# Spin-orbit coupling

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

 (36) (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 :
 (38)

with
 (39)

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
 (40)

where for and
 (41)

and for and
 (42)

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:
 (43)

 (44)

The non-local potentials and are non-zero within a certain radius. Then, the pseudopotential defined by Eq. (40) is written by
 (45) (46)

The non-local part is transformed by the Blochl projector into a separable form:
 (47)

with
 (48) (49) (50) (51)

where and is an orthonormal set defined by a norm , and is calculated by the following Gram-Schmidt orthogonalization:
 (52)

 (53)

Similarly,

 (54)

with
 (55) (56) (57) (58)

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
 (67) (68) (69) (70)

 (71) (72) (73) (74)

 (75) (76) (77) (78)

 (79) (80) (81) (82)

 (83) (84) (85) (86)

 (87) (88) (89) (90)

 (91) (92) (93) (94)

 (95) (96) (97) (98)

where the real spherical harmonics functions are denoted by , , and for -, -, and -orbitals, respectively.

# Non-collinear LDA+U

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

 (99)

with
 (100)

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:
 (101)

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:

 (102)

where and are the double counting corrections of LDA- and U-energies, respectively.

## Occupation number

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

 (103)

where, to count the occupation number , we define three occupation number operators given by
on-site
 (104)

full
 (105)

dual
 (106)

where is the dual orbital of a original non-orthogonal basis orbital , and is defined by
 (107)

with the overlap matrix between non-orthogonal basis orbitals. Then, the following bi-orthogonal relation is verified:
 (108)

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:
 (109)

where is the density matrix defined by
 (110)

with a density operator:
 (111)

For three definition of occupation number operators, on-site, full, and dual, the occupation numbers are given by
on-site

 (112)

full
 (113)

dual
 (114)

## Effective potential

The derivative of the total energy Eq. (99) with respect to LCAO coefficient is given by

 (115)

with
on-site
 (116)

full
 (117)

dual
 (118)

Substituting Eqs. (116)-(118) for the second term of Eq. (115), we see
on-site

 (119)

full

 (120)

dual

 (121)

Therefore, the effective projector potentials can be expressed by

on-site

 (122)

full

 (123)

dual

 (124)

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:

 (125)

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.

## Force on atom

The force on atom is evaluated by

 (126)

The first term can be calculated in the same way as in the collinear case. The second term is evaluated as follows:
 (127)

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.

## Enhancement of orbital polarization

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,...

## Density of states

Define

 (133)

Then, the density of states, , is given by
 (134)

 (135)

Also, the Mulliken populations, , are given by
 (136)

The local spin direction is determined by
 (137)

 (138)

## Zeeman term

The contribution to the total energy arising the Zeeman term is given by

 (139)

where
 (140)

 (141)

The vector components of the spin magnetic moment are given by
 (142) (143) (144)

with
 (145)

 (146)

where is given by
 (147)

After some alegebra we have
 (148) (149) (150)

The vector components of the orbital magnetic moment are given by
 (151)

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
 (157) (158) (159)

it can be shown that
 (160)

 (161)

 (162)

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.

## Bibliography

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