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

$$\vert \psi_{\nu} \rangle = \vert \varphi_{\nu}^{\alpha} \alpha \rangle + \vert \varphi_{\nu}^{\beta} \beta \rangle,$$ (1)

where $\vert \varphi_{\nu}^{\alpha} \alpha \rangle \equiv \vert \varphi_{\nu}^{\alpha} \rangle \vert \alpha \rangle$ with a spatial function $\vert \varphi_{\nu}^{\alpha} \rangle$ and a spin function $\vert \alpha \rangle$. 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

$$\hat{n} = \sum_{\nu}f_{\nu} \vert \psi_{\nu} \rangle \langle \psi_{\nu} \vert,$$

$$= \sum_{\nu}f_{\nu} \left( \vert \varphi_{\nu}^{\alpha}\alpha \rang... ..._{\nu}^{\alpha}\alpha \vert + \langle \varphi_{\nu}^{\beta}\beta \vert \right),$$ (2)

where $f_{\nu}$ should be a step function, but it is replaced by the Fermi function in the implementation of OpenMX. With the definition of density operator $\hat{n}$, a non-collinear electron density in real space is given by

$$n_{\sigma \sigma'} = \langle {\bf r}\sigma\vert \hat{n} \vert {\bf r} \sigma' \rangle,$$

$$= \sum_{\nu} f_{\nu} \varphi_{\nu}^{\sigma} \varphi_{\nu}^{\sigma',*},$$ (3)

where $\sigma,\sigma'=\alpha~{\rm or}~\beta$, and $\vert {\bf r}\rangle$ is a position eigenvector. The up- and down-spin densities $n'_{\uparrow}$, $n'_{\downarrow}$ at each point are defined by diagonalizing a matrix consisting of a non-collinear electron densities as follows:

$$\left( \begin{array}{cc} n'_{\uparrow} & 0 \\ 0 & n'_{\downarrow} \end{array}\right) = U n~U^{\dag },$$

$$= U \left( \begin{array}{cc} n_{\alpha \alpha} & n_{\alpha \beta} \\ n_{\beta \alpha} & n_{\beta \beta} \\ \end{array}\right) U^{\dag }.$$ (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

$$E_{\rm tot} = \sum_{\sigma=\alpha,\beta} \sum_{\nu} f_{\nu} \langle \varphi_{\n... ...\bf r}-{\bf r'} \vert} dv dv' + E_{\rm xc} \left\{ n_{\sigma \sigma'} \right\},$$ (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 $n'$ at each point is the sum of up- and down-spin densities $n'_{\uparrow}$, $n'_{\downarrow}$. Alternatively, the total energy $E_{\rm tot}$ can be expressed in terms of the Kohn-Sham eigenenergies ${\varepsilon_{\nu}}$ as follows:

$$E_{\rm tot} = E_{\rm band} - \frac{1}{2}\int n' V_{\rm H} dv - \int {\rm Tr}(V_{\rm xc}n) dv + E_{\rm xc},$$ (6)

where $V_{\rm xc}$ 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 $F$:

$$F = E_{\rm tot} + \sum_{\nu \nu'} \epsilon_{\nu \nu'} \left( \delta_{\nu \nu'} - \langle \psi_{\nu} \vert \psi_{\nu'} \rangle \right).$$ (7)

The variation of $F$ with respect to the spatial wave function $\varphi$ is found as:

$$\frac{\delta F}{\delta \varphi_{\mu}^{\sigma,*}} = \hat{T} \varphi_{\mu}^{\sigma} + \sum_{\sigma'}w_{\sigma \sigma'}... ... \varphi_{\mu}^{\sigma'} - \sum_{\nu} \epsilon_{\mu \nu} \varphi_{\nu}^{\sigma}$$ (8)

with

$$V_{\rm H} = \int \frac{d({\bf r})}{\vert {\bf r}-{\bf r'} \vert} dv,$$ (9)

$$V_{\rm xc}^{\sigma \sigma'} = \frac{\delta E_{\rm xc}}{\delta {\rm n}_{\sigma' \sigma}}.$$ (10)

By setting the variation of $F$ with respect to the spatial wave function $\varphi$ to zero, and considering a unitary transformation of $\varphi_{\nu}^{\sigma}$ so that $\epsilon_{\mu \nu}$ can be diagonalized, we can obtain the non-collinear Kohn-Sham equation as follows:

$$\left. \begin{array}{l} \frac{\delta F}{\delta \varphi_{\mu}^{\al... ...in{array}{l} \varphi_{\mu}^{\alpha}\\ \varphi_{\mu}^{\beta} \end{array}\right) = \varepsilon_{\mu} \left( \begin{array}{l} \varphi_{\mu}^{\alpha}\\ \varphi_{\mu}^{\beta} \end{array}\right).\quad\quad$$ (11)

We see that the off-diagonal potentials produce explicitly a direct interaction between $\alpha$ and $\beta$ spin components in this $\alpha$-$\beta$ coupled equation. The off-diagonal potentials consist of the exchange-correlation potential $V_{\rm xc}$ and the other contributions $w$ such as spin-orbit interactions.

The $U$-matrix in Eq. (4) which relates the non-collinear electron densities to the up- and down-spin densities is expressed by a rotation operator $D$ [4]:

$$D \equiv \exp \left( \frac{-i{\bf\hat{\sigma}} \cdot {\bf h}\phi}{2} \right)$$ (12)

with Pauli matrices

$$\sigma_{1} = \left( \begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array}\right), ~~~~... ...~~~~ \sigma_{3} = \left( \begin{array}{cc} 1 & 0 \\ 0 & -1 \end{array}\right),$$ (13)

where ${\bf h}$ is a unit vector along certain direction, and $\phi$ a rotational angle around ${\bf h}$. Then, consider the following two-step rotation of a unit vector (1,0) along the z-axis:

• First, rotate $\theta$ on the y-axis $\to \exp\left(-i \frac{\sigma_2\theta}{2} \right)$
• Second, rotate $\phi$ on the z-axis $\to \exp\left(-i \frac{\sigma_3\phi}{2} \right)$

The unit vector (1,0) along the z-axis is then transformed as follows:

$$\left( \begin{array}{c} 1 \\ 0 \end{array}\right) \Rightarrow \exp\left(-i \frac{\sigma_3\phi}{2} \right) \exp\left(-i \frac{\sigma_2\theta}{2} \right) \left( \begin{array}{c} 1 \\ 0 \end{array}\right)$$ (14)

where

$$\exp\left(-i \frac{\sigma_3\phi}{2} \right) \exp\left(-i \frac{\sigma_2\theta}{2} \right) = \left( \begin{array}{cc} \exp(-i\frac{\phi}{2}) & 0 \\ 0 & \exp(... ...eta}{2}) \\ \sin(\frac{\theta}{2}) & \cos(\frac{\theta}{2}) \end{array}\right)$$

$$= \left( \begin{array}{cc} \exp(-i\frac{\phi}{2})\cos(\frac{\theta}... ...ac{\theta}{2}) & \exp(i\frac{\phi}{2})\cos(\frac{\theta}{2}) \end{array}\right)$$ (15)

Thus, if the direction of the spin is specified by the Euler angle ($\theta, \phi$), the $U$-matrix in Eq. (4) is given by the conjugate transposed matrix of Eq. (15).

$$U = \left( \begin{array}{cc} \exp(i\frac{\phi}{2})\cos(\frac{\theta}{... ...c{\theta}{2}) & \exp(-i\frac{\phi}{2})\cos(\frac{\theta}{2}) \end{array}\right)$$ (16)

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

$$U n~U^{\dag } = U \left\{ \sum_{\nu} f_{\nu} \left( \begin{array}{c} \varphi_{\nu... ...\nu}^{\alpha,*} & \varphi_{\nu}^{\beta,*} \end{array}\right) \right\} U^{\dag }$$ (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 $\nu$. Since the exchange-correlation term is approximated by the LDA or GGA, once the non-collinear spin matrix $n$ is diagonalized, the diagonal up- and down-densities are used to evaluate the exchange-correlation potentials $\bar{V}_{\rm xc}$ within LDA or GGA:

$$\bar{V}_{\rm xc} = \left( \begin{array}{cc} V_{\rm xc}^{\uparrow} & 0 \\ 0 & V_{\rm xc}^{\downarrow} \\ \end{array}\right),$$

$$= \frac{1}{2}(V_{\rm xc}^{\uparrow}+V_{\rm xc}^{\downarrow}) I + \frac{1}{2}(V_{\rm xc}^{\uparrow}-V_{\rm xc}^{\downarrow}) \sigma_{3},$$

$$= V_{\rm xc}^{0} I + \Delta V_{\rm xc} \sigma_{3}.$$ (18)

Then, the potential $\bar{V}_{\rm xc}$ is transformed to the non-collinear exchange-correlation potential $V_{\rm xc}$ as follows:

$$V_{\rm xc} = U^{\dag } \bar{V}_{\rm xc} U,$$

$$= V_{\rm xc}^{0} I + \Delta V_{\rm xc} U^{\dag } \sigma_{3} U,$$

$$= V_{\rm xc}^{0} I + \Delta V_{\rm xc} \bar{\sigma}_{3},$$

$$= \left( \begin{array}{cc} V_{\rm xc}^{0} + \Delta V_{\rm xc} \cos(... ...n(\theta) & V_{\rm xc}^{0} - \Delta V_{\rm xc} \cos(\theta) \end{array}\right),$$ (19)

where

$$\bar{\sigma}_{3} = \left( \begin{array}{cc} \cos(\theta) & \exp(-i\phi)\sin(\theta)\\ \exp(i\phi)\sin(\theta) & -\cos(\theta) \end{array}\right).$$ (20)

The Euler angle ($\theta, \phi$) and the up- and down-spin densities ($n'_{\uparrow}$, $n'_{\downarrow}$) are determined from the non-collinear electron densities so that the following relation can be satisfied:

$$U n U^{\dag } = \left( \begin{array}{cc} n'_{\uparrow} & 0 \\ 0 & n'_{\downarrow} \\ \end{array}\right).$$ (21)

After some algebra, they are given by

$$\phi = -\arctan\left( \frac{{\rm Im}~n_{\alpha \beta}} {{\rm Re}~n_{\alpha \beta}} \right)$$ (22)

$$\theta = \arctan\left( \frac{2({\rm Re}~n_{\alpha \beta}\cos(\phi) -{\rm Im}~n_{\alpha \beta}\sin(\phi))} {n_{\alpha \alpha}-n_{\beta \beta}} \right)$$ (23)

$$n'_{\uparrow} = \frac{1}{2}(n_{\alpha \alpha} + n_{\beta \beta}) + \frac{1}{2}(n_... ...alpha \beta}\cos(\phi) -{\rm Im}~n_{\alpha \beta}\sin(\phi) \right)\sin(\theta)$$ (24)

$$n'_{\downarrow} = \frac{1}{2}(n_{\alpha \alpha} + n_{\beta \beta}) - \frac{1}{2}(n_... ...alpha \beta}\cos(\phi) -{\rm Im}~n_{\alpha \beta}\sin(\phi) \right)\sin(\theta)$$ (25)

Then, it is noted that the effective potential $V_{\rm eff}$ in Eq. (11) can be written in Pauli matrices as follows:

$$V_{\rm eff} = V_{0} \sigma_{0} + \Delta V_{\rm xc} \bar{\sigma}_{3} + W,$$

$$= V_{0} \sigma_{0} + {\bf b}\cdot \hat{\sigma} + W,$$ (26)

where

$$V_{\rm eff}^0 = V_{\rm H} + V_{\rm xc}^{0},$$ (27)

$$W = \left( \begin{array}{cc} w_{\alpha \alpha} & w_{\alpha \beta}\\ w_{\beta \alpha} & w_{\beta \beta}\\ \end{array}\right),$$ (28)

$$b_1 = \Delta V_{\rm xc} \sin(\theta)\cos(\phi),$$ (29)

$$b_2 = \Delta V_{\rm xc} \sin(\theta)\sin(\phi),$$ (30)

$$b_3 = \Delta V_{\rm xc} \cos(\theta),$$ (31)

$$\hat{\sigma} = (\sigma_1, \sigma_2, \sigma_3).$$ (32)

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

$$n({\bf r}) = \frac{1}{2} \left( {\rm N}({\bf r})\sigma_{0} + {\bf m}({\bf r})\cdot \sigma \right)$$ (33)

with

$${\rm N}({\bf r}) = \sum_{\nu}f_{\nu} \psi_{\nu}^{\dag }({\bf r}) \psi_{\nu}({\bf r}),$$ (34)

$${\bf m}({\bf r}) = \sum_{\nu}f_{\nu} \psi_{\nu}^{\dag }({\bf r}) \hat{\sigma} \psi_{\nu}({\bf r}),$$ (35)

where $\sigma_0$ is a $2\times 2$ 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

$$\frac{dG_{nlj}}{dr} + \frac{\kappa}{r} G_{nlj} - a \left[ \frac{2}{a^2} + \varepsilon_{nlj} - V(r) \right] F_{nlj} = 0,$$ (36)

$$\frac{dF_{nlj}}{dr} - \frac{\kappa}{r} F_{nlj} + a \left[ \varepsilon_{nlj} - V(r) \right] G_{nlj} = 0,$$ (37)

where $G$ and $L$ are the majority and minority components of the radial wave function. $a\equiv 1/c$ (1/137.036 in a.u.). $\kappa=l$ and $\kappa=-(l+1)$ for $j=l-\frac{1}{2}$ and $j=l+\frac{1}{2}$, respectively. Combining both Eqs. and eliminating $F$, we have the following equation for $G$:

$$\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$$ (38)

with

$$M(r) = 1 + \frac{a^2(\varepsilon_{nlj}-V)}{2}.$$ (39)

By solving numerically Eq. (38) and generating j-dependent pseudo potential $V_{j}^{\rm ps}$ by the Troullier and Martine (TM) scheme, we can define a general pseudopotential by

$$V_{\rm ps} = \sum_{lm} \left[ \vert \Phi_{J}^{M}\rangle V^{l+\frac{1}{2}}_{\rm... ...'}^{M'}\rangle V^{l-\frac{1}{2}}_{\rm ps} \langle \Phi_{J'}^{M'} \vert \right],$$ (40)

where for $J=l+\frac{1}{2}$ and $M=m+\frac{1}{2}$

$$\vert \Phi_{J}^{M}\rangle = \left( \frac{l+m+1}{2l+1} \right)^\frac{1}{2} \vert Y_{l}^{m} \ra... ...c{l-m}{2l+1} \right)^\frac{1}{2} \vert Y_{l}^{m+1} \rangle \vert \beta \rangle,$$ (41)

and for $J'=l-\frac{1}{2}$ and $M'=m-\frac{1}{2}$

$$\vert \Phi_{J'}^{M'}\rangle = \left( \frac{l-m+1}{2l+1} \right)^\frac{1}{2} \vert Y_{l}^{m-1} \... ...rac{l+m}{2l+1} \right)^\frac{1}{2} \vert Y_{l}^{m} \rangle \vert \beta \rangle.$$ (42)

The $\Phi_{J}^{M}$ and $\Phi_{J'}^{M'}$ are constituents of the eigenfunction of Dirac equation. Since $-J \leq M \leq J$ and $-J' \leq M' \leq J'$, the degeneracies of $J$ and $J'$ are $2(l+1)$ and $2l$, respectively. In the use of the pseudopotential defined by Eq. (40), it is transformed to a separable form. By introducing a local potential $V^{\rm L}$ which approaches $-\frac{Z_{\rm eff}}{r}$ as $r$ increases, the j-dependent pseudo potential is divided into two contributions:

$$V^{l+\frac{1}{2}}_{\rm ps} = V^{l+\frac{1}{2}}_{\rm NL} + V_{\rm L},$$ (43)

$$V^{l-\frac{1}{2}}_{\rm ps} = V^{l-\frac{1}{2}}_{\rm NL} + V_{\rm L}.$$ (44)

The non-local potentials $V^{l+\frac{1}{2}}_{\rm NL}$ and $V^{l-\frac{1}{2}}_{\rm NL}$ are non-zero within a certain radius. Then, the pseudopotential defined by Eq. (40) is written by

$$V_{\rm ps} = V_{\rm L} + \sum_{lm} \left[ \vert \Phi_{J}^{M}\rangle V^{l+\frac... ...'}^{M'}\rangle V^{l-\frac{1}{2}}_{\rm NL} \langle \Phi_{J'}^{M'} \vert \right],$$ (45)

$$= V_{\rm L} + \hat{V}^{l+\frac{1}{2}}_{\rm NL} + \hat{V}^{l-\frac{1}{2}}_{\rm NL}.$$ (46)

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

$$\hat{V}^{l+\frac{1}{2}}_{\rm NL} = \sum_{lm} \vert \Phi_{J}^{M}\rangle V^{l+\frac{1}{2}}_{\rm NL} \langle \Phi_{J}^{M} \vert,$$

$$= \sum_{lm} \sum_{\zeta} \vert V^{l+\frac{1}{2}}_{\rm NL} \bar{R}_{... ...\zeta}} \langle \bar{R}_{J\zeta} \Phi_{J}^{M} V^{l+\frac{1}{2}}_{\rm NL} \vert,$$

$$= \sum_{l\zeta} \frac{1}{c_{J\zeta}} \left[ \hat{P}_{\alpha \alpha}... ...ta} + \hat{P}_{\alpha \beta}^{J\zeta} + \hat{P}_{\beta \alpha}^{J\zeta} \right]$$ (47)

with

$$\hat{P}_{\alpha \alpha}^{J\zeta} = \sum_{m=-l-1}^{l} \left( \frac{l+m+1}{2l+1} \right) \vert C_{J\zeta} Y_{l}^{m} \alpha \rangle \langle C_{J\zeta} Y_{l}^{m} \alpha \vert,$$ (48)

$$\hat{P}_{\beta \beta}^{J\zeta} = \sum_{m=-l-1}^{l} \left( \frac{l-m}{2l+1} \right) \vert C_{J\zeta} Y_{l}^{m+1} \beta \rangle \langle C_{J\zeta} Y_{l}^{m+1} \beta \vert,$$ (49)

$$\hat{P}_{\alpha \beta}^{J\zeta} = \sum_{m=-l-1}^{l} \left( \frac{l+m+1}{2l+1} \right)^{\frac{1}{2}}... ...C_{J\zeta} Y_{l}^{m} \alpha \rangle \langle C_{J\zeta} Y_{l}^{m+1} \beta \vert,$$ (50)

$$\hat{P}_{\beta \alpha}^{J\zeta} = \sum_{m=-l-1}^{l} \left( \frac{l+m+1}{2l+1} \right)^{\frac{1}{2}}... ...C_{J\zeta} Y_{l}^{m+1} \beta \rangle \langle C_{J\zeta} Y_{l}^{m} \alpha \vert,$$ (51)

where $C_{J\zeta}\equiv \bar{R}_{J\zeta} V^{l+\frac{1}{2}}_{\rm NL}$ and $\bar{R}$ is an orthonormal set defined by a norm $\int r^2 dr R V^{l+\frac{1}{2}}_{\rm NL}R'$, and is calculated by the following Gram-Schmidt orthogonalization:

$$\bar{R}_{J\zeta} = R_{J\zeta} - \sum_{\eta}^{\zeta-1} \bar{R}_{J\eta} \frac{1}{c_{J\zeta}} \int r^2 dr \bar{R}_{J\eta} V^{l+\frac{1}{2}}_{\rm NL} R_{J\eta},$$ (52)

$$c_{J\zeta} = \int r^2 dr \bar{R}_{J\zeta} V^{l+\frac{1}{2}}_{\rm NL} \bar{R}_{J\zeta},$$ (53)

Similarly,

$$\hat{V}^{l-\frac{1}{2}}_{\rm NL} = \sum_{lm} \vert \Phi_{J'}^{M'}\rangle V^{l-\frac{1}{2}}_{\rm NL} \langle \Phi_{J'}^{M'} \vert,$$

$$= \sum_{\zeta} \vert V^{l-\frac{1}{2}}_{\rm NL} \bar{R}_{J'\zeta} \... ...ta}} \langle \bar{R}_{J'\zeta} \Phi_{J'}^{M'} V^{l-\frac{1}{2}}_{\rm NL} \vert,$$

$$= \sum_{l\zeta} \frac{1}{c_{J'\zeta}} \left[ \hat{P}_{\alpha \alpha... ...} - \hat{P}_{\alpha \beta}^{J'\zeta} - \hat{P}_{\beta \alpha}^{J'\zeta} \right]$$ (54)

with

$$\hat{P}_{\alpha \alpha}^{J'\zeta} = \sum_{m=-l+1}^{l} \left( \frac{l-m+1}{2l+1} \right) \vert C_{J'\zeta} Y_{l}^{m-1} \alpha \rangle \langle C_{J'\zeta} Y_{l}^{m-1} \alpha \vert,$$ (55)

$$\hat{P}_{\beta \beta}^{J'\zeta} = \sum_{m=-l+1}^{l} \left( \frac{l+m}{2l+1} \right) \vert C_{J'\zeta} Y_{l}^{m} \beta \rangle \langle C_{J'\zeta} Y_{l}^{m} \beta \vert,$$ (56)

$$\hat{P}_{\alpha \beta}^{J'\zeta} = \sum_{m=-l+1}^{l} \left( \frac{l-m+1}{2l+1} \right)^{\frac{1}{2}}... ...{J'\zeta} Y_{l}^{m-1} \alpha \rangle \langle C_{J'\zeta} Y_{l}^{m} \beta \vert,$$ (57)

$$\hat{P}_{\beta \alpha}^{J'\zeta} = \sum_{m=-l+1}^{l} \left( \frac{l-m+1}{2l+1} \right)^{\frac{1}{2}}... ...{J'\zeta} Y_{l}^{m} \beta \rangle \langle C_{J'\zeta} Y_{l}^{m-1} \alpha \vert.$$ (58)

Moreover, by unitary transforming the complex spherical harmonics functions $Y$ into the real spherical harmonics function $\bar{Y}$, we obtain the following expressions:

$$\hat{P}_{\alpha \alpha}^{J\zeta} = \sum_{mm'} \vert C_{J\zeta} \bar{Y}_{l}^{m} \alpha \rangle F^{0}_{l,mm'} \langle C_{J\zeta} \bar{Y}_{l}^{m'} \alpha \vert,$$ (59)

$$\hat{P}_{\beta \beta}^{J\zeta} = \sum_{mm'} \vert C_{J\zeta} \bar{Y}_{l}^{m} \beta \rangle F^{1}_{l,mm'} \langle C_{J\zeta} \bar{Y}_{l}^{m'} \beta \vert,$$ (60)

$$\hat{P}_{\alpha \beta}^{J\zeta} = \sum_{mm'} \vert C_{J\zeta} \bar{Y}_{l}^{m} \alpha \rangle F^{2}_{l,mm'} \langle C_{J\zeta} \bar{Y}_{l}^{m'} \beta \vert,$$ (61)

$$\hat{P}_{\beta \alpha}^{J\zeta} = \sum_{mm'} \vert C_{J\zeta} \bar{Y}_{l}^{m} \beta \rangle F^{3}_{l,mm'} \langle C_{J\zeta} \bar{Y}_{l}^{m'} \alpha \vert,$$ (62)

$$% \hat{P}_{\alpha \alpha}^{J'\zeta} = \sum_{mm'} \vert C_{J'\zeta} \bar{Y}_{l}^{m} \alpha \rangle G^{0}_{l,mm'} \langle C_{J'\zeta} \bar{Y}_{l}^{m'} \alpha \vert,$$ (63)

$$\hat{P}_{\beta \beta}^{J'\zeta} = \sum_{mm'} \vert C_{J'\zeta} \bar{Y}_{l}^{m} \beta \rangle G^{1}_{l,mm'} \langle C_{J'\zeta} \bar{Y}_{l}^{m'} \beta \vert,$$ (64)

$$\hat{P}_{\alpha \beta}^{J'\zeta} = \sum_{mm'} \vert C_{J'\zeta} \bar{Y}_{l}^{m} \alpha \rangle G^{2}_{l,mm'} \langle C_{J'\zeta} \bar{Y}_{l}^{m'} \beta \vert,$$ (65)

$$\hat{P}_{\beta \alpha}^{J'\zeta} = \sum_{mm'} \vert C_{J'\zeta} \bar{Y}_{l}^{m} \beta \rangle G^{3}_{l,mm'} \langle C_{J'\zeta} \bar{Y}_{l}^{m'} \alpha \vert$$ (66)

with

$$F^{0}_{0} = 1,$$ (67)

$$F^{0}_{1} = \frac{2}{3}I + \frac{i}{3} \left( \begin{array}{ccc} 0 & -1 & 0\\ 1 & 0 & 0\\ 0 & 0 & 0 \end{array}\right),$$ (68)

$$F^{0}_{2} = \frac{3}{5}I + \frac{i}{5} \left( \begin{array}{ccccc} 0 & 0 & 0 ... ...2 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & -1\\ 0 & 0 & 0 & 1 & 0\\ \end{array}\right),$$ (69)

$$F^{0}_{3} = \frac{4}{7}I + \frac{i}{7} \left( \begin{array}{ccccccc} 0 & 0 & ... ... 0 & 0 & 0 & 0 & 0 & 0 & -3\\ 0 & 0 & 0 & 0 & 0 & 3 & 0\\ \end{array}\right),$$ (70)

$$F^{1}_{0} = 1,$$ (71)

$$F^{1}_{1} = (F^{0}_{1})^*,$$ (72)

$$F^{1}_{2} = (F^{0}_{2})^*,$$ (73)

$$F^{1}_{3} = (F^{0}_{3})^*,$$ (74)

$$F^{2}_{0} = 0,$$ (75)

$$F^{2}_{1} = \frac{1}{3} \left( \begin{array}{ccc} 0 & 0 & 1\\ 0 & 0 & 0\\ -... ...ft( \begin{array}{ccc} 0 & 0 & 0\\ 0 & 0 & -1\\ 0 & 1 & 0 \end{array}\right),$$ (76)

$$F^{2}_{2} = \frac{1}{5} \left( \begin{array}{ccccc} 0 & 0 & 0 & -\sqrt{3} & 0... ...1 & 0\\ 0 & 0 & 1 & 0 & 0\\ -\sqrt{3} & -1 & 0 & 0 & 0\\ \end{array}\right),$$ (77)

$$F^{2}_{3} = \frac{1}{7} \left( \begin{array}{ccccccc} 0 & -\sqrt{6} & 0 & 0 &... ...& 0 & 0 & 0\\ 0 & 0 & 0 & 0 & \sqrt{\frac{3}{2}} & 0 & 0\\ \end{array}\right)$$

$$+ \frac{i}{7} \left( \begin{array}{ccccccc} 0 & 0 & \sqrt{6} & 0 & ... ...} & 0 & 0\\ 0 & 0 & 0 & -\sqrt{\frac{3}{2}} & 0 & 0 & 0\\ \end{array}\right),$$ (78)

$$F^{3}_{0} = 0,$$ (79)

$$F^{3}_{1} = (F^{2}_{1})^{\dag },$$ (80)

$$F^{3}_{2} = (F^{2}_{2})^{\dag },$$ (81)

$$F^{3}_{3} = (F^{2}_{3})^{\dag },$$ (82)

$$G^{0}_{0} = 0,$$ (83)

$$G^{0}_{1} = \frac{1}{3}I - \frac{i}{3} \left( \begin{array}{ccc} 0 & -1 & 0\\ 1 & 0 & 0\\ 0 & 0 & 0 \end{array}\right),$$ (84)

$$G^{0}_{2} = \frac{2}{5}I - \frac{i}{5} \left( \begin{array}{ccccc} 0 & 0 & 0 ... ...2 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & -1\\ 0 & 0 & 0 & 1 & 0\\ \end{array}\right),$$ (85)

$$G^{0}_{3} = \frac{3}{7}I - \frac{i}{7} \left( \begin{array}{ccccccc} 0 & 0 & ... ... 0 & 0 & 0 & 0 & 0 & 0 & -3\\ 0 & 0 & 0 & 0 & 0 & 3 & 0\\ \end{array}\right),$$ (86)

$$G^{1}_{0} = 0,$$ (87)

$$G^{1}_{1} = (G^{0}_{1})^*,$$ (88)

$$G^{1}_{2} = (G^{0}_{2})^*,$$ (89)

$$G^{1}_{3} = (G^{0}_{3})^*,$$ (90)

$$G^{2}_{0} = 0,$$ (91)

$$G^{2}_{1} = F^{2}_{1},$$ (92)

$$G^{2}_{2} = F^{2}_{2},$$ (93)

$$G^{2}_{3} = F^{2}_{3},$$ (94)

$$G^{3}_{0} = 0,$$ (95)

$$G^{3}_{1} = F^{3}_{1},$$ (96)

$$G^{3}_{2} = F^{3}_{2},$$ (97)

$$G^{3}_{3} = F^{3}_{3},$$ (98)

where the real spherical harmonics functions $\bar{Y}$ are denoted by $(x,y,z)$, $(3z^2-r^2,x^2-y^2,xy,xz,yz)$, and $(5z^2-3r^2,5xz^2-xr^2,5yz^2-yr^2,zx^2-zy^2,xyz,x^3-3xy^2,3yx^2-y^3)$ for $p$-, $d$-, and $f$-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

$$E_{\rm LDA+U} = E_{\rm LDA} + E_{\rm U}$$ (99)

with

$$E_{\rm U} = \frac{1}{2} \sum_{i} \sum_{p} \sum_{l} U_{ipl} \left[ {\rm Tr}(N_{ipl}) - {\rm Tr}(N_{ipl} N_{ipl}) \right],$$

$$= \frac{1}{2} \sum_{s} U_{s} \left[ {\rm Tr}(N_{s}) - {\rm Tr}(N_{s} N_{s}) \right],$$ (100)

where $i$ is a site index, $l$ an angular momentum quantum number, $p$ a multiplicity number of radial basis functions, and $s$ an organized index of $(ipl)$. $N$ is an diagonalized occupation matrix with the size of $2(2l+1) \times 2(2l+1)$. The $U$ is the effective Coulomb electron-electron interaction energy. Also, $E_{\rm LDA}$ 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 $s$, Eq. (100) can be transformed as follows:

$$E_{\rm U} = \frac{1}{2} \sum_{s} U_{s} \left[ {\rm Tr}(A_{s} N_{s} A_{s}^{\dag }) - {\rm Tr}(A_{s}N_{s} A_{s}^{\dag } A_{s}N_{s} A_{s}^{\dag }) \right],$$

$$= \frac{1}{2} \sum_{s} U_{s} \left[ {\rm Tr}(n_{s}) - {\rm Tr}(n_{s} n_{s}) \right],$$

$$= \frac{1}{2} \sum_{s} U_{s} \left[ \sum_{\sigma m}n_{s,mm}^{\sigma... ...\sigma m,\sigma'm'}n_{s,mm'}^{\sigma\sigma'} n_{s,m'm}^{\sigma'\sigma} \right],$$ (101)

where $\sigma, \sigma' =\alpha~{\rm and}~\beta$. In this Eq. (101), although off-diagonal occupation terms in each sub-shell $s$ 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 $s\equiv (ipl)$. 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 $E_{\rm LDA+U}$ can be expressed in terms of the Kohn-Sham eigenenergies ${\varepsilon_{\nu}}$ as follows:

$$E_{\rm LDA+U} = E_{\rm LDA} + E_{\rm U},$$

$$= E_{\rm band} +\left[ E_{\rm ee} + E_{\rm cc} + E_{\rm xc} -\sum_{... ...\nu} \langle \psi_{\nu} \vert \hat{v}_{\rm U} \vert \psi_{\nu} \rangle \right],$$

$$= E_{\rm band} + \Delta E_{\rm LDA} + \frac{1}{2} \sum_{s} U_{s} \sum_{\sigma m, \sigma'm'} n_{s,mm'}^{\sigma \sigma'} n_{s,m'm}^{\sigma'\sigma},$$

$$= E_{\rm band} + \Delta E_{\rm LDA} + \Delta E_{\rm U},$$ (102)

where $\Delta E_{\rm LDA}$ and $\Delta E_{\rm U}$ are the double counting corrections of LDA- and U-energies, respectively.

Occupation number

The occupation number $n$ (which is written by an italic font, while the electron density, appears in Sec. 1, is denoted by a roman font) is defined by

$$n_{smm'}^{\sigma\sigma'} = \sum_{\nu} f_{\nu} \langle \psi_{\nu} \vert \hat{n}_{s,m m'}^{\sigma\sigma'} \vert \psi_{\nu} \rangle,$$ (103)

where, to count the occupation number $n$, we define three occupation number operators given by
on-site

$$\hat{n}_{smm'}^{\sigma\sigma'} = \vert \tilde{sm\sigma} \rangle \langle \tilde{sm'\sigma'} \vert,$$ (104)

full

$$\hat{n}_{smm'}^{\sigma\sigma'} = \vert sm\sigma \rangle \langle sm'\sigma'\vert,$$ (105)

dual

$$\hat{n}_{smm'}^{\sigma\sigma'} = \frac{1}{2} \left( \vert \tilde{sm\sigma} \rangle \langle sm'\sigma' \vert + \vert sm\sigma \rangle \langle \tilde{sm'\sigma'} \vert \right),$$ (106)

where $\vert \tilde{sm\sigma} \rangle$ is the dual orbital of a original non-orthogonal basis orbital $\vert sm\sigma \rangle$, and is defined by

$$\vert \tilde{sm\sigma} \rangle = \sum_{s'm'} S_{sm,s'm'}^{-1} \vert s'm'\sigma \rangle$$ (107)

with the overlap matrix $S$ between non-orthogonal basis orbitals. Then, the following bi-orthogonal relation is verified:

$$\langle \tilde{sm\sigma} \vert sm'\sigma' \rangle = \delta_{sm\sigma ,s'm'\sigma'}.$$ (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:

$${\rm Tr}(n) = \frac{1}{2} \left\{ {\rm Tr}(\rho S) + {\rm Tr}(S\rho) \right\} = N_{\rm ele},$$ (109)

where $\rho$ is the density matrix defined by

$$\rho_{sm,s'm'}^{\sigma\sigma'} = \sum_{\nu} f_{\nu} \langle \psi_{\nu} \vert \hat{\rho}_{sm,s' m'}^{\sigma\sigma'} \vert \psi_{\nu} \rangle,$$

$$= \sum_{\nu} f_{\nu} c_{\nu,sm}^{\sigma,*} c_{\nu,s'm'}^{\sigma'}$$ (110)

with a density operator:

$$\hat{\rho}_{sm,s'm'}^{\sigma\sigma'} = \vert \tilde{sm\sigma} \rangle \langle \tilde{s'm'\sigma'} \vert.$$ (111)

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

$$n_{smm'}^{\sigma\sigma'} = \rho_{sm,sm'}^{\sigma\sigma'},$$ (112)

full

$$n_{smm'}^{\sigma\sigma'} = \sum_{tn,t'n'} \rho_{tn,t'n'}^{\sigma\sigma'} S_{tn,sm}S_{sm',t'n'},$$ (113)

dual

$$n_{smm'}^{\sigma\sigma'} = \frac{1}{2} \sum_{tn} \left\{ \rho_{sm,tn}^{\sigma\sigma'} S_{tn,sm'} + S_{sm,tn} \rho_{tn,sm'}^{\sigma\sigma'} \right\},$$ (114)

Effective potential

The derivative of the total energy Eq. (99) with respect to LCAO coefficient $c_{\nu,tn}^{\sigma}$ is given by

$$\frac{\partial E_{\rm LDA+U}} {\partial c_{\nu,tn}^{\sigma,*}} = \frac{\partial E_{\rm LDA}} {\partial c_{\nu,tn}^{\sigma,*}} + \frac{\partial E_{\rm U}} {\partial c_{\nu,tn}^{\sigma,*}},$$

$$= \frac{\partial E_{\rm LDA}} {\partial c_{\nu,tn}^{\sigma,*}} + \s... ...a'}} \frac{\partial n_{smm'}^{\sigma\sigma'}} {\partial c_{\nu,tn}^{\sigma,*}},$$

$$= \frac{\partial E_{\rm LDA}} {\partial c_{\nu,tn}^{\sigma,*}} + \s... ...a'}} \frac{\partial n_{smm'}^{\sigma\sigma'}} {\partial c_{\nu,tn}^{\sigma,*}},$$

$$=$$