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 1p0 : -0.45) ( 1 O 2p0 : -0.28) ( 2 O 1p0 : -0.45)
( 2 O 2p0 : -0.28)
Molecular Orbital No. 7:
--------------------------
Occupation: 2.000 Energy: -0.70486 a.u.
( 1 O 1p-1: 0.45) ( 1 O 2p-1: 0.28) ( 2 O 1p-1: 0.45)
( 2 O 2p-1: 0.28)
Molecular Orbital No. 8:
--------------------------
Occupation: 1.000 Energy: -0.20397 a.u.
( 1 O 1p0 : -0.53) ( 1 O 2p0 : -0.40) ( 2 O 1p0 : 0.53)
( 2 O 2p0 : 0.40)
Molecular Orbital No. 9:
--------------------------
Occupation: 1.000 Energy: -0.20397 a.u.
( 1 O 1p-1: -0.53) ( 1 O 2p-1: -0.40) ( 2 O 1p-1: 0.53)
( 2 O 2p-1: 0.40)
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.57) ( 1 O 1p0 : -0.37) ( 1 O 2p-1: -0.52)
( 1 O 2p0 : 0.34) ( 2 O 1p-1: 0.57) ( 2 O 1p0 : -0.37)
( 2 O 2p-1: -0.52) ( 2 O 2p0 : 0.34)
Molecular Orbital No. 13:
--------------------------
Occupation: 0.000 Energy: 1.10988 a.u.
( 1 O 1p-1: -0.37) ( 1 O 1p0 : -0.57) ( 1 O 2p-1: 0.34)
( 1 O 2p0 : 0.52) ( 2 O 1p-1: -0.37) ( 2 O 1p0 : -0.57)
( 2 O 2p-1: 0.34) ( 2 O 2p0 : 0.52)
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_drv.get_energy() : 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.
print(casscf_drv.CIVecs[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\).
casscf_drv.get_natural_occupations()
array([1.96299851, 1.96299851, 1.96098176, 1.03637598, 1.03637598,
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(molecule)
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: 1.999 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: 1.997 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: 1.994 Energy: -0.70203 a.u.
( 1 O 1p-1: -0.15) ( 1 O 1p0 : -0.42) ( 1 O 2p0 : -0.27)
( 2 O 1p-1: -0.15) ( 2 O 1p0 : -0.42) ( 2 O 2p0 : -0.27)
Molecular Orbital No. 7:
--------------------------
Occupation: 1.994 Energy: -0.70203 a.u.
( 1 O 1p-1: 0.42) ( 1 O 1p0 : -0.15) ( 1 O 2p-1: 0.27)
( 2 O 1p-1: 0.42) ( 2 O 1p0 : -0.15) ( 2 O 2p-1: 0.27)
Molecular Orbital No. 8:
--------------------------
Occupation: 1.000 Energy: -0.20095 a.u.
( 1 O 1p-1: -0.40) ( 1 O 1p0 : 0.37) ( 1 O 2p-1: -0.28)
( 1 O 2p0 : 0.26) ( 2 O 1p-1: 0.40) ( 2 O 1p0 : -0.37)
( 2 O 2p-1: 0.28) ( 2 O 2p0 : -0.26)
Molecular Orbital No. 9:
--------------------------
Occupation: 1.000 Energy: -0.20095 a.u.
( 1 O 1p-1: 0.37) ( 1 O 1p0 : 0.40) ( 1 O 2p-1: 0.26)
( 1 O 2p0 : 0.28) ( 2 O 1p-1: -0.37) ( 2 O 1p0 : -0.40)
( 2 O 2p-1: -0.26) ( 2 O 2p0 : -0.28)
Molecular Orbital No. 10:
--------------------------
Occupation: 0.006 Energy: 1.21925 a.u.
( 1 O 1p-1: -0.64) ( 1 O 1p0 : 0.22) ( 1 O 2p-1: 0.56)
( 1 O 2p0 : -0.19) ( 2 O 1p-1: -0.64) ( 2 O 1p0 : 0.22)
( 2 O 2p-1: 0.56) ( 2 O 2p0 : -0.19)
Molecular Orbital No. 11:
--------------------------
Occupation: 0.006 Energy: 1.21925 a.u.
( 1 O 1p-1: 0.22) ( 1 O 1p0 : 0.64) ( 1 O 2p-1: -0.19)
( 1 O 2p0 : -0.56) ( 2 O 1p-1: 0.22) ( 2 O 1p0 : 0.64)
( 2 O 2p-1: -0.19) ( 2 O 2p0 : -0.56)
Molecular Orbital No. 12:
--------------------------
Occupation: 0.003 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.001 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)