Preliminary notes

Note

Here and below sybol \( {\boldsymbol \sigma} \) is reserved for the vector of Pauli matrices. Spin components are labeled with \( \alpha \) or \( \beta \).

Wave-function of spin-1/2 particle is a two-component spinor:

\[ {\bf \Psi}({\bf r}) = \left( \begin{array}{c} \psi^{\uparrow}({\bf r}) \\ \psi^{\downarrow}({\bf r}) \end{array} \right) \]

Operator of spin:

\[ {\bf \hat S}=\frac{\hbar}{2}{\boldsymbol \sigma}, \]

Pauli matrices:

\[ \sigma_x=\left( \begin{array}{cc} 0 & 1 \\ 1 & 0 \\ \end{array} \right) \, \sigma_y=\left( \begin{array}{cc} 0 & -i \\ i & 0 \\ \end{array} \right) \, \sigma_z=\left( \begin{array}{cc} 1 & 0 \\ 0 & -1 \\ \end{array} \right) \]

Spin moment of an electron in quantum state \( | \Psi \rangle \):

\[ {\bf S}=\langle \Psi | {\bf \hat S} | \Psi \rangle = \frac{\hbar}{2} \langle \Psi | {\boldsymbol \sigma} | \Psi \rangle \]

Spin magnetic moment of electron:

\[ {\bf \mu}_e=\gamma_e {\bf S}, \]

where \( \gamma_e \) is the gyromagnetic ratio for the electron.

\[ \gamma_e=-\frac{g_e \mu_B}{\hbar} \;\;\; \mu_B=\frac{e\hbar}{2m_ec} \]

Here \( g_e \) is a g-factor for electron which is ~2, and \( \mu_B \) - Bohr magneton (defined as positive constant). Finally, magnetic moment of electron:

\[ {\bf \mu}_e=-{\bf \mu}_B \langle \Psi | {\boldsymbol \sigma} | \Psi \rangle \]

Potential energy of magnetic dipole in magnetic field:

\[ U=-{\bf B}{\bf \mu}={\bf \mu}_B {\bf B} \langle \Psi | {\boldsymbol \sigma} | \Psi \rangle \]

Density and magnetization

In magnetic calculations we have charge density \( \rho({\bf r}) \) (scalar function) and magnetization density \( {\bf m}({\bf r}) \) (vector function). Density is defined as:

\[ \rho({\bf r}) = \sum_{j}^{occ} \Psi_{j}^{*}({\bf r}){\bf I} \Psi_{j}({\bf r}) = \sum_{j}^{occ} \psi_{j}^{\uparrow *}({\bf r}) \psi_{j}^{\uparrow}({\bf r}) + \psi_{j}^{\downarrow *}({\bf r}) \psi_{j}^{\downarrow}({\bf r}) \]

Magnetization is defined as:

\[ {\bf m}({\bf r}) = \sum_{j}^{occ} \Psi_{j}^{*}({\bf r}) {\boldsymbol \sigma} \Psi_{j}({\bf r}) \]

\[ m_x({\bf r}) = \sum_{j}^{occ} \psi_{j}^{\uparrow *}({\bf r}) \psi_{j}^{\downarrow}({\bf r}) + \psi_{j}^{\downarrow *}({\bf r}) \psi_{j}^{\uparrow}({\bf r}) \]

\[ m_y({\bf r}) = \sum_{j}^{occ} -i \psi_{j}^{\uparrow *}({\bf r}) \psi_{j}^{\downarrow}({\bf r}) + i \psi_{j}^{\downarrow *}({\bf r}) \psi_{j}^{\uparrow}({\bf r}) \]

\[ m_z({\bf r}) = \sum_{j}^{occ} \psi_{j}^{\uparrow *}({\bf r}) \psi_{j}^{\uparrow}({\bf r}) - \psi_{j}^{\downarrow *}({\bf r}) \psi_{j}^{\downarrow}({\bf r}) \]

Density and magnetization can be grouped into a \( 2 \times 2 \) density matrix, which is defined as:

