Navigation

Hartree-Fock (HF) Method

Derivation

The interacting Hamiltonian for many body Schrödinger equation is (see the general QFT notes for derivation):

i\hbar\partial_t\ket{\Psi(t)} = \hat H\ket{\Psi(t)} \hat H = \hat T + \hat V = \sum_{ij} c_i^\dag\braket{i|T|j}c_j + \half \sum_{ijkl} c_i^\dag c_j^\dag\braket{ij|V|kl}c_l c_k

where \ket{i} are spin orbitals (thus the integration over \omega below) and:

\braket{i|T|j} = \int \chi_i^*({\bf x}) \left( -\half\nabla^2 - \sum_n {Z_n\over | {\bf x} -{\bf R}_n | }\right) \chi_j({\bf x})\d^3 x\, \d\omega \braket{ij|V|kl} = \int \chi_i^*({\bf x}) \chi_j^*({\bf y}) {1\over | {\bf x} - {\bf y} | } \chi_k({\bf x}) \chi_l({\bf y}) \d^3 x\, \d\omega_x\, \d^3 y\, \d\omega_y

We would like to minimize the energy E = \braket{\Psi | \hat H | \Psi } using the following basis for Z electrons:

\ket{\Psi} = c_1^\dag c_2^\dag \cdots c_Z^\dag\ket{0}

We express the energy E in this basis:

E = \braket{\Psi | \hat H | \Psi } = = \bra{0}c_Z\cdots c_2 c_1 \hat H c_1^\dag c_2^\dag \cdots c_Z^\dag\ket{0} = = \bra{0}c_Z\cdots c_2 c_1 \left(\sum_{i,j=1}^Z c_i^\dag\braket{i|T|j}c_j + \half \sum_{i,j,k,l=1}^Z c_i^\dag c_j^\dag\braket{ij|V|kl}c_l c_k\right) c_1^\dag c_2^\dag \cdots c_Z^\dag\ket{0} = = \sum_{i=1}^Z \braket{i|T|i} + \half \sum_{i,j=1}^Z \left(\braket{ij|V|ij}-\braket{ij|V|ji}\right)

We minimize it with the constrain \braket{i|j} = \delta_{ij}:

\delta \left(E - \sum_{i,j=1}^Z \epsilon_{ij} \braket{i|j}\right) = 0

We obtain:

(1)T \ket{i} + \sum_{j=1}^Z \left(\braket{j|V|ij}-\braket{j|V|ji}\right) = \epsilon_i \ket{i}

in the x-representation:

\braket{{\bf x} | T | i} + \sum_{j=1}^Z \left(\bra{{\bf x}}\braket{j|V|ij}- \bra{{\bf x}}\braket{j|V|ji}\right) = \epsilon_i \braket{{\bf x} | i} \braket{{\bf x} | T | i} + \sum_{j=1}^Z \left(\braket{j{\bf x}|V|ij}- \braket{j{\bf x}|V|ji}\right) = \epsilon_i \braket{{\bf x} | i}

And writing the individual terms explicitly (in this section, all orbitals are spin orbitals):

\braket{{\bf x} | i} = \psi_i({\bf x}) \braket{{\bf x} | T | i} = \left(-\half \nabla^2 -\sum_n {Z_n\over | {\bf x} -{\bf R}_n | } \right)\psi_i({\bf x}) \braket{j{\bf x}|V|ij} = \int \psi_j^*({\bf y}){1\over|{\bf x}-{\bf y}|} \psi_i({\bf x})\psi_j({\bf y}) \d^3 y = \int {|\psi_j({\bf y})|^2\over|{\bf x}-{\bf y}|} \d^3 y\,\,\psi_i({\bf x}) \braket{j{\bf x}|V|ji} = \int \psi_j^*({\bf y}){1\over|{\bf x}-{\bf y}|} \psi_j({\bf x})\psi_i({\bf y}) \d^3 y = \int {\psi_i({\bf y})\psi_j^*({\bf y})\over|{\bf x}-{\bf y}|} \d^3 y\,\,\psi_j({\bf x})

we get the Hartree-Fock equations:

(2)\left(-\half \nabla^2 -{Z\over |{\bf x}|} + \int {\sum_{j=1}^Z|\psi_j({\bf y})|^2\over|{\bf x}-{\bf y}|} \d^3 y\right)\psi_i({\bf x}) - \sum_{j=1}^Z\int {\psi_i({\bf y})\psi_j^*({\bf y})\over|{\bf x}-{\bf y}|} \d^3 y\,\,\psi_j({\bf x}) = \epsilon_i \psi_i({\bf x})

Let’s introduce the number density n({\bf x}), Hartree potential V_H({\bf x}) and nonlocal exchange potential V_x with its kernel U({\bf x}, {\bf y}):

n({\bf x}) = \sum_{j=1}^Z|\psi_j({\bf y})|^2 V_H({\bf x}) = \int {\sum_{j=1}^Z|\psi_j({\bf y})|^2\over|{\bf x}-{\bf y}|} \d^3 y = \int {n({\bf y})\over|{\bf x}-{\bf y}|} \d^3 y \hat V_x f({\bf x}) = - \sum_{j=1}^Z\int {f({\bf y})\psi_j^*({\bf y})\over|{\bf x}-{\bf y}|} \d^3 y\,\,\psi_j({\bf x}) = \int U({\bf x}, {\bf y}) f({\bf y}) \d^3 y U({\bf x}, {\bf y}) = - \sum_{j=1}^Z {\psi_j({\bf x})\psi_j^*({\bf y})\over|{\bf x}-{\bf y}|}

then we can write the HF equations as:

\left(-\half \nabla^2 -{Z\over |{\bf x}|} + V_H({\bf x}) + \hat V_x \right)\psi_i({\bf x}) = \epsilon_i \psi_i({\bf x}) \left(-\half \nabla^2 -{Z\over |{\bf x}|} + V_H({\bf x}) \right)\psi_i({\bf x}) + \int U({\bf x}, {\bf y}) \psi_i({\bf y}) \d^3 y = \epsilon_i \psi_i({\bf x})

The Hartree potential can be calculated by solving the Poisson equation:

\nabla^2V_H({\bf x}) = -4\pi n({\bf x})

where:

n({\bf x}) = \sum_{i=1}^Z |\psi_i({\bf x})|^2

The application of the exchange potential \hat V_x on any function f({\bf x}) can be calculated by:

\hat V_x f({\bf x}) = - \sum_{j=1}^Z W_{fj}({\bf x})\psi_j({\bf x}) W_{fj}({\bf x}) = \int {f({\bf y})\psi_j^*({\bf y})\over|{\bf x}-{\bf y}|} \d^3 y \nabla^2 W_{fj}({\bf x}) = -4\pi f({\bf x})\psi_j^*({\bf x})

Roothaan Equations For Closed Shell Systems

Starting from (1) and integrating over spins we get (here i, k are spatial orbitals, not spin orbitals):

(3)T \ket{i} + \sum_{k=1}^{N/2} \left(2\braket{k|V|ik}-\braket{k|V|ki}\right) = \epsilon_i \ket{i}

We introduce basis functions \ket{\mu} by (below the greek letters are basis functions, latin letters are spatial orbitals):

\ket{i} = \sum_\nu C_{\nu i} \ket{\nu}

substitute into (3) and multiply by \bra{\mu} from the left:

(4)\sum_\nu \braket{\mu | T | \nu} C_{\nu i} + \sum_\nu\sum_{k=1}^{N/2} \left(2\braket{\mu k|V|\nu k} -\braket{\mu k|V|k\nu}\right) C_{\nu i} = \epsilon_i \sum_\nu \braket{\mu | \nu} C_{\nu i}

Now we expand the functions \ket{k}:

(5)\sum_\nu \braket{\mu | T | \nu} C_{\nu i} + \sum_\nu \sum_{\alpha\beta} \left(2\sum_{k=1}^{N/2} C_{\alpha k} C_{\beta k}^* \right) \left(\braket{\mu \beta|V|\nu \alpha} -\half\braket{\mu \beta|V|\alpha \nu}\right) C_{\nu i} = \epsilon_i \sum_\nu \braket{\mu | \nu} C_{\nu i}

we introduce the density matrix:

\hat\rho = 2 \sum_{k=1}^{N/2} \ket{k}\bra{k} = \sum_{\alpha\beta} \ket{\alpha} 2 \sum_{k=1}^{N/2} C_{\alpha k} C_{\beta k}^*\bra{\beta} = \sum_{\alpha\beta} \ket{\alpha}P_{\alpha\beta}\bra{\beta} P_{\alpha\beta} = 2 \sum_{k=1}^{N/2} C_{\alpha k} C_{\beta k}^*

and get:

(6)\sum_\nu \left( \braket{\mu | T | \nu} + \sum_{\alpha\beta} P_{\alpha\beta} \left(\braket{\mu \beta|V|\nu \alpha} -\half \braket{\mu \beta|V|\alpha \nu}\right) \right) C_{\nu i} = \epsilon_i \sum_\nu \braket{\mu | \nu} C_{\nu i}

introducing:

F_{\mu\nu} = H_{\mu\nu}^{\mbox{core}} + G_{\mu\nu} H_{\mu\nu}^{\mbox{core}} = \braket{\mu | T | \nu} G_{\mu\nu} = \sum_{\alpha\beta} P_{\alpha\beta} \left(\braket{\mu \beta|V|\nu \alpha} -\half \braket{\mu \beta|V|\alpha \nu}\right) S_{\mu\nu} = \braket{\mu | \nu}

the equation (6) is:

(7)\sum_\nu F_{\mu\nu} C_{\nu i} = \epsilon_i \sum_\nu S_{\mu\nu} C_{\nu i}

These are the Roothaan equations. It is a generalized eigenvalue problem.

Total energy is given by (the i,j in the first equation are spin orbitals, in the other equations i,j are spatial orbitals):

E = \sum_{i} \braket{i|T|i} + \half \sum_{i,j} \left(\braket{ij|V|ij}-\braket{ij|V|ji}\right) = = \sum_{i} 2 I(i) + \sum_{i,j} \left(2 J(i, j) - K(i, j)\right) = = \sum_{i} 2 \braket{i|T|i} + \sum_{i,j} \left(2\braket{ij|V|ij}- \braket{ij|V|ji}\right) = = \sum_{i} 2 \braket{i|T|i} + 2\sum_{i,j} \left(\braket{ij|V|ij}-\half \braket{ij|V|ji}\right) = = \sum_{\mu\nu} \sum_{i} 2 \braket{\mu|T|\nu} C_{\nu i} C_{\mu i}^* + 2\sum_{\mu\nu}\sum_{\alpha\beta} \sum_{i,j} C_{\nu i} C_{\mu i}^* C_{\alpha j} C_{\beta j}^* \left(\braket{\mu\beta|V|\nu\alpha}- \half \braket{\mu\beta|V|\alpha\nu}\right) = = \sum_{\mu\nu} \braket{\mu|T|\nu} P_{\nu\mu} + \half \sum_{\mu\nu}\sum_{\alpha\beta} P_{\nu\mu} P_{\alpha\beta} \left(\braket{\mu\beta|V|\nu\alpha}- \half \braket{\mu\beta|V|\alpha\nu}\right) = = \sum_{\mu\nu} P_{\nu\mu} (H_{\mu\nu}^{\mbox{core}} +\half G_{\mu\nu}) = = \sum_{\mu\nu} P_{\nu\mu} (\half H_{\mu\nu}^{\mbox{core}} +\half (H_{\mu\nu}^{\mbox{core}} + G_{\mu\nu})) = = \half \sum_{\mu\nu} P_{\nu\mu} (H_{\mu\nu}^{\mbox{core}} +F_{\mu\nu})

The same thing can be derived in x-representation starting from (2) and introducing spatial orbitals:

(8)\left(-\half \nabla^2 -{Z\over |{\bf x}|} + \int {2\sum_{k=1}^{N/2}|\psi_k({\bf y})|^2\over|{\bf x}-{\bf y}|} \d^3 y\right)\psi_i({\bf x}) - \sum_{k=1}^{N/2}\int {\psi_i({\bf y})\psi_k^*({\bf y})\over|{\bf x}-{\bf y}|} \d^3 y\,\,\psi_k({\bf x}) = \epsilon_i \psi_i({\bf x})

We introduce basis functions \phi_\mu:

\psi_i({\bf x}) = \sum_\nu C_{\nu i} \phi_\nu({\bf x})

substitute into (8) and also multiply the whole equation by \phi_\mu^* and integrate over {\bf x}:

(9)\sum_\nu \int \phi_\mu^*({\bf x}) \left(-\half \nabla^2 -{Z\over |{\bf x}|} + \int {2\sum_{k=1}^{N/2}|\psi_k({\bf y})|^2\over|{\bf x}-{\bf y}|} \d^3 y\right)\phi_\nu({\bf x}) \d^3 x\, C_{\nu i} -\sum_\nu \int \phi_\mu^*({\bf x}) \sum_{k=1}^{N/2}\int {\phi_\nu({\bf y})\psi_k^*({\bf y})\over|{\bf x}-{\bf y}|} \d^3 y\,\,\psi_k({\bf x})\d^3 x\, C_{\nu i} = \epsilon_i \sum_\nu \int \phi_\mu^*({\bf x}) \phi_\nu({\bf x})\d^3 x \, C_{\nu i}

This can be written as:

\sum_\nu F_{\mu\nu} C_{\nu i} = \epsilon_i \sum_\nu S_{\mu\nu} C_{\nu i} F_{\mu\nu} = H_{\mu\nu}^{\mbox{core}} + G_{\mu\nu} = T_{\mu\nu} + V_{\mu\nu} + G_{\mu\nu}

where:

T_{\mu\nu} = \int \phi_\mu^*({\bf x}) \left(-\half \nabla^2 \right) \phi_\nu({\bf x}) \d^3 x = \half \int \nabla \phi_\mu^*({\bf x}) \cdot \nabla \phi_\nu({\bf x}) \d^3 x V_{\mu\nu} = \int \phi_\mu^*({\bf x}) \left(-{Z\over |{\bf x}|}\right) \phi_\nu({\bf x}) \d^3 x G_{\mu\nu} = \int \phi_\mu^*({\bf x}) \left( \int {2\sum_{k=1}^{N/2}|\psi_k({\bf y})|^2\over|{\bf x}-{\bf y}|} \d^3 y\right)\phi_\nu({\bf x}) \d^3 x -\int \phi_\mu^*({\bf x}) \sum_{k=1}^{N/2}\int {\phi_\nu({\bf y})\psi_k^*({\bf y})\over|{\bf x}-{\bf y}|} \d^3 y\,\,\psi_k({\bf x})\d^3 x S_{\mu\nu} = \int \phi_\mu^*({\bf x}) \phi_\nu({\bf x})\d^3 x

Introducing the density matrix and density:

\rho({\bf x}, {\bf y}) = \braket{{\bf x} | \hat \rho | {\bf y}} = \sum_{\alpha\beta} \braket{{\bf x}|\alpha}P_{\alpha\beta} \braket{\beta|{\bf y}} = \sum_{\alpha\beta} \phi_\alpha({\bf x}) P_{\alpha\beta} \phi_\beta^*({\bf y}) P_{\alpha\beta} = 2 \sum_{k=1}^{N/2} C_{\alpha k} C_{\beta k}^* \rho({\bf x}) = 2 \sum_{k=1}^{N/2} | \psi_k({\bf x})|^2 = 2 \sum_{k=1}^{N/2} | \braket{{\bf x}|k}|^2 = 2 \sum_{k=1}^{N/2} \braket{{\bf x}|k}\braket{k|{\bf x}} = \braket{{\bf x}|\hat \rho|{\bf x}} = \sum_{\alpha\beta} \phi_\alpha({\bf x}) P_{\alpha\beta} \phi_\beta^*({\bf x})

Expanding the \psi_k functions and using the density matrix we get for G_{\mu\nu}:

G_{\mu\nu} = \sum_{\alpha\beta} P_{\alpha\beta} \int \phi_\mu^*({\bf x}) \left( \int {\phi_\beta^*({\bf y})\phi_\alpha({\bf y})\over|{\bf x}-{\bf y}|} \d^3 y\right)\phi_\nu({\bf x}) \d^3 x -\half \sum_{\alpha\beta} P_{\alpha\beta} \int \phi_\mu^*({\bf x}) \int {\phi_\nu({\bf y})\phi_\beta^*({\bf y})\over|{\bf x}-{\bf y}|} \d^3 y\,\,\phi_\alpha({\bf x})\d^3 x

or

G_{\mu\nu} = \sum_{\alpha\beta} P_{\alpha\beta} \int {\phi_\mu^*({\bf x}) \phi_\nu({\bf x}) \phi_\beta^*({\bf y}) \phi_\alpha({\bf y}) -\half \phi_\mu^*({\bf x}) \phi_\alpha({\bf x}) \phi_\beta^*({\bf y}) \phi_\nu({\bf y}) \over | {\bf x}-{\bf y}|} \d^3 x\, \d^3 y \equiv \sum_{\alpha\beta} P_{\alpha\beta} \left(\braket{\mu \beta|{1\over r_{12}}|\nu \alpha} -\half \braket{\mu \beta|{1\over r_{12}}|\alpha \nu}\right)

In physical and chemistry notation this is written as:

G_{\mu\nu} = \sum_{\alpha\beta} P_{\alpha\beta} \left(\braket{\mu \beta|\nu \alpha} -\half \braket{\mu \beta|\alpha \nu}\right) = \sum_{\alpha\beta} P_{\alpha\beta} \left((\mu \nu|\beta \alpha) -\half (\mu \alpha|\beta \nu)\right)

Note that this notation implicitly assumes the {1\over r_{12}} factor, so for example \braket{\mu \beta|\nu \alpha} actually means \braket{\mu \beta|{1\over r_{12}}|\nu \alpha} and one has to understand this from the context.

