Electronic structure calculations
The PyQInt package provides two flavors of Hartree-Fock electronic structure calculations: restricted Hartree-Fock (RHF) and unrestricted Hartree-Fock (UHF). The restricted formulation is suitable for closed-shell systems, whereas the unrestricted formulation allows for the treatment of open-shell systems with unpaired electrons.
Both approaches expose detailed intermediate results, allowing the user to inspect orbital coefficients, density matrices, Fock matrices, and energy contributions throughout the self-consistent field (SCF) procedure.
Restricted Hartree-Fock
Below, an example of an restricted Hartree-Fock calculation is given.
from pyqint import PyQInt, MoleculeBuilder, HF
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
def main():
# calculate sto3g coefficients for h2o
cgfs, coeff = calculate_co()
# visualize orbitals
fig, ax = plt.subplots(2,5, figsize=(12, 5), dpi=144)
sz = 3
for i in range(0,2):
for j in range(0,5):
dens = plot_wavefunction(cgfs, coeff[:,i*5+j], sz=sz)
limit = max(abs(np.min(dens)), abs(np.max(dens)) )
im = ax[i,j].contourf(dens, origin='lower',
extent=[-sz, sz, -sz, sz], cmap='PiYG', vmin=-limit, vmax=limit,
levels=11)
im = ax[i,j].contour(dens, origin='lower', colors='black',
extent=[-sz, sz, -sz, sz], vmin=-limit, vmax=limit,
levels=11)
ax[i,j].set_xlabel('x [Bohr]')
ax[i,j].set_ylabel('z [Bohr]')
ax[i,j].set_aspect('equal', adjustable='box')
ax[i,j].set_xticks(np.linspace(-3,3, 7))
ax[i,j].set_yticks(np.linspace(-3,3, 7))
ax[i,j].grid(linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()
def calculate_co():
mol = MoleculeBuilder.from_name('CO')
result = HF(mol, 'sto3g').rhf()
return result['cgfs'], result['orbc']
def plot_wavefunction(cgfs, coeff, sz=3.5):
# build integrator
integrator = PyQInt()
# build grid
x = np.linspace(-sz, sz, 150)
z = np.linspace(-sz, sz, 150)
xx, zz = np.meshgrid(x,z)
yy = np.zeros(len(x) * len(z))
grid = np.vstack([xx.flatten(), yy, zz.flatten()]).reshape(3,-1).T
res = integrator.plot_wavefunction(grid, coeff, cgfs).reshape((len(z), len(x)))
return res
if __name__ == '__main__':
main()
Canonical molecular orbitals of CO visualized using contour plots.
Hint
The manual construction of contour plots shown above provides full control
over grid generation and visualization. For convenience, PyQInt
also provides the ContourPlotter helper class (see
ContourPlotter helper class) which encapsulates this workflow and allows
rapid generation of grids of molecular orbital contour plots with minimal
boilerplate code.
Result dictionary (RHF)
The result of a Hartree-Fock calculation is captured inside a dictionary object. This dictionary objects contains the following keys
Key |
Description |
|---|---|
|
Final energy of the electronic structure calculation |
|
Whether the convergence criteria were met |
|
Number of SCF iterations |
|
List of elements and their position in Bohr units |
|
List of contracted Gaussian functional objects |
|
List of energies during the self-convergence procedure |
|
Orbital energies (converged) (array of M elements) |
|
Orbital coefficients (converged) (matrix of N x M elements) |
|
Density matrix \(\mathbf{P}\) |
|
Fock matrix \(\mathbf{F}\) |
|
Unitary transformation matrix \(\mathbf{X}\) |
|
Overlap matrix \(\mathbf{S}\) |
|
Kinetic energy matrix \(\mathbf{T}\) |
|
Nuclear attraction matrix \(\mathbf{V}\) |
|
Core Hamiltonian matrix \(\mathbf{H_\textrm{core}}\) |
|
Two-electron tensor object \((i,j,k,l)\) |
|
Time statistics object |
|
Sum of kinetic and nuclear attraction energy |
|
Total kinetic energy |
|
Total nuclear attraction energy |
|
Total electron-electron repulsion energy |
|
Total exchange energy |
|
Electrostatic repulsion energy of the nuclei |
|
Total number of electrons |
|
Molecule class |
|
Forces on the atoms (if calculated, else |
To provide an example how one can use the above data, let us consider the situation wherein the user wants to decompose the individual components of the total energy as given by
Via the script below, one can easily verify that the above equation holds and that the total energy is indeed the sum of the kinetic, nuclear attraction, electron-electron repulsion, exchange and nuclear repulsion energies within a Hartree-Fock calculation.
from pyqint import MoleculeBuilder,HF
mol = MoleculeBuilder.from_name('ch4')
mol.name = 'CH4'
res = HF(mol, 'sto3g').rhf()
print()
print('Kinetic energy: ', res['ekin'])
print('Nuclear attraction energy: ', res['enuc'])
print('Electron-electron repulsion: ', res['erep'])
print('Exchange energy: ', res['ex'])
print('Repulsion between nuclei: ', res['enucrep'])
print()
print('Total energy: ', res['energy'])
print('Sum of the individual terms: ',
res['ekin'] + res['enuc'] + res['erep'] + res['ex'] + res['enucrep'])
The output of the above script yields:
Kinetic energy: 39.426137798761935
Nuclear attraction energy: -118.63789212479075
Electron-electron repulsion: 32.73242708923475
Exchange energy: -6.609004693243019
Repulsion between nuclei: 13.362026647057352
Total energy: -39.726305282979084
Sum of the individual terms: -39.72630528297975
Convergence
The self-consistent field procedure is considered converged when, between two
subsequent iterations, the total energy changes by less than tolerance
(default \(10^{-9}\) Ht) and the largest element of the commutator
\(\mathbf{FPS} - \mathbf{SPF}\), expressed in the orthogonalized basis,
is smaller than tolerance_grad (default
\(\sqrt{\texttt{tolerance}}\)). This commutator vanishes only when the
density matrix is self-consistent with the Fock matrix it produces, and it is
also the error vector used by the DIIS procedure.
When these criteria are not met within itermax iterations (default 100),
a RuntimeWarning is issued. The converged entry of the result
dictionary reports the outcome and niter the number of iterations used.
res = HF(mol, 'sto3g').rhf()
if not res['converged']:
raise RuntimeError('SCF did not converge')
By default, canonical orthogonalization is used. It removes those linear
combinations of basis functions whose overlap eigenvalue is smaller than
\(10^{-7}\), as these make the calculation numerically unstable. In that
case, there are fewer molecular orbitals \(M\) than basis functions
\(N\). Symmetric orthogonalization (ortho='symmetric') keeps all
basis functions and warns when the basis set is nearly linearly dependent.
The rhf routine requires an even number of electrons; open-shell
systems are treated with uhf.
Unrestricted Hartree-Fock (UHF)
For systems containing unpaired electrons, such as radicals or atoms with open
shells, the restricted Hartree-Fock approximation is no longer appropriate.
In such cases, PyQInt provides an unrestricted Hartree-Fock (UHF)
implementation via the uhf() method.
In UHF, separate spatial orbitals are used for spin-up (\(\alpha\)) and spin-down (\(\beta\)) electrons. This leads to distinct α and β density matrices and Fock matrices, allowing the electronic structure to exhibit spin polarization.
As a representative example, consider the methyl radical (CH₃), which has one unpaired electron and therefore requires an unrestricted treatment. The geometry used below places the carbon atom at the origin, with three hydrogen atoms arranged in a planar configuration with H-C-H angles of 120 degrees and a C-H bond length of 2.039 Bohr.
from pyqint import Molecule, HF
import numpy as np
R = 2.039
sqrt3 = np.sqrt(3.0)
mol = Molecule()
mol.add_atom('C', 0.0, 0.0, 0.0)
mol.add_atom('H', R, 0.0, 0.0)
mol.add_atom('H', -0.5 * R, 0.5 * sqrt3 * R, 0.0)
mol.add_atom('H', -0.5 * R, -0.5 * sqrt3 * R, 0.0)
# CH3 is a doublet (multiplicity = 2)
res = HF(mol, 'sto3g').uhf(multiplicity=2)
print('Total energy:', res['energy'])
Result dictionary (UHF)
The result of an unrestricted Hartree-Fock (UHF) calculation is likewise captured inside a dictionary object. Compared to restricted Hartree-Fock, the UHF result dictionary contains spin-resolved quantities, reflecting the fact that different spatial orbitals are used for spin-up (\(\alpha\)) and spin-down (\(\beta\)) electrons.
Below, the most important entries of the UHF result dictionary are listed.
Key |
Description |
|---|---|
|
Final total electronic energy |
|
Whether the convergence criteria were met |
|
Number of SCF iterations |
|
List of elements and their position in Bohr units |
|
List of contracted Gaussian functional objects |
|
List of total energies during the SCF convergence procedure |
|
Converged orbital energies for spin-up (\(\alpha\)) electrons |
|
Converged orbital energies for spin-down (\(\beta\)) electrons |
|
Orbital coefficient matrix for \(\alpha\) electrons |
|
Orbital coefficient matrix for \(\beta\) electrons |
|
Density matrix \(\mathbf{P}^{\alpha}\) for spin-up electrons |
|
Density matrix \(\mathbf{P}^{\beta}\) for spin-down electrons |
|
Total density matrix \(\mathbf{P} = \mathbf{P}^{\alpha} + \mathbf{P}^{\beta}\) |
|
Fock matrix for \(\alpha\) electrons |
|
Fock matrix for \(\beta\) electrons |
|
Unitary transformation matrix \(\mathbf{X}\) |
|
Overlap matrix \(\mathbf{S}\) |
|
Kinetic energy matrix \(\mathbf{T}\) |
|
Nuclear attraction matrix \(\mathbf{V}\) |
|
Core Hamiltonian matrix \(\mathbf{H_{\textrm{core}}}\) |
|
Two-electron integral tensor \((i,j,k,l)\) |
|
Time statistics object |
In unrestricted Hartree-Fock, the exchange contribution is spin-dependent. Accordingly, the exchange energy is split into separate \(\alpha\) and \(\beta\) components.
Key |
Description |
|---|---|
|
Sum of kinetic and nuclear attraction energy |
|
Coulomb (Hartree) energy |
|
Exchange energy contribution from \(\alpha\) electrons |
|
Exchange energy contribution from \(\beta\) electrons |
|
Electrostatic repulsion energy of the nuclei |
In addition to the electronic structure data, the UHF result dictionary contains explicit information about the spin state of the system.
Key |
Description |
|---|---|
|
Total number of electrons |
|
Number of spin-up (\(\alpha\)) electrons |
|
Number of spin-down (\(\beta\)) electrons |
|
Spin multiplicity \((2S + 1)\) |
|
Expectation value \(\langle \hat{S}^{2} \rangle\) of the UHF wavefunction |
|
Exact value \(S(S+1)\) for the requested spin state |
Note
Unlike restricted Hartree-Fock, unrestricted Hartree-Fock allows for spin
polarization. As a consequence, the resulting wavefunction is generally not
an exact eigenfunction of the total spin operator, which may lead to spin
contamination. The amount of spin contamination is given by the
difference between s2 and s2_exact, where
with \(S_z = (N_{\alpha} - N_{\beta})/2\). For the methyl radical above, \(\langle \hat{S}^{2} \rangle = 0.7652\), slightly above the exact value of 0.75 for a doublet.
The number of electrons may be set explicitly via the nelec argument,
e.g. to model ions. The multiplicity must be compatible with the number of
electrons: an even number of electrons requires an odd multiplicity and vice
versa. Otherwise, a ValueError is raised.
Warning
The UHF procedure starts from the core Hamiltonian guess. For some
open-shell systems, e.g. triplet O2, this guess may lead to a
self-consistent solution which is not the ground state. Comparing the
result with a calculation started from different initial orbitals (via the
orbc_init argument) is advised in such cases.
Custom basis sets
Besides the basis sets offered by PyQInt, one can also use a custom
basis set defined by the user. The rhf routine accepts either a basis set
for its basis argument, or alternatively a list of cgf objects.
In the example code shown below, the latter is done.
from pyqint import Molecule, HF, CGF
mol = Molecule()
mol.add_atom('H', 0.0000, 0.0000, 0.3561150187, unit='angstrom')
mol.add_atom('H', 0.0000, 0.0000, -0.3561150187, unit='angstrom')
nuclei = mol.get_nuclei()
cgfs = []
for n in nuclei:
cgf = CGF(n[0])
cgf.add_gto(0.154329, 3.425251, 0, 0, 0)
cgf.add_gto(0.535328, 0.623914, 0, 0, 0)
cgf.add_gto(0.444635, 0.168855, 0, 0, 0)
cgfs.append(cgf)
res = HF(mol, basis=cgfs).rhf(verbose=True)
Hint
A nice website to find a large collection of Gaussian Type basis set coefficients is https://www.basissetexchange.org/.