LDA+U Method: Ver. 1.0

1 Total energy

In conjunction with on-site terms of the unrestricted Hartree-Fock theory, the total energy of a LDA+U method [1] within the collinear spin treatment could be defined by

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

with

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

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

where $i$ is a site index, $l$ an angular momemtum quantum number, $p$ a multiplicity number of radial basis functions, $\sigma$ a spin index, and $s$ an organized index of $(ipl)$. $N$ is an diagonalized occupation matrix. $U$ is the effective Coulomb electron-electron interaction energy. Considering the rotational invariance of total energy with respect to each subshell $s$, Eq. (2) can be transformed as follows:

$E_{\mathrm{U}}$ $=$ ${\displaystyle \frac{1}{2}}{\displaystyle \sum _{\sigma }}{\displaystyle \sum _s}{U}_s\left[\mathrm{Tr}(A_s{N}_s^{\sigma }{A}_s^{†})-\mathrm{Tr}(A_s{N}_s^{\sigma }{A}_s^{†}{A}_s{N}_s^{\sigma }{A}_s^{†})\right],$ (3)

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

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

In the Eq. (3), although off-diagonal occupation terms in each subshell $s$ are taken into account, however, those between subshells are neglected. This treatment is consistent with their rotational invariant functional by Dudarev et al. [2], 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 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 $E_{\mathrm{LDA}+\mathrm{U}}$ can be expressed in terms of the Kohn-Sham eigenenergies ${\epsilon }_{\nu }^{\sigma }$ as follows:

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

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

$=$ $E_{\mathrm{band}}+\mathrm{\Delta }{E}_{\mathrm{LDA}}+{\displaystyle \frac{1}{2}}{\displaystyle \sum _{\sigma }}{\displaystyle \sum _s}{U}_s{\displaystyle \sum _{m,m^{\prime }}}{n}_{sm{m}^{\prime }}^{\sigma }{n}_{s{m}^{\prime }m}^{\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 couting corrections of LDA- and U-energies, respectively.

2 Occupation number

The occupation number $n$ may be defined by

$n_{sm{m}^{\prime }}^{\sigma }$ $=$ ${\displaystyle \sum _{\nu }}{f}_{\nu }⟨{\psi }_{\nu }^{\sigma }|{\widehat{n}}_{sm{m}^{\prime }}^{\sigma }|{\psi }_{\nu }^{\sigma }⟩,$ (5)

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

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

full

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

dual

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

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 ⟩$ (9)

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

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

The on-site and full occupation number operators have been proposed by Eschrig et al. [3] and Pickett et al. [4], 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:

${\displaystyle \sum _{\sigma }}\mathrm{Tr}(n^{\sigma })={\displaystyle \sum _{\sigma }}{\displaystyle \frac{1}{2}}\left\{\mathrm{Tr}(S{\rho }^{\sigma })+\mathrm{Tr}({\rho }^{\sigma }S)\right\}=N_{\mathrm{ele}},$ (11)

where ${\rho }^{\sigma }$ is the density matrix defined by

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

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

with a density operator:

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

The notes limit the discussion to non-Bloch wave functions for simplicity, but the extension is straightforward. 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 }={\rho }_{sm,s{m}^{\prime }}^{\sigma },$ (14)

full

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

dual

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

3 Effective potential

The derivative of the total energy Eq. (1) 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 ,*}}},$ (17)

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

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

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

with
on-site

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

full

${\displaystyle \frac{\partial n_{sm{m}^{\prime }}^{\sigma }}{\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 }$ (19)

dual

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

Substituting Eqs. (18)-(20) for the second term of Eq. (17), we see
on-site

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

full

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

dual

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

Therefore, the effective projector potentials ${\widehat{v}}_{\mathrm{U}}^{\sigma }$ can be expressed by
on-site

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

full

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

dual

${\widehat{v}}_{\mathrm{U}}^{\sigma }={\displaystyle \frac{1}{2}}{\displaystyle \sum _{sm{m}^{\prime }}}\left[|\widetilde{sm\sigma }⟩v_{\mathrm{U},sm{m}^{\prime }}^{\sigma }⟨s{m}^{\prime }\sigma |+|sm\sigma ⟩v_{\mathrm{U},sm{m}^{\prime }}^{\sigma }⟨\widetilde{s{m}^{\prime }\sigma }|\right].$ (26)

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}}^{\sigma }|t^{\prime }{n}^{\prime }\sigma ⟩$ $=$ ${\displaystyle \frac{1}{2}}{\displaystyle \sum _{m^{\prime }}}{v}_{\mathrm{U},tn{m}^{\prime }}^{\sigma }{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 },$ (27)

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

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 by Eqs. (25) and (26) because of the overlap.

4 Force on atom

The derivative of the total energy with respect to atomic coordinates ${\tau }_k$ consists of two contributions:

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

The first term can be evaluated in the same way as in the LDA. The second term is given by

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

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

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

Considering $H{c}_{\nu }={\epsilon }_{\nu }S{c}_{\nu }$ and $c^{†}Sc=\mathrm{I}$, the first and second terms in Eq. (30) can be transformed into derivatives of the overlap matrix. The third term in Eq. (30) is analytically differentiated, since it contains just two-center integrals.

5 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^{\sigma }=V^{†}{n}_s^{\sigma }V\mathit{\hspace{1em}\hspace{0.25em}}{d}_s^{\sigma }:\mathrm{ascending\; order}$ (31)

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

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

$d_{2l}^{\prime }=1,$

$\mathrm{\dots },$

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

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

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

$\mathrm{back}\mathrm{trasform}$ $n^{\prime }{}_s{}^{\sigma }=Vd^{\prime }{}_m{}^{\sigma }{V}^{†}$ (35)