Two Particle Matrix Element

The two particle matrix element is:

(10)(ij|kl) = \int {\psi_i^*({\bf x})\psi_j({\bf x})\psi_k^*({\bf x}')\psi_l({\bf x}') \over | {\bf x} - {\bf x}' |} \d^3 x \d^3 x' = = \braket{ik|jl} = \int {\psi_i^*({\bf x})\psi_k^*({\bf x}')\psi_j({\bf x})\psi_l({\bf x}') \over | {\bf x} - {\bf x}' |} \d^3 x \d^3 x'

The \braket{ik|jl} is called the physicists’ notation because the \ket{jl} and \ket{ik} kets are:

\ket{jl}=\psi_j({\bf x})\psi_l({\bf x}') \ket{ik}=\psi_i({\bf x})\psi_k({\bf x}')

The (ij|kl) is called the chemists’ notation. From (10) there are two types of symmetries — interchanging of the dummy variables:

(ij|kl) = \braket{ik|jl} = \int {\psi_i^*({\bf x})\psi_j({\bf x})\psi_k^*({\bf x}')\psi_l({\bf x}') \over | {\bf x} - {\bf x}' |} \d^3 x \d^3 x' = = \int {\psi_i^*({\bf x}')\psi_j({\bf x}') \psi_k^*({\bf x})\psi_l({\bf x}) \over | {\bf x}' - {\bf x} |} \d^3 x' \d^3 x = (kl|ij) = \braket{ki|lj}

and taking complex conjugate:

(ij|kl)^* = \braket{ik|jl}^* = \left( \int {\psi_i^*({\bf x})\psi_k^*({\bf x}')\psi_j({\bf x})\psi_l({\bf x}') \over | {\bf x} - {\bf x}' |} \d^3 x \d^3 x'\right)^* = = \int {\psi_i({\bf x})\psi_k({\bf x}')\psi_j^*({\bf x})\psi_l^*({\bf x}') \over | {\bf x} - {\bf x}' |} \d^3 x \d^3 x' = \braket{jl|ik} = (ji|lk)

If the matrix elements are real, then:

(ij|kl) = \braket{ik|jl} = \braket{jl|ik} = (ji|lk)

In general those are the only symmetries (4 total).

If however, the \psi_i({\bf x}) functions are real, then there are additional symmetries: an exchange i \leftrightarrow j and k \leftrightarrow l. The symmetries of (ij|kl) are exchange of i with j or k with l or the ij and kl pairs (8 total):

(ij|kl) = (ji|kl) = (ij|lk) = (ji|lk) = = (kl|ij) = (lk|ij) = (kl|ji) = (lk|ji)

So if we view (ij|kl) as two boxes (\cdot | \cdot ) then we can permute the labels in the given box “\cdot”, as well as exchange the boxes (the only thing we cannot do is to take one particle from one box and put it into the other). As such the box “\cdot” is a pair of two electrons (in any order) and the two electron integral assigns a unique number to a pair of such boxes (in any order). The symmetries of the \braket{ik|jl} symbol are:

\braket{ik|jl} = \braket{jk|il} = \braket{il|jk} = \braket{jl|ik} = = \braket{ki|lj} = \braket{li|kj} = \braket{kj|li} = \braket{lj|ki}

Example I: the Slater integral R^k(i, j, k, l) has all 8 symmetries (Slater integral uses physical notation).

Example II: In spherical symmetry, the 2-particle integrals can be written as (see below for derivation):

\braket{\alpha \beta | \gamma \delta} = \braket{n_\alpha l_\alpha m_\alpha \ n_\beta l_\beta m_\beta | n_\gamma l_\gamma m_\gamma \ n_\delta l_\delta m_\delta} = (n_\alpha l_\alpha m_\alpha \ n_\gamma l_\gamma m_\gamma | n_\beta l_\beta m_\beta \ n_\delta l_\delta m_\delta) = = \sum_{k=\max(| l_\alpha-l_\gamma| ,| l_\beta-l_\delta| , | m_\alpha-m_\gamma| )}^{ \min(l_\alpha+l_\gamma, l_\beta+l_\delta) }\!\!\!\!\!\!\!\!\!\!\!\! c^k(l_\alpha, m_\alpha, l_\gamma, m_\gamma) c^k(l_\delta, m_\delta, l_\beta, m_\beta) \delta_{m_\alpha+m_\beta- m_\gamma-m_\delta, 0} R^k(n_\alpha l_\alpha, n_\beta l_\beta, n_\gamma l_\gamma, n_\delta l_\delta)

They only have 4 symmetries, because spherical harmonics are complex. In particular:

\braket{\beta \alpha | \delta \gamma} = \sum_{k=\max(| l_\beta-l_\delta| ,| l_\alpha-l_\gamma| , | m_\beta-m_\delta| )}^{ \min(l_\beta+l_\delta, l_\alpha+l_\gamma) }\!\!\!\!\!\!\!\!\!\!\!\! c^k(l_\beta, m_\beta, l_\delta, m_\delta) c^k(l_\gamma, m_\gamma, l_\alpha, m_\alpha) \delta_{m_\beta+m_\alpha- m_\delta-m_\gamma, 0} R^k(n_\beta l_\beta, n_\alpha l_\alpha, n_\delta l_\delta, n_\gamma l_\gamma) = = \sum_{k=\max(| l_\beta-l_\delta|, | l_\alpha-l_\gamma|, | m_\beta-m_\delta| )}^{\min(l_\beta+l_\delta, l_\alpha+l_\gamma) }\!\!\!\!\!\!\!\!\!\!\!\! c^k(l_\alpha, m_\alpha, l_\gamma, m_\gamma) c^k(l_\delta, m_\delta, l_\beta, m_\beta) (-1)^{m_\alpha-m_\gamma + m_\beta-m_\delta} \delta_{m_\beta+m_\alpha- m_\delta-m_\gamma, 0} R^k(n_\alpha l_\alpha, n_\beta l_\beta, n_\gamma l_\gamma, n_\delta l_\delta) = \braket{\alpha \beta | \gamma \delta}

and

\braket{\gamma \delta | \alpha \beta} = \braket{n_\gamma l_\gamma m_\gamma \ n_\delta l_\delta m_\delta | n_\alpha l_\alpha m_\alpha \ n_\beta l_\beta m_\beta} = (n_\gamma l_\gamma m_\gamma \ n_\alpha l_\alpha m_\alpha | n_\delta l_\delta m_\delta \ n_\beta l_\beta m_\beta) = = \sum_{k=\max(| l_\gamma-l_\alpha| ,| l_\delta-l_\beta| , | m_\gamma-m_\alpha| )}^{ \min(l_\gamma+l_\alpha, l_\delta+l_\beta) }\!\!\!\!\!\!\!\!\!\!\!\! c^k(l_\gamma, m_\gamma, l_\alpha, m_\alpha) c^k(l_\beta, m_\beta, l_\delta, m_\delta) \delta_{m_\gamma+m_\delta- m_\alpha-m_\beta, 0} R^k(n_\gamma l_\gamma, n_\delta l_\delta, n_\alpha l_\alpha, n_\beta l_\beta) = \braket{\alpha \beta | \gamma \delta}

We used the symmetries of the Slater integrals as well as the c^k coefficients that change a sign, but thanks to the \delta_{m_\gamma+m_\delta- m_\alpha-m_\beta, 0}, the overall sign does not change. The other two symmetries are missing, i.e. \braket{\gamma \beta | \alpha \delta} \ne \braket{\alpha \beta | \gamma \delta} and \braket{\alpha \delta | \gamma \beta} \ne \braket{\alpha \beta | \gamma \delta}.

General Matrix Elements in Spherical Symmetry

Spherical symmetry is this particular choice of a basis:

\phi_\mu({\bf x}) = \phi_{n_\mu l_\mu m_\mu}({\bf x}) = {\phi_{n_\mu l_\mu}(r)\over r} Y_{l_\mu m_\mu}(\Omega)

It can be shown that the solutions are of the form:

\psi_i({\bf x}) = \psi_{nlm}({\bf x}) = {P_{n l}(r)\over r} Y_{l m}(\Omega)

We can now write:

\psi_i({\bf x}) = \sum_\nu C_{\nu i} \phi_\nu({\bf x}) \int \phi_\mu({\bf x})\psi_i({\bf x}) \d^3 x = \sum_\nu S_{\mu\nu} C_{\nu i} \sum_{\mu} S_{\mu\nu}^{-1} \int \phi_\mu({\bf x})\psi_i({\bf x}) \d^3 x = C_{\nu i} \sum_{\mu} \delta_{l_\mu l_\nu} \delta_{m_\mu m_\nu} \left(S_{n_\mu n_\nu}^{l_\mu}\right)^{-1} \int \phi_{n_\mu l_\mu m_\mu}({\bf x}) \psi_{nlm}({\bf x}) \d^3 x = C_{n_\nu l_\nu m_\nu; n l m} \sum_{\mu} \delta_{l_\mu l_\nu} \delta_{m_\mu m_\nu} \left(S_{n_\mu n_\nu}^{l_\mu}\right)^{-1} \int {\phi_{n_\mu l_\mu}(r)\over r} Y_{l_\mu m_\mu}(\Omega) {P_{n l}(r)\over r} Y_{l m}(\Omega) r^2 \d r \d \Omega = C_{n_\nu l_\nu m_\nu; n l m} \sum_{\mu} \delta_{l_\mu l_\nu} \delta_{m_\mu m_\nu} \left(S_{n_\mu n_\nu}^{l_\mu}\right)^{-1} \int \phi_{n_\mu l_\mu}(r) P_{n l}(r) \d r \delta_{l l_\mu} \delta_{m m_\mu} = C_{n_\nu l_\nu m_\nu; n l m} \delta_{l l_\nu} \delta_{m m_\nu} \sum_{n_\mu} \left(S_{n_\mu n_\nu}^{l}\right)^{-1} \int \phi_{n_\mu l}(r) P_{n l}(r) \d r = C_{n_\nu l_\nu m_\nu; n l m} \delta_{l l_\nu} \delta_{m m_\nu} C_{n_\nu n}^l = C_{n_\nu l_\nu m_\nu; n l m}

where

C_{n_\nu n}^l = \sum_{n_\mu} \left(S_{n_\mu n_\nu}^{l}\right)^{-1} \int \phi_{n_\mu l}(r) P_{n l}(r) \d r

Also we get:

\psi_{nlm}({\bf x}) = \sum_\nu C_{n_\nu l_\nu m_\nu; n l m} \phi_{n_\nu l_\nu m_\nu}({\bf x}) {P_{n l}(r)\over r} Y_{l m}(\Omega) = \sum_\nu C_{n_\nu l_\nu m_\nu; n l m} {\phi_{n_\nu l_\nu}(r)\over r} Y_{l_\nu m_\nu}(\Omega) = = \sum_\nu \delta_{l l_\nu} \delta_{m m_\nu} C_{n_\nu n}^l {\phi_{n_\nu l_\nu}(r)\over r} Y_{l_\nu m_\nu}(\Omega) = = \sum_{n_\nu} C_{n_\nu n}^l {\phi_{n_\nu l}(r)\over r} Y_{l m}(\Omega)

From which it follows:

P_{n l}(r) = \sum_{n_\nu} C_{n_\nu n}^l \phi_{n_\nu l}(r)

The \mu index runs over all combinations of n l m. In particular, here is an example of one possible way to index the basis of 12 radial functions for each l:

\begin{array}{r|rrr} \mu & n_\mu & l_\mu & m_\mu \\ \hline 1 & 1 & 0 & 0 \\ 2 & 2 & 0 & 0 \\ 3 & 3 & 0 & 0 \\ \cdots & \cdots & & \\ 12 & 12 & 0 & 0 \\ 13 & 1 & 1 & -1 \\ 14 & 2 & 1 & -1 \\ 15 & 3 & 1 & -1 \\ \cdots & \cdots & & \\ 24 & 12 & 1 & -1 \\ 25 & 1 & 1 & 0 \\ 26 & 2 & 1 & 0 \\ 27 & 3 & 1 & 0 \\ \cdots & \cdots & & \\ 36 & 12 & 1 & 0 \\ 37 & 1 & 1 & 1 \\ 38 & 2 & 1 & 1 \\ 39 & 3 & 1 & 1 \\ \cdots & \cdots & & \\ 48 & 12 & 1 & 1 \\ 49 & 1 & 2 & -2 \\ 50 & 2 & 2 & -2 \\ 51 & 3 & 2 & -2 \\ \cdots & \cdots & & \\ \end{array}

So the radial index n_\mu always starts from 1 for each l_\mu.

Overlap

The overlap matrix element

S_{\mu\nu} = \int \phi_\mu^*({\bf x}) \phi_\nu({\bf x})\d^3 x

becomes

S_{\mu\nu} = S_{n_\mu l_\mu m_\mu n_\nu l_\nu m_\nu} = \int {\phi_{n_\mu l_\mu}(r)\over r} Y_{l_\mu m_\mu}^*(\Omega) {\phi_{n_\nu l_\nu}(r)\over r} Y_{l_\nu m_\nu}(\Omega) r^2 \d r \d \Omega= = \delta_{l_\mu l_\nu} \delta_{m_\mu m_\nu} \int_0^\infty \phi_{n_\mu l_\mu}(r) \phi_{n_\nu l_\nu}(r) \d r = = \delta_{l_\mu l_\nu} \delta_{m_\mu m_\nu} \int_0^\infty \phi_{n_\mu l_\mu}(r) \phi_{n_\nu l_\mu}(r) \d r = = \delta_{l_\mu l_\nu} \delta_{m_\mu m_\nu} S^{l_\mu}_{n_\mu n_\nu}

where

S^l_{n_\mu n_\nu} = \int_0^\infty \phi_{n_\mu l}(r) \phi_{n_\nu l}(r) \d r

Potential

The potential matrix element

V_{\mu\nu} = \int \phi_\mu^*({\bf x}) \left(-{Z\over |{\bf x}|}\right) \phi_\nu({\bf x})\d^3 x

becomes

V_{\mu\nu} = V_{n_\mu l_\mu m_\mu n_\nu l_\nu m_\nu} = \int {\phi_{n_\mu l_\mu}(r) \over r} Y_{l_\mu m_\mu}^*(\Omega) \left(-{Z\over r}\right) {\phi_{n_\nu l_\nu}(r)\over r} Y_{l_\nu m_\nu}(\Omega) r^2 \d r \d \Omega= = \delta_{l_\mu l_\nu} \delta_{m_\mu m_\nu} \int_0^\infty \phi_{n_\mu l_\mu}(r) \left(-{Z\over r}\right) \phi_{n_\nu l_\nu}(r) \d r = = \delta_{l_\mu l_\nu} \delta_{m_\mu m_\nu} \int_0^\infty \phi_{n_\mu l_\mu}(r) \left(-{Z\over r}\right) \phi_{n_\nu l_\mu}(r) \d r = = \delta_{l_\mu l_\nu} \delta_{m_\mu m_\nu} V^{l_\mu}_{n_\mu n_\nu}

where

V^l_{n_\mu n_\nu} = \int_0^\infty \phi_{n_\mu l}(r) \left(-{Z\over r}\right) \phi_{n_\nu l}(r) \d r

Kinetic

The kinetic matrix element

T_{\mu\nu} = \int \phi_\mu^*({\bf x}) \left(-\half \nabla^2 \right) \phi_\nu({\bf x}) \d^3 x

becomes

T_{\mu\nu} = T_{n_\mu l_\mu m_\mu n_\nu l_\nu m_\nu} = \int {\phi_{n_\mu l_\mu}(r)\over r} Y_{l_\mu m_\mu}^*(\Omega) \left(\left(-\half \nabla^2\right) {\phi_{n_\nu l_\nu}(r)\over r} Y_{l_\nu m_\nu}(\Omega) \right) r^2 \d r \d \Omega = = \int {\phi_{n_\mu l_\mu}(r)\over r} Y_{l_\mu m_\mu}^*(\Omega) \left(\left(-\half {\partial^2\over\partial r^2} -{1\over r}{\partial\over\partial r} +{l_\nu (l_\nu+1)\over 2r^2}\right) {\phi_{n_\nu l_\nu}(r)\over r} Y_{l_\nu m_\nu}(\Omega) \right) r^2 \d r \d \Omega = = \delta_{l_\mu l_\nu} \delta_{m_\mu m_\nu} \int_0^\infty {\phi_{n_\mu l_\mu}(r)\over r} \left(\left(-\half {\partial^2\over\partial r^2} -{1\over r}{\partial\over\partial r} +{l_\mu (l_\mu+1)\over 2r^2}\right) {\phi_{n_\nu l_\nu}(r)\over r} \right) r^2 \d r = = \delta_{l_\mu l_\nu} \delta_{m_\mu m_\nu} T^{l_\mu}_{\mu\nu}

where

T^l_{n_\mu n_\nu} = \int_0^\infty {\phi_{n_\mu l}(r)\over r} \left(\left(-\half {\partial^2\over\partial r^2} -{1\over r}{\partial\over\partial r} +{l (l+1)\over 2r^2}\right) {\phi_{n_\nu l}(r)\over r} \right) r^2 \d r = = \int_0^\infty {\phi_{n_\mu l}(r)\over r} \left(\left(-{1\over 2r} {\partial^2\over\partial r^2}r +{l (l+1)\over 2 r^2}\right) {\phi_{n_\nu l}(r)\over r} \right) r^2 \d r = = \int_0^\infty \phi_{n_\mu l}(r) \left(-\half {\partial^2\over\partial r^2} +{l (l+1)\over 2 r^2}\right) \phi_{n_\nu l}(r) \d r = = \int_0^\infty \left( -\half \phi_{n_\mu l}(r) \phi_{n_\nu l}''(r) +\phi_{n_\mu l}(r) {l (l+1)\over 2 r^2} \phi_{n_\nu l}(r)\right) \d r = = \int_0^\infty \left( \half \phi_{n_\mu l}'(r) \phi_{n_\nu l}'(r) +\phi_{n_\mu l}(r) {l (l+1)\over 2 r^2} \phi_{n_\nu l}(r)\right) \d r

G Matrix

The G matrix is given by:

G_{\mu\nu} = \sum_{\alpha\beta} P_{\alpha\beta} \left(\braket{\mu \beta|\nu \alpha} -\half \braket{\mu \beta|\alpha \nu}\right) = \sum_{\alpha\beta} P_{\alpha\beta} \left((\mu \nu|\beta \alpha) -\half (\mu \alpha|\beta \nu)\right)

Now we use:

\braket{\alpha \beta | \gamma \delta} = \braket{n_\alpha l_\alpha m_\alpha \ n_\beta l_\beta m_\beta | n_\gamma l_\gamma m_\gamma \ n_\delta l_\delta m_\delta} = (n_\alpha l_\alpha m_\alpha \ n_\gamma l_\gamma m_\gamma | n_\beta l_\beta m_\beta \ n_\delta l_\delta m_\delta) = = \sum_{k=\max(| l_\alpha-l_\gamma| ,| l_\beta-l_\delta| , | m_\alpha-m_\gamma| )}^{ \min(l_\alpha+l_\gamma, l_\beta+l_\delta) }\!\!\!\!\!\!\!\!\!\!\!\! c^k(l_\alpha, m_\alpha, l_\gamma, m_\gamma) c^k(l_\delta, m_\delta, l_\beta, m_\beta) \delta_{m_\alpha+m_\beta- m_\gamma-m_\delta, 0} R^k(n_\alpha l_\alpha, n_\beta l_\beta, n_\gamma l_\gamma, n_\delta l_\delta)

and

P_{\alpha\beta} = 2 \sum_{i} C_{\alpha i} C_{\beta i} P_{n_\alpha l_\alpha m_\alpha; n_\beta l_\beta m_\beta} = 2 \sum_{i} C_{n_\alpha l_\alpha m_\alpha; n_i l_i m_i} C_{n_\beta l_\beta m_\beta; n_i l_i m_i} = = 2 \sum_{i} \delta_{l_i l_\alpha} \delta_{m_i m_\alpha} \delta_{l_i l_\beta} \delta_{m_i m_\beta} C_{n_\alpha n_i}^{l_i} C_{n_\beta n_i}^{l_i} = = 2 \delta_{l_\alpha l_\beta}\delta_{m_\alpha m_\beta} \sum_{n_i} C_{n_\alpha n_i}^{l_\alpha} C_{n_\beta n_i}^{l_\alpha} = = \delta_{l_\alpha l_\beta}\delta_{m_\alpha m_\beta} P^{l_\alpha}_{n_\alpha n_\beta}

where

P^l_{n_\alpha n_\beta} = 2 \sum_{n_i} C_{n_\alpha n_i}^l C_{n_\beta n_i}^l

Where the sum over all occupied orbitals i can be written as:

\sum_i = \sum_{n_i l_i m_i} = \sum_{l=0}^\infty \sum_{m_l=-l}^l \sum_{n_l=1}^\infty

Where the sum over l is the outer sum, both m=m_l and n=n_l depend on l. We get:

G_{\mu\nu} = G_{n_\mu l_\mu m_\mu n_\nu l_\nu m_\nu} = \sum_{n_\alpha l_\alpha m_\alpha n_\beta l_\beta m_\beta} P_{n_\alpha l_\alpha m_\alpha n_\beta l_\beta m_\beta} \left(\braket{\mu \beta|\nu \alpha} -\half \braket{\mu \beta|\alpha \nu}\right) = J_{\mu\nu} - K_{\mu\nu}

The first part is the direct term, the second part the exchange term. Let’s treat the direct term first:

J_{\mu\nu} = \sum_{n_\alpha l_\alpha m_\alpha n_\beta l_\beta m_\beta} P_{n_\alpha l_\alpha m_\alpha n_\beta l_\beta m_\beta} \braket{\mu \beta|\nu \alpha} = = \sum_{n_\alpha l_\alpha m_\alpha n_\beta l_\beta m_\beta} P_{n_\alpha l_\alpha m_\alpha n_\beta l_\beta m_\beta} \sum_{k=\max(| l_\mu-l_\nu| ,| l_\beta-l_\alpha|, | m_\mu-m_\nu| )}^{ \min(l_\mu+l_\nu, l_\beta+l_\alpha) }\!\!\!\!\!\!\!\!\!\!\!\! c^k(l_\mu, m_\mu, l_\nu, m_\nu) c^k(l_\alpha, m_\alpha, l_\beta, m_\beta) \delta_{m_\mu+m_\beta- m_\nu-m_\alpha, 0} R^k(n_\mu l_\mu, n_\beta l_\beta, n_\nu l_\nu, n_\alpha l_\alpha) = = \sum_{n_\alpha l_\alpha m_\alpha n_\beta l_\beta m_\beta} \delta_{l_\alpha l_\beta} \delta_{m_\alpha m_\beta} P^{l_\alpha}_{n_\alpha n_\beta} \sum_{k=\max(| l_\mu-l_\nu| ,| l_\beta-l_\alpha|, | m_\mu-m_\nu| )}^{ \min(l_\mu+l_\nu, l_\beta+l_\alpha) }\!\!\!\!\!\!\!\!\!\!\!\! c^k(l_\mu, m_\mu, l_\nu, m_\nu) c^k(l_\alpha, m_\alpha, l_\beta, m_\beta) \delta_{m_\mu+m_\beta- m_\nu-m_\alpha, 0} R^k(n_\mu l_\mu, n_\beta l_\beta, n_\nu l_\nu, n_\alpha l_\alpha) = = \sum_{n_\alpha n_\beta} \sum_{l m} P^{l}_{n_\alpha n_\beta} \sum_{k=\max(| l_\mu-l_\nu| , | m_\mu-m_\nu| )}^{ \min(l_\mu+l_\nu, 2l) }\!\!\!\!\!\!\!\!\!\!\!\! c^k(l_\mu, m_\mu, l_\nu, m_\nu) c^k(l, m, l, m) \delta_{m_\mu- m_\nu, 0} R^k(n_\mu l_\mu, n_\beta l, n_\nu l_\nu, n_\alpha l) = = \delta_{m_\mu m_\nu} \sum_{n_\alpha n_\beta} \sum_{l} P^{l}_{n_\alpha n_\beta} c^0(l_\mu, m_\mu, l_\nu, m_\nu) (2l+1) R^0(n_\mu l_\mu, n_\beta l, n_\nu l_\nu, n_\alpha l) = = \delta_{l_\mu l_\nu} \delta_{m_\mu m_\nu} \sum_{l} (2l+1)\sum_{n_\alpha n_\beta} P^{l}_{n_\alpha n_\beta} R^0(n_\mu l_\mu, n_\beta l, n_\nu l_\mu, n_\alpha l)

For the exchange term we get:

K_{\mu\nu} = \half \sum_{n_\alpha l_\alpha m_\alpha n_\beta l_\beta m_\beta} P_{n_\alpha l_\alpha m_\alpha n_\beta l_\beta m_\beta} \braket{\mu \beta|\alpha \nu} = = \half \sum_{n_\alpha l_\alpha m_\alpha n_\beta l_\beta m_\beta} P_{n_\alpha l_\alpha m_\alpha n_\beta l_\beta m_\beta} \sum_{k=\max(| l_\mu-l_\alpha| ,| l_\beta-l_\nu|, | m_\mu-m_\alpha| )}^{ \min(l_\mu+l_\alpha, l_\beta+l_\nu) }\!\!\!\!\!\!\!\!\!\!\!\! c^k(l_\mu, m_\mu, l_\alpha, m_\alpha) c^k(l_\nu, m_\nu, l_\beta, m_\beta) \delta_{m_\mu+m_\beta- m_\alpha-m_\nu, 0} R^k(n_\mu l_\mu, n_\beta l_\beta, n_\alpha l_\alpha, n_\nu l_\nu) = = \half \sum_{n_\alpha l_\alpha m_\alpha n_\beta l_\beta m_\beta} \delta_{l_\alpha l_\beta} \delta_{m_\alpha m_\beta} P^{l_\alpha}_{n_\alpha n_\beta} \sum_{k=\max(| l_\mu-l_\alpha| ,| l_\beta-l_\nu|, | m_\mu-m_\alpha| )}^{ \min(l_\mu+l_\alpha, l_\beta+l_\nu) }\!\!\!\!\!\!\!\!\!\!\!\! c^k(l_\mu, m_\mu, l_\alpha, m_\alpha) c^k(l_\nu, m_\nu, l_\beta, m_\beta) \delta_{m_\mu+m_\beta- m_\alpha-m_\nu, 0} R^k(n_\mu l_\mu, n_\beta l_\beta, n_\alpha l_\alpha, n_\nu l_\nu) = = \half \delta_{m_\mu m_\nu} \sum_{n_\alpha n_\beta} \sum_{l m} P^l_{n_\alpha n_\beta} \sum_{k=\max(| l_\mu-l| ,| l-l_\nu|, | m_\mu-m| )}^{ \min(l_\mu+l, l+l_\nu) }\!\!\!\!\!\!\!\!\!\!\!\! c^k(l_\mu, m_\mu, l, m) c^k(l_\nu, m_\nu, l, m) R^k(n_\mu l_\mu, n_\beta l, n_\alpha l, n_\nu l_\nu) = = \half \delta_{m_\mu m_\nu} \sum_{n_\alpha n_\beta} \sum_{l m} P^l_{n_\alpha n_\beta} \sum_{k=\max(| l_\mu-l| ,| l-l_\nu|, | m_\mu-m| )}^{ \min(l_\mu+l, l+l_\nu) }\!\!\!\!\!\!\!\!\!\!\!\! c^k(l_\mu, m_\mu, l, m) c^k(l_\nu, m_\mu, l, m) R^k(n_\mu l_\mu, n_\beta l, n_\alpha l, n_\nu l_\nu)

For l_\mu = l_\nu this can be written as:

K_{\mu\nu} = \half \delta_{l_\mu l_\nu} \delta_{m_\mu m_\nu} \sum_{n_\alpha n_\beta} \sum_{l m} P^l_{n_\alpha n_\beta} \sum_{k=\max(| l_\mu-l|, | m_\mu-m| )}^{l_\mu+l}\!\!\!\!\!\!\!\!\!\!\!\! c^k(l_\mu, m_\mu, l, m) c^k(l_\mu, m_\mu, l, m) R^k(n_\mu l_\mu, n_\beta l, n_\alpha l, n_\nu l_\mu) = = \half \delta_{l_\mu l_\nu} \delta_{m_\mu m_\nu} \sum_{n_\alpha n_\beta} \sum_{l} P^l_{n_\alpha n_\beta} \sum_{k=| l_\mu-l|}^{l_\mu+l} \sqrt{2l+1 \over 2l_\mu +1} c^k(l_\mu, 0, l, 0) R^k(n_\mu l_\mu, n_\beta l, n_\alpha l, n_\nu l_\mu) = = \delta_{l_\mu l_\nu} \delta_{m_\mu m_\nu} \sum_{l} (2l+1) \sum_{n_\alpha n_\beta} P^l_{n_\alpha n_\beta} \sum_{k=| l_\mu-l|}^{l_\mu+l} \half \begin{pmatrix} l_\mu & k & l \\ 0 & 0 & 0 \end{pmatrix}^2 R^k(n_\mu l_\mu, n_\beta l, n_\alpha l, n_\nu l_\mu)

All together we get:

G_{\mu\nu} = \delta_{l_\mu l_\nu} \delta_{m_\mu m_\nu} \sum_{l} (2l+1)\sum_{n_\alpha n_\beta} P^{l}_{n_\alpha n_\beta} \left( R^0(n_\mu l_\mu, n_\beta l, n_\nu l_\mu, n_\alpha l) - \sum_{k=| l_\mu-l|}^{l_\mu+l} \half \begin{pmatrix} l_\mu & k & l \\ 0 & 0 & 0 \end{pmatrix}^2 R^k(n_\mu l_\mu, n_\beta l, n_\alpha l, n_\nu l_\mu) \right) = = \delta_{l_\mu l_\nu} \delta_{m_\mu m_\nu} G^{l_\mu}_{n_\mu n_\nu}

where

G^{l}_{n_\mu n_\nu} = \sum_{l'} (2l'+1)\sum_{n_\alpha n_\beta} P^{l'}_{n_\alpha n_\beta} \left( R^0(n_\mu l, n_\beta l', n_\nu l, n_\alpha l') - \sum_{k=| l-l'| }^{l+l'} \half \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 R^k(n_\mu l, n_\beta l', n_\alpha l', n_\nu l) \right)

Note: performing the sum over n_\alpha and n_\beta we get:

G^{l}_{n_\mu n_\nu} = \sum_{l'} 2(2l'+1)\sum_{n'}\sum_{n_\alpha n_\beta} C^{l'}_{n_\alpha n'}C^{l'}_{n_\beta n'} \left( R^0(n_\mu l, n_\beta l', n_\nu l, n_\alpha l') - \sum_{k=| l-l'| }^{l+l'} \half \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 R^k(n_\mu l, n_\beta l', n_\alpha l', n_\nu l) \right) = = \sum_{l'} 2(2l'+1)\sum_{n'} \left( R^0(n_\mu l, n' l', n_\nu l, n' l') - \sum_{k=| l-l'| }^{l+l'} \half \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 R^k(n_\mu l, n' l', n' l', n_\nu l) \right)

Radial Roothaan Equations

The Roothaan equations are:

\sum_\nu (T + V + G)_{\mu\nu} C_{\nu i} = \epsilon_i \sum_\nu S_{\mu\nu} C_{\nu i} \sum_\nu (T + V + G)_{n_\mu l_\mu m_\mu n_\nu l_\nu m_\nu} C_{n_\nu l_\nu m_\nu n_i l_i m_i} = \epsilon_{n_i l_i m_i} \sum_\nu S_{n_\mu l_\mu m_\mu n_\nu l_\nu m_\nu} C_{n_\nu l_\nu m_\nu n_i l_i m_i} \sum_\nu (T + V + G)_{n_\mu l_\mu m_\mu n_\nu l_\nu m_\nu} \delta_{l_i l_\nu} \delta_{m_i m_\nu} C_{n_\nu n_i}^{l_i} = \epsilon_{n_i l_i m_i} \sum_\nu S_{n_\mu l_\mu m_\mu n_\nu l_\nu m_\nu} \delta_{l_i l_\nu} \delta_{m_i m_\nu} C_{n_\nu n_i}^{l_i} \sum_{n_\nu} (T + V + G)_{n_\mu l_\mu m_\mu n_\nu l_i m_i} C_{n_\nu n_i}^{l_i} = \epsilon_{n_i l_i m_i} \sum_{n_\nu} S_{n_\mu l_\mu m_\mu n_\nu l_i m_i} C_{n_\nu n_i}^{l_i} \delta_{l_\mu l_i} \delta_{m_\mu m_i} \sum_{n_\nu} (T + V + G)_{n_\mu n_\nu}^{l_i} C_{n_\nu n_i}^{l_i} = \delta_{l_\mu l_i} \delta_{m_\mu m_i} \epsilon_{n_i l_i m_i} \sum_{n_\nu} S_{n_\mu n_\nu}^{l_i} C_{n_\nu n_i}^{l_i}

The eigenvalues will be degenerate with respect to m_i and so the radial Roothaan equations are to be solved for each l:

\sum_{n_\nu} (T + V + G)_{n_\mu n_\nu}^l C_{n_\nu n_i}^l = \epsilon_{n_i l} \sum_{n_\nu} S_{n_\mu n_\nu}^l C_{n_\nu n_i}^l

The total energy is:

E = \half \sum_{\mu\nu} P_{\nu\mu} (H_{\mu\nu}^{\mbox{core}} +F_{\mu\nu}) = = \sum_{\mu\nu} P_{\nu\mu} (F_{\mu\nu} - \half G_{\mu\nu}) = = \sum_{\mu\nu} \delta_{l_\mu l_\nu} \delta_{m_\mu m_\nu} P_{n_\mu n_\nu}^{l_\mu} \delta_{l_\mu l_\nu} \delta_{m_\mu m_\nu} (F - \half G)^{l_\mu}_{n_\mu n_\nu} = = \sum_{l_\mu}\sum_{m_\mu}\sum_{n_\mu n_\nu} P_{n_\mu n_\nu}^{l_\mu} (F - \half G)^{l_\mu}_{n_\mu n_\nu} = = \sum_{l_\mu} \sum_{n_\mu n_\nu} (2l_\mu+1) P_{n_\mu n_\nu}^{l_\mu} (F - \half G)^{l_\mu}_{n_\mu n_\nu} = = \sum_l \sum_{n_\mu n_\nu} (2l+1) P_{n_\mu n_\nu}^l (F - \half G)^l_{n_\mu n_\nu} = = \sum_l \sum_{n_\mu n_\nu} (2l+1) \sum_n 2 C^l_{n_\mu n} C^l_{n_\nu n} (F^l_{n_\mu n_\nu} - \half G^l_{n_\mu n_\nu}) = = \sum_l \sum_{n_\mu n_\nu} (2l+1) \sum_n 2 C^l_{n_\mu n} C^l_{n_\nu n} \left(\epsilon_{n l} S^l_{n_\mu n_\nu} - \half \sum_{l'} 2(2l'+1)\sum_{n'} \right. \left. \left( R^0(n_\mu l, n' l', n_\nu l, n' l') - \sum_{k=| l-l'| }^{l+l'} \half \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 R^k(n_\mu l, n' l', n' l', n_\nu l) \right) \right) = = \sum_l \sum_n 2(2l+1) \left( \epsilon_{n l} -\sum_{l'} \sum_{n'} (2l'+1) \left( R^0(n l, n' l', n l, n' l') -\half \sum_{k=|l-l'|}^{k=l+l'} \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 R^k(n l, n' l', n' l', n l) \right) \right)

Two particle

We use the following functions for \psi:

\psi_i({\bf x}) = {P_{n_1l_1}(r)\over r} Y_{l_1m_1}(\Omega) \psi_j({\bf x}) = {P_{n_1'l_1'}(r)\over r} Y_{l_1'm_1'}(\Omega) \psi_k({\bf x}) = {P_{n_2l_2}(r)\over r} Y_{l_2m_2}(\Omega) \psi_l({\bf x}) = {P_{n_2'l_2'}(r)\over r} Y_{l_2'm_2'}(\Omega)

And the multipole expansion:

{1\over |{\bf x}-{\bf x}'|} = \sum_{k,q}{r_{<}^k\over r_{>}^{k+1}} {4\pi\over 2k+1}Y_{kq}(\Omega)Y_{kq}^*(\Omega')

And we get:

(ij|kl) = \braket{ik|jl} = (1 1' | 2 2') = \braket{1 2 | 1' 2'} = = \braket{l_1 m_1 l_2 m_2 | {1\over |{\bf x} - {\bf x}'|} | l_1' m_1' l_2' m_2'} = =\int {\psi_i^*({\bf x})\psi_j({\bf x})\psi_k^*({\bf x}')\psi_l({\bf x}') \over |{\bf x} - {\bf x}'|} \d^3 x \d^3 x' = = \int {P_{n_1l_1}(r)\over r} Y_{l_1m_1}^*(\Omega) {P_{n_1'l_1'}(r)\over r} Y_{l_1'm_1'}(\Omega) {P_{n_2l_2}(r')\over r'} Y_{l_2m_2}^*(\Omega') {P_{n_2'l_2'}(r')\over r'} Y_{l_2'm_2'}(\Omega') \sum_{k,q}{r_{<}^k\over r_{>}^{k+1}} {4\pi\over 2k+1}Y_{kq}(\Omega)Y_{kq}^*(\Omega') r^2 r'^2 \d r \d r' \d \Omega \d \Omega' = = \sum_{k,q} \int Y_{l_1m_1}^*(\Omega) Y_{l_1'm_1'}(\Omega) Y_{kq}(\Omega) \d \Omega \int Y_{l_2m_2}^*(\Omega') Y_{l_2'm_2'}(\Omega') Y_{kq}^*(\Omega') \d \Omega' \int {r_{<}^k\over r_{>}^{k+1}} {4\pi\over 2k+1} P_{n_1l_1}(r) P_{n_1'l_1'}(r) P_{n_2l_2}(r') P_{n_2'l_2'}(r') \d r \d r' = = \sum_{k,q} (-1)^{m_1+m_2+q} \int Y_{l_1,-m_1}(\Omega) Y_{l_1'm_1'}(\Omega) Y_{kq}(\Omega) \d \Omega \int Y_{l_2,-m_2}(\Omega') Y_{l_2'm_2'}(\Omega') Y_{k,-q}(\Omega') \d \Omega' \int {r_{<}^k\over r_{>}^{k+1}} {4\pi\over 2k+1} P_{n_1l_1}(r) P_{n_1'l_1'}(r) P_{n_2l_2}(r') P_{n_2'l_2'}(r') \d r \d r' = = \sum_{k,q} (-1)^{m_1+m_2+q} \sqrt{(2l_1+1)(2l_1'+1)(2k+1)\over 4\pi} \begin{pmatrix} l_1 & l_1' & k \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_1 & l_1' & k \\ -m_1 & m_1' & q \end{pmatrix} \sqrt{(2l_2+1)(2l_2'+1)(2k+1)\over 4\pi} \begin{pmatrix} l_2 & l_2' & k \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_2 & l_2' & k \\ -m_2 & m_2' & -q \end{pmatrix} \int {r_{<}^k\over r_{>}^{k+1}} {4\pi\over 2k+1} P_{n_1l_1}(r) P_{n_1'l_1'}(r) P_{n_2l_2}(r') P_{n_2'l_2'}(r') \d r \d r' = = \sum_k \sqrt{(2l_1+1)(2l_1'+1)(2l_2+1)(2l_2'+1)} \begin{pmatrix} l_1 & l_1' & k \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_2 & l_2' & k \\ 0 & 0 & 0 \end{pmatrix} \sum_{q=-k}^k (-1)^{m_1+m_2+q} \begin{pmatrix} l_1 & l_1' & k \\ -m_1 & m_1' & q \end{pmatrix} \begin{pmatrix} l_2 & l_2' & k \\ -m_2 & m_2' & -q \end{pmatrix} \int {r_{<}^k\over r_{>}^{k+1}} P_{n_1l_1}(r) P_{n_1'l_1'}(r) P_{n_2l_2}(r') P_{n_2'l_2'}(r') \d r \d r' = = \sum_{k=\max(| l_1-l_1'| ,| l_2-l_2'| , | m_1-m_1'| )}^{ \min(l_1+l_1', l_2+l_2') }\!\!\!\!\!\!\!\!\!\!\!\! \sqrt{(2l_1+1)(2l_1'+1)(2l_2+1)(2l_2'+1)} (-1)^{m_1+m_2'} \delta_{m_1+m_2- m_1'-m_2', 0} \begin{pmatrix} l_1 & l_1' & k \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_2 & l_2' & k \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_1 & l_1' & k \\ -m_1 & m_1' & m_1-m_1' \end{pmatrix} \begin{pmatrix} l_2 & l_2' & k \\ -m_2 & m_2' & m_2-m_2' \end{pmatrix} \int {r_{<}^k\over r_{>}^{k+1}} P_{n_1l_1}(r) P_{n_1'l_1'}(r) P_{n_2l_2}(r') P_{n_2'l_2'}(r') \d r \d r'

In the last step we used the fact that the 3j symbols are zero unless -m_1+m_1'+q=0 and -m_2+m_2'-q=0, from which it follows that q=m_1-m_1'=-m_2+m_2' and so one of the 3j symbols is zero unless m_1+m_2-m_1'-m_2'=0, which is expressed by \delta_{m_1+m_2- m_1'-m_2', 0}. Given this condition, the sum over q must be such that one q is equal to m_1-m_1'=-m_2+m_2', which means that k \ge |m_1 - m_1'| = |m_2 - m_2'| otherwise the 3j symbols will be zero. Finally, k must also satisfy the conditions |l_1-l_1'| \le k \le l_1+l_1' and |l_2-l_2'| \le k \le l_2+l_2'. The sign factor (-1)^{m_1+m_2+q} = (-1)^{m_1+m_2+m_1-m_1'} =(-1)^{m_1+m_2-m_2+m_2'} is equal to both (-1)^{m_1+m_2'} and (-1)^{m_2-m_1'} so we just used the former.

We can write this using the c^k symbols as:

(ij|kl) = \sum_{k=\max(| l_1-l_1'| ,| l_2-l_2'| , | m_1-m_1'| )}^{ \min(l_1+l_1', l_2+l_2') }\!\!\!\!\!\!\!\!\!\!\!\! \sqrt{(2l_1+1)(2l_1'+1)(2l_2+1)(2l_2'+1)} (-1)^{m_2-m_2'} (-1)^{-m_1} (-1)^{-m_2} \delta_{m_1+m_2- m_1'-m_2', 0} \begin{pmatrix} l_1 & l_1' & k \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_2 & l_2' & k \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_1 & l_1' & k \\ -m_1 & m_1' & m_1-m_1' \end{pmatrix} \begin{pmatrix} l_2 & l_2' & k \\ -m_2 & m_2' & m_2-m_2' \end{pmatrix} \int {r_{<}^k\over r_{>}^{k+1}} P_{n_1l_1}(r) P_{n_1'l_1'}(r) P_{n_2l_2}(r') P_{n_2'l_2'}(r') \d r \d r' = = \sum_{k=\max(| l_1-l_1'| ,| l_2-l_2'| , | m_1-m_1'| )}^{ \min(l_1+l_1', l_2+l_2') }\!\!\!\!\!\!\!\!\!\!\!\! c^k(l_1, m_1, l_1', m_1') c^k(l_2, m_2, l_2', m_2') (-1)^{m_2-m_2'} \delta_{m_1+m_2- m_1'-m_2', 0} \int {r_{<}^k\over r_{>}^{k+1}} P_{n_1l_1}(r) P_{n_1'l_1'}(r) P_{n_2l_2}(r') P_{n_2'l_2'}(r') \d r \d r' = = \sum_{k=\max(| l_1-l_1'| ,| l_2-l_2'| , | m_1-m_1'| )}^{ \min(l_1+l_1', l_2+l_2') }\!\!\!\!\!\!\!\!\!\!\!\! c^k(l_1, m_1, l_1', m_1') c^k(l_2', m_2', l_2, m_2) \delta_{m_1+m_2- m_1'-m_2', 0} \int {r_{<}^k\over r_{>}^{k+1}} P_{n_1l_1}(r) P_{n_1'l_1'}(r) P_{n_2l_2}(r') P_{n_2'l_2'}(r') \d r \d r' = = \sum_{k=\max(| l_1-l_1'| ,| l_2-l_2'| , | m_1-m_1'| )}^{ \min(l_1+l_1', l_2+l_2') }\!\!\!\!\!\!\!\!\!\!\!\! c^k(l_1, m_1, l_1', m_1') c^k(l_2', m_2', l_2, m_2) \delta_{m_1+m_2- m_1'-m_2', 0} R^k(n_1 l_1, n_2 l_2, n_1' l_1', n_2' l_2')

We can also couple the angular momenta as follows:

\ket{l_1 l_2 LM} = \sum_{m_1 m_2} (l_1 m_1 l_2 m_2 | LM) \ket{l_1 m_1} \ket{l_2 m_2} = = \sum_{m_1 m_2} (-1)^{l_1-l_2+M}\sqrt{2L+1} \begin{pmatrix} l_1 & l_2 & L \\ m_1 & m_2 & -M \end{pmatrix} \ket{l_1 m_1} \ket{l_2 m_2}

and we get for the matrix elements:

\braket{l_1 l_2 LM | {1\over |{\bf x} - {\bf x}'|} | l_1' l_2' L' M'} = = \sum_{m_1 m_2} \sum_{m_1' m_2'} (-1)^{l_1-l_2+l_1'-l_2'+M+M'} \sqrt{(2L+1)(2L'+1)} \begin{pmatrix} l_1 & l_2 & L \\ m_1 & m_2 & -M \end{pmatrix} \begin{pmatrix} l_1' & l_2' & L' \\ m_1' & m_2' & -M' \end{pmatrix} \bra{l_1 m_1} \bra{l_2 m_2} {1\over |{\bf x} - {\bf x}'|} \ket{l_1' m_1'} \ket{l_2 m_2'} = = \sum_{m_1 m_2} \sum_{m_1' m_2'} (-1)^{l_1-l_2+l_1'-l_2'+M+M'} \sqrt{(2L+1)(2L'+1)} \begin{pmatrix} l_1 & l_2 & L \\ m_1 & m_2 & -M \end{pmatrix} \begin{pmatrix} l_1' & l_2' & L' \\ m_1' & m_2' & -M' \end{pmatrix} (ij|kl) = = \sum_{m_1 m_2} \sum_{m_1' m_2'} (-1)^{l_1-l_2+l_1'-l_2'+M+M'} \sqrt{(2L+1)(2L'+1)} \begin{pmatrix} l_1 & l_2 & L \\ m_1 & m_2 & -M \end{pmatrix} \begin{pmatrix} l_1' & l_2' & L' \\ m_1' & m_2' & -M' \end{pmatrix} \sum_k \sqrt{(2l_1+1)(2l_1'+1)(2l_2+1)(2l_2'+1)} \begin{pmatrix} l_1 & l_1' & k \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_2 & l_2' & k \\ 0 & 0 & 0 \end{pmatrix} \sum_{q=-k}^k (-1)^{m_1+m_2+q} \begin{pmatrix} l_1 & l_1' & k \\ -m_1 & m_1' & q \end{pmatrix} \begin{pmatrix} l_2 & l_2' & k \\ -m_2 & m_2' & -q \end{pmatrix} \int {r_{<}^k\over r_{>}^{k+1}} P_{n_1l_1}(r) P_{n_1'l_1'}(r) P_{n_2l_2}(r') P_{n_2'l_2'}(r') \d r \d r' = = \sum_{m_1 m_2} \sum_{m_1' m_2'} (-1)^{l_1-l_2+l_1'-l_2'} (2L+1) \delta_{MM'}\delta_{LL'} \begin{pmatrix} l_1 & l_2 & L \\ m_1 & m_2 & -M \end{pmatrix} \begin{pmatrix} l_1' & l_2' & L \\ m_1' & m_2' & -M \end{pmatrix} \sum_k \sqrt{(2l_1+1)(2l_1'+1)(2l_2+1)(2l_2'+1)} \begin{pmatrix} l_1 & l_1' & k \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_2 & l_2' & k \\ 0 & 0 & 0 \end{pmatrix} \sum_{q=-k}^k (-1)^{m_1+m_2+q} \begin{pmatrix} l_1 & l_1' & k \\ -m_1 & m_1' & q \end{pmatrix} \begin{pmatrix} l_2 & l_2' & k \\ -m_2 & m_2' & -q \end{pmatrix} \int {r_{<}^k\over r_{>}^{k+1}} P_{n_1l_1}(r) P_{n_1'l_1'}(r) P_{n_2l_2}(r') P_{n_2'l_2'}(r') \d r \d r' = = (-1)^{l_1-l_2+l_1'-l_2'} (2L+1) \sum_k \sqrt{(2l_1+1)(2l_1'+1)(2l_2+1)(2l_2'+1)} \begin{pmatrix} l_1 & l_1' & k \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_2 & l_2' & k \\ 0 & 0 & 0 \end{pmatrix} \delta_{MM'}\delta_{LL'} (-1)^{l_1+l_1'+L} \begin{Bmatrix} l_1 & l_2 & L \\ l_2' & l_1' & k \end{Bmatrix} \int {r_{<}^k\over r_{>}^{k+1}} P_{n_1l_1}(r) P_{n_1'l_1'}(r) P_{n_2l_2}(r') P_{n_2'l_2'}(r') \d r \d r' = = \sum_k \int {r_{<}^k\over r_{>}^{k+1}} P_{n_1l_1}(r) P_{n_1'l_1'}(r) P_{n_2l_2}(r') P_{n_2'l_2'}(r') \d r \d r' (-1)^{L-l_2-l_2'} (2L+1) \delta_{MM'}\delta_{LL'}\sqrt{(2l_1+1)(2l_1'+1)(2l_2+1)(2l_2'+1)} \begin{pmatrix} l_1 & l_1' & k \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_2 & l_2' & k \\ 0 & 0 & 0 \end{pmatrix} \begin{Bmatrix} l_1 & l_2 & L \\ l_2' & l_1' & k \end{Bmatrix} = = \sum_k \int {r_{<}^k\over r_{>}^{k+1}} P_{n_1l_1}(r) P_{n_1'l_1'}(r) P_{n_2l_2}(r') P_{n_2'l_2'}(r') \d r \d r' (-1)^{l_1 + l_1' + L} (2L+1) \delta_{LL'}\delta_{MM'}\sqrt{(2l_1+1)(2l_1'+1)(2l_2+1)(2l_2'+1)} \begin{pmatrix} l_1 & k & l_1' \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_2 & k & l_2' \\ 0 & 0 & 0 \end{pmatrix} \begin{Bmatrix} l_1 & l_2 & L \\ l_2' & l_1' & k \end{Bmatrix}

Where we used the 6j symbol:

\begin{Bmatrix} l_1 & l_2 & L \\ l_2' & l_1' & k \end{Bmatrix} =\sum_{m_1 m_2 m_1' m_2' M q} (-1)^{l_1+l_2+l_1'+l_2'+L+k -m_1-m_2-m_1'-m_2'-M-q} \begin{pmatrix} l_1 & l_2 & L \\ m_1 & m_2 & -M \end{pmatrix} \begin{pmatrix} l_1 & l_1' & k \\ -m_1 & m_1' & q \end{pmatrix} \begin{pmatrix} l_2' & l_1' & L \\ m_2' & -m_1' & M \end{pmatrix} \begin{pmatrix} l_2' & l_2 & k \\ -m_2' & -m_2 & -q \end{pmatrix} = =\sum_{m_1 m_2 m_1' m_2' M q} (-1)^{l_1+l_2+l_1'+l_2'+L+k -m_1-m_2-m_1'-m_2'-M-q} \begin{pmatrix} l_1 & l_2 & L \\ m_1 & m_2 & -M \end{pmatrix} \begin{pmatrix} l_1 & l_1' & k \\ -m_1 & m_1' & q \end{pmatrix} \begin{pmatrix} l_1' & l_2' & L \\ m_1' & -m_2' & -M \end{pmatrix} (-1)^{l_2+l_2'+k} \begin{pmatrix} l_2 & l_2' & k \\ -m_2 & -m_2' & -q \end{pmatrix} = =\sum_{m_1 m_2 m_1' m_2' q} \delta_{M, m_1' + m_2'} (-1)^{l_1+l_1'+L} (-1)^{m_1+m_2+q} \begin{pmatrix} l_1 & l_2 & L \\ m_1 & m_2 & -M \end{pmatrix} \begin{pmatrix} l_1 & l_1' & k \\ -m_1 & m_1' & q \end{pmatrix} \begin{pmatrix} l_1' & l_2' & L \\ m_1' & +m_2' & -M \end{pmatrix} \begin{pmatrix} l_2 & l_2' & k \\ -m_2 & +m_2' & -q \end{pmatrix}

Where we have renamed -m_2' to m_2'.

Slater Type Orbitals (STO)

In this section we express the matrix elements in the STO basis. It turns out that all integrals that we need can be expressed in terms of the following simple integral (where n,\zeta \ge 0):

(11)\int_0^\infty r^n e^{-\zeta r} \d r =\int_0^\infty \left(x\over\zeta\right)^n e^{-x} {\d x\over\zeta} ={1\over\zeta^{n+1}} \int_0^\infty x^n e^{-x} \d x ={\Gamma(n+1)\over\zeta^{n+1}} ={n!\over\zeta^{n+1}}

The STO basis function for the radial Schrödinger equation for P(r) is:

(12)P_{n\zeta}(r) = N_{n\zeta} r^n e^{-\zeta r}

Where the normalization constant N_{n\zeta} is such that the STO orbital is normalized as the radial wavefunction P(r):

1 = \int_0^\infty P_{n\zeta}^2(r) \d r = N_{n\zeta}^2 \int_0^\infty r^{2n} e^{-2\zeta r} \d r = N_{n\zeta}^2 {(2n)!\over (2\zeta)^{2n+1}}

from which we get:

N_{n\zeta} = \sqrt{(2\zeta)^{2n+1}\over (2n)!}

Note that for R(r)={P(r)\over r} we get the following STO basis function:

(13)R_{n\zeta}(r) = {P_{n\zeta}(r)\over r} = N_{n\zeta} r^{n-1} e^{-\zeta r}

One uses either (12) or (13) depending on whether one solves the radial Schrödinger equation for P or for R={P\over r}.

Overlap

S_{ij} = \int P_{n_i \zeta_i}(r) P_{n_j \zeta_j}(r) \d r = = \int N_{n_i\zeta_i} r^{n_i} e^{-\zeta_i r} N_{n_j\zeta_j} r^{n_j} e^{-\zeta_j r} \d r = = N_{n_i\zeta_i} N_{n_j\zeta_j} \int r^{n_i + n_j} e^{-(\zeta_i+\zeta_j) r} \d r = = N_{n_i\zeta_i} N_{n_j\zeta_j} {(n_i + n_j)! \over (\zeta_i+\zeta_j)^{n_i+n_j+1}}

Potential

V_{ij} = \int P_{n_i\zeta_i}(r) \left(-{Z\over r}\right) P_{n_j\zeta_j}(r) \d r = = \int N_{n_i\zeta_i} r^{n_i} e^{-\zeta_i r} \left(-{Z\over r}\right) N_{n_j\zeta_j} r^{n_j} e^{-\zeta_j r} \d r = = -Z N_{n_i\zeta_i} N_{n_j\zeta_j} \int r^{n_i+n_j-1} e^{-(\zeta_i+\zeta_j) r} \d r = = -Z N_{n_i\zeta_i} N_{n_j\zeta_j} {(n_i + n_j - 1)! \over (\zeta_i+\zeta_j)^{n_i+n_j}}

Kinetic

T_{ij} = \int \left( \half P_{n_i\zeta_i}'(r) P_{n_j\zeta_j}'(r) +P_{n_i\zeta_i}(r) {l (l+1)\over 2 r^2} P_{n_j\zeta_j}(r) \right) \d r = = \half N_{n_i\zeta_i}N_{n_j\zeta_j}\int \left( {\d\over\d r}(r^{n_i} e^{-\zeta_i r}) {\d\over\d r}(r^{n_j} e^{-\zeta_j r}) +r^{n_i} e^{-\zeta_i r} {l (l+1)\over r^2} r^{n_j} e^{-\zeta_j r} \right) \d r = = \half N_{n_i\zeta_i}N_{n_j\zeta_j}\int \left( (n_i r^{n_i-1} e^{-\zeta_i r} -\zeta_i r^{n_i} e^{-\zeta_i r}) (n_j r^{n_j-1} e^{-\zeta_j r} -\zeta_j r^{n_j} e^{-\zeta_j r}) +l (l+1) r^{n_i+n_j-2} e^{-(\zeta_i+\zeta_j) r} \right) \d r = = \half N_{n_i\zeta_i}N_{n_j\zeta_j}\int \left( \left( {n_i n_j\over r^2} -{n_i \zeta_j+n_j\zeta_i \over r} +\zeta_i \zeta_j \right) r^{n_i+n_j} e^{-(\zeta_i+\zeta_j) r} +l (l+1) r^{n_i+n_j-2} e^{-(\zeta_i+\zeta_j) r} \right) \d r = = \half N_{n_i\zeta_i}N_{n_j\zeta_j}\int \left( (n_i n_j + l(l+1)) r^{n_i+n_j-2} e^{-(\zeta_i+\zeta_j) r} -(n_i \zeta_j+n_j\zeta_i) r^{n_i+n_j-1} e^{-(\zeta_i+\zeta_j) r} +\zeta_i \zeta_j r^{n_i+n_j } e^{-(\zeta_i+\zeta_j) r} \right) \d r = = \half N_{n_i\zeta_i}N_{n_j\zeta_j}\left( (n_i n_j + l(l+1)) {(n_i + n_j - 2)! \over (\zeta_i+\zeta_j)^{n_i +n_j - 1}} -(n_i \zeta_j+n_j\zeta_i) {(n_i + n_j - 1)! \over (\zeta_i+\zeta_j)^{n_i +n_j }} +\zeta_i \zeta_j {(n_i + n_j )! \over (\zeta_i+\zeta_j)^{n_i +n_j + 1}} \right)

Two particle

In this section we also need the following integral:

\int_u^\infty r^n e^{-\zeta r} \d r = {1\over\zeta^{n+1}} \int_{\zeta u}^\infty x^n e^{-x} \d x = {\Gamma(n+1, \zeta u)\over\zeta^{n+1}} = {n!\over \zeta^{n+1}} e^{-\zeta u} \sum_{\nu=0}^n {u^\nu \zeta^\nu \over \nu!}

where

\Gamma(n, x) = \int_x^\infty t^{n-1} e^{-t}\d t = (n-1)! e^{-x} \sum_{\nu=0}^{n-1} {x^\nu\over \nu!}

is the incomplete gamma function.

The Slater integral is

R^k(i, j, k, l) = \int_0^\infty \int_0^\infty {r_<^k\over r_>^{k+1}} P_i(r) P_k(r) P_j(r') P_l(r') \d r \d r' = = \int_0^\infty \d r \int_0^r \d r' {r'^k\over r^{k+1}} P_i(r) P_k(r) P_j(r') P_l(r') + + \int_0^\infty \d r \int_r^\infty \d r' {r^k\over r'^{k+1}} P_i(r) P_k(r) P_j(r') P_l(r') = = \int_0^\infty \d r \int_0^r \d r' {r'^k\over r^{k+1}} P_i(r) P_k(r) P_j(r') P_l(r') + + \int_0^\infty \d r' \int_0^{r'} \d r {r^k\over r'^{k+1}} P_i(r) P_k(r) P_j(r') P_l(r') = = \int_0^\infty \d r \int_0^r \d r' {r'^k\over r^{k+1}} P_i(r) P_k(r) P_j(r') P_l(r') + + \int_0^\infty \d r \int_0^r \d r' {r'^k\over r^{k+1}} P_i(r') P_k(r') P_j(r) P_l(r) = = R^k_\triangle(i, j, k, l) + R^k_\triangle(j, i, l, k)

where R^k_\triangle(i, j, k, l) is the integral over the lower triangle (assuming r is the x-axis and r' is the y-axis), that is r>r':

R^k_\triangle(i, j, k, l) = \int_0^\infty \d r \int_0^{r} \d r' {r'^k\over r^{k+1}} P_i(r) P_k(r) P_j(r') P_l(r') = = \int_0^\infty \d r' \int_{r'}^\infty \d r {r'^k\over r^{k+1}} P_i(r) P_k(r) P_j(r') P_l(r') = = N_{n_i\zeta_i}N_{n_j\zeta_j}N_{n_k\zeta_k}N_{n_l\zeta_l} \int_0^\infty \d r' \int_{r'}^\infty \d r {r'^k\over r^{k+1}} r^{n_i + n_k} r'^{n_j+n_l} e^{-(\zeta_i+\zeta_k)r} e^{-(\zeta_j+\zeta_l)r'} = = N_{n_i\zeta_i}N_{n_j\zeta_j}N_{n_k\zeta_k}N_{n_l\zeta_l} \int_0^\infty \d r' r'^k r'^{n_j+n_l} e^{-(\zeta_j+\zeta_l)r'} {(n_i+n_k-k-1)!\over (\zeta_i+\zeta_k)^{n_i+n_k-k}} e^{-(\zeta_i+\zeta_k) r'} \sum_{\nu=0}^{n_i+n_k-k-1} {r'^\nu (\zeta_i+\zeta_k)^\nu \over \nu!} = = N_{n_i\zeta_i}N_{n_j\zeta_j}N_{n_k\zeta_k}N_{n_l\zeta_l} {(n_i+n_k-k-1)!\over (\zeta_i+\zeta_k)^{n_i+n_k-k}} \sum_{\nu=0}^{n_i+n_k-k-1} {(n_j+n_l+k+\nu)! (\zeta_i+\zeta_k)^\nu \over \nu! (\zeta_i+\zeta_j+\zeta_k+\zeta_l)^{n_j+n_l+k+\nu+1} } = = N_{n_i\zeta_i}N_{n_j\zeta_j}N_{n_k\zeta_k}N_{n_l\zeta_l} {(n_i+n_k-k-1)!\over (\zeta_i+\zeta_k)^{n_i+n_k-k}} H^k_{ijkl}

where:

H^k_{ijkl} = \sum_{\nu=0}^{n_i+n_k-k-1}{ (n_j+n_l+k+\nu)! (\zeta_i+\zeta_k)^\nu\over \nu!(\zeta_i+\zeta_j+\zeta_k+\zeta_l)^{n_j+n_l+k+\nu+1}}

Much more tedious method is the following:

\int_0^u r^n e^{-\zeta r} \d r = \int_0^\infty r^n e^{-\zeta r} \d r -\int_u^\infty r^n e^{-\zeta r} \d r = {n!\over \zeta^{n+1}}\left(1- e^{-\zeta u} \sum_{\nu=0}^n {u^\nu \zeta^\nu \over \nu!}\right)

The Slater integral is given by:

R^k(i, j, k, l) = \int_0^\infty {r_<^k\over r_>^{k+1}} P_i(r) P_k(r) P_j(r') P_l(r') \d r \d r' = = \int_0^\infty P_i(r) P_k(r) {Y^k(P_j P_l, r)\over r} \d r

where

Y^k(f(r), r) = r \int_0^\infty {r_<^k\over r_>^{k+1}} f(r') \d r' = {1\over r^k} \int_0^r r'^k f(r') \d r' + r^{k+1} \int_r^\infty {1\over r'^{k+1}} f(r') \d r'

and we get:

Y^k(P_j(r) P_l(r), r) = {1\over r^k} \int_0^r r'^k P_j(r') P_l(r') \d r' + r^{k+1} \int_r^\infty {1\over r'^{k+1}} P_j(r') P_l(r') \d r' = = {N_{n_j \zeta_j}N_{n_l \zeta_l}\over r^k} \int_0^r r'^k r'^{n_j}e^{-\zeta_j r'} r'^{n_l}e^{-\zeta_l r'} \d r' + N_{n_j \zeta_j}N_{n_l \zeta_l} r^{k+1} \int_r^\infty {1\over r'^{k+1}} r'^{n_j}e^{-\zeta_j r'} r'^{n_l}e^{-\zeta_l r'} \d r' = = {N_{n_j \zeta_j}N_{n_l \zeta_l}\over r^k} \int_0^r r'^{n_j+n_l+k}e^{-(\zeta_j+\zeta_l) r'} \d r' + N_{n_j \zeta_j}N_{n_l \zeta_l} r^{k+1} \int_r^\infty r'^{n_j+n_l-k-1}e^{-(\zeta_j+\zeta_l) r'} \d r' = = {N_{n_j \zeta_j}N_{n_l \zeta_l}\over r^k} {(n_j+n_l+k)!\over(\zeta_j+\zeta_l)^{n_j+n_l+k+1}} \left( 1 - e^{-(\zeta_j+\zeta_l)r} \sum_{\nu=0}^{n_j+n_l+k} { r^\nu(\zeta_j+\zeta_l)^\nu \over \nu! } \right) + + N_{n_j \zeta_j}N_{n_l \zeta_l} r^{k+1} {(n_j+n_l-k-1)!\over (\zeta_j+\zeta_l)^{n_j+n_l-k}} e^{-(\zeta_j+\zeta_l)r} \sum_{\nu=0}^{n_j+n_l-k-1}{r^\nu(\zeta_j+\zeta_l)^\nu\over \nu!}

Putting everything together we get:

R^k(i, j, k, l) = N_{n_i \zeta_i}N_{n_j \zeta_j}N_{n_k \zeta_k}N_{n_l \zeta_l} \int_0^\infty r^{n_i}e^{-\zeta_i r} r^{n_k}e^{-\zeta_k r} {1\over r} \left( {1\over r^k} {(n_j+n_l+k)!\over(\zeta_j+\zeta_l)^{n_j+n_l+k+1}} \left( 1 - e^{-(\zeta_j+\zeta_l)r} \sum_{\nu=0}^{n_j+n_l+k} { r^\nu(\zeta_j+\zeta_l)^\nu \over \nu! } \right) + \right. \left. + r^{k+1} {(n_j+n_l-k-1)!\over (\zeta_j+\zeta_l)^{n_j+n_l-k}} e^{-(\zeta_j+\zeta_l)r} \sum_{\nu=0}^{n_j+n_l-k-1}{r^\nu(\zeta_j+\zeta_l)^\nu\over \nu!} \right) \d r = = N_{n_i \zeta_i}N_{n_j \zeta_j}N_{n_k \zeta_k}N_{n_l \zeta_l} \left( {(n_j+n_l+k)!\over(\zeta_j+\zeta_l)^{n_j+n_l+k+1}} \int_0^\infty r^{n_i+n_k-k-1}e^{-(\zeta_i+\zeta_k) r} \left( 1 - e^{-(\zeta_j+\zeta_l)r} \sum_{\nu=0}^{n_j+n_l+k} { r^\nu(\zeta_j+\zeta_l)^\nu \over \nu! } \right) \d r + \right. \left. + \int_0^\infty e^{-(\zeta_i+\zeta_j+\zeta_k+\zeta_l) r} {(n_j+n_l-k-1)!\over (\zeta_j+\zeta_l)^{n_j+n_l-k}} \sum_{\nu=0}^{n_j+n_l-k-1}{ r^{n_i+n_k+k+\nu}(\zeta_j+\zeta_l)^\nu\over \nu!} \d r \right) = = N_{n_i \zeta_i}N_{n_j \zeta_j}N_{n_k \zeta_k}N_{n_l \zeta_l} \left( {(n_j+n_l+k)!\over(\zeta_j+\zeta_l)^{n_j+n_l+k+1}} \left( {(n_i+n_k-k-1)!\over (\zeta_i+\zeta_k)^{n_i+n_k-k}} - \sum_{\nu=0}^{n_j+n_l+k} {(\zeta_j+\zeta_l)^\nu (n_i+n_k-k+\nu-1)! \over \nu! (\zeta_i+\zeta_j+\zeta_k+\zeta_l)^{n_i+n_k-k+\nu}} \right) + \right. \left. + {(n_j+n_l-k-1)!\over (\zeta_j+\zeta_l)^{n_j+n_l-k}} \sum_{\nu=0}^{n_j+n_l-k-1}{ (n_i+n_k+k+\nu)! (\zeta_j+\zeta_l)^\nu\over \nu!(\zeta_i+\zeta_j+\zeta_k+\zeta_l)^{n_i+n_k+k+\nu+1}} \right) = = N_{n_i \zeta_i}N_{n_j \zeta_j}N_{n_k \zeta_k}N_{n_l \zeta_l} \left( {(n_j+n_l+k)!\over(\zeta_j+\zeta_l)^{n_j+n_l+k+1}} \left( {(n_i+n_k-k-1)!\over (\zeta_i+\zeta_k)^{n_i+n_k-k}} - H^{-k-1}_{jilk} \right) + \right. \left. + {(n_j+n_l-k-1)!\over (\zeta_j+\zeta_l)^{n_j+n_l-k}} H^k_{jilk} \right) = = N_{n_i \zeta_i}N_{n_j \zeta_j}N_{n_k \zeta_k}N_{n_l \zeta_l} \left( {(n_i+n_k-k-1)!\over (\zeta_i+\zeta_k)^{n_i+n_k-k}} H^k_{ijkl} + {(n_j+n_l-k-1)!\over (\zeta_j+\zeta_l)^{n_j+n_l-k}} H^k_{jilk} \right)

where:

H^k_{ijkl} = \sum_{\nu=0}^{n_i+n_k-k-1}{ (n_j+n_l+k+\nu)! (\zeta_i+\zeta_k)^\nu\over \nu!(\zeta_i+\zeta_j+\zeta_k+\zeta_l)^{n_j+n_l+k+\nu+1}}

Gaussian Type Orbitals (GTO)

In this section we express the matrix elements in the GTO basis. It turns out that all integrals that we need can be expressed in terms of the following simple integral (where n,\zeta \ge 0):

(14)\int_0^\infty r^n e^{-\zeta r^2} \d r = {1\over\sqrt{\zeta^{n+1}}} \int_0^\infty x^n e^{-x^2} \d x = {1\over2\sqrt{\zeta^{n+1}}} \int_0^\infty t^{n-1\over 2} e^{-t} \d t = {\Gamma({n+1\over2})\over2 \zeta^{n+1\over2}} = = \begin{cases} (n-1)!! \sqrt{\pi\over 2 (2\zeta)^{n+1}} & \text{for even $n$} \\ {\left(n-1\over2\right)!\over 2 \sqrt{\zeta^{n+1}}} & \text{for odd $n$} \\ \end{cases}

The GTO basis function for the radial Schrödinger equation for P(r) is:

(15)P_{n\zeta}(r) = N_{n\zeta} r^n e^{-\zeta r^2}

However, unlike STO, the GTO functions must satisfy the condition n = l + 2i + 1 where i=0, 1, 2, \dots (this condition is later used in (14) to determine whether n is even or odd). The normalization constant N_{n\zeta} is such that the GTO orbital is normalized as the radial wavefunction P(r):

1 = \int_0^\infty P_{n\zeta}^2(r) \d r = N_{n\zeta}^2 \int_0^\infty r^{2n} e^{-2\zeta r^2} \d r = N_{n\zeta}^2 (2n-1)!! \sqrt{\pi\over 2 (4\zeta)^{2n+1}}

from which we get:

N_{n\zeta} = \sqrt{{1\over (2n-1)!!}\sqrt{2 (4\zeta)^{2n+1}\over\pi}}

Note that for R(r)={P(r)\over r} we get the following GTO basis function:

(16)R_{n\zeta}(r) = {P_{n\zeta}(r)\over r} = N_{n\zeta} r^{n-1} e^{-\zeta r^2}

One uses either (12) or (13) depending on whether one solves the radial Schrödinger equation for P or for R={P\over r}.

Overlap

S_{ij} = \int P_{n_i \zeta_i}(r) P_{n_j \zeta_j}(r) \d r = = \int N_{n_i\zeta_i} r^{n_i} e^{-\zeta_i r^2} N_{n_j\zeta_j} r^{n_j} e^{-\zeta_j r^2} \d r = = N_{n_i\zeta_i} N_{n_j\zeta_j} \int r^{n_i + n_j} e^{-(\zeta_i+\zeta_j) r^2} \d r = = N_{n_i\zeta_i} N_{n_j\zeta_j} {\Gamma({n_i+n_j+1\over 2}) \over 2(\zeta_i+\zeta_j)^{n_i+n_j+1\over 2}} = = N_{n_i\zeta_i} N_{n_j\zeta_j} (n_i + n_j -1)!! \sqrt{\pi\over 2(2\zeta_i+2\zeta_j)^{n_i+n_j+1}}

Potential

V_{ij} = \int P_{n_i\zeta_i}(r) \left(-{Z\over r}\right) P_{n_j\zeta_j}(r) \d r = = \int N_{n_i\zeta_i} r^{n_i} e^{-\zeta_i r^2} \left(-{Z\over r}\right) N_{n_j\zeta_j} r^{n_j} e^{-\zeta_j r^2} \d r = = -Z N_{n_i\zeta_i} N_{n_j\zeta_j} \int r^{n_i+n_j-1} e^{-(\zeta_i+\zeta_j) r^2} \d r = = -Z N_{n_i\zeta_i} N_{n_j\zeta_j} {\Gamma({n_i+n_j\over 2}) \over 2(\zeta_i+\zeta_j)^{n_i+n_j\over 2}} = = -Z N_{n_i\zeta_i} N_{n_j\zeta_j} {\left(n_i + n_j - 2\over 2\right)! \over 2\sqrt{(\zeta_i+\zeta_j)^{n_i+n_j}}}

Kinetic

T_{ij} = \int \left( \half P_{n_i\zeta_i}'(r) P_{n_j\zeta_j}'(r) +P_{n_i\zeta_i}(r) {l (l+1)\over 2 r^2} P_{n_j\zeta_j}(r) \right) \d r = = \half N_{n_i\zeta_i}N_{n_j\zeta_j}\int \left( {\d\over\d r}(r^{n_i} e^{-\zeta_i r^2}) {\d\over\d r}(r^{n_j} e^{-\zeta_j r^2}) +r^{n_i} e^{-\zeta_i r^2} {l (l+1)\over r^2} r^{n_j} e^{-\zeta_j r^2} \right) \d r = = \half N_{n_i\zeta_i}N_{n_j\zeta_j}\int \left( (n_i r^{n_i-1} e^{-\zeta_i r^2} -2\zeta_i r^{n_i+1} e^{-\zeta_i r^2}) (n_j r^{n_j-1} e^{-\zeta_j r^2} -2\zeta_j r^{n_j+1} e^{-\zeta_j r^2}) +l (l+1) r^{n_i+n_j-2} e^{-(\zeta_i+\zeta_j) r^2} \right) \d r = = \half N_{n_i\zeta_i}N_{n_j\zeta_j}\int \left( \left( {n_i n_j\over r^2} -2(n_i \zeta_j+n_j\zeta_i) +4\zeta_i \zeta_j r^2 \right) r^{n_i+n_j} e^{-(\zeta_i+\zeta_j) r^2} +l (l+1) r^{n_i+n_j-2} e^{-(\zeta_i+\zeta_j) r^2} \right) \d r = = \half N_{n_i\zeta_i}N_{n_j\zeta_j}\int \left( (n_i n_j + l(l+1)) r^{n_i+n_j-2} e^{-(\zeta_i+\zeta_j) r^2} -2(n_i \zeta_j+n_j\zeta_i) r^{n_i+n_j} e^{-(\zeta_i+\zeta_j) r^2} +4\zeta_i \zeta_j r^{n_i+n_j+2} e^{-(\zeta_i+\zeta_j) r^2} \right) \d r = = \half N_{n_i\zeta_i}N_{n_j\zeta_j} \left( (n_i n_j + l(l+1)) {\Gamma({n_i+n_j-1\over 2}) \over 2(\zeta_i+\zeta_j)^{n_i+n_j-1\over 2}} -2(n_i \zeta_j+n_j\zeta_i) {\Gamma({n_i+n_j+1\over 2}) \over 2(\zeta_i+\zeta_j)^{n_i+n_j+1\over 2}} +4\zeta_i \zeta_j {\Gamma({n_i+n_j+3\over 2}) \over 2(\zeta_i+\zeta_j)^{n_i+n_j+3\over 2}} \right) \d r = = \half N_{n_i\zeta_i}N_{n_j\zeta_j}\Bigg( (n_i n_j + l(l+1)) {(n_i + n_j - 3)!! \sqrt{\pi\over 2(2\zeta_i+2\zeta_j)^{n_i + n_j - 1}}} + -2(n_i \zeta_j+n_j\zeta_i) {(n_i + n_j - 1)!! \sqrt{\pi\over 2(2\zeta_i+2\zeta_j)^{n_i + n_j + 1}}} + +4\zeta_i \zeta_j {(n_i + n_j + 1)!! \sqrt{\pi\over 2(2\zeta_i+2\zeta_j)^{n_i + n_j + 3}}} \Bigg)

Two particle

We will need the following integral:

\int_u^\infty r^n e^{-\zeta r^2} \d r = {1\over2\zeta^{n+1\over2}} \int_{\zeta u^2}^\infty x^{n-1\over2} e^{-x} \d x = {\Gamma({n+1\over2}, \zeta u^2)\over2\zeta^{n+1\over2}}

Just like for STO, we get:

R^k(i, j, k, l) = R^k_\triangle(i, j, k, l) + R^k_\triangle(j, i, l, k)

where R^k_\triangle(i, j, k, l) is the integral over the lower triangle (assuming r is the x-axis and r' is the y-axis), that is r>r':

R^k_\triangle(i, j, k, l) = \int_0^\infty \d r \int_0^{r} \d r' {r'^k\over r^{k+1}} P_i(r) P_k(r) P_j(r') P_l(r') = = \int_0^\infty \d r' \int_{r'}^\infty \d r {r'^k\over r^{k+1}} P_i(r) P_k(r) P_j(r') P_l(r') = = N_{n_i\zeta_i}N_{n_j\zeta_j}N_{n_k\zeta_k}N_{n_l\zeta_l} \int_0^\infty \d r' \int_{r'}^\infty \d r {r'^k\over r^{k+1}} r^{n_i + n_k} r'^{n_j+n_l} e^{-(\zeta_i+\zeta_k)r^2} e^{-(\zeta_j+\zeta_l)r'^2} = = N_{n_i\zeta_i}N_{n_j\zeta_j}N_{n_k\zeta_k}N_{n_l\zeta_l} \int_0^\infty \d r' r'^k r'^{n_j+n_l} e^{-(\zeta_j+\zeta_l)r'^2} {\Gamma({n_i+n_k-k\over2}, (\zeta_i+\zeta_k)r'^2)\over 2(\zeta_i+\zeta_k)^{n_i+n_k-k\over2}} = = N_{n_i\zeta_i}N_{n_j\zeta_j}N_{n_k\zeta_k}N_{n_l\zeta_l} \int_0^\infty \d r' r'^k r'^{n_j+n_l} e^{-(\zeta_i+\zeta_j+\zeta_k+\zeta_l)r'^2} {\left({n_i+n_k-k\over2}-1\right)! \over 2(\zeta_i+\zeta_k)^{n_i+n_k-k\over2}} \sum_{\nu=0}^{n_i+n_k-k-2\over2}{(\zeta_i+\zeta_k)^\nu r'^{2\nu} \over \nu!} = = N_{n_i\zeta_i}N_{n_j\zeta_j}N_{n_k\zeta_k}N_{n_l\zeta_l} {\left({n_i+n_k-k\over2}-1\right)! \over 2(\zeta_i+\zeta_k)^{n_i+n_k-k\over2}} \sum_{\nu=0}^{n_i+n_k-k-2\over2}{(\zeta_i+\zeta_k)^\nu \over \nu!} {\Gamma({n_j+n_l+k+2\nu+1\over2})\over 2 (\zeta_i+\zeta_j+\zeta_k+\zeta_l)^{n_j+n_l+k+2\nu+1\over2}} = = N_{n_i\zeta_i}N_{n_j\zeta_j}N_{n_k\zeta_k}N_{n_l\zeta_l} {\left({n_i+n_k-k\over2}-1\right)! \over 2(\zeta_i+\zeta_k)^{n_i+n_k-k\over2}} \sum_{\nu=0}^{n_i+n_k-k-2\over2}{(\zeta_i+\zeta_k)^\nu \over \nu!} {(n_j+n_l+k+2\nu-1)!!\sqrt{\pi}\over 2^{n_j+n_l+k+2\nu+2\over 2} (\zeta_i+\zeta_j+\zeta_k+\zeta_l)^{n_j+n_l+k+2\nu+1\over2}} = = N_{n_i\zeta_i}N_{n_j\zeta_j}N_{n_k\zeta_k}N_{n_l\zeta_l} {\sqrt{\pi}\over 2^{p+2}\sqrt{(\zeta_i+\zeta_j+\zeta_k +\zeta_l)^{2p+1}}} H^k_{ijkl}

where p={n_i+n_j+n_k+n_l-2\over2} and:

H^k_{ijkl} = \sum_{\nu=0}^{n_i+n_k-k\over2}{ (2p-2\nu-1)!! \left(n_i+n_k-k\over 2\right)! 2^\nu (\zeta_i+\zeta_j+\zeta_k+\zeta_l)^\nu \over (n_i+n_k-\nu)! (\lambda_i + \lambda_k)^{1+\nu}}

Exchange Integral in Spherical Symmetry

Let’s calculate the exchange integral

\int {\psi_i^*({\bf x})\psi_j({\bf x})\psi_j^*({\bf x}')\psi_i({\bf x}') \over |{\bf x} - {\bf x}'|} \d^3 x \d^3 x'

for the particular choice of the functions \psi:

\psi_i({\bf x}) = {P_{nl}(r)\over r} Y_{lm}(\Omega) \psi_j({\bf x}) = {P_{n'l'}(r)\over r} Y_{l'm'}(\Omega)

We use multipole expansion:

{1\over |{\bf x}-{\bf x}'|} = \sum_{k,q}{r_{<}^k\over r_{>}^{k+1}} {4\pi\over 2k+1}Y_{kq}(\Omega)Y_{kq}^*(\Omega')

And we get:

\int {\psi_i^*({\bf x})\psi_j({\bf x})\psi_j^*({\bf x}')\psi_i({\bf x}') \over |{\bf x} - {\bf x}'|} \d^3 x \d^3 x' = = \int {P_{nl}(r)\over r} Y_{lm}^*(\Omega) {P_{n'l'}(r)\over r} Y_{l'm'}(\Omega) {P_{n'l'}(r')\over r'} Y_{l'm'}^*(\Omega') {P_{nl}(r')\over r'} Y_{lm}(\Omega') \sum_{k,q}{r_{<}^k\over r_{>}^{k+1}} {4\pi\over 2k+1}Y_{kq}(\Omega)Y_{kq}^*(\Omega') r^2 r'^2 \d r \d r' \d \Omega \d \Omega' = = \sum_{k,q} \int Y_{lm}^*(\Omega) Y_{l'm'}(\Omega) Y_{kq}(\Omega) \d \Omega \int Y_{l'm'}^*(\Omega') Y_{lm}(\Omega') Y_{kq}^*(\Omega') \d \Omega' \int {r_{<}^k\over r_{>}^{k+1}} {4\pi\over 2k+1} P_{nl}(r) P_{n'l'}(r) P_{n'l'}(r') P_{nl}(r') \d r \d r' = = \sum_{k,q} \int Y_{lm}^*(\Omega) Y_{l'm'}(\Omega) Y_{kq}(\Omega) \d \Omega (-1)^{m+m'+q} \int Y_{l',-m'}(\Omega') Y_{l,-m}^*(\Omega') Y_{k,-q}(\Omega') \d \Omega' \int {r_{<}^k\over r_{>}^{k+1}} {4\pi\over 2k+1} P_{nl}(r) P_{n'l'}(r) P_{n'l'}(r') P_{nl}(r') \d r \d r' = = \sum_{k} c^k(l, m, l', m') \sqrt{2k+1\over 4\pi} (-1)^{m+m'+m-m'} c^k(l, -m, l', -m') \sqrt{2k+1\over 4\pi} \int {r_{<}^k\over r_{>}^{k+1}} {4\pi\over 2k+1} P_{nl}(r) P_{n'l'}(r) P_{n'l'}(r') P_{nl}(r') \d r \d r' = = \sum_{k} c^k(l, m, l', m') c^k(l, -m, l', -m') \int {r_{<}^k\over r_{>}^{k+1}} P_{nl}(r) P_{n'l'}(r) P_{n'l'}(r') P_{nl}(r') \d r \d r' = = \sum_{k} c^k(l, m, l', m') c^k(l, m, l', m') \int {r_{<}^k\over r_{>}^{k+1}} P_{nl}(r) P_{n'l'}(r) P_{n'l'}(r') P_{nl}(r') \d r \d r' = = \sum_{k=|l-l'|}^{l+l'} \left( c^k(l, m, l', m') \right)^2 \int {r_{<}^k\over r_{>}^{k+1}} P_{nl}(r) P_{n'l'}(r) P_{n'l'}(r') P_{nl}(r') \d r \d r'

Occupation Numbers

We have a sum over N electron states like this:

\sum_{i=1}^N A_i({\bf x}) = \sum_{nlms} A_{nlms}({\bf x})

where A_{nlms}({\bf x}) are some functions that depend on the state numbers (for example squares of the wavefunctions). Then there are two options — either there is a way to sum over the m and s degrees of freedom, so that the sum can be written exactly as:

\sum_{nlms} A_{nlms}({\bf x}) = \sum_{nlms} B_{nl}({\bf x})

where B_{nl} (that don’t depend on m and s) will in general be different to A_{nlms}, but the sum will be the same. Or we have to approximate the sum (for example by averaging over the angles, or in some other way) as:

\sum_{nlms} A_{nlms}({\bf x}) \to \sum_{nlms} B_{nl}({\bf x})

In either case, the occupation numbers f_{nl} are simply the number of times the functions B_{nl}({\bf x}) appear in the sum for the given n and l:

\sum_{nlms} B_{nl}({\bf x}) = \sum_{nl} f_{nl} B_{nl}({\bf x})

So for closed shells atoms, it is always:

f_{nl} = 2(2l+1)

because there are two spins, and 2l+1 possibilities for m, for open shell atoms, f_{nl} is anything between 0 and 2l+1.

Example I

As an example, let’s say that after some calculation for closed shell systems we get exactly:

\sum_{nlms} A_{nlms}({\bf x}) = \sum_{nl} 2(2l+1) B_{nl}({\bf x})

Then because there are exactly 2(2l+1) states in the nl shell, we write the above as:

\sum_{nlms} A_{nlms}({\bf x}) = \sum_{nl} 2(2l+1) B_{nl}({\bf x}) = \sum_{nl} f_{nl} B_{nl}({\bf x})

Then we do similar calculation for the open shell system, and we have to use some approximations to get the following formula, where the B_{nl}({\bf x}) happen to be exactly the same as for the closed shell system:

\sum_{nlms} A_{nlms}({\bf x}) \to \sum_{nlm} 2 B_{nl}({\bf x})

Then we denote by f_{nl} the number of electrons in the shell nl (at least one of them will be open, for which nl we have f_{nl} < 2(2l+1)), and we can write the above as:

\sum_{nlms} A_{nlms}({\bf x}) \to \sum_{nlm} 2 B_{nl}({\bf x}) = \sum_{nl} f_{nl} B_{nl}({\bf x})

Example II

The usual chemical occupation numbers for the Uranium atom are:

f_{1l} & = 2 (2l+1) \\ f_{2l} & = 2 (2l+1) \\ f_{3l} & = 2 (2l+1) \\ f_{4l} & = 2 (2l+1) \\ f_{5l} & = 2 (2l+1)\quad\quad\mbox{for $l\le2$} \\ f_{53} & = 3 \\ f_{60} & = 2 \\ f_{61} & = 6 \\ f_{62} & = 1 \\ f_{70} & = 2 \\

So the n=5, l=3 and n=6, l=2 shells are open, all others are closed. By summing all these f_{nl}, we get 92 states as expected:

\sum_{nl} f_{nl} = 2 + (2+6) + (2+6+10) + (2+6+10+14) + (2+6+10) + + 3 + 2 + 6 + 1 + 2 = 92

Code:

def f_nl(n, l):
    if n < 5 or (n == 5 and l <= 2):
        return 2*(2*l+1)
    else:
        d = {
            (5, 3): 3,
            (6, 0): 2,
            (6, 1): 6,
            (6, 2): 1,
            (7, 0): 2,
            }
        if (n, l) in d:
            return d[n, l]
        else:
            return 0

print "Sum f_nl =", sum([f_nl(n, l) for n in range(8) for l in range(n)])

prints:

Sum f_nl = 92

Hartree Screening Functions

Hartree screening function Y^k(f, r) is defined as:

Y^k(f, r) = r \int_0^\infty {r_{<}^k\over r_{>}^{k+1}} f(r') \d r'

and it occurs in many formulas in the Hartree Fock theory, so this section shows how to calculate it. It depends on k and a function f(r).

We first do the integral:

Y^k(f, r) = r \int_0^\infty {r_{<}^k\over r_{>}^{k+1}} f(r') \d r' = r \int_0^r {r'^k\over r^{k+1}} f(r') \d r' + r \int_r^\infty {r^k\over r'^{k+1}} f(r') \d r' = = {1\over r^k} \int_0^r {x^k} f(x) \d x + r^{k+1} \int_r^\infty {1\over x^{k+1}} f(x) \d x =Z^k(r) + r^{k+1} \int_r^\infty {1\over x^{k+1}} f(x) \d x

where:

Z^k(r) = {1\over r^k} \int_0^r {x^k} f(x) \d x {\d Z^k(r) \over \d r}= -{k\over r} Z^k(r) + f(r) Z^k(0) = 0

Now we differentiate Y^k(r):

{\d Y^k(r) \over \d r} = {\d Z^k(r) \over \d r} + {k+1\over r} r^{k+1} \int_r^\infty {1\over x^{k+1}} f(x) \d x -f(r) = = -{k\over r} Z^k(r) + f(r) + {k+1\over r} r^{k+1} \int_r^\infty {1\over x^{k+1}} f(x) \d x -f(r) = = -{k\over r} Z^k(r) + {k+1\over r} r^{k+1} \int_r^\infty {1\over x^{k+1}} f(x) \d x = = -{k\over r} Z^k(r) + {k+1\over r} (Y^k(r) - Z^k(r)) = = -{2k+1\over r} Z^k(r) + {k+1\over r} Y^k(r)

Also Y^k(\infty) = Z^k(\infty), so we get the following set of first order differential equations with boundary conditions:

\left({\d\over\d r} - {k+1\over r}\right) Y^k(r) = -{2k+1\over r} Z^k(r) \left({\d\over\d r} + {k\over r}\right) Z^k(r) = f(r) Y^k(\infty) = Z^k(\infty) Z^k(0) = 0

One way to calculate the Hartree screening function is to integrate the second equation from the left using the boundary condition Z^k(0) = 0 and then integrate the first equation from the right, using the boundary condition Y^k(\infty) = Z^k(\infty).

Another way is to obtain one second order equation. Expressing Z^k from the first equation:

Z^k(r) = -{r\over 2k+1}\left({\d\over\d r} - {k+1\over r}\right) Y^k(r) = =-{r\over 2k+1}{\d Y^k(r)\over \d r} + {k+1\over 2k+1} Y^k(r)

and substituting into the second equation we get:

-\left({\d\over\d r} + {k\over r}\right) \left({r\over 2k+1}{\d Y^k(r)\over \d r} + {k+1\over 2k+1} Y^k\right) = f(r) -{r\over 2k+1}\left({\d^2\over\d r^2} - {k(k+1)\over r^2}\right) Y^k(r) = f(r) \left(-{\d^2\over\d r^2} + {k(k+1)\over r^2}\right) Y^k(r) = {2k+1\over r} f(r)

With boundary condition on the left:

Z^k(0) = {k+1\over 2k+1} Y^k(0) = 0 Y^k(0) = 0

and on the right:

Z^k(r) =-{r\over 2k+1}{\d Y^k(r)\over \d r} + {k+1\over 2k+1} Y^k(r) = Y^k(r) -{r\over 2k+1}{\d Y^k(r)\over \d r} - {k\over 2k+1} Y^k(r) = 0 {\d Y^k(r)\over \d r} + {k\over r} Y^k(r) = 0

which for r\to\infty becomes:

\left.{\d Y^k(r)\over \d r}\right|_{r=\infty} = 0

but in practise, it’s better to use the former Newton (Robin) boundary condition. We have obtained one second order equation for Y^k(r)

\left(-{\d^2\over\d r^2} + {k(k+1)\over r^2}\right) Y^k(r) = {2k+1\over r} f(r)

with boundary conditions:

Y^k(0) = 0 {\d Y^k(r)\over \d r} + {k\over r} Y^k(r) = 0

The weak formulation is:

\int_0^{r_{max}} Y^k{}'(r) v'(r) + {k(k+1)\over r^2} Y^k(r) v(r) \d r -[Y^k{}'(r)v(r)]_0^{r_{max}} = \int_0^{r_{max}} {2k+1\over r}f(r)v(r) \d r

The boundary term can be simplified using the boundary conditions as:

-[Y^k{}'(r)v(r)]_0^{r_{max}} = -Y^k{}'(r_{max})v(r_{max}) + Y^k{}'(0) v(0) = -Y^k{}'(r_{max})v(r_{max}) = {k\over r_{max}} Y^k(r_{max})v(r_{max})

so we get

\int_0^{r_{max}} Y^k{}'(r) v'(r) + {k(k+1)\over r^2} Y^k(r) v(r) \d r + {k\over r_{max}} Y^k(r_{max})v(r_{max}) = \int_0^{r_{max}} {2k+1\over r}f(r)v(r) \d r

where the test functions v(r) have the constrain v(0)=0 on the left boundary and no constrain on the right.

Hartree Potential in Spherical Symmetry

For both open and closed shell atoms we get exactly:

V_H({\bf x}) = \int {n({\bf y})\over|{\bf x}-{\bf y}|} \d^3 y = \int {\sum_{j=1}^Z|\psi_j({\bf y})|^2\over|{\bf x}-{\bf y}|} \d^3 y = = 2\sum_{nlm}\int {|Y_{lm}(\Omega')|^2 P_{nl}^2(r')\over |{\bf x}-{\bf y}|} \d\Omega' \d r' = = 2\sum_{nlm}\sum_{l'm'}\int {r_<^{l'}\over r_>^{l'+1}} {4\pi\over 2l'+1} Y_{lm}^*(\Omega')Y_{lm}(\Omega') Y_{l'm'}^*(\Omega)Y_{l'm'}(\Omega') P_{nl}^2(r') \d\Omega' \d r' = = 2\sum_{nlm}\sum_{l'}\int {r_<^{l'}\over r_>^{l'+1}} {4\pi\over 2l'+1} Y_{l'0}^*(\Omega) \sqrt{2l'+1\over 4\pi} c^{l'}(l, m, l, m) P_{nl}^2(r') \d r' = = 2\sum_{nl}\sum_{l'=0}^{2l}\int {r_<^{l'}\over r_>^{l'+1}} \sqrt{4\pi\over 2l'+1} Y_{l'0}^*(\Omega) \sum_m c^{l'}(l, m, l, m) P_{nl}^2(r') \d r' = = 2\sum_{nl}\int {1\over r_>} \sqrt{4\pi} Y_{00}^*(\Omega) \sum_m c^0(l, m, l, m) P_{nl}^2(r') \d r' + + 2\sum_{nl}\sum_{l'=1}^{2l}\int {r_<^{l'}\over r_>^{l'+1}} \sqrt{4\pi\over 2l'+1} Y_{l'0}^*(\Omega) \sum_m c^{l'}(l, m, l, m) P_{nl}^2(r') \d r' = = \sum_{nl}2\sum_m c^0(l, m, l, m) \int {1\over r_>} P_{nl}^2(r') \d r' + + 2\sum_{nl}\sum_{l'=1}^{2l}\int {r_<^{l'}\over r_>^{l'+1}} \sqrt{4\pi\over 2l'+1} Y_{l'0}^*(\Omega) \sum_m c^{l'}(l, m, l, m) P_{nl}^2(r') \d r' = = \sum_{nl}2\sum_m \int {1\over r_>} P_{nl}^2(r') \d r' + + 2\sum_{nl}\sum_{l'=1}^{2l}\int {r_<^{l'}\over r_>^{l'+1}} \sqrt{4\pi\over 2l'+1} Y_{l'0}^*(\Omega) \sum_m c^{l'}(l, m, l, m) P_{nl}^2(r') \d r' = = \sum_{nl} f_{nl} \int {1\over r_>} P_{nl}^2(r') \d r' + + 2\sum_{nl}\sum_{l'=1}^{2l} \sqrt{4\pi\over 2l'+1} Y_{l'0}^*(\Omega)\sum_m c^{l'}(l, m, l, m) \int {r_<^{l'}\over r_>^{l'+1}} P_{nl}^2(r') \d r'

For closed shell atoms we use the fact, that

\sum_{m=-l}^l c^{l'}(l, m, l, m) = (2l+1) \delta_{l' 0}

and the second term disappears, and for open shell atoms we have to use the central field approximation: we average the integral for V_H over the angles:

V_H({\bf x}) \to V_H(r) = {1\over 4\pi} \int V_H({\bf x})\, \d \Omega

and using the fact, that

\int Y_{l'0}^*(\Omega)\, \d \Omega = \sqrt{4\pi} \delta_{l' 0}

the second term disappears as well. We got the same expression for both open shell (with central field approximation) and closed shell (no approximation) atoms. The radial charge density is:

n(r) = {1\over 4\pi} \sum_{nl} f_{nl} \left(P_{nl}(r)\over r\right)^2

So we got:

V_H(r) = \sum_{nl} f_{nl} \int {1\over r_>} P_{nl}^2(r') \d r' =\int {4\pi n(r') r'^2 \over r_>} \d r' = {Y^0(4\pi n(r) r^2, r) \over r}

The Hartree screening function Y^0(4\pi n(r) r^2, r) is given by the equation:

-{\d^2\over\d r^2} Y^0(r) = {1\over r} 4\pi n(r) r^2

So V_H(r) satisfies the radial Poisson equation:

(V_H(r) r)'' = -{1\over r} 4\pi n(r) r^2 V_H''(r) r + 2V_H'(r) = - 4\pi n(r) r V_H''(r) + {2\over r}V_H'(r) = -4\pi n(r)

Nonlocal Exchange Potential in Spherical Symmetry

Similarly, we calculate:

\sum_{j=1}^Z\int {\psi_i({\bf x'})\psi_j^*({\bf x'})\over|{\bf x}-{\bf x'}|} \d^3 x'\,\,\psi_j({\bf x}) = = \sum_{n'l'm'}\sum_{k,q}\int {P_{nl}(r')\over r'} Y_{lm}(\Omega') {P_{n'l'}(r')\over r'} Y_{l'm'}^*(\Omega') {P_{n'l'}(r)\over r} Y_{l'm'}(\Omega) {r_{<}^k\over r_{>}^{k+1}} {4\pi\over 2k+1}Y_{kq}(\Omega)Y_{kq}^*(\Omega') r'^2 \d r' \d \Omega' = = \sum_{n'l'm'}\sum_{k,q} {P_{n'l'}(r)\over r} {4\pi\over 2k+1} \int Y_{lm}(\Omega') Y_{l'm'}^*(\Omega') Y_{kq}^*(\Omega') Y_{l'm'}(\Omega) Y_{kq}(\Omega) \d \Omega' \int {r_{<}^k\over r_{>}^{k+1}} P_{nl}(r') P_{n'l'}(r') \d r' = = \sum_{n'l'}\sum_{k} {P_{n'l'}(r)\over r} {4\pi\over 2k+1} {2k+1\over 4\pi} \sqrt{2l'+1\over 2l+1} c^k(l', 0, l, 0) Y_{lm}(\Omega) \int {r_{<}^k\over r_{>}^{k+1}} P_{nl}(r') P_{n'l'}(r') \d r' = = {Y_{lm}(\Omega)\over r} \sum_{n'l'}\sum_{k=|l-l'|}^{k=l+l'} \sqrt{2l'+1\over 2l+1} c^k(l', 0, l, 0) \int {r_{<}^k\over r_{>}^{k+1}} P_{nl}(r') P_{n'l'}(r') \d r'\, P_{n'l'}(r) = = {Y_{lm}(\Omega)\over r} \sum_{n'l'}\sum_{k=|l-l'|}^{k=l+l'} (2l'+1) \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 \int {r_{<}^k\over r_{>}^{k+1}} P_{nl}(r') P_{n'l'}(r') \d r'\, P_{n'l'}(r) = = {Y_{lm}(\Omega)\over r} \sum_{n'l'} f_{n'l'} \sum_{k=|l-l'|}^{k=l+l'} \half \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 \int {r_{<}^k\over r_{>}^{k+1}} P_{nl}(r') P_{n'l'}(r') \d r'\, P_{n'l'}(r)

Functions with different spins don’t contribute to the sum, so there is no multiplication by 2. We assumed closed shells atoms (we summed over all m' in the above). We used the result of the integral in Example VI and also:

(17)\sqrt{2l'+1\over 2l+1} c^k(l', 0, l, 0) = \sqrt{2l'+1\over 2l+1} \sqrt{(2l'+1)(2l+1)} \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 = (2l'+1) \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2

Radial Hartree-Fock Equations

Using the above integrals, the HF equations become:

-\half P_{nl}''(r) + \left({l(l+1)\over 2r^2} -{Z\over r} + V_H(r)\right)P_{nl}(r) + -\sum_{n'l'} f_{n'l'} \sum_{k=|l-l'|}^{k=l+l'} \half \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 \int {r_{<}^k\over r_{>}^{k+1}} P_{nl}(r') P_{n'l'}(r') \d r'\, P_{n'l'}(r) = \epsilon_{nl} P_{nl}(r)

with:

V_H(r) = \sum_{nl} f_{nl} \int {1\over r_>} P_{nl}^2(r') \d r'

Using the Hartree screening functions, the HF equations are:

-\half P_{nl}''(r) + \left({l(l+1)\over 2r^2} -{Z\over r} + V_H(r)\right)P_{nl}(r) + -\sum_{n'l'} f_{n'l'} \sum_{k=|l-l'|}^{k=l+l'} \half \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 {Y^k(P_{nl}(r) P_{n'l'}(r), r) \over r} P_{n'l'}(r) = \epsilon_{nl} P_{nl}(r)

with:

V_H(r) = \sum_{nl} f_{nl} {Y^0(P_{nl}^2(r), r) \over r} = {Y^0(4\pi n(r) r^2, r) \over r}

Total Energy

The total energy is:

E = \sum_a 2(2l_a+1) \left( \epsilon_a -\sum_b (2l_b+1) \left( R_0(a, b, a, b)-\sum_l \half \begin{pmatrix} l_a & l & l_b \\ 0 & 0 & 0 \end{pmatrix}^2 R_l(a, b, b, a)\right) \right) = = \sum_{nl} f_{nl} \left( \epsilon_{nl} -\sum_{n'l'} \half f_{n'l'} \left( R_0(nl, n'l', nl, n'l')-\sum_k \half \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 R_l(nl, n'l', n'l', nl)\right) \right) = = \sum_{nl} f_{nl} \epsilon_{nl} -\sum_{nl}\sum_{n'l'} \half f_{nl} f_{n'l'} \left( \int_0^\infty P_{nl}^2(r) {Y^0(P_{n'l'}^2(r), r)\over r} \d r -\sum_k \half \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 \int_0^\infty P_{nl}(r) P_{n'l'}(r) {Y^l(P_{nl}(r) P_{n'l'}(r), r)\over r} \d r \right) = = \sum_{nl} f_{nl} \epsilon_{nl} -\half \int_0^\infty 4\pi n(r) r^2 {Y^0(4\pi n(r) r^2, r)\over r} \d r + +\sum_{nl}\sum_{n'l'} \half f_{nl} f_{n'l'} \sum_k \half \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 \int_0^\infty P_{nl}(r) P_{n'l'}(r) {Y^l(P_{nl}(r) P_{n'l'}(r), r)\over r} \d r = = \sum_{nl} f_{nl} \epsilon_{nl} -2\pi \int_0^\infty V_H(r) n(r) r^2 \d r +\sum_{nl}\sum_{n'l'} \half f_{nl} f_{n'l'} \sum_k \half \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 \int_0^\infty P_{nl}(r) P_{n'l'}(r) {Y^l(P_{nl}(r) P_{n'l'}(r), r)\over r} \d r

where:

R_l(a, b, c, d) = \int_0^\infty P_a(r) P_c(r) {Y^l(P_b(r) P_d(r), r)\over r} \d r

Example: Helium

For Helium atom, the only nonzero occupation numbers are:

f_{10} = 2

and the sum over n'l' simplifies to:

\sum_{n'l'} f_{n'l'} \sum_{k=|l-l'|}^{k=l+l'} \half \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 = f_{10} \half \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}^2 = f_{10} \half = 1

so we only need to solve for the 1s state and we get:

-\half P_{10}''(r) + \left(-{Z\over r} + V_H(r)\right)P_{10}(r) -{Y^0(P_{10}(r) P_{10}(r), r) \over r} P_{10}(r) = \epsilon_{10} P_{10}(r)

with:

V_H(r) = 2 {Y^0(P_{10}^2(r), r) \over r} = {Y^0(4\pi n(r) r^2, r) \over r}

We can combine the equations:

-\half P_{10}''(r) + \left(-{Z\over r} + 2 {Y^0(P_{10}^2(r), r) \over r}\right)P_{10}(r) -{Y^0(P_{10}^2(r), r) \over r} P_{10}(r) = \epsilon_{10} P_{10}(r)

and we obtain:

-\half P_{10}''(r) + \left(-{Z\over r} + {Y^0(P_{10}^2(r), r) \over r}\right)P_{10}(r) = \epsilon_{10} P_{10}(r)

FEM

The weak formulation is (u(r) = P_{nl}(r)):

\int_0^\infty \left( \half u'(r) v'(r) + \left({l(l+1)\over 2r^2} -{Z\over r} + V_H(r)\right)u(r)v(r) \right) \d r+ -\sum_{n'l'} f_{n'l'} \sum_{k=|l-l'|}^{k=l+l'} \half \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 \int_0^\infty v(r) P_{n'l'}(r) {Y^k(u(r) P_{n'l'}(r), r)\over r} \d r = \epsilon \int_0^\infty u(r)v(r)\d r

for closed shell atoms:

\int_0^\infty \left( \half u'(r) v'(r) + \left({l(l+1)\over 2r^2} -{Z\over r} + V_H(r)\right)u(r)v(r) \right) \d r+ -\sum_{n'l'} 2(2l'+1) \sum_{k=|l-l'|}^{k=l+l'} \half \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 \int_0^\infty v(r) P_{n'l'}(r) {Y^k(u(r) P_{n'l'}(r), r)\over r} \d r = \epsilon \int_0^\infty u(r)v(r)\d r

or (here we use the il index to label all functions u_i for the given l)

\int_0^\infty \left( \half u_{il}'(r) v'(r) + \left({l(l+1)\over 2r^2} -{Z\over r} + V_H(r)\right)u_{il}(r)v(r) \right) \d r+ -\sum_{n'l'} \sum_{k=|l-l'|}^{k=l+l'} \sqrt{2l'+1\over 2l+1}c^k(l, 0, l', 0) \int_0^\infty v(r) P_{n'l'}(r) {Y^k(u_{il}(r) P_{n'l'}(r), r)\over r} \d r = \epsilon \int_0^\infty u_{il}(r)v(r)\d r

Introducing radial basis \phi_{\mu l}(r) (where \mu, \nu labels all basis functions for the given l):

u_{il}(r) = \sum_\nu C_{\nu i l} \phi_{\nu l}(r) v(r) = \phi_{\mu l}(r)

we get (here i is again restricted for the subset corresponding to the given l):

\sum_\nu \int_0^\infty \left( \half \phi_{\mu l}'(r) \phi_{\nu l}'(r) + \left({l(l+1)\over 2r^2} -{Z\over r} + V_H(r)\right) \phi_{\mu l}(r) \phi_{\nu l}(r) \right) \d r\,C_{\nu i l}+ -\sum_\nu \sum_{n'l'} \sum_{k=|l-l'|}^{k=l+l'} \sqrt{2l'+1\over 2l+1}c^k(l, 0, l', 0) \int_0^\infty \phi_{\mu l}(r) P_{n'l'}(r) {Y^k(\phi_{\nu l}(r) P_{n'l'}(r), r)\over r} \d r\,C_{\nu i l} = \epsilon \sum_\nu \int_0^\infty \phi_{\mu l}(r)\phi_{\nu l}(r) \d r\,C_{\nu i l}

This can be written as

\sum_\nu F_{\mu\nu}^l C_{\nu i l} = \epsilon_i \sum_\nu S_{\mu\nu}^l C_{\nu i l} F_{\mu\nu}^l = H_{\mu\nu}^{l \mbox{core}} + G_{\mu\nu}^l = T_{\mu\nu}^l + V_{\mu\nu}^l + G_{\mu\nu}^l

where

T_{\mu\nu}^l = \int_0^\infty \half \phi_{\mu l}'(r) \phi_{\nu l}'(r) + \phi_{\mu l}(r){l(l+1)\over 2r^2} \phi_{\nu l}(r) \d r V_{\mu\nu}^l = \int_0^\infty \phi_{\mu l}(r)\left(-{Z\over r}\right) \phi_{\nu l}(r) \d r G_{\mu\nu}^l = \int_0^\infty \phi_{\mu l}(r) V_H(r) \phi_{\nu l}(r) \d r + -\sum_{n'l'} \sum_{k=|l-l'|}^{k=l+l'} \sqrt{2l'+1\over 2l+1}c^k(l, 0, l', 0) \int_0^\infty \phi_{\mu l}(r) P_{n'l'}(r) {Y^k(\phi_{\nu l}(r) P_{n'l'}(r), r)\over r} \d r S_{\mu\nu}^l = \int_0^\infty \phi_{\mu l}(r)\phi_{\nu l}(r) \d r V_H(r) = \sum_{n'l'} 2(2l'+1) {Y^0(P_{n'l'}^2(r), r) \over r}

The indices n' l' go over all occupied orbitals P_{n' l'}. Introducing density:

n({\bf x}) = 2 \sum_{k=1}^{N/2} | \psi_k({\bf x})|^2 = 2 \sum_{nlm} | \psi_{nlm}({\bf x})|^2 = 2 \sum_{nlm} {P_{nl}^2(r)\over r^2} |Y_{lm}(\Omega)|^2 = = 2 \sum_{nl} {P_{nl}^2(r)\over r^2} {2l+1\over 4\pi} = {1\over 4\pi} \sum_{nl} 2(2l+1) {P_{nl}^2(r)\over r^2} = n(r)

We introduce the density matrix P_{\alpha\beta}^l (where as before \alpha, \beta run over basis functions for the given l only):

P_{\alpha\beta}^l = 2 \sum_{nm} C_{\alpha nlm} C_{\beta nlm} = \sum_{n} 2(2l+1) C_{\alpha nl} C_{\beta nl}

where the C_{\alpha nl} coefficients are the same for all m corresponding to the given l. The index n runs over all occupied states for the given l. We can write n(r) as

P_{nl}(r) = \sum_\alpha C_{\alpha nl} \phi_{\alpha l}(r) n(r) = {1\over 4\pi} \sum_{nl} 2(2l+1) {P_{nl}^2(r)\over r^2} = {1\over 4\pi} \sum_{nl} 2(2l+1) \sum_{\alpha\beta} C_{\alpha nl}C_{\beta nl} {\phi_{\alpha l}(r)\phi_{\beta l}(r)\over r^2} = = {1\over 4\pi} \sum_l \sum_{\alpha\beta} {\phi_{\alpha l}(r) P_{\alpha\beta}^l \phi_{\beta l}(r)\over r^2}

Finally we get:

V_H(r) = \sum_{nl} 2(2l+1) {Y^0(P_{nl}^2(r), r) \over r} = {Y^0(4\pi n(r) r^2, r) \over r} = \sum_l \sum_{\alpha\beta} P_{\alpha\beta}^l {Y^0(\phi_{\alpha l}(r) \phi_{\beta l}(r), r) \over r} \int_0^\infty \phi_{\mu l}(r) V_H(r) \phi_{\nu l}(r) \d r = \sum_{l'} \sum_{\alpha\beta} P_{\alpha\beta}^{l'} \int_0^\infty \phi_{\mu l}(r) \phi_{\nu l}(r) {Y^0(\phi_{\alpha l'}(r) \phi_{\beta l'}(r), r) \over r}\d r = \sum_{l'} \sum_{\alpha\beta} P_{\alpha\beta}^{l'} R^0(\mu l, \beta l', \nu l, \alpha l')

and

-\sum_{n'l'} \sum_{k=|l-l'|}^{k=l+l'} \sqrt{2l'+1\over 2l+1}c^k(l, 0, l', 0) \int_0^\infty \phi_{\mu l}(r) P_{n'l'}(r) {Y^k(\phi_{\nu l}(r) P_{n'l'}(r), r)\over r} \d r = = -\sum_{n'l'} \sum_{k=|l-l'|}^{k=l+l'} 2(2l'+1) \half \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 \sum_{\alpha\beta} C_{\alpha n'l'} C_{\beta n'l'} \int_0^\infty \phi_{\mu l}(r) \phi_{\alpha l'}(r) {Y^k(\phi_{\nu l}(r) \phi_{\beta l'}(r), r)\over r} \d r = = -\half \sum_{l'} \sum_{\alpha\beta} P^{l'}_{\alpha\beta} \sum_{k=|l-l'|}^{k=l+l'} \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 R^k(\mu l, \beta l', \alpha l', \nu l)

So we get:

G_{\mu\nu}^l = \sum_{l'} \sum_{\alpha\beta} P_{\alpha\beta}^{l'} R^0(\mu l, \beta l', \nu l, \alpha l') -\half \sum_{l'} \sum_{\alpha\beta} P^{l'}_{\alpha\beta} \sum_{k=|l-l'|}^{k=l+l'} \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 R^k(\mu l, \beta l', \alpha l', \nu l) = = \sum_{l'=0}^\infty \sum_{\alpha\beta} P^{l'}_{\alpha\beta} \left( R^0(\mu l, \beta l', \nu l, \alpha l') -\half \sum_{k=|l-l'|}^{k=l+l'} \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 R^k(\mu l, \beta l', \alpha l', \nu l) \right)

The density matrix is zero if there are no occupied orbitals for the given l'.

The total energy is:

E = \half \sum_{\mu\nu} P_{\mu\nu} (H_{\mu\nu}^{\mbox{core}} + F_{\mu\nu}) = E = \half \sum_l \sum_{\mu\nu} P^l_{\mu\nu} (H^{l\mbox{core}}_{\mu\nu} + F^l_{\mu\nu}) = = \sum_l \sum_{\mu\nu} P^l_{\mu\nu} (F^l_{\mu\nu} - \half G^l_{\mu\nu}) = = \sum_l \sum_n 2(2l+1) \left( \epsilon_{n l} -\sum_{l'} \sum_{n'} (2l'+1) \left( R^0(n l, n' l', n l, n' l') -\half \sum_{k=|l-l'|}^{k=l+l'} \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 R^k(n l, n' l', n' l', n l) \right) \right)

where we used:

\sum_l \sum_{\mu\nu} P^l_{\mu\nu} F^l_{\mu\nu} = = \sum_l \sum_{\mu\nu} \sum_n 2(2l+1) C_{\mu nl}C_{\nu nl} F^l_{\mu\nu} = = \sum_l \sum_{\mu\nu} \sum_n 2(2l+1) C_{\mu nl}\epsilon_{n l} S^l_{\mu\nu} C_{\nu n l} = = \sum_l \sum_n 2(2l+1) \epsilon_{n l}

and

\half \sum_l \sum_{\mu\nu} P^l_{\mu\nu} G^l_{\mu\nu} = = \half \sum_l \sum_{\mu\nu} P^l_{\mu\nu} \sum_{l'} \sum_{\alpha\beta} P^{l'}_{\alpha\beta} \left( R^0(\mu l, \beta l', \nu l, \alpha l') -\half \sum_{k=|l-l'|}^{k=l+l'} \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 R^k(\mu l, \beta l', \alpha l', \nu l) \right) = = \half \sum_l \sum_{\mu\nu} \sum_n 2(2l+1) C_{\mu nl}C_{\nu nl} \sum_{l'} \sum_{\alpha\beta} \sum_{n'} 2(2l'+1) C_{\alpha n'l'}C_{\beta n'l'} \left( R^0(\mu l, \beta l', \nu l, \alpha l') -\half \sum_{k=|l-l'|}^{k=l+l'} \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 R^k(\mu l, \beta l', \alpha l', \nu l) \right) = = \sum_l \sum_n 2(2l+1) \sum_{l'} \sum_{n'} (2l'+1) \left( R^0(n l, n' l', n l, n' l') -\half \sum_{k=|l-l'|}^{k=l+l'} \begin{pmatrix} l & k & l' \\ 0 & 0 & 0 \end{pmatrix}^2 R^k(n l, n' l', n' l', n l) \right)

4-Index Transformation

The 4-index transformation is a way to convert the two particle integrals over basis functions (\alpha \beta | \gamma \delta) into two particle integrals over atomic (or molecular) orbitals (ij|kl):

(ij | kl) = \sum_{\alpha \beta \gamma \delta} C_{\alpha i} C_{\beta j} C_{\gamma k} C_{\delta l} (\alpha \beta | \gamma \delta)

Green’s Functions

The self energy up to a second order is given by:

\Sigma_{ij}(E) = \sum_{ars} {(ir|as) (2(jr|as)-(js|ar))\over E + \epsilon_a -\epsilon_r - \epsilon_s} + \sum_{abr} {(ia|br) (2(ja|br)-(jb|ar))\over E + \epsilon_r -\epsilon_a - \epsilon_b}

The a, b are occupied orbitals, r, s are virtual orbitals. The Dyson equation is:

(18){\bf G}(E) = {\bf G}_0(E) + {\bf G}_0(E) {\boldsymbol\Sigma}(E) {\bf G}(E)

where

{\bf G}_0(E) = {1\over E - {\boldsymbol \epsilon}} = \sum_k {\ket{k}\bra{k}\over E - \epsilon_k} {\bf G}(E) = {1\over E - {\bf E}} = \sum_k {\ket{k_{int}}\bra{k_{int}}\over E - E_k}

The \epsilon_k and \ket{k} are HF energies and eigenvectors, E_k and \ket{k_{int}} are interacting energies and eigenvectors. The matrix \boldsymbol\epsilon is a diagonal matrix of the \epsilon_k eigenvalues, \bf E is a diagonal matrix of the E_k eigenvalues. Any Green’s function {\bf G}(E) (interacting or not) can be written using the spectral density function {\bf A}(z) as follows:

{\bf G}(E) = {1\over E - {\boldsymbol \epsilon}} = \sum_k {\ket{k}\bra{k}\over E - \epsilon_k} = = \sum_k \int_{-\infty}^{\infty} {\ket{k}\bra{k}\over E - z} \delta(z-\epsilon_k) \d z = = \int_{-\infty}^{\infty} {{\bf A}(z)\over E-z} \d z

where

{\bf A}(z) = \sum_k \ket{k}\bra{k} \delta(z-\epsilon_k) = = \sum_k \ket{k}\bra{k} \lim_{\eta\to0} {\eta\over\pi} {1\over (z-\epsilon_k)^2 + \eta^2} = = \sum_k \ket{k}\bra{k} \lim_{\eta\to0} {i\over 2\pi} \left({1\over z-\epsilon_k+i\eta} -{1\over z-\epsilon_k-i\eta}\right) = = \lim_{\eta\to0} {i\over 2\pi} \left({\bf G}(z+i\eta) - {\bf G}(z-i\eta)\right)

From (18) we get:

{1\over {\bf G}_0(E)} {\bf G}(E) = 1 + {\boldsymbol \Sigma}(E) {\bf G}(E) {1\over {\bf G}_0(E)} = {1\over {\bf G}(E)} + {\boldsymbol\Sigma}(E) {\bf G}(E) = {1\over {1\over {\bf G}_0(E)} - {\boldsymbol \Sigma}(E)} = {1\over E - {\boldsymbol \epsilon} - {\boldsymbol \Sigma}(E)}

The poles E_k of the Green’s function {\bf G}(E) are then given by:

D(E) = \det (E - {\boldsymbol \epsilon} - {\boldsymbol \Sigma}(E)) = 0

or equivalently:

({\boldsymbol \Sigma}(E_k) + {\boldsymbol \epsilon})\ket{v} = E_k\ket{v}

and from the theory of matrices:

\Tr {1\over E-{\bf h}} = \Tr {{\partial\over\partial E} (E-{\bf h})\over E-{\bf h}} = \Tr {\partial\over\partial E}\log(E-{\bf h}) = {\partial\over\partial E} \Tr \log(E-{\bf h}) = = {\partial\over\partial E} \log | \det(E-{\bf h}) | = {1\over\det(E-{\bf h})} {\partial \det(E-{\bf h}) \over\partial E}

one obtains that

G_{ij}(E) = (-1)^{i+j+1} {\partial \log D(E) \over \partial ({\boldsymbol \Sigma}(E_k) + {\boldsymbol \epsilon})_{ji}}

and

\Tr {\bf G}(E) = \sum_{k} G_{kk}(E) = {\partial \log | D(E) | \over\partial E}

The number of particles N can be calculated as follows (a are occupied orbitals, k are all orbitals):

N = \sum_a | \braket{{\bf r} | a}|^2 = = \sum_a \braket{{\bf r} | a}\braket{a | {\bf r}} = = \Tr \sum_a \ket{a}\bra{a} = = {1\over 2\pi i} \Tr \int_C \sum_k {\ket{k}\bra{k}\over E-E_k}\d E = = {1\over 2\pi i} \int_C \Tr {\bf G}(E) \,\d E = = {1\over 2\pi i} \int_C {\partial \log | D(E) | \over\partial E} \,\d E

The countour C only encloses poles E_a corresponding to occupied orbitals a. Similarly for the total energy E_{tot}:

E_{tot} = \sum_a E_a | \braket{{\bf r} | a}|^2 = = \sum_a \braket{{\bf r} | a} E_a \braket{a | {\bf r}} = = \Tr \sum_a \ket{a} E_a \bra{a} = = {1\over 2\pi i} \Tr \int_C \sum_k {\ket{k}\bra{k}\over E-E_k} E\,\d E = = {1\over 2\pi i} \int_C \Tr {\bf G}(E) E\,\d E = = {1\over 2\pi i} \int_C {\partial \log | D(E) | \over\partial E} E\,\d E

For doubly filled orbitals we multiply the expressions by 2.

Table Of Contents

Previous topic

Density Functional Theory (DFT)

Next topic

Projector Augmented-Wave Method (PAW)

This Page

Navigation

© Copyright 2009-2011, Ondřej Čertík. Created using Sphinx 1.2b1.