\[ {\boldsymbol \rho}({\bf r}) = \frac{1}{2} \Big( {\bf I}\rho({\bf r}) + {\boldsymbol \sigma} {\bf m}({\bf r})\Big) = \frac{1}{2} \left( \begin{array}{cc} \rho({\bf r}) + m_z({\bf r}) & m_x({\bf r}) - i m_y({\bf r}) \\ m_x({\bf r}) + i m_y({\bf r}) & \rho({\bf r}) - m_z({\bf r}) \end{array} \right) = \sum_{j}^{occ} \left( \begin{array}{cc} \psi_{j}^{\uparrow *}({\bf r}) \psi_{j}^{\uparrow}({\bf r}) & \psi_{j}^{\downarrow *}({\bf r}) \psi_{j}^{\uparrow}({\bf r}) \\ \psi_{j}^{\uparrow *}({\bf r}) \psi_{j}^{\downarrow}({\bf r}) & \psi_{j}^{\downarrow *}({\bf r}) \psi_{j}^{\downarrow}({\bf r}) \end{array} \right) \]

or simply

\[ \rho_{\alpha \beta}({\bf r}) = \sum_{j}^{occ} \psi_{j}^{\beta *}({\bf r})\psi_{j}^{\alpha}({\bf r}) \]

Pay attention to the order of spin indices in the \( 2 \times 2 \) density matrix. External potential \( v^{ext}({\bf r}) \) and external magnetic field \( {\bf B}^{ext}({\bf r}) \) can also be grouped into a \( 2 \times 2 \) matrix:

\[ V_{\alpha\beta}^{ext}({\bf r})=\Big({\bf I}v^{ext}({\bf r})+\mu_{B}{\boldsymbol \sigma}{\bf B}^{ext}({\bf r}) \Big) = \left( \begin{array}{cc} v^{ext}({\bf r})+\mu_{B}B_z^{ext}({\bf r}) & \mu_{B} \Big( B_x^{ext}({\bf r})-iB_y^{exp}({\bf r}) \Big) \\ \mu_{B} \Big( B_x^{ext}({\bf r})+iB_y^{ext}({\bf r}) \Big) & v^{ext}({\bf r})-\mu_{B}B_z^{ext}({\bf r}) \end{array} \right) \]

Let's check that potential energy in external fields can be written in the following way:

\[ E^{ext}=\int Tr \Big( {\boldsymbol \rho}({\bf r}) {\bf V}^{ext}({\bf r}) \Big) d{\bf r} = \sum_{\alpha\beta} \int \rho_{\alpha\beta}({\bf r})V_{\beta\alpha}^{ext}({\bf r}) d{\bf r} \]

(below \( {\bf r} \), \( ext \) and \( \int d{\bf r} \) are dropped for simplicity)

\[ \begin{eqnarray} \rho_{11}V_{11} &= \frac{1}{2}(\rho+m_z)(v+\mu_{B}B_z) = \frac{1}{2}(\rho v +\mu_{B} \rho B_z + m_z v + \mu_{B}m_zB_z) \\ \rho_{22}V_{22} &= \frac{1}{2}(\rho-m_z)(v-\mu_{B}B_z) = \frac{1}{2}(\rho v -\mu_{B} \rho B_z - m_z v + \mu_{B}m_zB_z) \\ \rho_{12}V_{21} &= \frac{1}{2}(m_x-im_y)\Big( \mu_{B}( B_x+iB_y) \Big) = \frac{\mu_B}{2}(m_xB_x+im_xB_y-im_yB_x+m_yB_y) \\ \rho_{21}V_{12} &= \frac{1}{2}(m_x+im_y)\Big( \mu_{B}( B_x-iB_y) \Big) = \frac{\mu_B}{2}(m_xB_x-im_xB_y+im_yB_x+m_yB_y) \end{eqnarray} \]

The sum of this four terms will give exactly \( \int \rho({\bf r}) v^{ext}({\bf r}) + \mu_{B}{\bf m}({\bf r}) {\bf B}^{ext}({\bf r}) d{\bf r}\)

Total energy variation

To derive Kohn-Sham equations we need to write total energy functional of density matrix \( \rho_{\alpha\beta}({\bf r}) \):

\[ E^{tot}[\rho_{\alpha\beta}] = E^{kin}[\rho_{\alpha\beta}] + E^{H}[\rho_{\alpha\beta}] + E^{ext}[\rho_{\alpha\beta}] + E^{XC}[\rho_{\alpha\beta}] \]