After diagonalizing each subshell matrix consisting of occupation numbers, we introduce a polarized redistribution scheme given by Eq. (33) while keeping Eq. (34). Then, by a back transformation Eq. (35), we can obtain a polarized occupation matrix for each subshell. 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 subshell, and any orbitals: p,d,f,…

6 Orbital optimization within LDA+U

In the orbital optimization within LDA+U, let us assume that the effective U-potential in the LDA+U method is applied to the primitive basis orbital $\chi$ instead of the optimized basis orbital $\varphi$, which is more natural in a physical sense than the opposite assumption.

A Kohn-Sham (KS) orbital ${\psi }_{\mu }$ in the orbital optimization method is expressed by a linear combination of primitive orbitals $\chi$:

$|{\psi }_{\nu }^{\sigma }⟩$ $=$ ${\displaystyle \sum _{i\alpha }}{c}_{\mu ,i\alpha }^{\sigma }|{\varphi }_{i\alpha }⟩,$ (36)

$=$ ${\displaystyle \sum _{i\alpha }}{c}_{\mu ,i\alpha }^{\sigma }\left\{{\displaystyle \sum _q}{a}_{i\alpha q}|{\chi }_{i\eta }⟩\right\},$

$=$ ${\displaystyle \sum _{i\alpha }}{\displaystyle \sum _q}{c}_{\mu ,i\alpha }^{\sigma }{a}_{i\alpha q}|{\chi }_{i\eta }⟩,$

$=$ ${\displaystyle \sum _{i\eta }}\left\{{\displaystyle \sum _p}{c}_{\mu ,i\alpha }^{\sigma }{a}_{i\alpha q}\right\}|{\chi }_{i\eta }⟩,$

$=$ ${\displaystyle \sum _{i\eta }}{b}_{\mu ,i\eta }^{\sigma }|{\chi }_{i\eta }⟩,$

where $\alpha \equiv (plm)$, $\eta \equiv (qlm)$, $c$ and $b$ are LCAO coefficients for contracted and primitive orbitals, respectively, and $a$ contraction coefficients. For simplicity we consider an non-Bloch expression of the one-particle wave functions, but the extention of the below description to Bloch wave functions is straightforward. Assuming that the occupation number operators defined by Eqs. (6)-(8) are constructed by the primitive orbitals, we have the occupation numbers for the on-site, full, and dual given by
on-site

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

full

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

dual

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

where ${\varrho }^{\sigma }$ is the primitive density matrix defined by

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

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

with a primitive density operator:

${\widehat{\varrho }}_{sm,s^{\prime }{m}^{\prime }}^{\sigma }$ $=$ $|{\tilde{\chi }}_{sm\sigma }⟩⟨{\tilde{\chi }}_{s^{\prime }{m}^{\prime }\sigma }|.$ (41)

Moreover, by defining a contracted density operator:

${\widehat{\rho }}_{sm,s^{\prime }{m}^{\prime }}^{\sigma }$ $=$ $|{\tilde{\varphi }}_{sm\sigma }⟩⟨{\tilde{\varphi }}_{s^{\prime }{m}^{\prime }\sigma }|,$ (42)

we have the contracted density matrix ${\rho }^{\sigma }$ given by

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

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

Then, the primitive density matrix $\varrho$ is written by the contracted density matrix $\rho$ as follows:

${\varrho }_{iqlm,i^{\prime }{q}^{\prime }{l}^{\prime }{m}^{\prime }}^{\sigma }$ $=$ ${\displaystyle \sum _{p,p^{\prime }}}{a}_{iplmq}{a}_{i^{\prime }{p}^{\prime }{l}^{\prime }{m}^{\prime }{q}^{\prime }}{\rho }_{iplm,i^{\prime }{p}^{\prime }{l}^{\prime }{m}^{\prime }}^{\sigma }.$ (44)

Considering the variation of the total energy Eq. (1) with respect to $b$, we find the effective potentials of the LDA+U method with respect to the primitive basis orbital. They are given by the same expression as Eqs. (24)-(26), while the occupation number is given by Eqs. (37)-(39). After the Hamiltonian matrix with respect to the primitive basis orbital $\chi$ is constructed, it is transformed to that of the optimized basis orbital $\varphi$ as follows:

$⟨{\varphi }_{iplm}|\widehat{H}|{\varphi }_{i^{\prime }{p}^{\prime }{l}^{\prime }{m}^{\prime }}⟩$ $=$ ${\displaystyle \sum _{q,q^{\prime }}}{a}_{iplmq}{a}_{i^{\prime }{p}^{\prime }{l}^{\prime }{m}^{\prime }{q}^{\prime }}⟨{\chi }_{iqlm}|\widehat{H}|{\chi }_{i^{\prime }{q}^{\prime }{l}^{\prime }{m}^{\prime }}⟩.$ (45)

The Hamiltonian matrix with respect to the contracted basis orbital is diagonalized. The procedure is summarized as follows:

1. diagonalize the contracted Hamiltonian $⟨{\varphi }_{iplm}|\widehat{H}|{\varphi }_{i^{\prime }{p}^{\prime }{l}^{\prime }{m}^{\prime }}⟩$

2. calculate the contracted density matrix by Eq. (43)

3. calculate the primitive density matrix by Eq. (44)

4. calculate the occupation number by Eq. (37), (38), or (39)

5. construct the Hamitonian by Eq. (24), (25), or (26)

6. contract the Hamitonian by Eq. (45)

7. return 1

Although the optimization procedure of the contracted coefficients $a$ is not discussed here, it can be easily verified that the same procedure as in the LDA method is derived. Thus, the orbital optimization can be performed within the LDA+U method as well as the LDA method.