# Orbitals and basis sets#

## Electronic wave functions and molecular orbitals#

State vectors representing single-electron systems reside in a Hilbert space formed as a direct product

$\mathcal{V} = \mathcal{V}^\mathrm{o} \otimes \mathcal{V}^\mathrm{c} \otimes \mathcal{V}^\mathrm{s}$

The underlying vector spaces are associated with orbital, charge, and spin degrees of freedom. The elements of $$\mathcal{V}$$ can be expressed in terms of

$\begin{split} \psi^{\mathrm{L}\alpha}(\mathbf{r},t) \otimes \begin{pmatrix} 1\\0 \end{pmatrix} \otimes \begin{pmatrix} 1\\0 \end{pmatrix}; \quad \psi^{\mathrm{L}\beta}(\mathbf{r},t) \otimes \begin{pmatrix} 1\\0 \end{pmatrix} \otimes \begin{pmatrix} 0\\1 \end{pmatrix} \end{split}$
$\begin{split} \psi^{\mathrm{S}\alpha}(\mathbf{r},t) \otimes \begin{pmatrix} 0\\1 \end{pmatrix} \otimes \begin{pmatrix} 1\\0 \end{pmatrix} ; \quad \psi^{\mathrm{S}\beta}(\mathbf{r},t) \otimes \begin{pmatrix} 0\\1 \end{pmatrix} \otimes \begin{pmatrix} 0\\1 \end{pmatrix} \end{split}$

with a sum that equals a general wave function for a single-electron system, also known as a spinor,

$\begin{split} \psi(\mathbf{r},t) = \begin{pmatrix} \psi^{\mathrm{L}\alpha} \\ \psi^{\mathrm{L}\beta} \\ \psi^{\mathrm{S}\alpha} \\ \psi^{\mathrm{S}\beta} \end{pmatrix} \end{split}$

In this general case, the particle density becomes

$n(\mathbf{r},t) = \psi^\dagger \psi = \sum_\mathrm{c}^{\mathrm{L,S}} \sum_\sigma^{\alpha,\beta} \left| \psi^{\mathrm{c}\sigma}(\mathbf{r},t) \right|^2$

Spinor components $$\psi^{\mathrm{S}\alpha}$$ and $$\psi^{\mathrm{S}\beta}$$ are small and sometimes ignored, reducing the electronic wave function to two-component form

$\begin{split} \psi(\mathbf{r},t) = \begin{pmatrix} \psi^{\mathrm{L}\alpha} \\ \psi^{\mathrm{L}\beta} \end{pmatrix} \end{split}$

For a system described by a spin-independent Hamiltonian, the spatial and spin degrees of freedom are separable and we arrive at (suppressing the time variable) single-electron wave functions, or spin orbitals, of the form

$\begin{split} \psi(\mathbf{r}) = \phi(\mathbf{r}) \sigma ; \quad \sigma = c_\alpha \begin{pmatrix} 1 \\ 0 \end{pmatrix} + c_\beta \begin{pmatrix} 0 \\ 1 \end{pmatrix} \end{split}$

In the collinear approximation, all spin orbitals are eigenfunctions of the spin operator $$\hat{s}_z$$ such that every spatial function (or molecular orbital), $$\phi(\mathrm{r})$$, gives rise to two spin orbitals

$\begin{split} \psi_{a}(\mathbf{r}) = \phi_a(\mathbf{r}) \begin{pmatrix} 1 \\ 0 \end{pmatrix} ; \quad \psi_{\bar{a}}(\mathbf{r}) = \phi_a(\mathbf{r}) \begin{pmatrix} 0 \\ 1 \end{pmatrix} \end{split}$

referred to as $$\alpha$$- and $$\beta$$-spin orbitals, respectively.

### Linear combination of atomic orbitals#

In quantum chemistry, molecular orbitals (MOs) are normally expanded in a set of atom-centered basis functions, or localized atomic orbitals (AOs)

