Natural orbitals#
The one-particle density and the associated spin-orbital density matrix, \(D^\mathrm{SO}\), are related as follows
A nonzero contribution to the density requires the spin of both orbitals to be identical and the summation can therefore be separated by introducing the \(\alpha\)- and \(\beta\)-spin densities and instead run over molecular orbitals
Natural orbitals is a single set of spatial orbitals in which the one-particle density is expanded
and the elements of the diagonal density matrix, \(D_{pp}^\mathrm{NO}\), in this basis are referred to as occupation numbers.
Restricted orbitals#
In restricted calculations, we have a common set of MOs for the \(\alpha\)- and \(\beta\)-spin orbitals and can write
The density matrix is Hermitian and can be diagonalized by a unitary transformation. We get
In restricted closed (RHF) and open-shell Hartree–Fock (ROHF), the natural orbitals equal the HF orbitals and the occupation numbers equal to 0, 1, and 2 for unoccuppied, singly occupied, and doubly occupied MOs, respectively.
For the case of a multi-configurational wave function, we study the oxygen molecule as an example.
import multipsi as mtp
import numpy as np
import veloxchem as vlx
mol_str = """2
O -0.6 0.0 0.0
O 0.6 0.0 0.0
"""
molecule = vlx.Molecule.read_xyz_string(mol_str)
basis = vlx.MolecularBasis.read(molecule, "cc-pvdz", ostream=None)
molecule.set_multiplicity(3)
rohf_drv = vlx.ScfRestrictedOpenDriver()
rohf_drv.ostream.mute()
rohf_results = rohf_drv.compute(molecule, basis)
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
print(f" RHF energy: {rohf_drv.get_scf_energy() : 14.8f}")
RHF energy: -149.60946112
The MOs are occupied by nine \(\alpha\)- and seven \(\beta\)-spin electrons.
rohf_drv.molecular_orbitals.print_orbitals(molecule, basis)
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 4:
--------------------------
Occupation: 2.000 Energy: -1.09511 a.u.
( 1 O 2s : -0.36) ( 1 O 3s : -0.45) ( 2 O 2s : 0.36)
( 2 O 3s : 0.45)
Molecular Orbital No. 5:
--------------------------
Occupation: 2.000 Energy: -0.73155 a.u.
( 1 O 3s : -0.27) ( 1 O 1p+1: 0.46) ( 1 O 2p+1: 0.24)
( 2 O 3s : -0.27) ( 2 O 1p+1: -0.46) ( 2 O 2p+1: -0.24)
Molecular Orbital No. 6:
--------------------------
Occupation: 2.000 Energy: -0.70486 a.u.
( 1 O 1p-1: 0.44) ( 1 O 2p-1: 0.28) ( 2 O 1p-1: 0.44)
( 2 O 2p-1: 0.28)
Molecular Orbital No. 7:
--------------------------
Occupation: 2.000 Energy: -0.70486 a.u.
( 1 O 1p0 : 0.44) ( 1 O 2p0 : 0.28) ( 2 O 1p0 : 0.44)
( 2 O 2p0 : 0.28)
Molecular Orbital No. 8:
--------------------------
Occupation: 1.000 Energy: -0.20397 a.u.
( 1 O 1p-1: 0.51) ( 1 O 1p0 : 0.16) ( 1 O 2p-1: 0.38)
( 2 O 1p-1: -0.51) ( 2 O 1p0 : -0.16) ( 2 O 2p-1: -0.38)
Molecular Orbital No. 9:
--------------------------
Occupation: 1.000 Energy: -0.20397 a.u.
( 1 O 1p-1: 0.16) ( 1 O 1p0 : -0.51) ( 1 O 2p0 : -0.38)
( 2 O 1p-1: -0.16) ( 2 O 1p0 : 0.51) ( 2 O 2p0 : 0.38)
Molecular Orbital No. 10:
--------------------------
Occupation: 0.000 Energy: 0.47356 a.u.
( 1 O 2s : 0.19) ( 1 O 3s : 1.28) ( 1 O 1p+1: 0.40)
( 1 O 2p+1: 1.23) ( 2 O 2s : -0.19) ( 2 O 3s : -1.28)
( 2 O 1p+1: 0.40) ( 2 O 2p+1: 1.23)
Molecular Orbital No. 11:
--------------------------
Occupation: 0.000 Energy: 1.08237 a.u.
( 1 O 3s : -0.99) ( 1 O 1p+1: 0.71) ( 1 O 2p+1: -1.33)
( 2 O 3s : 0.99) ( 2 O 1p+1: 0.71) ( 2 O 2p+1: -1.33)
Molecular Orbital No. 12:
--------------------------
Occupation: 0.000 Energy: 1.10988 a.u.
( 1 O 1p-1: -0.29) ( 1 O 1p0 : -0.61) ( 1 O 2p-1: 0.27)
( 1 O 2p0 : 0.56) ( 2 O 1p-1: -0.29) ( 2 O 1p0 : -0.61)
( 2 O 2p-1: 0.27) ( 2 O 2p0 : 0.56)
Molecular Orbital No. 13:
--------------------------
Occupation: 0.000 Energy: 1.10988 a.u.
( 1 O 1p-1: -0.61) ( 1 O 1p0 : 0.29) ( 1 O 2p-1: 0.56)
( 1 O 2p0 : -0.27) ( 2 O 1p-1: -0.61) ( 2 O 1p0 : 0.29)
( 2 O 2p-1: 0.56) ( 2 O 2p0 : -0.27)
Main AO character of MOs
Active orbitals (8 electrons)
\(3\sigma_u = p_z^1 + p_z^2\)
\(1\pi_g = p_{x,y}^1 - p_{x,y}^2\)
\(1\pi_u = p_{x,y}^1 + p_{x,y}^2\)
\(3\sigma_g = p_z^1 - p_z^2\)
Inactive orbitals (8 electrons)
\(2\sigma_u = 2s^1 - 2s^2\)
\(2\sigma_g = 2s^1 + 2s^2\)
\(1\sigma_u = 1s^1 - 1s^2\)
\(1\sigma_g = 1s^1 + 1s^2\)
casscf_drv = mtp.McscfDriver()
orbital_space = mtp.OrbSpace(molecule, rohf_drv.mol_orbs)
orbital_space.cas(8, 6)
casscf_results = casscf_drv.compute(molecule, basis, orbital_space)
print(f"CAS(8,6) energy: {casscf_results["energies"][0] : 12.8f}")
CAS(8,6) energy: -149.70836718
The wave function is predominantly a sum of two determinants with \(\alpha\)-spin occupation of the \(\pi_u\) and \(\pi_g\) molecular orbitals, respectively.
casscf_results["ci_vectors"].print(0)
Determinant coef. weight
222aa0 -0.972 0.945
aa2220 0.155 0.024
The occupation numbers are (i) equal to 2 for the four inactive orbitals, (ii) close to 2 for \(3\sigma_g\) and \(1\pi_u\), (iii) close to 1 for \(1\pi_g\), and (iv) close to 0 for \(3\sigma_u\).
print(casscf_results["natural_occupations"][0])
[1.96299851 1.96299851 1.96098176 1.03637599 1.03637599 0.04026925]
Unrestricted orbitals#
In unrestricted calculations, we have separate sets of MOs for the \(\alpha\)- and \(\beta\)-spin orbitals. However, one set can be expressed in terms of the other according to
where \(S^{\alpha\beta}\) is the overlap between the two sets of orbitals that, in turn, can be expressed in terms of the overlap matrix in the AO basis
The one-particle density becomes
and the natural orbitals and occupation numbers are obtained from a diagonalization of this total density matrix.
Let us consider the unrestricted Hartree–Fock state of the oxygen molecule.
uhf_drv = vlx.ScfUnrestrictedDriver()
uhf_drv.ostream.mute()
uhf_results = uhf_drv.compute(molecule, basis)
print(f"UDFT energy: {uhf_drv.get_scf_energy() : 14.8f}")
UDFT energy: -149.62899231
norb = basis.get_dimensions_of_basis()
n_alpha_electrons = 9
n_beta_electrons = 7
S = uhf_results["S"]
C_alpha = uhf_results["C_alpha"]
C_beta = uhf_results["C_beta"]
D_alpha = np.zeros((norb, norb))
for i in range(n_alpha_electrons):
D_alpha[i, i] = 1
D_beta = np.zeros((norb, norb))
for i in range(n_beta_electrons):
D_beta[i, i] = 1
S_ab = np.einsum("ap, bq, ab -> pq", C_alpha, C_beta, S)
D = D_alpha + np.einsum("pr, qs, rs -> pq", S_ab, S_ab, D_beta)
occ_numbers, U = np.linalg.eigh(D)
print("Orbital Occ. number")
for orb, occ in enumerate(list(occ_numbers)[::-1]):
if occ > 1e-3:
print(f"{orb + 1 : >2} {occ : 12.3f}")
Orbital Occ. number
1 2.000
2 2.000
3 2.000
4 1.999
5 1.997
6 1.994
7 1.994
8 1.000
9 1.000
10 0.006
11 0.006
12 0.003
We can also get the natural orbitals and occupation numbers from the natural_orbitals method in VeloxChem.
natural_orbs = uhf_drv.natural_orbitals(uhf_results)
natural_orbs.print_orbitals(molecule, basis)
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 4:
--------------------------
Occupation: 3.998 Energy: -1.34641 a.u.
( 1 O 2s : -0.33) ( 1 O 3s : -0.37) ( 1 O 1p+1: 0.16)
( 2 O 2s : -0.33) ( 2 O 3s : -0.37) ( 2 O 1p+1: -0.16)
Molecular Orbital No. 5:
--------------------------
Occupation: 3.994 Energy: -1.09919 a.u.
( 1 O 2s : -0.36) ( 1 O 3s : -0.45) ( 2 O 2s : 0.36)
( 2 O 3s : 0.45)
Molecular Orbital No. 6:
--------------------------
Occupation: 3.988 Energy: -0.70203 a.u.
( 1 O 1p-1: 0.44) ( 1 O 2p-1: 0.28) ( 2 O 1p-1: 0.44)
( 2 O 2p-1: 0.28)
Molecular Orbital No. 7:
--------------------------
Occupation: 3.988 Energy: -0.70203 a.u.
( 1 O 1p0 : -0.44) ( 1 O 2p0 : -0.28) ( 2 O 1p0 : -0.44)
( 2 O 2p0 : -0.28)
Molecular Orbital No. 8:
--------------------------
Occupation: 2.000 Energy: -0.20095 a.u.
( 1 O 1p-1: -0.52) ( 1 O 1p0 : 0.16) ( 1 O 2p-1: -0.37)
( 2 O 1p-1: 0.52) ( 2 O 1p0 : -0.16) ( 2 O 2p-1: 0.37)
Molecular Orbital No. 9:
--------------------------
Occupation: 2.000 Energy: -0.20095 a.u.
( 1 O 1p-1: 0.16) ( 1 O 1p0 : 0.52) ( 1 O 2p0 : 0.37)
( 2 O 1p-1: -0.16) ( 2 O 1p0 : -0.52) ( 2 O 2p0 : -0.37)
Molecular Orbital No. 10:
--------------------------
Occupation: 0.012 Energy: 1.21925 a.u.
( 1 O 1p0 : 0.67) ( 1 O 2p0 : -0.58) ( 2 O 1p0 : 0.67)
( 2 O 2p0 : -0.58)
Molecular Orbital No. 11:
--------------------------
Occupation: 0.012 Energy: 1.21925 a.u.
( 1 O 1p-1: -0.67) ( 1 O 2p-1: 0.58) ( 2 O 1p-1: -0.67)
( 2 O 2p-1: 0.58)
Molecular Orbital No. 12:
--------------------------
Occupation: 0.006 Energy: 0.65627 a.u.
( 1 O 2s : 0.41) ( 1 O 3s : 0.41) ( 1 O 1p+1: 0.32)
( 1 O 2p+1: 0.74) ( 1 O 1d+2: -0.20) ( 2 O 2s : -0.41)
( 2 O 3s : -0.41) ( 2 O 1p+1: 0.32) ( 2 O 2p+1: 0.74)
( 2 O 1d+2: 0.20)
Molecular Orbital No. 13:
--------------------------
Occupation: 0.002 Energy: 2.07548 a.u.
( 1 O 1s : 0.38) ( 1 O 2s : 0.94) ( 1 O 3s : -0.80)
( 1 O 2p+1: 0.18) ( 1 O 1d+2: -0.35) ( 1 O 1d0 : 0.20)
( 2 O 1s : 0.38) ( 2 O 2s : 0.94) ( 2 O 3s : -0.80)
( 2 O 2p+1: -0.18) ( 2 O 1d+2: -0.35) ( 2 O 1d0 : 0.20)