Kinetic energy of non-interacting electrons is written in the following way:

\[ E^{kin}[\rho_{\alpha\beta}] \equiv E^{kin}[\Psi[\rho_{\alpha\beta}]] = -\frac{1}{2} \sum_{i}^{occ}\sum_{\alpha} \int \psi_{i}^{\alpha*}({\bf r}) \nabla^{2} \psi_{i}^{\alpha}({\bf r})d^3{\bf r} \]

Hartree energy:

\[ E^{H}[\rho_{\alpha\beta}]= \frac{1}{2} \iint \frac{\rho({\bf r})\rho({\bf r'})}{|{\bf r}-{\bf r'}|} d{\bf r} d{\bf r'} = \frac{1}{2} \iint \sum_{\alpha\beta}\delta_{\alpha\beta} \frac{\rho_{\alpha\beta}({\bf r}) \rho({\bf r'})}{|{\bf r}-{\bf r'}|} d{\bf r} d{\bf r'} \]

where \( \rho({\bf r}) = Tr \rho_{\alpha\beta}({\bf r}) \).

Exchange-correlation energy:

\[ E^{XC}[\rho_{\alpha\beta}({\bf r})] \equiv E^{XC}[\rho({\bf r}),|{\bf m}({\bf r})|] = \int \rho({\bf r}) \eta_{XC}(\rho({\bf r}), m({\bf r})) d{\bf r} = \int \rho({\bf r}) \eta_{XC}(\rho^{\uparrow}({\bf r}), \rho_{\downarrow}({\bf r})) d{\bf r} \]

Now we can write the total energy variation over auxiliary orbitals with constrain of orbital normalization:

\[ \frac{\delta \Big( E^{tot}+\varepsilon_i \big( 1-\sum_{\alpha} \int \psi^{\alpha*}_{i}({\bf r}) \psi^{\alpha}_{i}({\bf r})d{\bf r} \big) \Big)} {\delta \psi_{i}^{\gamma*}({\bf r})} = 0 \]

We will use the following chain rule:

\[ \frac{\delta F[\rho_{\alpha\beta}]}{\delta \psi_{i}^{\gamma *}({\bf r})} = \sum_{\alpha' \beta'} \frac{\delta F[\rho_{\alpha\beta}]}{\delta \rho_{\alpha'\beta'}({\bf r})} \frac{\delta \rho_{\alpha'\beta'}({\bf r})}{\delta \psi_{i}^{\gamma *}({\bf r})} = \sum_{\alpha'\beta'}\frac{\delta F[\rho_{\alpha\beta}]}{\delta \rho_{\alpha'\beta'}({\bf r})} \psi_{i}^{\alpha'}({\bf r}) \delta_{\beta'\gamma} = \sum_{\alpha'}\frac{\delta F[\rho_{\alpha\beta}]}{\delta \rho_{\alpha'\gamma}({\bf r})}\psi_{i}^{\alpha'}({\bf r}) \]

Variation of the normalization integral:

\[ \frac{\delta \sum_{\alpha} \int \psi_{i}^{\alpha*}({\bf r}) \psi_{i}^{\alpha}({\bf r}) d {\bf r} }{\delta \psi_{i}^{\gamma *}({\bf r})} = \psi_{i}^{\gamma}({\bf r}) \]

Variation of the kinetic energy functional:

\[ \frac{\delta E^{kin}}{\delta \psi_{i}^{\gamma*}({\bf r})} = -\frac{1}{2} \sum_{\alpha} \nabla^{2} \psi_{i}^{\alpha}({\bf r}) \delta_{\alpha\gamma} = -\frac{1}{2}\nabla^{2}\psi_{i}^{\gamma}({\bf r}) \]

Variation of the Hartree energy functional:

\[ \frac{\delta E^{H}[\rho_{\alpha\beta}]}{\delta \psi_{i}^{\gamma *}({\bf r})} = \sum_{\alpha'} \sum_{\alpha\beta} \delta_{\alpha\beta} \frac{1}{2} \int \frac{ \rho({\bf r'})}{|{\bf r}-{\bf r'}|} d{\bf r'} \delta_{\alpha\alpha'}\delta_{\beta\gamma} \psi_{i}^{\alpha'}({\bf r}) = v^{H}({\bf r}) \psi_{i}^{\gamma}({\bf r}) \]