$\phi_a(\mathbf{r}) = \sum_\alpha c_{\alpha a} \chi_\alpha(\mathbf{r} - \mathbf{R}_\alpha)$

where $$\mathbf{R}_\alpha$$ denotes the atomic position center of basis function $$\chi_\alpha$$, and the expansion coefficients $$c_{\alpha a}$$ are known as molecular orbital (MO) coefficients. Using the Dirac notation, this linear combination of atomic orbitals (LCAO) expansion takes the form

$| \phi_a \rangle = \sum_\alpha | \chi_\alpha \rangle c_{\alpha a}$

It can be convenient to collect the set of MOs (and AOs) as a row vector and introduce a notation

$| \overline{\phi} \rangle = (\ldots, | \phi_a \rangle, \ldots)$

The LCAO expansion for the entire set of MOs can then be compactly written as

$| \overline{\phi} \rangle = | \overline{\chi} \rangle C$

where $$C$$ is the MO coefficient matrix.

import veloxchem as vlx

mol_str = """
C        0.00000000    0.00000000    0.00000000
O        0.00000000    0.00000000    1.43
"""
molecule = vlx.Molecule.read_str(mol_str, units='angstrom')
basis = vlx.MolecularBasis.read(molecule, "sto-3g")

scf_drv = vlx.ScfRestrictedDriver()
scf_drv.compute(molecule, basis)


### Ordering of atomic orbitals#

The VeloxChem program orders AO basis functions based on

1. Angular momentum quantum number, $$l$$

2. Angular momentum quantum number, $$m_l$$

3. User-defined order of atoms

In our present example of carbon monoxide using a mimal basis set, the order of basis functions become

$$\chi_{1s}^{C}$$, $$\chi_{2s}^{C}$$, $$\chi_{1s}^{O}$$, $$\chi_{2s}^{O}$$ $$\chi_{2py}^{C}$$, $$\chi_{2py}^{O}$$, $$\chi_{2pz}^{C}$$, $$\chi_{2pz}^{O}$$, $$\chi_{2px}^{C}$$, $$\chi_{2px}^{O}$$

import numpy as np
np.set_printoptions(precision = 3, suppress = True, linewidth = 170)

C = scf_drv.scf_tensors['C']
print('MO coefficient matrix:\n', C)

MO coefficient matrix:
[[ 0.001  0.994  0.091  0.228  0.109  0.     0.     0.     0.     0.105]
[-0.004  0.024 -0.26  -0.836 -0.47  -0.    -0.    -0.    -0.    -0.536]
[ 0.995  0.     0.231 -0.114  0.036 -0.    -0.     0.    -0.    -0.089]
[ 0.023 -0.005 -0.856  0.526 -0.196  0.     0.    -0.     0.     0.547]
[ 0.     0.     0.    -0.     0.    -0.095 -0.42  -0.034  0.912 -0.   ]
[-0.    -0.    -0.     0.     0.    -0.187 -0.826  0.02  -0.548  0.   ]
[-0.004  0.004 -0.141  0.003  0.621  0.     0.     0.    -0.    -0.942]
[-0.003  0.002  0.112  0.38  -0.578 -0.    -0.     0.    -0.    -0.852]
[-0.     0.    -0.    -0.     0.    -0.42   0.095 -0.912 -0.034 -0.   ]
[-0.    -0.    -0.    -0.     0.    -0.826  0.187  0.548  0.02   0.   ]]


### Overlap matrix#

The overlap between two basis functions is

$S_{\alpha\beta} = \langle \chi_\alpha | \chi_\beta \rangle$

and the overlap matrix can be written as

$S = \langle \overline{\chi} | \overline{\chi} \rangle$

where

$\langle \overline{\chi} | = (\ldots, \langle \chi_\alpha |, \ldots)^T$
S = scf_drv.scf_tensors['S']
print('Overlap matrix:\n', S)

