Non-Collinear Spin Density Functional: Ver. 1.0

1 Non-collinear spin density functional

A two component spinor wave function is defined by

$|{\psi }_{\nu }⟩$ $=$ $|{\phi }_{\nu }^{\alpha }\alpha ⟩+|{\phi }_{\nu }^{\beta }\beta ⟩,$ (1)

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

$\widehat{n}$ $=$ ${\displaystyle \sum _{\nu }}{f}_{\nu }|{\psi }_{\nu }⟩⟨{\psi }_{\nu }|,$ (2)

$=$ ${\displaystyle \sum _{\nu }}{f}_{\nu }\left(|{\phi }_{\nu }^{\alpha }\alpha ⟩+|{\phi }_{\nu }^{\beta }\beta ⟩\right)\left(⟨{\phi }_{\nu }^{\alpha }\alpha |+⟨{\phi }_{\nu }^{\beta }\beta |\right),$

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 $\widehat{n}$, a non-collinear electron density in real space is given by

$n_{\sigma {\sigma }^{\prime }}$ $=$ $⟨𝐫\sigma |\widehat{n}|𝐫{\sigma }^{\prime }⟩,$ (3)

$=$ ${\displaystyle \sum _{\nu }}{f}_{\nu }{\phi }_{\nu }^{\sigma }{\phi }_{\nu }^{{\sigma }^{\prime },*},$

where $\sigma ,{\sigma }^{\prime }=\alpha$ $\mathrm{or}$ $\beta$, and $|𝐫⟩$ is a position eigenvector. The up- and down-spin densities $n_{\uparrow }^{\prime }$, $n_{\downarrow }^{\prime }$ at each point are defined by diagonalizing a matrix consisting of a non-collinear electron densities as follows:

$=$ $Un{U}^{†},$ (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_{\mathrm{tot}}$ $=$ ${\displaystyle \sum _{\sigma =\alpha ,\beta }}{\displaystyle \sum _{\nu }}{f}_{\nu }⟨{\phi }_{\nu }^{\sigma }|\widehat{T}|{\phi }_{\nu }^{\sigma }⟩+{\displaystyle \sum _{\sigma {\sigma }^{\prime }}}{\displaystyle \int w_{\sigma {\sigma }^{\prime }}{n}_{{\sigma }^{\prime }\sigma }}+{\displaystyle \frac{1}{2}}{\displaystyle \int \int \frac{n^{\prime }(𝐫)n^{\prime }({𝐫}^{\prime })}{|𝐫-{𝐫}^{\prime }|}𝑑v𝑑v^{\prime }}+E_{\mathrm{xc}}\left\{n_{\sigma {\sigma }^{\prime }}\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^{\prime }$ at each point is the sum of up- and down-spin densities $n_{\uparrow }^{\prime }$, $n_{\downarrow }^{\prime }$. Alternatively, the total energy $E_{\mathrm{tot}}$ can be expressed in terms of the Kohn-Sham eigenenergies ${\epsilon }_{\nu }$ as follows:

$E_{\mathrm{tot}}$ $=$ $E_{\mathrm{band}}-{\displaystyle \frac{1}{2}}{\displaystyle \int n^{\prime }{V}_{\mathrm{H}}𝑑v}-{\displaystyle \int \mathrm{Tr}(V_{\mathrm{xc}}n)𝑑v}+E_{\mathrm{xc}},$ (6)

where $V_{\mathrm{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_{\mathrm{tot}}+{\displaystyle \sum _{\nu {\nu }^{\prime }}}{ϵ}_{\nu {\nu }^{\prime }}\left({\delta }_{\nu {\nu }^{\prime }}-⟨{\psi }_{\nu }|{\psi }_{{\nu }^{\prime }}⟩\right).$ (7)

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

${\displaystyle \frac{\delta F}{\delta {\phi }_{\mu }^{\sigma ,*}}}$ $=$ $\widehat{T}{\phi }_{\mu }^{\sigma }+{\displaystyle \sum _{{\sigma }^{\prime }}}{w}_{\sigma {\sigma }^{\prime }}{\phi }_{\mu }^{{\sigma }^{\prime }}+V_{\mathrm{H}}{\phi }_{\mu }^{\sigma }+{\displaystyle \sum _{{\sigma }^{\prime }}}{V}_{\mathrm{xc}}^{\sigma {\sigma }^{\prime }}{\phi }_{\mu }^{{\sigma }^{\prime }}-{\displaystyle \sum _{\nu }}{ϵ}_{\mu \nu }{\phi }_{\nu }^{\sigma }$ (8)

with

$V_{\mathrm{H}}$ $=$ ${\displaystyle \int \frac{d(𝐫)}{|𝐫-{𝐫}^{\prime }|}𝑑v},$ (9)

$V_{\mathrm{xc}}^{\sigma {\sigma }^{\prime }}$ $=$ ${\displaystyle \frac{\delta E_{\mathrm{xc}}}{\delta {\mathrm{n}}_{{\sigma }^{\prime }\sigma }}}.$ (10)

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

$=$ ${\epsilon }_{\mu }\left(\begin{array}{c}{\phi }_{\mu }^{\alpha }\hfill \\ {\phi }_{\mu }^{\beta }\hfill \end{array}\right).$ (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_{\mathrm{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({\displaystyle \frac{-i\widehat{\sigma }\cdot 𝐡\varphi }{2}}\right)$ (12)

with Pauli matrices

${\sigma }_1$ $=$ (13)

where $𝐡$ is a unit vector along certain direction, and $\varphi$ a rotational angle around $𝐡$. 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 $\varphi$ on the z-axis $\to \exp \left(-i\frac{{\sigma }_3\varphi }{2}\right)$

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

$\left(\begin{array}{c}\hfill 1\hfill \\ \hfill 0\hfill \end{array}\right)$ $\Rightarrow$ $\exp \left(-i{\displaystyle \frac{{\sigma }_3\varphi }{2}}\right)\exp \left(-i{\displaystyle \frac{{\sigma }_2\theta }{2}}\right)\left(\begin{array}{c}\hfill 1\hfill \\ \hfill 0\hfill \end{array}\right),$ (14)

where

$\exp \left(-i{\displaystyle \frac{{\sigma }_3\varphi }{2}}\right)\exp \left(-i{\displaystyle \frac{{\sigma }_2\theta }{2}}\right)$ $=$ (15)

$=$

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

$U$ $=$ (16)

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

$Un{U}^{†}$ $=$ (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 ${\overline{V}}_{\mathrm{xc}}$ within LDA or GGA:

${\overline{V}}_{\mathrm{xc}}$ $=$ (18)

$=$ ${\displaystyle \frac{1}{2}}(V_{\mathrm{xc}}^{\uparrow }+V_{\mathrm{xc}}^{\downarrow })I+{\displaystyle \frac{1}{2}}(V_{\mathrm{xc}}^{\uparrow }-V_{\mathrm{xc}}^{\downarrow }){\sigma }_3,$

$=$ $V_{\mathrm{xc}}^{0}I+\mathrm{\Delta }{V}_{\mathrm{xc}}{\sigma }_3.$

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

$V_{\mathrm{xc}}$ $=$ $U^{†}{\overline{V}}_{\mathrm{xc}}U,$ (19)

$=$ $V_{\mathrm{xc}}^{0}I+\mathrm{\Delta }{V}_{\mathrm{xc}}{U}^{†}{\sigma }_{3}U,$

$=$ $V_{\mathrm{xc}}^{0}I+\mathrm{\Delta }{V}_{\mathrm{xc}}{\overline{\sigma }}_3,$

$=$

where

${\overline{\sigma }}_3$ $=$ (20)

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

$Un{U}^{†}$ $=$ (21)

After some algebra, they are given by

$\varphi$ $=$ $-\arctan \left({\displaystyle \frac{\mathrm{Im}{n}_{\alpha \beta }}{\mathrm{Re}{n}_{\alpha \beta }}}\right).$ (22)

$\theta$ $=$ $\arctan \left({\displaystyle \frac{2(\mathrm{Re}{n}_{\alpha \beta }\cos (\varphi )-\mathrm{Im}{n}_{\alpha \beta }\sin (\varphi ))}{n_{\alpha \alpha }-n_{\beta \beta }}}\right).$ (23)

$n_{\uparrow }^{\prime }$ $=$ ${\displaystyle \frac{1}{2}}(n_{\alpha \alpha }+n_{\beta \beta })+{\displaystyle \frac{1}{2}}(n_{\alpha \alpha }-n_{\beta \beta })\cos (\theta )+\left(\mathrm{Re}{n}_{\alpha \beta }\cos (\varphi )-\mathrm{Im}{n}_{\alpha \beta }\sin (\varphi )\right)\sin (\theta ).$ (24)

$n_{\downarrow }^{\prime }$ $=$ ${\displaystyle \frac{1}{2}}(n_{\alpha \alpha }+n_{\beta \beta })-{\displaystyle \frac{1}{2}}(n_{\alpha \alpha }-n_{\beta \beta })\cos (\theta )-\left(\mathrm{Re}{n}_{\alpha \beta }\cos (\varphi )-\mathrm{Im}{n}_{\alpha \beta }\sin (\varphi )\right)\sin (\theta ).$ (25)

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

$V_{\mathrm{eff}}$ $=$ $V_0{\sigma }_0+\mathrm{\Delta }{V}_{\mathrm{xc}}{\overline{\sigma }}_3+W,$ (26)

$=$ $V_0{\sigma }_0+𝐛\cdot \widehat{\sigma }+W,$

where

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

$W$ $=$ (28)

$b_1$ $=$ $\mathrm{\Delta }{V}_{\mathrm{xc}}\sin (\theta )\cos (\varphi ),$ (29)

$b_2$ $=$ $\mathrm{\Delta }{V}_{\mathrm{xc}}\sin (\theta )\sin (\varphi ),$ (30)

$b_3$ $=$ $\mathrm{\Delta }{V}_{\mathrm{xc}}\cos (\theta ),$ (31)

$\widehat{\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(𝐫)$ $=$ ${\displaystyle \frac{1}{2}}\left(\mathrm{N}(𝐫){\sigma }_0+𝐦(𝐫)\cdot \sigma \right)$ (33)

with

$\mathrm{N}(𝐫)$ $=$ ${\displaystyle \sum _{\nu }}{f}_{\nu }{\psi }_{\nu }^{†}(𝐫){\psi }_{\nu }(𝐫),$ (34)

$𝐦(𝐫)$ $=$ ${\displaystyle \sum _{\nu }}{f}_{\nu }{\psi }_{\nu }^{†}(𝐫)\widehat{\sigma }{\psi }_{\nu }(𝐫),$ (35)

where ${\sigma }_0$ is a $2\times 2$ unit matrix.

2 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

${\displaystyle \frac{d{G}_{nlj}}{dr}}+{\displaystyle \frac{\kappa }{r}}{G}_{nlj}-a\left[{\displaystyle \frac{2}{a^2}}+{\epsilon }_{nlj}-V(r)\right]F_{nlj}=0,$ (36)

${\displaystyle \frac{d{F}_{nlj}}{dr}}-{\displaystyle \frac{\kappa }{r}}{F}_{nlj}+a\left[{\epsilon }_{nlj}-V(r)\right]G_{nlj}=0,$ (37)

where $G$ and $F$ 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[{\displaystyle \frac{1}{2M(r)}}\left({\displaystyle \frac{d^2}{d{r}^2}}+{\displaystyle \frac{a^2}{2M(r)}}{\displaystyle \frac{dV}{dr}}{\displaystyle \frac{d}{dr}}+{\displaystyle \frac{a^2}{2M(r)}}{\displaystyle \frac{\kappa }{r}}{\displaystyle \frac{dV}{dr}}-{\displaystyle \frac{\kappa (\kappa +1)}{r^2}}\right)+{\epsilon }_{nlj}-V\right]G_{nlj}=0$ (38)

with

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

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

$V_{\mathrm{ps}}$ $=$ ${\displaystyle \sum _{lm}}\left[|{\mathrm{\Phi }}_J^M⟩V_{\mathrm{ps}}^{l+\frac{1}{2}}⟨{\mathrm{\Phi }}_J^M|+|{\mathrm{\Phi }}_{J^{\prime }}^{M^{\prime }}⟩V_{\mathrm{ps}}^{l-\frac{1}{2}}⟨{\mathrm{\Phi }}_{J^{\prime }}^{M^{\prime }}|\right],$ (40)

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

$|{\mathrm{\Phi }}_J^M⟩$ $=$ ${\left({\displaystyle \frac{l+m+1}{2l+1}}\right)}^{\frac{1}{2}}|Y_l^m⟩|\alpha ⟩+{\left({\displaystyle \frac{l-m}{2l+1}}\right)}^{\frac{1}{2}}|Y_l^{m+1}⟩|\beta ⟩,$ (41)

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

$|{\mathrm{\Phi }}_{J^{\prime }}^{M^{\prime }}⟩$ $=$ ${\left({\displaystyle \frac{l-m+1}{2l+1}}\right)}^{\frac{1}{2}}|Y_l^{m-1}⟩|\alpha ⟩-{\left({\displaystyle \frac{l+m}{2l+1}}\right)}^{\frac{1}{2}}|Y_l^m⟩|\beta ⟩.$ (42)

The ${\mathrm{\Phi }}_J^M$ and ${\mathrm{\Phi }}_{J^{\prime }}^{M^{\prime }}$ are constituents of the eigenfunction of Dirac equation. Since $-J\le M\le J$ and $-J^{\prime }\le M^{\prime }\le J^{\prime }$, the degeneracies of $J$ and $J^{\prime }$ 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^{\mathrm{L}}$ which approaches $-\frac{Z_{\mathrm{eff}}}{r}$ as $r$ increases, the j-dependent pseudo potential is divided into two contributions:

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

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

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

$V_{\mathrm{ps}}$ $=$ $V_{\mathrm{L}}+{\displaystyle \sum _{lm}}\left[|{\mathrm{\Phi }}_J^M⟩V_{\mathrm{NL}}^{l+\frac{1}{2}}⟨{\mathrm{\Phi }}_J^M|+|{\mathrm{\Phi }}_{J^{\prime }}^{M^{\prime }}⟩V_{\mathrm{NL}}^{l-\frac{1}{2}}⟨{\mathrm{\Phi }}_{J^{\prime }}^{M^{\prime }}|\right],$ (45)

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

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

${\widehat{V}}_{\mathrm{NL}}^{l+\frac{1}{2}}$ $=$ ${\displaystyle \sum _{lm}}|{\mathrm{\Phi }}_J^M⟩V_{\mathrm{NL}}^{l+\frac{1}{2}}⟨{\mathrm{\Phi }}_J^M|,$ (47)

$=$ ${\displaystyle \sum _{lm}}{\displaystyle \sum _{\zeta }}|V_{\mathrm{NL}}^{l+\frac{1}{2}}{\overline{R}}_{J\zeta }{\mathrm{\Phi }}_J^M⟩{\displaystyle \frac{1}{c_{J\zeta }}}⟨{\overline{R}}_{J\zeta }{\mathrm{\Phi }}_J^M{V}_{\mathrm{NL}}^{l+\frac{1}{2}}|,$

$=$ ${\displaystyle \sum _{l\zeta }}{\displaystyle \frac{1}{c_{J\zeta }}}\left[{\widehat{P}}_{\alpha \alpha }^{J\zeta }+{\widehat{P}}_{\beta \beta }^{J\zeta }+{\widehat{P}}_{\alpha \beta }^{J\zeta }+{\widehat{P}}_{\beta \alpha }^{J\zeta }\right]$

with

${\widehat{P}}_{\alpha \alpha }^{J\zeta }$ $=$ ${\displaystyle \sum _{m=-l-1}^l}\left({\displaystyle \frac{l+m+1}{2l+1}}\right)|C_{J\zeta }{Y}_l^m\alpha ⟩⟨C_{J\zeta }{Y}_l^m\alpha |,$ (48)

${\widehat{P}}_{\beta \beta }^{J\zeta }$ $=$ ${\displaystyle \sum _{m=-l-1}^l}\left({\displaystyle \frac{l-m}{2l+1}}\right)|C_{J\zeta }{Y}_l^{m+1}\beta ⟩⟨C_{J\zeta }{Y}_l^{m+1}\beta |,$ (49)

${\widehat{P}}_{\alpha \beta }^{J\zeta }$ $=$ ${\displaystyle \sum _{m=-l-1}^l}{\left({\displaystyle \frac{l+m+1}{2l+1}}\right)}^{\frac{1}{2}}{\left({\displaystyle \frac{l-m}{2l+1}}\right)}^{\frac{1}{2}}|C_{J\zeta }{Y}_l^m\alpha ⟩⟨C_{J\zeta }{Y}_l^{m+1}\beta |,$ (50)

${\widehat{P}}_{\beta \alpha }^{J\zeta }$ $=$ ${\displaystyle \sum _{m=-l-1}^l}{\left({\displaystyle \frac{l+m+1}{2l+1}}\right)}^{\frac{1}{2}}{\left({\displaystyle \frac{l-m}{2l+1}}\right)}^{\frac{1}{2}}|C_{J\zeta }{Y}_l^{m+1}\beta ⟩⟨C_{J\zeta }{Y}_l^m\alpha |,$ (51)

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

${\overline{R}}_{J\zeta }$ $=$ $R_{J\zeta }-{\displaystyle \sum _{\eta }^{\zeta -1}}{\overline{R}}_{J\eta }{\displaystyle \frac{1}{c_{J\zeta }}}{\displaystyle \int r^2𝑑r{\overline{R}}_{J\eta }{V}_{\mathrm{NL}}^{l+\frac{1}{2}}{R}_{J\eta }},$ (52)

$c_{J\zeta }$ $=$ ${\displaystyle \int r^2𝑑r{\overline{R}}_{J\zeta }{V}_{\mathrm{NL}}^{l+\frac{1}{2}}{\overline{R}}_{J\zeta }},$ (53)

Similarly,

${\widehat{V}}_{\mathrm{NL}}^{l-\frac{1}{2}}$ $=$ ${\displaystyle \sum _{lm}}|{\mathrm{\Phi }}_{J^{\prime }}^{M^{\prime }}⟩V_{\mathrm{NL}}^{l-\frac{1}{2}}⟨{\mathrm{\Phi }}_{J^{\prime }}^{M^{\prime }}|,$ (54)

$=$ ${\displaystyle \sum _{\zeta }}|V_{\mathrm{NL}}^{l-\frac{1}{2}}{\overline{R}}_{J^{\prime }\zeta }{\mathrm{\Phi }}_{J^{\prime }}^{M^{\prime }}⟩{\displaystyle \frac{1}{c_{J^{\prime }\zeta }}}⟨{\overline{R}}_{J^{\prime }\zeta }{\mathrm{\Phi }}_{J^{\prime }}^{M^{\prime }}{V}_{\mathrm{NL}}^{l-\frac{1}{2}}|,$

$=$ ${\displaystyle \sum _{l\zeta }}{\displaystyle \frac{1}{c_{J^{\prime }\zeta }}}\left[{\widehat{P}}_{\alpha \alpha }^{J^{\prime }\zeta }+{\widehat{P}}_{\beta \beta }^{J^{\prime }\zeta }-{\widehat{P}}_{\alpha \beta }^{J^{\prime }\zeta }-{\widehat{P}}_{\beta \alpha }^{J^{\prime }\zeta }\right]$

with

${\widehat{P}}_{\alpha \alpha }^{J^{\prime }\zeta }$ $=$ ${\displaystyle \sum _{m=-l+1}^l}\left({\displaystyle \frac{l-m+1}{2l+1}}\right)|C_{J^{\prime }\zeta }{Y}_l^{m-1}\alpha ⟩⟨C_{J^{\prime }\zeta }{Y}_l^{m-1}\alpha |,$ (55)

${\widehat{P}}_{\beta \beta }^{J^{\prime }\zeta }$ $=$ ${\displaystyle \sum _{m=-l+1}^l}\left({\displaystyle \frac{l+m}{2l+1}}\right)|C_{J^{\prime }\zeta }{Y}_l^m\beta ⟩⟨C_{J^{\prime }\zeta }{Y}_l^m\beta |,$ (56)

${\widehat{P}}_{\alpha \beta }^{J^{\prime }\zeta }$ $=$ ${\displaystyle \sum _{m=-l+1}^l}{\left({\displaystyle \frac{l-m+1}{2l+1}}\right)}^{\frac{1}{2}}{\left({\displaystyle \frac{l+m}{2l+1}}\right)}^{\frac{1}{2}}|C_{J^{\prime }\zeta }{Y}_l^{m-1}\alpha ⟩⟨C_{J^{\prime }\zeta }{Y}_l^m\beta |,$ (57)

${\widehat{P}}_{\beta \alpha }^{J^{\prime }\zeta }$ $=$ ${\displaystyle \sum _{m=-l+1}^l}{\left({\displaystyle \frac{l-m+1}{2l+1}}\right)}^{\frac{1}{2}}{\left({\displaystyle \frac{l+m}{2l+1}}\right)}^{\frac{1}{2}}|C_{J^{\prime }\zeta }{Y}_l^m\beta ⟩⟨C_{J^{\prime }\zeta }{Y}_l^{m-1}\alpha |.$ (58)

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

${\widehat{P}}_{\alpha \alpha }^{J\zeta }$ $=$ ${\displaystyle \sum _{m{m}^{\prime }}}|C_{J\zeta }{\overline{Y}}_l^m\alpha ⟩F_{l,m{m}^{\prime }}^0⟨C_{J\zeta }{\overline{Y}}_l^{m^{\prime }}\alpha |,$ (59)

${\widehat{P}}_{\beta \beta }^{J\zeta }$ $=$ ${\displaystyle \sum _{m{m}^{\prime }}}|C_{J\zeta }{\overline{Y}}_l^m\beta ⟩F_{l,m{m}^{\prime }}^1⟨C_{J\zeta }{\overline{Y}}_l^{m^{\prime }}\beta |,$ (60)

${\widehat{P}}_{\alpha \beta }^{J\zeta }$ $=$ ${\displaystyle \sum _{m{m}^{\prime }}}|C_{J\zeta }{\overline{Y}}_l^m\alpha ⟩F_{l,m{m}^{\prime }}^2⟨C_{J\zeta }{\overline{Y}}_l^{m^{\prime }}\beta |,$ (61)

${\widehat{P}}_{\beta \alpha }^{J\zeta }$ $=$ ${\displaystyle \sum _{m{m}^{\prime }}}|C_{J\zeta }{\overline{Y}}_l^m\beta ⟩F_{l,m{m}^{\prime }}^3⟨C_{J\zeta }{\overline{Y}}_l^{m^{\prime }}\alpha |,$ (62)

${\widehat{P}}_{\alpha \alpha }^{J^{\prime }\zeta }$ $=$ ${\displaystyle \sum _{m{m}^{\prime }}}|C_{J^{\prime }\zeta }{\overline{Y}}_l^m\alpha ⟩G_{l,m{m}^{\prime }}^0⟨C_{J^{\prime }\zeta }{\overline{Y}}_l^{m^{\prime }}\alpha |,$ (63)

${\widehat{P}}_{\beta \beta }^{J^{\prime }\zeta }$ $=$ ${\displaystyle \sum _{m{m}^{\prime }}}|C_{J^{\prime }\zeta }{\overline{Y}}_l^m\beta ⟩G_{l,m{m}^{\prime }}^1⟨C_{J^{\prime }\zeta }{\overline{Y}}_l^{m^{\prime }}\beta |,$ (64)

${\widehat{P}}_{\alpha \beta }^{J^{\prime }\zeta }$ $=$ ${\displaystyle \sum _{m{m}^{\prime }}}|C_{J^{\prime }\zeta }{\overline{Y}}_l^m\alpha ⟩G_{l,m{m}^{\prime }}^2⟨C_{J^{\prime }\zeta }{\overline{Y}}_l^{m^{\prime }}\beta |,$ (65)

${\widehat{P}}_{\beta \alpha }^{J^{\prime }\zeta }$ $=$ ${\displaystyle \sum _{m{m}^{\prime }}}|C_{J^{\prime }\zeta }{\overline{Y}}_l^m\beta ⟩G_{l,m{m}^{\prime }}^3⟨C_{J^{\prime }\zeta }{\overline{Y}}_l^{m^{\prime }}\alpha |$ (66)

with

$F_0^0$ $=$ $1,$ (67)

$F_1^0$ $=$ (68)

$F_2^0$ $=$ (69)

$F_3^0$ $=$ (70)

$F_0^1=1,$ (71)

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

$F_2^1={(F_2^0)}^{*},$ (73)

$F_3^1={(F_3^0)}^{*},$ (74)

$F_0^2$ $=$ $0,$ (75)

$F_1^2$ $=$ (76)

$F_2^2$ $=$ (77)

$F_3^2$ $=$ (78)

$+$

$F_0^3=0,$ (79)

$F_1^3={(F_1^2)}^{†},$ (80)

$F_2^3={(F_2^2)}^{†},$ (81)

$F_3^3={(F_3^2)}^{†},$ (82)

$G_0^0$ $=$ $0,$ (83)

$G_1^0$ $=$ (84)

$G_2^0$ $=$ (85)

$G_3^0$ $=$ (86)

$G_0^1=0,$ (87)

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

$G_2^1={(G_2^0)}^{*},$ (89)

$G_3^1={(G_3^0)}^{*},$ (90)

$G_0^2=0,$ (91)

$G_1^2=F_1^2,$ (92)

$G_2^2=F_2^2,$ (93)

$G_3^2=F_3^2,$ (94)

$G_0^3=0,$ (95)

$G_1^3=F_1^3,$ (96)

$G_2^3=F_2^3,$ (97)

$G_3^3=F_3^3,$ (98)

where the real spherical harmonics functions $\overline{Y}$ are denoted by $(x,y,z)$, $(3z^2-r^2,x^2-y^2,xy,xz,yz)$, and $(5z^2-3r^2,5x{z}^2-x{r}^2,5y{z}^2-y{r}^2,z{x}^2-z{y}^2,xyz,x^3-3x{y}^2,3y{x}^2-y^3)$ for $p$-, $d$-, and $f$-orbitals, respectively.

3 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_{\mathrm{LDA}+\mathrm{U}}$ $=$ $E_{\mathrm{LDA}}+E_{\mathrm{U}}$ (99)

with

$E_{\mathrm{U}}$ $=$ ${\displaystyle \frac{1}{2}}{\displaystyle \sum _i}{\displaystyle \sum _p}{\displaystyle \sum _l}{U}_{ipl}\left[\mathrm{Tr}(N_{ipl})-\mathrm{Tr}(N_{ipl}{N}_{ipl})\right],$ (100)

$=$ ${\displaystyle \frac{1}{2}}{\displaystyle \sum _s}{U}_s\left[\mathrm{Tr}(N_s)-\mathrm{Tr}(N_s{N}_s)\right],$

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_{\mathrm{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_{\mathrm{U}}$ $=$ ${\displaystyle \frac{1}{2}}{\displaystyle \sum _s}{U}_s\left[\mathrm{Tr}(A_s{N}_s{A}_s^{†})-\mathrm{Tr}(A_s{N}_s{A}_s^{†}{A}_s{N}_s{A}_s^{†})\right],$ (101)

$=$ ${\displaystyle \frac{1}{2}}{\displaystyle \sum _s}{U}_s\left[\mathrm{Tr}(n_s)-\mathrm{Tr}(n_s{n}_s)\right],$

$=$ ${\displaystyle \frac{1}{2}}{\displaystyle \sum _s}{U}_s\left[{\displaystyle \sum _{\sigma m}}{n}_{s,mm}^{\sigma \sigma }-{\displaystyle \sum _{\sigma m,{\sigma }^{\prime }{m}^{\prime }}}{n}_{s,m{m}^{\prime }}^{\sigma {\sigma }^{\prime }}{n}_{s,m^{\prime }m}^{{\sigma }^{\prime }\sigma }\right],$

where $\sigma ,{\sigma }^{\prime }=\alpha$ $\mathrm{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_{\mathrm{LDA}+\mathrm{U}}$ can be expressed in terms of the Kohn-Sham eigenenergies ${\epsilon }_{\nu }$ as follows:

$E_{\mathrm{LDA}+\mathrm{U}}$ $=$ $E_{\mathrm{LDA}}+E_{\mathrm{U}},$ (102)

$=$ $E_{\mathrm{band}}+\left[E_{\mathrm{ee}}+E_{\mathrm{cc}}+E_{\mathrm{xc}}-{\displaystyle \sum _{\nu }}⟨{\psi }_{\nu }|{\widehat{v}}_{\mathrm{LDA}}|{\psi }_{\nu }⟩\right]+\left[E_{\mathrm{U}}-{\displaystyle \sum _{\nu }}⟨{\psi }_{\nu }|{\widehat{v}}_{\mathrm{U}}|{\psi }_{\nu }⟩\right],$

$=$ $E_{\mathrm{band}}+\mathrm{\Delta }{E}_{\mathrm{LDA}}+{\displaystyle \frac{1}{2}}{\displaystyle \sum _s}{U}_s{\displaystyle \sum _{\sigma m,{\sigma }^{\prime }{m}^{\prime }}}{n}_{s,m{m}^{\prime }}^{\sigma {\sigma }^{\prime }}{n}_{s,m^{\prime }m}^{{\sigma }^{\prime }\sigma },$

$=$ $E_{\mathrm{band}}+\mathrm{\Delta }{E}_{\mathrm{LDA}}+\mathrm{\Delta }{E}_{\mathrm{U}},$

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

3.1 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_{sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}$ $=$ ${\displaystyle \sum _{\nu }}{f}_{\nu }⟨{\psi }_{\nu }|{\widehat{n}}_{s,m{m}^{\prime }}^{\sigma {\sigma }^{\prime }}|{\psi }_{\nu }⟩,$ (103)

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

${\widehat{n}}_{sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}$ $=$ $|\widetilde{sm\sigma }⟩⟨\widetilde{s{m}^{\prime }{\sigma }^{\prime }}|,$ (104)

full

${\widehat{n}}_{sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}$ $=$ $|sm\sigma ⟩⟨s{m}^{\prime }{\sigma }^{\prime }|,$ (105)

dual

${\widehat{n}}_{sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}$ $=$ ${\displaystyle \frac{1}{2}}\left(|\widetilde{sm\sigma }⟩⟨s{m}^{\prime }{\sigma }^{\prime }|+|sm\sigma ⟩⟨\widetilde{s{m}^{\prime }{\sigma }^{\prime }}|\right),$ (106)

where $|\widetilde{sm\sigma }⟩$ is the dual orbital of a original non-orthogonal basis orbital $|sm\sigma ⟩$, and is defined by

$|\widetilde{sm\sigma }⟩$ $=$ ${\displaystyle \sum _{s^{\prime }{m}^{\prime }}}{S}_{sm,s^{\prime }{m}^{\prime }}^{-1}|s^{\prime }{m}^{\prime }\sigma ⟩$ (107)

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

$⟨\widetilde{sm\sigma }|s{m}^{\prime }{\sigma }^{\prime }⟩$ $=$ ${\delta }_{sm\sigma ,s^{\prime }{m}^{\prime }{\sigma }^{\prime }}.$ (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:

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

where $\rho$ is the density matrix defined by

${\rho }_{sm,s^{\prime }{m}^{\prime }}^{\sigma {\sigma }^{\prime }}$ $=$ ${\displaystyle \sum _{\nu }}{f}_{\nu }⟨{\psi }_{\nu }|{\widehat{\rho }}_{sm,s^{\prime }{m}^{\prime }}^{\sigma {\sigma }^{\prime }}|{\psi }_{\nu }⟩,$ (110)

$=$ ${\displaystyle \sum _{\nu }}{f}_{\nu }{c}_{\nu ,sm}^{\sigma ,*}{c}_{\nu ,s^{\prime }{m}^{\prime }}^{{\sigma }^{\prime }}$

with a density operator:

${\widehat{\rho }}_{sm,s^{\prime }{m}^{\prime }}^{\sigma {\sigma }^{\prime }}$ $=$ $|\widetilde{sm\sigma }⟩⟨\widetilde{s^{\prime }{m}^{\prime }{\sigma }^{\prime }}|.$ (111)

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

$n_{sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}$ $=$ ${\rho }_{sm,s{m}^{\prime }}^{\sigma {\sigma }^{\prime }},$ (112)

full

$n_{sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}$ $=$ ${\displaystyle \sum _{tn,t^{\prime }{n}^{\prime }}}{\rho }_{tn,t^{\prime }{n}^{\prime }}^{\sigma {\sigma }^{\prime }}{S}_{tn,sm}{S}_{s{m}^{\prime },t^{\prime }{n}^{\prime }},$ (113)

dual

$n_{sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}$ $=$ ${\displaystyle \frac{1}{2}}{\displaystyle \sum _{tn}}\left\{{\rho }_{sm,tn}^{\sigma {\sigma }^{\prime }}{S}_{tn,s{m}^{\prime }}+S_{sm,tn}{\rho }_{tn,s{m}^{\prime }}^{\sigma {\sigma }^{\prime }}\right\},$ (114)

3.2 Effective potential

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

${\displaystyle \frac{\partial E_{\mathrm{LDA}+\mathrm{U}}}{\partial c_{\nu ,tn}^{\sigma ,*}}}$ $=$ ${\displaystyle \frac{\partial E_{\mathrm{LDA}}}{\partial c_{\nu ,tn}^{\sigma ,*}}}+{\displaystyle \frac{\partial E_{\mathrm{U}}}{\partial c_{\nu ,tn}^{\sigma ,*}}},$ (115)

$=$ ${\displaystyle \frac{\partial E_{\mathrm{LDA}}}{\partial c_{\nu ,tn}^{\sigma ,*}}}+{\displaystyle \sum _{{\sigma }^{\prime }}}{\displaystyle \sum _{sm{m}^{\prime }}}{\displaystyle \frac{\partial E_{\mathrm{U}}}{\partial n_{sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}}}{\displaystyle \frac{\partial n_{sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}}{\partial c_{\nu ,tn}^{\sigma ,*}}},$

$=$ ${\displaystyle \frac{\partial E_{\mathrm{LDA}}}{\partial c_{\nu ,tn}^{\sigma ,*}}}+{\displaystyle \sum _{{\sigma }^{\prime }}}{\displaystyle \sum _{sm{m}^{\prime }}}{U}_s({\displaystyle \frac{1}{2}}{\delta }_{\sigma {\sigma }^{\prime }}{\delta }_{m{m}^{\prime }}-n_{s{m}^{\prime }m}^{{\sigma }^{\prime }\sigma }){\displaystyle \frac{\partial n_{sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}}{\partial c_{\nu ,tn}^{\sigma ,*}}},$

$=$ ${\displaystyle \frac{\partial E_{\mathrm{LDA}}}{\partial c_{\nu ,tn}^{\sigma ,*}}}+{\displaystyle \sum _{{\sigma }^{\prime }}}{\displaystyle \sum _{sm{m}^{\prime }}}{v}_{\mathrm{U},sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}{\displaystyle \frac{\partial n_{sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}}{\partial c_{\nu ,tn}^{\sigma ,*}}}$

with
on-site

${\displaystyle \frac{\partial n_{sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}}{\partial c_{\nu ,tn}^{\sigma ,*}}}$ $=$ ${\delta }_{st}{\delta }_{mn}{c}_{\nu ,s{m}^{\prime }}^{{\sigma }^{\prime }},$ (116)

full

${\displaystyle \frac{\partial n_{sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}}{\partial c_{\nu ,tn}^{\sigma ,*}}}$ $=$ ${\displaystyle \sum _{t^{\prime }{n}^{\prime }}}{S}_{tn,sm}{S}_{s{m}^{\prime },t^{\prime }{n}^{\prime }}{c}_{\nu ,t^{\prime }{n}^{\prime }}^{{\sigma }^{\prime }},$ (117)

dual

${\displaystyle \frac{\partial n_{sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}}{\partial c_{\nu ,tn}^{\sigma ,*}}}$ $=$ ${\displaystyle \frac{1}{2}}\left\{{\delta }_{st}{\delta }_{mn}{\displaystyle \sum _{t^{\prime }{n}^{\prime }}}{c}_{\nu ,t^{\prime }{n}^{\prime }}^{{\sigma }^{\prime }}{S}_{t^{\prime }{n}^{\prime },s{m}^{\prime }}+S_{sm,tn}{c}_{\nu ,s{m}^{\prime }}^{{\sigma }^{\prime }}\right\}.$ (118)

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

${\displaystyle \sum _{sm{m}^{\prime }}}{v}_{\mathrm{U},sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}{\displaystyle \frac{\partial n_{sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}}{\partial c_{\nu ,tn}^{\sigma ,*}}}$ $=$ ${\displaystyle \sum _{{\sigma }^{\prime \prime \prime }}}{\displaystyle \sum _{t^{\prime }{n}^{\prime }}}⟨tn\sigma |\left[{\displaystyle \sum _{{\sigma }^{\prime }{\sigma }^{\prime \prime }}}{\displaystyle \sum _{sm{m}^{\prime }}}|\widetilde{sm{\sigma }^{\prime }}⟩v_{\mathrm{U},sm{m}^{\prime }}^{{\sigma }^{\prime }{\sigma }^{\prime \prime }}⟨\widetilde{s{m}^{\prime }{\sigma }^{\prime \prime }}|\right]|t^{\prime }{n}^{\prime }{\sigma }^{\prime \prime \prime }⟩c_{\nu ,t^{\prime }{n}^{\prime }}^{{\sigma }^{\prime \prime \prime }},$ (119)

full

${\displaystyle \sum _{sm{m}^{\prime }}}{v}_{\mathrm{U},sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}{\displaystyle \frac{\partial n_{sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}}{\partial c_{\nu ,tn}^{\sigma ,*}}}$ $=$ ${\displaystyle \sum _{{\sigma }^{\prime \prime \prime }}}{\displaystyle \sum _{t^{\prime }{n}^{\prime }}}⟨tn\sigma |\left[{\displaystyle \sum _{{\sigma }^{\prime }{\sigma }^{\prime \prime }}}{\displaystyle \sum _{sm{m}^{\prime }}}|sm{\sigma }^{\prime }⟩v_{\mathrm{U},sm{m}^{\prime }}^{{\sigma }^{\prime }{\sigma }^{\prime \prime }}⟨s{m}^{\prime }{\sigma }^{\prime \prime }|\right]|t^{\prime }{n}^{\prime }{\sigma }^{\prime \prime \prime }⟩c_{\nu ,t^{\prime }{n}^{\prime }}^{{\sigma }^{\prime \prime \prime }},$ (120)

dual

${\displaystyle \sum _{sm{m}^{\prime }}}{v}_{\mathrm{U},sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}{\displaystyle \frac{\partial n_{sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}}{\partial c_{\nu ,tn}^{\sigma ,*}}}$ $=$ ${\displaystyle \sum _{{\sigma }^{\prime \prime \prime }}}{\displaystyle \sum _{t^{\prime }{n}^{\prime }}}⟨tn\sigma |{\displaystyle \frac{1}{2}}{\displaystyle \sum _{{\sigma }^{\prime }{\sigma }^{\prime \prime }}}{\displaystyle \sum _{sm{m}^{\prime }}}\left[|\widetilde{sm{\sigma }^{\prime }}⟩v_{\mathrm{U},sm{m}^{\prime }}^{{\sigma }^{\prime }{\sigma }^{\prime \prime }}⟨s{m}^{\prime }{\sigma }^{\prime \prime }|+|sm{\sigma }^{\prime }⟩v_{\mathrm{U},sm{m}^{\prime }}^{{\sigma }^{\prime }{\sigma }^{\prime \prime }}⟨\widetilde{s{m}^{\prime }{\sigma }^{\prime \prime }}|\right]|t^{\prime }{n}^{\prime }{\sigma }^{\prime \prime \prime }⟩c_{\nu ,t^{\prime }{n}^{\prime }}^{{\sigma }^{\prime \prime \prime }}.$ (121)

Therefore, the effective projector potentials ${\widehat{v}}_{\mathrm{U}}$ can be expressed by

on-site

${\widehat{v}}_{\mathrm{U}}$ $=$ ${\displaystyle \sum _{\sigma {\sigma }^{\prime }}}{\displaystyle \sum _{sm{m}^{\prime }}}|\widetilde{sm\sigma }⟩v_{\mathrm{U},sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}⟨\widetilde{s{m}^{\prime }{\sigma }^{\prime }}|,$ (122)

full

${\widehat{v}}_{\mathrm{U}}$ $=$ ${\displaystyle \sum _{\sigma {\sigma }^{\prime }}}{\displaystyle \sum _{sm{m}^{\prime }}}|sm\sigma ⟩v_{\mathrm{U},sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}⟨s{m}^{\prime }{\sigma }^{\prime }|,$ (123)

dual

${\widehat{v}}_{\mathrm{U}}$ $=$ ${\displaystyle \frac{1}{2}}{\displaystyle \sum _{\sigma {\sigma }^{\prime }}}{\displaystyle \sum _{sm{m}^{\prime }}}\left[|\widetilde{sm\sigma }⟩v_{\mathrm{U},sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}⟨s{m}^{\prime }{\sigma }^{\prime }|+|sm\sigma ⟩v_{\mathrm{U},sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}⟨\widetilde{s{m}^{\prime }{\sigma }^{\prime }}|\right].$ (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:

$⟨tn\sigma |{\widehat{v}}_{\mathrm{U}}|t^{\prime }{n}^{\prime }{\sigma }^{\prime }⟩$ $=$ ${\displaystyle \frac{1}{2}}{\displaystyle \sum _{m^{\prime }}}{v}_{\mathrm{U},tn{m}^{\prime }}^{\sigma {\sigma }^{\prime }}{S}_{t{m}^{\prime },t^{\prime }{n}^{\prime }}+{\displaystyle \frac{1}{2}}{\displaystyle \sum _m}{S}_{tn,t^{\prime }m}{v}_{\mathrm{U},t^{\prime }m{n}^{\prime }}^{\sigma {\sigma }^{\prime }},$ (125)

$=$ ${(⟨t^{\prime }{n}^{\prime }{\sigma }^{\prime }|{\widehat{v}}_{\mathrm{U}}|tn\sigma ⟩)}^{*}.$

It should be noted that in the full and dual the $v_U^{\sigma }$ of the site $i$ can affect the different sites by the projector potentials Eqs. (123) and (124) because of the overlap.

3.3 Force on atom

The force on atom is evaluated by

${\displaystyle \frac{\partial E_{\mathrm{LDA}+\mathrm{U}}}{\partial {𝐑}_k}}$ $=$ ${\displaystyle \frac{\partial E_{\mathrm{LDA}}}{\partial {𝐑}_k}}+{\displaystyle \frac{\partial E_{\mathrm{U}}}{\partial {𝐑}_k}}.$ (126)

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

${\displaystyle \frac{\partial E_{\mathrm{U}}}{\partial {𝐑}_k}}$ $=$ ${\displaystyle \sum _{\sigma {\sigma }^{\prime }}}{\displaystyle \sum _{sm{m}^{\prime }}}{\displaystyle \frac{\partial E_{\mathrm{U}}}{\partial n_{sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}}}{\displaystyle \frac{\partial n_{sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}}{\partial {𝐑}_k}},$ (127)

$=$ ${\displaystyle \sum _{\sigma {\sigma }^{\prime }}}{\displaystyle \sum _{sm{m}^{\prime }}}{v}_{\mathrm{U},sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}{\displaystyle \frac{\partial n_{sm{m}^{\prime }}^{\sigma {\sigma }^{\prime }}}{\partial {𝐑}_k}},$

$=$ ${\displaystyle \sum _{\nu }}{f}_{\nu }{\displaystyle \sum _{\sigma {\sigma }^{\prime }}}{\displaystyle \sum _{tn,t^{\prime }{n}^{\prime }}}\{{\displaystyle \frac{\partial c_{\nu ,tn}^{\sigma ,*}}{\partial {𝐑}_k}}⟨tn\sigma |{\widehat{v}}_{\mathrm{U}}|t^{\prime }{n}^{\prime }{\sigma }^{\prime }⟩c_{\nu ,t^{\prime }{n}^{\prime }}^{{\sigma }^{\prime }}+c_{\nu ,tn}^{\sigma ,*}⟨tn\sigma |{\widehat{v}}_{\mathrm{U}}|t^{\prime }{n}^{\prime }{\sigma }^{\prime }⟩{\displaystyle \frac{\partial c_{\nu ,t^{\prime }{n}^{\prime }}^{{\sigma }^{\prime }}}{\partial {𝐑}_k}}$

$\mathrm{ }\mathit{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.5em}}+c_{\nu ,tn}^{\sigma ,*}{c}_{\nu ,t^{\prime }{n}^{\prime }}^{{\sigma }^{\prime }}{\displaystyle \frac{\partial ⟨tn\sigma |{\widehat{v}}_{\mathrm{U}}|t^{\prime }{n}^{\prime }{\sigma }^{\prime }⟩}{\partial {𝐑}_k}}\}.$

Considering $H{c}_{\nu }={\epsilon }_{\nu }S{c}_{\nu }$ and $C^{†}SC=\mathrm{I}$, 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.

3.4 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 $t_{2g}$ and $e_g$ states, and the five of seven d-electrons are occupied in $t_{2g}$ and $e_g$ states of the majority spin, and remaining two d-electrons are occupied in the $t_{2g}$ state of the minority spin. Then, it depends on the initial occupancy ratios for the $t_{2g}$ states of the minority spin how the remaining two d-electrons are occupied in three $t_{2g}$ 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:

$\mathrm{diagonalize}$ $d_s=V^{†}{n}_{s}V\mathit{\hspace{1em}\hspace{1em}}{d}_s:\mathrm{ascending\; order}$ (128)

$\mathrm{summation}$ $D={\displaystyle \sum _{m=1}^{2(2l+1)}}{d}_{sm}$ (129)

$\mathrm{redistribution}$ $d_{4l+2}^{\prime }=1,$ (130)

$d_{4l+1}^{\prime }=1,$

$\mathrm{\dots },$

$d_m^{\prime }=D-(4l+2-m),$

$d_{m-1}^{\prime }=0,\mathrm{\dots }.,$

$\mathrm{where\; D}$ $=$ ${\displaystyle \sum _m}{d}_m^{\prime }$ (131)

$\mathrm{backtrasform}$ $n^{\prime }{}_s=Vd^{\prime }{}_m{V}^{†}.$ (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,…

3.5 Density of states

Define

$P_{i\kappa }^{\sigma {\sigma }^{\prime }}(E)$ $=$ ${\displaystyle \sum _{𝐤}}{\displaystyle \sum _{\nu }}{c}_{𝐤\nu ,i\kappa }^{\sigma ,*}{c}_{𝐤\nu ,i^{\prime }{\kappa }^{\prime }}^{{\sigma }^{\prime }}{S}_{i^{\prime }{\kappa }^{\prime },i\kappa }\delta (E-{\epsilon }_{\nu 𝐤}).$ (133)

Then, the density of states, $D$, is given by

$D_{i\kappa }^{\uparrow }(E)$ $=$ ${\displaystyle \frac{1}{2}}(P_{i\kappa }^{\alpha \alpha }+P_{i\kappa }^{\beta \beta })+{\displaystyle \frac{1}{2}}(P_{i\kappa }^{\alpha \alpha }-P_{i\kappa }^{\beta \beta })\cos ({\theta }_i)+\left(\mathrm{Re}{P}_{i\kappa }^{\alpha \beta }\cos ({\varphi }_i)-\mathrm{Im}{P}_{i\kappa }^{\alpha \beta }\sin ({\varphi }_i)\right)\sin ({\theta }_i),$ (134)

$D_{i\kappa }^{\downarrow }(E)$ $=$ ${\displaystyle \frac{1}{2}}(P_{i\kappa }^{\alpha \alpha }+P_{i\kappa }^{\beta \beta })-{\displaystyle \frac{1}{2}}(P_{i\kappa }^{\alpha \alpha }-P_{i\kappa }^{\beta \beta })\cos ({\theta }_i)-\left(\mathrm{Re}{P}_{i\kappa }^{\alpha \beta }\cos ({\varphi }_i)-\mathrm{Im}{P}_{i\kappa }^{\alpha \beta }\sin ({\varphi }_i)\right)\sin ({\theta }_i).$ (135)

Also, the Mulliken populations, $Q$, are given by

$Q_i^{\sigma {\sigma }^{\prime }}$ $=$ ${\displaystyle \sum _{\kappa }}{\displaystyle \int 𝑑Ef(E)P_{i\kappa }^{\sigma {\sigma }^{\prime }}(E)}.$ (136)

The local spin direction is determined by

${\varphi }_i$ $=$ $-\arctan ({\displaystyle \frac{\mathrm{Im}{Q}_i^{\alpha \beta }}{\mathrm{Re}{Q}_i^{\alpha \beta }}}),$ (137)

${\theta }_i$ $=$ $\arctan ({\displaystyle \frac{2(\mathrm{Re}{Q}_i^{\alpha \beta }\cos ({\varphi }_i)-\mathrm{Im}{Q}_i^{\alpha \beta }\sin ({\varphi }_i))}{Q_i^{\alpha \alpha }-Q_i^{\beta \beta }}}.)$ (138)

3.6 Zeeman term

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

$E_z$ $=$ $E_{zs}+E_{zo},$ (139)

where

$E_{zs}$ $=$ ${\displaystyle \sum _i}{𝐁}_i^s\cdot {𝐬}_i={\displaystyle \sum _i}(B_{ix}^s{s}_{ix}+B_{iy}^s{s}_{iy}+B_{iz}^s{s}_{iz}),$ (140)

$E_{zo}$ $=$ ${\displaystyle \frac{1}{2}}{\displaystyle \sum _i}{𝐁}_i^o\cdot {𝐥}_i={\displaystyle \frac{1}{2}}{\displaystyle \sum _i}(B_{ix}^o{l}_{ix}+B_{iy}^o{l}_{iy}+B_{iz}^o{l}_{iz}).$ (141)

The vector components of the spin magnetic moment are given by

$s_{ix}$ $=$ ${\displaystyle \frac{1}{2}}(N^{\uparrow }-N^{\downarrow })\sin ({\theta }_i)\cos ({\varphi }_i),$ (142)

$s_{iy}$ $=$ ${\displaystyle \frac{1}{2}}(N^{\uparrow }-N^{\downarrow })\sin ({\theta }_i)\sin ({\varphi }_i),$ (143)

$s_{iz}$ $=$ ${\displaystyle \frac{1}{2}}(N^{\uparrow }-N^{\downarrow })\cos ({\theta }_i)$ (144)

with

$N_{i\uparrow }^{\prime }$ $=$ ${\displaystyle \frac{1}{2}}(N_{i\alpha \alpha }+N_{i\beta \beta })+{\displaystyle \frac{1}{2}}(N_{i\alpha \alpha }-N_{i\beta \beta })\cos ({\theta }_i)+\left(\mathrm{Re}{N}_{i\alpha \beta }\cos ({\varphi }_i)-\mathrm{Im}{N}_{i\alpha \beta }\sin ({\varphi }_i)\right)\sin ({\theta }_i),$ (145)

$N_{i\downarrow }^{\prime }$ $=$ ${\displaystyle \frac{1}{2}}(N_{i\alpha \alpha }+N_{i\beta \beta })-{\displaystyle \frac{1}{2}}(N_{i\alpha \alpha }-N_{i\beta \beta })\cos ({\theta }_i)-\left(\mathrm{Re}{N}_{i\alpha \beta }\cos ({\varphi }_i)-\mathrm{Im}{N}_{i\alpha \beta }\sin ({\varphi }_i)\right)\sin ({\theta }_i),$ (146)

where $N_{i\sigma {\sigma }^{\prime }}$ is given by

$N_{i\sigma {\sigma }^{\prime }}$ $=$ $\mathrm{Tr}(n_i^{\sigma {\sigma }^{\prime }}).$ (147)

After some alegebra we have

$s_{ix}$ $=$ ${\displaystyle \frac{1}{2}}(N_{i\alpha \beta }+N_{i\beta \alpha }),$ (148)

$s_{iy}$ $=$ ${\displaystyle \frac{i}{2}}(N_{i\alpha \beta }-N_{i\beta \alpha }),$ (149)

$s_{iz}$ $=$ ${\displaystyle \frac{1}{2}}(N_{i\alpha \alpha }-N_{i\beta \beta }).$ (150)

The vector components of the orbital magnetic moment are given by

$l_{iv}$ $=$ ${\displaystyle \int 𝑑E\sum _{𝐤}\sum _{\nu }f(E)⟨{\psi }_{𝐤\nu }|{\widehat{l}}_v|{\psi }_{𝐤\nu }⟩\delta (E-{\epsilon }_{𝐤\nu })},$ (151)

$=$ ${\displaystyle \int 𝑑E\sum _{𝐤}\sum _{\nu }f(E)\left[⟨{\phi }_{𝐤\nu }^{\alpha }|{\widehat{l}}_v|{\phi }_{𝐤\nu }^{\alpha }⟩+⟨{\phi }_{𝐤\nu }^{\beta }|{\widehat{l}}_v|{\phi }_{𝐤\nu }^{\beta }⟩\right]\delta (E-{\epsilon }_{𝐤\nu })},$

$=$ ${\displaystyle \sum _{𝐤}}{\displaystyle \sum _{\nu }}f({\epsilon }_{𝐤\nu })\left[{\displaystyle \sum _{\kappa ,{\kappa }^{\prime }}}{c}_{𝐤\nu ,i\kappa }^{\alpha ,*}{c}_{𝐤\nu ,i{\kappa }^{\prime }}^{\alpha }⟨{\varphi }_{i\kappa }^{\alpha }|{\widehat{l}}_v|{\varphi }_{i{\kappa }^{\prime }}^{\alpha }⟩+c_{𝐤\nu ,i\kappa }^{\beta ,*}{c}_{𝐤\nu ,i{\kappa }^{\prime }}^{\beta }⟨{\varphi }_{i\kappa }^{\beta }|{\widehat{l}}_v|{\varphi }_{i{\kappa }^{\prime }}^{\beta }⟩\right],$

$=$ ${\displaystyle \sum _{\kappa ,{\kappa }^{\prime }}}{\rho }_{i\kappa ,i{\kappa }^{\prime }}^{\alpha \alpha }⟨{\varphi }_{i\kappa }^{\alpha }|{\widehat{l}}_v|{\varphi }_{i{\kappa }^{\prime }}^{\alpha }⟩+{\rho }_{i\kappa ,i{\kappa }^{\prime }}^{\beta \beta }⟨{\varphi }_{i\kappa }^{\beta }|{\widehat{l}}_v|{\varphi }_{i{\kappa }^{\prime }}^{\beta }⟩,$

where $v=x,y,\mathrm{or}z$. Noting that

${\widehat{l}}_x$ $=$ ${\displaystyle \frac{1}{2}}({\widehat{l}}_{+}+{\widehat{l}}_{-}),$ (152)

${\widehat{l}}_y$ $=$ ${\displaystyle \frac{1}{2i}}({\widehat{l}}_{+}-{\widehat{l}}_{-}),$ (153)

${\widehat{l}}_z{Y}_l^m$ $=$ $m{Y}_l^m,$ (154)

${\widehat{l}}_{+}{Y}_l^m$ $=$ $\sqrt{(l-m)(l+m+1)}{Y}_l^{m+1},$ (155)

${\widehat{l}}_{-}{Y}_l^m$ $=$ $\sqrt{(l+m)(l-m+1)}{Y}_l^{m-1},$ (156)

and considering a unitary transformation of the spherical harmonic functions into a set of real harmonic functions defined by

$Y_{px}$ $=$ ${\displaystyle \frac{1}{\sqrt{2}}}(-Y_1^{-1}+Y_1^1),$ (157)

$Y_{py}$ $=$ ${\displaystyle \frac{1}{i\sqrt{2}}}(Y_1^{-1}+Y_1^1),$ (158)

$Y_{pz}$ $=$ $-Y_1^0,$ (159)

it can be shown that

$\left(\begin{array}{c}\hfill ⟨Y_{px}|\hfill \\ \hfill ⟨Y_{py}|\hfill \\ \hfill ⟨Y_{pz}|\hfill \end{array}\right){\widehat{l}}_x(|Y_{px}⟩,|Y_{py}⟩,|Y_{pz}⟩)$ $=$ (160)

$\left(\begin{array}{c}\hfill ⟨Y_{px}|\hfill \\ \hfill ⟨Y_{py}|\hfill \\ \hfill ⟨Y_{pz}|\hfill \end{array}\right){\widehat{l}}_y(|Y_{px}⟩,|Y_{py}⟩,|Y_{pz}⟩)$ $=$ (161)

$\left(\begin{array}{c}\hfill ⟨Y_{px}|\hfill \\ \hfill ⟨Y_{py}|\hfill \\ \hfill ⟨Y_{pz}|\hfill \end{array}\right){\widehat{l}}_z(|Y_{px}⟩,|Y_{py}⟩,|Y_{pz}⟩)$ $=$ (162)

It is noted that the expectation values of ${\widehat{l}}_v$ in terms of the real harmonic functions are purely imaginary numbers. The unitary transformation for the other $L$-components can be found in a subroutine ’Set_Comp2Real()’ in ’SetPara_DFT.c’. Thus, one can obtain the matrix representation for ${\widehat{l}}_v$ in terms of the real harmonic functions.