Variation of the external energy functional:

\[ \frac{\delta E^{ext}[\rho_{\alpha\beta}]}{\delta \psi_{i}^{\gamma*}({\bf r}) } = \sum_{\alpha'} \sum_{\alpha\beta} V_{\beta\alpha}^{ext}({\bf r}) \delta_{\alpha\alpha'} \delta_{\beta\gamma} \psi_{i}^{\alpha'}({\bf r})= \sum_{\alpha} V_{\gamma\alpha}^{ext}({\bf r}) \psi_{i}^{\alpha}({\bf r}) \]

Variation of the exchange-correlation functional:

\[ \frac{\delta E^{XC}[\rho_{\alpha\beta}]}{ \delta \psi_{i}^{\gamma*}({\bf r}) } = \frac{\delta E^{XC}[\rho_{\alpha\beta}]}{ \delta \rho({\bf r})} \frac{\delta \rho({\bf r})}{\delta \psi_{i}^{\gamma*}({\bf r})} + \frac{\delta E^{XC}[\rho_{\alpha\beta}]}{ \delta m({\bf r})} \sum_{p=x,y,z} \frac{\delta m({\bf r})}{ \delta m_p({\bf r})} \frac{\delta m_p({\bf r})}{\delta \psi_{i}^{\gamma*}({\bf r})} \]

where \( m({\bf r}) = |{\bf m}({\bf r})|\) is the length of magnetization vector.

First term:

\[ \frac{\delta E^{XC}[\rho_{\alpha\beta}]}{ \delta \rho({\bf r})} \frac{\delta \rho({\bf r})}{\delta \psi_{i}^{\gamma*}({\bf r})} = v^{XC}({\bf r}) \psi_{i}^{\gamma}({\bf r}) \]

Second term:

\[ \frac{\delta E^{XC}[\rho_{\alpha\beta}]}{ \delta m({\bf r})} \sum_{p=x,y,z} \frac{\delta m({\bf r})}{ \delta m_p({\bf r})} \frac{\delta m_p({\bf r})}{\delta \psi_{i}^{\gamma*}({\bf r})} = B^{XC}({\bf r}) \hat {\bf m} \sum_{\beta} {\boldsymbol \sigma}_{\gamma \beta} \psi_{i}^{\beta}({\bf r}) \]

where \( B^{XC}({\bf r}) = \frac{\delta E^{XC}[\rho_{\alpha\beta}]}{ \delta m({\bf r})} \), \( \hat {\bf m}({\bf r}) = \frac{\delta m({\bf r})}{ \delta {\bf m}({\bf r})} \) is the unit vector, parallel to \( {\bf m}({\bf r}) \) and \( {\bf m}({\bf r}) = \sum_{i} \sum_{\alpha \beta} \psi_{i}^{\alpha*}({\bf r}) {\boldsymbol \sigma}_{\alpha \beta} \psi_i^{\beta}({\bf r}) \)

Similarly to external potential, exchange-correlation potential can be grouped into \( 2 \times 2 \) matrix::

\[ \frac{\delta E^{XC}[\rho_{\alpha\beta}]}{\delta \rho_{\alpha'\beta'}({\bf r})} \equiv V^{XC}_{\beta'\alpha'}({\bf r}) = \Big( {\bf I}v^{XC}({\bf r}) + {\bf B}^{XC}({\bf r}) {\boldsymbol \sigma} \Big)_{\beta'\alpha'} \]

where \({\bf B}^{XC}({\bf r}) = \hat {\bf m}({\bf r})B^{XC}({\bf r}) \) – exchange-correlation magnetic field, parallel to \( {\bf m}({\bf r}) \) at each point in space. We can now collect \( v^{H}({\bf r}) \), \( V_{\alpha\beta}^{ext}({\bf r}) \) and \(V_{\alpha\beta}^{XC}({\bf r}) \) to one effective potential:

\[ V^{eff}_{\alpha\beta}({\bf r}) = v^{H}({\bf r})\delta_{\alpha\beta} + V_{\alpha\beta}^{ext}({\bf r}) + V_{\alpha\beta}^{XC}({\bf r}) = \Big({\bf I}\big(v^{H}({\bf r})+v^{ext}({\bf r})+v^{XC}({\bf r})\big) + {\boldsymbol \sigma}\big( \mu_{B}{\bf B}^{ext}({\bf r}) + {\bf B}^{XC}({\bf r})\big)\Big)_{\alpha\beta} \]

and finally, we arrive to the following Kohn-Sham equation for each component \( \gamma \) of spinor wave-functions:

\[ -\frac{1}{2}\nabla^{2}\psi_{i}^{\gamma}({\bf r}) + \sum_{\alpha} V_{\gamma\alpha}^{eff}({\bf r}) \psi_{i}^{\alpha}({\bf r}) = \varepsilon_i \psi_{i}^{\gamma}({\bf r}) \]

or in matrix form

\[ \left( \begin{array}{cc} -\frac{1}{2}\nabla^2+V^{eff}_{\uparrow \uparrow} & V^{eff}_{\uparrow \downarrow} \\ V^{eff}_{\downarrow \uparrow} & -\frac{1}{2}\nabla^2+V^{eff}_{\downarrow \downarrow} \end{array}\right) \left(\begin{array}{c} \psi_{i}^{\uparrow}({\bf r}) \\ \psi_{i}^{\downarrow} ({\bf r}) \end{array} \right) = \varepsilon_i \left(\begin{array}{c} \psi_{i}^{\uparrow}({\bf r}) \\ \psi_{i}^{\downarrow}({\bf r}) \end{array} \right) \]

Second-variational approach

Suppose that we know first \( N_{fv} \) solutions of the following equation (so-called first variational equation):

\[ \Big(-\frac{1}{2}\nabla^2+v^{H}({\bf r})+v^{ext}({\bf r})+v^{XC}({\bf r}) \Big)\phi_{i}({\bf r}) = \epsilon_i \phi_{i}({\bf r}) \]

We can write expansion of the components of spinor wave-functions \( \psi^{\alpha}_j({\bf r}) \) in terms of first-variational states \( \phi_i({\bf r}) \):

\[ \psi_{j}^{\alpha}({\bf r}) = \sum_{i}^{N_{fv}}C_{ij}^{\alpha}\phi_{i}({\bf r}) \]

Next, we switch to matrix equation:

\[ \langle \Psi_{j'}| \hat H | \Psi_{j} \rangle = \varepsilon_j \delta_{j'j} \\ \sum_{\alpha \alpha'} \sum_{ii'} C_{i'j'}^{\alpha'*} C_{ij}^{\alpha} \langle \phi_{i'} | \hat H_{\alpha' \alpha} | \phi_{i} \rangle = \sum_{\alpha \alpha'} \sum_{ii'} C_{i'j'}^{\alpha'*} C_{ij}^{\alpha} H_{\alpha'i', \alpha i} = \varepsilon_j \delta_{j'j} \]

We can combine indices \( \{i,\alpha\} \) into one 'flat' index \( \nu \). If we also assume that the number of spinor wave-functions is equal to \( 2 N_{fv} \) then we arrive to the well-known eigen decomposition:

\[ \sum_{\nu'\nu} C_{\nu' j'}^{*} H_{\nu'\nu} C_{\nu j} = \varepsilon_j \delta_{j'j} \]

The expression for second-variational Hamiltonian is simple:

\[ \langle \phi_{i'}|\hat H_{\alpha'\alpha} |\phi_{i} \rangle = \langle \phi_{i'} | \Big(-\frac{1}{2}\nabla^2 + v^{H}({\bf r}) + v^{ext}({\bf r}) + v^{XC}({\bf r}) \Big) \delta_{\alpha\alpha'}|\phi_{i}\rangle + \langle \phi_{i'} | {\boldsymbol \sigma}_{\alpha\alpha'} \Big( \mu_{B}{\bf B}^{ext}({\bf r})+ {\bf B}^{XC}({\bf r})\Big) | \phi_{i}\rangle =\\ \epsilon_{i}\delta_{i'i}\delta_{\alpha\alpha'} + {\boldsymbol \sigma}_{\alpha\alpha'} \langle \phi_{i'} | \Big( \mu_{B}{\bf B}^{ext}({\bf r}) + {\bf B}^{XC}({\bf r})\Big) | \phi_{i}\rangle \]