Natural orbitals

Natural orbitals#

The one-particle density and the associated spin-orbital density matrix, \(D^\mathrm{SO}\), are related as follows

\[ n(\mathbf{r}) \equiv \sum_{p,q} \psi^\dagger_p(\mathbf{r}) D_{pq}^\mathrm{SO} \psi_q(\mathbf{r}) \]

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

\[\begin{align*} n(\mathbf{r}) & = \sum_{p,q} \Big( \big[ \phi^\alpha_p(\mathbf{r}) \big]^\ast D_{pq}^\alpha \phi_q^\alpha(\mathbf{r}) + \big[ \phi^\beta_p(\mathbf{r}) \big]^\ast D_{pq}^\beta \phi_q^\beta(\mathbf{r}) \Big) \end{align*}\]

Natural orbitals is a single set of spatial orbitals in which the one-particle density is expanded

\[ n(\mathbf{r}) = \sum_{p} \big[ \phi_p^\mathrm{NO}(\mathbf{r}) \big]^\ast D_{pp}^\mathrm{NO} \phi_p^\mathrm{NO}(\mathbf{r}) \]

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

\[\begin{align*} n(\mathbf{r}) & = \sum_{p,q} \phi_p^\ast(\mathbf{r}) D_{pq} \phi_q(\mathbf{r}) \\ D_{pq} & = D_{pq}^\alpha + D_{pq}^\beta \end{align*}\]

The density matrix is Hermitian and can be diagonalized by a unitary transformation. We get

\[\begin{align*} \mathbf{D}^\mathrm{NO} & = \mathbf{U}^\dagger \mathbf{D} \mathbf{U} \\ \phi_p^\mathrm{NO}(\mathbf{r}) & = \sum_q \phi_q(\mathbf{r}) U_{qp} \end{align*}\]

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)

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

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.32) (   1 O   1p0 :    -0.31) (   1 O   2p-1:    -0.20)                              
               (   1 O   2p0 :    -0.20) (   2 O   1p-1:    -0.32) (   2 O   1p0 :    -0.31)                              
               (   2 O   2p-1:    -0.20) (   2 O   2p0 :    -0.20)                                                        
                                                                                                                          
               Molecular Orbital No.   7:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:   -0.70486 a.u.                                                                  
               (   1 O   1p-1:     0.31) (   1 O   1p0 :    -0.32) (   1 O   2p-1:     0.20)                              
               (   1 O   2p0 :    -0.20) (   2 O   1p-1:     0.31) (   2 O   1p0 :    -0.32)                              
               (   2 O   2p-1:     0.20) (   2 O   2p0 :    -0.20)                                                        
                                                                                                                          
               Molecular Orbital No.   8:                                                                                 
               --------------------------                                                                                 
               Occupation: 1.000 Energy:   -0.20397 a.u.                                                                  
               (   1 O   1p-1:    -0.47) (   1 O   1p0 :     0.25) (   1 O   2p-1:    -0.35)                              
               (   1 O   2p0 :     0.18) (   2 O   1p-1:     0.47) (   2 O   1p0 :    -0.25)                              
               (   2 O   2p-1:     0.35) (   2 O   2p0 :    -0.18)                                                        
                                                                                                                          
               Molecular Orbital No.   9:                                                                                 
               --------------------------                                                                                 
               Occupation: 1.000 Energy:   -0.20397 a.u.                                                                  
               (   1 O   1p-1:     0.25) (   1 O   1p0 :     0.47) (   1 O   2p-1:     0.18)                              
               (   1 O   2p0 :     0.35) (   2 O   1p-1:    -0.25) (   2 O   1p0 :    -0.47)                              
               (   2 O   2p-1:    -0.18) (   2 O   2p0 :    -0.35)                                                        
                                                                                                                          
               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.44) (   1 O   1p0 :    -0.52) (   1 O   2p-1:     0.40)                              
               (   1 O   2p0 :     0.47) (   2 O   1p-1:    -0.44) (   2 O   1p0 :    -0.52)                              
               (   2 O   2p-1:     0.40) (   2 O   2p0 :     0.47)                                                        
                                                                                                                          
               Molecular Orbital No.  13:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    1.10988 a.u.                                                                  
               (   1 O   1p-1:     0.52) (   1 O   1p0 :    -0.44) (   1 O   2p-1:    -0.47)                              
               (   1 O   2p0 :     0.40) (   2 O   1p-1:     0.52) (   2 O   1p0 :    -0.44)                              
               (   2 O   2p-1:    -0.47) (   2 O   2p0 :     0.40)                                                        
                                                                                                                          

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

\[ \phi^\beta_q(\mathbf{r}) = \sum_p \phi^\alpha_p(\mathbf{r}) S_{pq}^{\alpha\beta} \]

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

\[ S_{pq}^{\alpha\beta} = \sum_{\delta,\gamma} \big[c_{\delta p}^\alpha\big]^\ast c^\beta_{\gamma q} S_{\delta \gamma} \]

The one-particle density becomes

\[\begin{align*} n(\mathbf{r}) & = \sum_{p,q} \big[ \phi^\alpha_p(\mathbf{r}) \big]^\ast D_{pq} \phi_q^\alpha(\mathbf{r}) \\ D_{pq} & = D_{pq}^\alpha + \sum_{r,s} \big[ S^{\alpha\beta}_{pr} \big]^\ast D_{rs}^\beta S^{\alpha\beta}_{qs} \end{align*}\]

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.43) (   1 O   2p-1:    -0.28) (   2 O   1p-1:    -0.43)                              
               (   2 O   2p-1:    -0.28)                                                                                  
                                                                                                                          
               Molecular Orbital No.   7:                                                                                 
               --------------------------                                                                                 
               Occupation: 1.994 Energy:   -0.70203 a.u.                                                                  
               (   1 O   1p0 :     0.43) (   1 O   2p0 :     0.28) (   2 O   1p0 :     0.43)                              
               (   2 O   2p0 :     0.28)                                                                                  
                                                                                                                          
               Molecular Orbital No.   8:                                                                                 
               --------------------------                                                                                 
               Occupation: 1.000 Energy:   -0.20095 a.u.                                                                  
               (   1 O   1p-1:    -0.53) (   1 O   2p-1:    -0.38) (   2 O   1p-1:     0.53)                              
               (   2 O   2p-1:     0.38)                                                                                  
                                                                                                                          
               Molecular Orbital No.   9:                                                                                 
               --------------------------                                                                                 
               Occupation: 1.000 Energy:   -0.20095 a.u.                                                                  
               (   1 O   1p0 :    -0.53) (   1 O   2p0 :    -0.38) (   2 O   1p0 :     0.53)                              
               (   2 O   2p0 :     0.38)                                                                                  
                                                                                                                          
               Molecular Orbital No.  10:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.006 Energy:    1.21925 a.u.                                                                  
               (   1 O   1p0 :     0.68) (   1 O   2p0 :    -0.59) (   2 O   1p0 :     0.68)                              
               (   2 O   2p0 :    -0.59)                                                                                  
                                                                                                                          
               Molecular Orbital No.  11:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.006 Energy:    1.21925 a.u.                                                                  
               (   1 O   1p-1:     0.68) (   1 O   2p-1:    -0.59) (   2 O   1p-1:     0.68)                              
               (   2 O   2p-1:    -0.59)                                                                                  
                                                                                                                          
               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)