Overlap matrix:
[[ 1.     0.248  0.     0.017  0.     0.     0.    -0.031  0.     0.   ]
[ 0.248  1.     0.022  0.255  0.     0.     0.    -0.251  0.     0.   ]
[ 0.     0.022  1.     0.237  0.     0.     0.036  0.     0.     0.   ]
[ 0.017  0.255  0.237  1.     0.     0.     0.337  0.     0.     0.   ]
[ 0.     0.     0.     0.     1.     0.132  0.     0.     0.     0.   ]
[ 0.     0.     0.     0.     0.132  1.     0.     0.     0.     0.   ]
[ 0.     0.     0.036  0.337  0.     0.     1.    -0.285  0.     0.   ]
[-0.031 -0.251  0.     0.     0.     0.    -0.285  1.     0.     0.   ]
[ 0.     0.     0.     0.     0.     0.     0.     0.     1.     0.132]
[ 0.     0.     0.     0.     0.     0.     0.     0.     0.132  1.   ]]


## Basis sets#

Our software comes distibuted with several of the standard basis sets and families of basis sets. For a given element, the available basis sets in the VeloxChem library can be provided from the MolecularBasis class object. Should a basis set not be included in the software distribution, it can be readily supplied by the user in a separate file in the working directory.

print('Available basis sets for carbon:\n', vlx.MolecularBasis.get_avail_basis('C'))

Available basis sets for carbon:
['6-31++G', '6-31+G', '6-31+G**', '6-311++G', '6-311++G(2D,2P)', '6-311++G**', '6-311+G', '6-311+G(2D,P)', '6-311G', '6-311G(2DF,2PD)', '6-31G', '6-31G(2DF,P)', 'AUG-CC-PCVDZ', 'AUG-CC-PCVQZ', 'AUG-CC-PCVTZ', 'AUG-CC-PVDZ', 'AUG-CC-PVTZ', 'CC-PCVDZ', 'CC-PVDZ', 'CC-PVDZ-RI', 'CC-PVTZ', 'CC-PVTZ-RI', 'D-AUG-CC-PVDZ', 'D-AUG-CC-PVQZ', 'D-AUG-CC-PVTZ', 'DEF2-QZVP', 'DEF2-QZVPD', 'DEF2-QZVPP', 'DEF2-QZVPPD', 'DEF2-RI-J', 'DEF2-SV(P)', 'DEF2-SVP', 'DEF2-SVPD', 'DEF2-TZVP', 'DEF2-TZVPP', 'DEF2-TZVPPD', 'MIN-CC-PVDZ', 'SADLEJ-PVTZ', 'STO-3G', 'STO-6G']


### Basis set assignment#

A basis set can be assigned to a molecule with the read method.

basis = vlx.MolecularBasis.read(molecule, 'cc-pVDZ')
print(basis.get_string('cc-pVDZ', molecule))

Molecular Basis (cc-pVDZ)
===========================

Basis: CC-PVDZ

Atom Contracted GTOs          Primitive GTOs

O   (3S,2P,1D)               (17S,4P,1D)
C   (3S,2P,1D)               (17S,4P,1D)

Contracted Basis Functions : 28
Primitive Basis Functions  : 68


The moleular basis set label is available with use of the get_label method.

vlx.MolecularBasis.get_label(basis)

'CC-PVDZ'


### Number of basis functions#

The number of contracted and primitive basis functions are available from the get_dimensions_of_basis and get_dimensions_of_primitive_basis methods, respectively.

nbas = vlx.MolecularBasis.get_dimensions_of_basis(basis, molecule)
nprim = vlx.MolecularBasis.get_dimensions_of_primitive_basis(basis, molecule)

print('Number of contracted basis functions:', nbas)
print('Number of primitive basis functions:', nprim)

Number of contracted basis functions: 28
Number of primitive basis functions: 68