PySCF: eigenvalues, non-integer occupation, CISD calculations, and frozen lists

Question

I have been performing CISD calculations with PySCF and I have a query surrounding the occupation of the orbitals. Here, the occupancies are calculated from the eigenvalues of the 1-electron restricted density matrix.

---

Taking the hydrogen fluoride molecule as an example, this is how the calculation is set up:

import numpy as np
from numpy import linalg as LA
from pyscf import gto, scf, lib, ci
lib.num_threads(1)

geometry = '''
      H       0.0     0.0     0.0
      F       0.0     0.0     1.1
'''
molecule = gto.Mole()
molecule.atom = geometry
molecule.basis = '3-21g'
molecule.build()

occupied_space = [occ for occ in range(int(molecule.nelectron * .5))]
print(occupied_space)

mf = scf.RHF(molecule)
mf.scf()

my_ci = ci.CISD(mf=mf,
                frozen=None,
                mo_coeff=None,
                mo_occ=None)
my_ci.kernel()
rdm1 = my_ci.make_rdm1()

evals = LA.eigvalsh(rdm1)[::-1]
print(evals)

When the CISD argument frozen = None is set, the following eigenvalues are obtained:

[1.99994686e+00 1.99455156e+00 1.99008540e+00 1.99008540e+00
 1.96069255e+00 3.79033668e-02 9.52687575e-03 9.52687575e-03
 5.33169546e-03 1.84641520e-03 5.03007156e-04]

This shows clear non-integer occupation of all the occupied/core orbitals.

If one sets frozen = [0, 1, 2, 3]; the eigenvalues/occupancies are:

[2.00000000e+00 2.00000000e+00 2.00000000e+00 2.00000000e+00
 1.99475055e+00 4.28307558e-03 7.83866124e-04 1.23523177e-04
 4.09351267e-05 1.78677364e-05 1.79726837e-07]

This shows that the HOMO is allowed to relax, smearing occupation into the virtual orbitals.

The really confusing part is when frozen = [1, 2, 3, 4]; which yields the following eigenvalues:

[2.00000000e+00 2.00000000e+00 2.00000000e+00 2.00000000e+00
 1.99998097e+00 7.30194999e-06 4.82147944e-06 3.45239124e-06
 3.45239124e-06 9.06692496e-12 1.69884163e-12]

Here, the HOMO (orbital 4, which should be frozen) has non-integer occupation, and the first unfrozen occupied orbital (orbital 0) retains complete integer occupation.

---

My question has two components. What is causing the occupation smearing in this CISD calculation (as this is not a DFT calculation), and why is the HOMO relaxing when frozen = [1, 2, 3, 4]?

Answer 1

One should not confuse the "smearing" used to facilitate self-consistent field convergence, especially for density functional calculations in periodic systems, and the role of fractional natural orbital occupation numbers that arise naturally from the many-particle nature of the exact wave function.

Natural orbitals are typically ordered in decreasing magnitude. In fact, I would bet a beer on PySCF reordering the natural orbitals so that the 2nd through 5th frozen orbital are reordered to the bottom. The occupation numbers are clearly different between the calculations, and you see that there's a significant difference in the strength of the correlation between the HOMO and the deepest core orbital: when you only correlate the HOMO, the occupation number for the lowest unoccupied natural orbital (LUNO) is $4.28\times10^{-3}$, whereas when you only correlate the F 1s orbital the LUNO is $7.30\times10^{-6}$.