Restricted versus unrestricted#
High spin states#
Single-reference methods#
By a high-spin state, we refer to an open-shell system where all unpaired electrons are chosen to have \(\alpha\)-spin. As detailed in section Electron spin, this case is particularly simple since the system can be described by a single Slater determinant.
There are two main ways to proceed:
Keep the molecular orbitals of the \(\alpha\) and \(\beta\) spin-orbitals identical but populating the \(\alpha\) spin-orbitals with more electrons. This is known as a restricted open-shell calculation.
Allow for different molecular orbitals for the \(\alpha\) and \(\beta\) spin-orbitals. This is known as an unrestricted calculation.
Let us illustrate these two approaches on the most important open-shell molecule on earth namely the oxygen molecule.
import multipsi as mtp
import veloxchem as vlx
au2kcal = 627.51
mol_str = """
O 0.0 0.0 -0.6
O 0.0 0.0 0.6
"""
molecule = vlx.Molecule.read_molecule_string(mol_str, units="angstrom")
basis = vlx.MolecularBasis.read(molecule, "cc-pvdz")
We determine the triplet ground state with the RODFT and UDFT approaches.
molecule.set_multiplicity(3)
triplet_rodft_drv = vlx.ScfRestrictedOpenDriver()
triplet_rodft_drv.xcfun = "B3LYP"
triplet_rodft_results = triplet_rodft_drv.compute(molecule, basis)
triplet_udft_drv = vlx.ScfUnrestrictedDriver()
triplet_udft_drv.xcfun = "B3LYP"
triplet_udft_results = triplet_udft_drv.compute(molecule, basis)
Owing to the larger flexibility in the reference state, the UDFT energy is slightly lower than the RODFT one.
print(f"ROB3LYP energy: {triplet_rodft_drv.get_scf_energy() : 14.8f}")
print(f" UB3LYP energy: {triplet_udft_drv.get_scf_energy() : 14.8f}")
ROB3LYP energy: -150.33014560
UB3LYP energy: -150.33392791
In the RODFT case, we have one set of MOs with the first seven being doubly occupied and MOs eight and nine being singly occupied.
triplet_rodft_drv.molecular_orbitals.print_orbitals(molecule, basis)
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 4:
--------------------------
Occupation: 2.000 Energy: -0.78650 a.u.
( 1 O 2s : -0.36) ( 1 O 3s : -0.46) ( 2 O 2s : 0.36)
( 2 O 3s : 0.46)
Molecular Orbital No. 5:
--------------------------
Occupation: 2.000 Energy: -0.52369 a.u.
( 1 O 3s : 0.28) ( 1 O 1p0 : -0.47) ( 1 O 2p0 : -0.23)
( 2 O 3s : 0.28) ( 2 O 1p0 : 0.47) ( 2 O 2p0 : 0.23)
Molecular Orbital No. 6:
--------------------------
Occupation: 2.000 Energy: -0.50972 a.u.
( 1 O 1p+1: 0.44) ( 1 O 2p+1: 0.26) ( 2 O 1p+1: 0.44)
( 2 O 2p+1: 0.26)
Molecular Orbital No. 7:
--------------------------
Occupation: 2.000 Energy: -0.50972 a.u.
( 1 O 1p-1: 0.44) ( 1 O 2p-1: 0.26) ( 2 O 1p-1: 0.44)
( 2 O 2p-1: 0.26)
Molecular Orbital No. 8:
--------------------------
Occupation: 1.000 Energy: -0.20092 a.u.
( 1 O 1p+1: -0.46) ( 1 O 1p-1: 0.28) ( 1 O 2p+1: -0.34)
( 1 O 2p-1: 0.20) ( 2 O 1p+1: 0.46) ( 2 O 1p-1: -0.28)
( 2 O 2p+1: 0.34) ( 2 O 2p-1: -0.20)
Molecular Orbital No. 9:
--------------------------
Occupation: 1.000 Energy: -0.20092 a.u.
( 1 O 1p+1: -0.28) ( 1 O 1p-1: -0.46) ( 1 O 2p+1: -0.20)
( 1 O 2p-1: -0.34) ( 2 O 1p+1: 0.28) ( 2 O 1p-1: 0.46)
( 2 O 2p+1: 0.20) ( 2 O 2p-1: 0.34)
Molecular Orbital No. 10:
--------------------------
Occupation: 0.000 Energy: 0.23772 a.u.
( 1 O 2s : 0.21) ( 1 O 3s : 1.06) ( 1 O 1p0 : 0.50)
( 1 O 2p0 : 1.00) ( 2 O 2s : -0.21) ( 2 O 3s : -1.06)
( 2 O 1p0 : 0.50) ( 2 O 2p0 : 1.00)
Molecular Orbital No. 11:
--------------------------
Occupation: 0.000 Energy: 0.79595 a.u.
( 1 O 3s : 1.12) ( 1 O 1p0 : -0.65) ( 1 O 2p0 : 1.46)
( 2 O 3s : -1.12) ( 2 O 1p0 : -0.65) ( 2 O 2p0 : 1.46)
Molecular Orbital No. 12:
--------------------------
Occupation: 0.000 Energy: 0.83473 a.u.
( 1 O 1p+1: 0.21) ( 1 O 1p-1: -0.64) ( 1 O 2p+1: -0.20)
( 1 O 2p-1: 0.60) ( 2 O 1p+1: 0.21) ( 2 O 1p-1: -0.64)
( 2 O 2p+1: -0.20) ( 2 O 2p-1: 0.60)
Molecular Orbital No. 13:
--------------------------
Occupation: 0.000 Energy: 0.83473 a.u.
( 1 O 1p+1: 0.64) ( 1 O 1p-1: 0.21) ( 1 O 2p+1: -0.60)
( 1 O 2p-1: -0.20) ( 2 O 1p+1: 0.64) ( 2 O 1p-1: 0.21)
( 2 O 2p+1: -0.60) ( 2 O 2p-1: -0.20)
In the UDFT case, there are two different sets of MOs with nine occupied \(\alpha\)-spin orbitals and seven occupied \(\beta\)-spin orbitals.
triplet_udft_drv.molecular_orbitals.print_orbitals(molecule, basis)
Spin Unrestricted Alpha Orbitals
--------------------------------
Molecular Orbital No. 5:
--------------------------
Occupation: 1.000 Energy: -0.56129 a.u.
( 1 O 1p-1: -0.46) ( 1 O 2p-1: -0.25) ( 2 O 1p-1: -0.46)
( 2 O 2p-1: -0.25)
Molecular Orbital No. 6:
--------------------------
Occupation: 1.000 Energy: -0.56129 a.u.
( 1 O 1p+1: 0.46) ( 1 O 2p+1: 0.25) ( 2 O 1p+1: 0.46)
( 2 O 2p+1: 0.25)
Molecular Orbital No. 7:
--------------------------
Occupation: 1.000 Energy: -0.54325 a.u.
( 1 O 3s : 0.28) ( 1 O 1p0 : -0.47) ( 1 O 2p0 : -0.22)
( 2 O 3s : 0.28) ( 2 O 1p0 : 0.47) ( 2 O 2p0 : 0.22)
Molecular Orbital No. 8:
--------------------------
Occupation: 1.000 Energy: -0.29937 a.u.
( 1 O 1p+1: -0.28) ( 1 O 1p-1: 0.46) ( 1 O 2p+1: -0.21)
( 1 O 2p-1: 0.33) ( 2 O 1p+1: 0.28) ( 2 O 1p-1: -0.46)
( 2 O 2p+1: 0.21) ( 2 O 2p-1: -0.33)
Molecular Orbital No. 9:
--------------------------
Occupation: 1.000 Energy: -0.29937 a.u.
( 1 O 1p+1: 0.46) ( 1 O 1p-1: 0.28) ( 1 O 2p+1: 0.33)
( 1 O 2p-1: 0.21) ( 2 O 1p+1: -0.46) ( 2 O 1p-1: -0.28)
( 2 O 2p+1: -0.33) ( 2 O 2p-1: -0.21)
Molecular Orbital No. 10:
--------------------------
Occupation: 0.000 Energy: 0.21614 a.u.
( 1 O 2s : 0.21) ( 1 O 3s : 1.02) ( 1 O 1p0 : 0.51)
( 1 O 2p0 : 0.98) ( 2 O 2s : -0.21) ( 2 O 3s : -1.02)
( 2 O 1p0 : 0.51) ( 2 O 2p0 : 0.98)
Molecular Orbital No. 11:
--------------------------
Occupation: 0.000 Energy: 0.78644 a.u.
( 1 O 3s : 1.11) ( 1 O 1p0 : -0.65) ( 1 O 2p0 : 1.46)
( 2 O 3s : -1.11) ( 2 O 1p0 : -0.65) ( 2 O 2p0 : 1.46)
Molecular Orbital No. 12:
--------------------------
Occupation: 0.000 Energy: 0.80682 a.u.
( 1 O 1p+1: -0.30) ( 1 O 1p-1: 0.59) ( 1 O 2p+1: 0.28)
( 1 O 2p-1: -0.57) ( 2 O 1p+1: -0.30) ( 2 O 1p-1: 0.59)
( 2 O 2p+1: 0.28) ( 2 O 2p-1: -0.57)
Molecular Orbital No. 13:
--------------------------
Occupation: 0.000 Energy: 0.80682 a.u.
( 1 O 1p+1: 0.59) ( 1 O 1p-1: 0.30) ( 1 O 2p+1: -0.57)
( 1 O 2p-1: -0.28) ( 2 O 1p+1: 0.59) ( 2 O 1p-1: 0.30)
( 2 O 2p+1: -0.57) ( 2 O 2p-1: -0.28)
Molecular Orbital No. 14:
--------------------------
Occupation: 0.000 Energy: 0.89114 a.u.
( 1 O 3s : 0.33) ( 1 O 1p0 : 0.58) ( 1 O 2p0 : -0.82)
( 2 O 3s : 0.33) ( 2 O 1p0 : -0.58) ( 2 O 2p0 : 0.82)
Spin Unrestricted Beta Orbitals
-------------------------------
Molecular Orbital No. 3:
--------------------------
Occupation: 1.000 Energy: -1.25526 a.u.
( 1 O 2s : -0.32) ( 1 O 3s : -0.26) ( 1 O 1p0 : -0.19)
( 2 O 2s : -0.32) ( 2 O 3s : -0.26) ( 2 O 1p0 : 0.19)
Molecular Orbital No. 4:
--------------------------
Occupation: 1.000 Energy: -0.74211 a.u.
( 1 O 2s : -0.35) ( 1 O 3s : -0.45) ( 2 O 2s : 0.35)
( 2 O 3s : 0.45)
Molecular Orbital No. 5:
--------------------------
Occupation: 1.000 Energy: -0.50336 a.u.
( 1 O 3s : 0.29) ( 1 O 1p0 : -0.46) ( 1 O 2p0 : -0.23)
( 2 O 3s : 0.29) ( 2 O 1p0 : 0.46) ( 2 O 2p0 : 0.23)
Molecular Orbital No. 6:
--------------------------
Occupation: 1.000 Energy: -0.45864 a.u.
( 1 O 1p-1: 0.42) ( 1 O 2p-1: 0.28) ( 2 O 1p-1: 0.42)
( 2 O 2p-1: 0.28)
Molecular Orbital No. 7:
--------------------------
Occupation: 1.000 Energy: -0.45864 a.u.
( 1 O 1p+1: 0.42) ( 1 O 2p+1: 0.28) ( 2 O 1p+1: 0.42)
( 2 O 2p+1: 0.28)
Molecular Orbital No. 8:
--------------------------
Occupation: 0.000 Energy: -0.10410 a.u.
( 1 O 1p+1: 0.21) ( 1 O 1p-1: 0.46) ( 1 O 2p+1: 0.18)
( 1 O 2p-1: 0.41) ( 2 O 1p+1: -0.21) ( 2 O 1p-1: -0.46)
( 2 O 2p+1: -0.18) ( 2 O 2p-1: -0.41)
Molecular Orbital No. 9:
--------------------------
Occupation: 0.000 Energy: -0.10410 a.u.
( 1 O 1p+1: -0.46) ( 1 O 1p-1: 0.21) ( 1 O 2p+1: -0.41)
( 1 O 2p-1: 0.18) ( 2 O 1p+1: 0.46) ( 2 O 1p-1: -0.21)
( 2 O 2p+1: 0.41) ( 2 O 2p-1: -0.18)
Molecular Orbital No. 10:
--------------------------
Occupation: 0.000 Energy: 0.26034 a.u.
( 1 O 2s : 0.21) ( 1 O 3s : 1.10) ( 1 O 1p0 : 0.49)
( 1 O 2p0 : 1.03) ( 2 O 2s : -0.21) ( 2 O 3s : -1.10)
( 2 O 1p0 : 0.49) ( 2 O 2p0 : 1.03)
Molecular Orbital No. 11:
--------------------------
Occupation: 0.000 Energy: 0.80590 a.u.
( 1 O 3s : 1.12) ( 1 O 1p0 : -0.66) ( 1 O 2p0 : 1.45)
( 2 O 3s : -1.12) ( 2 O 1p0 : -0.66) ( 2 O 2p0 : 1.45)
Molecular Orbital No. 12:
--------------------------
Occupation: 0.000 Energy: 0.86451 a.u.
( 1 O 1p+1: 0.46) ( 1 O 1p-1: -0.51) ( 1 O 2p+1: -0.42)
( 1 O 2p-1: 0.46) ( 2 O 1p+1: 0.46) ( 2 O 1p-1: -0.51)
( 2 O 2p+1: -0.42) ( 2 O 2p-1: 0.46)
It is a common procedure to determine a Kohn–Sham DFT value estimating the expectation value of the total spin operator in the same manner as would be the case for a Hartree–Fock wave function.
Note
Since the two-particle density is not directly available in DFT, the expectation value of a two-electron operator such as e.g. the total spin operator is also not available.
print(
f"<S^2> RODFT = {triplet_rodft_drv.compute_s2(molecule, triplet_rodft_results) : 8.6f}"
)
print(
f"<S^2> UDFT = {triplet_udft_drv.compute_s2(molecule, triplet_udft_results) : 8.6f}"
)
<S^2> RODFT = 2.000000
<S^2> UDFT = 2.006227
The ground state of oxygen is a triplet with spin quantum number \(S = 1\) and
This is exactly what we get with a restricted open-shell description, but not in the unrestricted case. In general, an unrestricted state is not a eigenfunction of \(\hat{S}^2\) and one refers to this situation by saying that the state is spin contaminated.
The relative simplicity of the unrestricted scheme together with its rather good performance makes it more widely used than the restricted open-shell form.
Note
The additional flexibility of unrestricted molecular orbitals does not lower the energy of closed-shell restricted states.
Multi-reference methods#
A multi-configurational method such as CASSCF can also be used and we note that a CASSCF wave function with only the high-spin open-shell MOs in the active space is equivalent to restricted open-shell Hartree–Fock (ROHF).
triplet_rohf_drv = vlx.ScfRestrictedOpenDriver()
triplet_rohf_results = triplet_rohf_drv.compute(molecule, basis)
# By default, OrbSpace includes all open-shells in the active space
space = mtp.OrbSpace(molecule, triplet_rohf_drv.mol_orbs)
triplet_mcscf_drv = mtp.McscfDriver()
mcscf_results = triplet_mcscf_drv.compute(molecule, basis, space)
print(f" ROHF energy: {triplet_rohf_drv.get_scf_energy() : 14.8f}")
print(f"CASSCF(2,2) energy: {triplet_mcscf_drv.get_energy() : 14.8f}")
ROHF energy: -149.60946112
CASSCF(2,2) energy: -149.60946112
Low spin states#
Single-reference methods: broken symmetry#
To address low spin open-shell systems, we can study singlet states of the oxygen molecule. The two lowest singlet states are \(a^1\Delta_g\) and \(b^1\Sigma^+_g\) and experimentally they are found to be 22.5 and 37.5 kcal/mol above the ground triplet state, \(X^3\Sigma_g^-\), respectively.
Let us first compute the lowest closed-shell singlet state using DFT.
molecule.set_multiplicity(1)
singlet_rdft_drv = vlx.ScfRestrictedDriver()
singlet_rdft_drv.xcfun = "B3LYP"
singlet_rdft_results = singlet_rdft_drv.compute(molecule, basis)
The \(\pi\)-orbitals in the closed-shell state are no longer doubly degenerate.
singlet_rdft_drv.molecular_orbitals.print_orbitals(molecule, basis)
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 4:
--------------------------
Occupation: 2.000 Energy: -0.79052 a.u.
( 1 O 2s : -0.36) ( 1 O 3s : -0.46) ( 2 O 2s : 0.36)
( 2 O 3s : 0.46)
Molecular Orbital No. 5:
--------------------------
Occupation: 2.000 Energy: -0.52746 a.u.
( 1 O 3s : 0.28) ( 1 O 1p0 : -0.47) ( 1 O 2p0 : -0.23)
( 2 O 3s : 0.28) ( 2 O 1p0 : 0.47) ( 2 O 2p0 : 0.23)
Molecular Orbital No. 6:
--------------------------
Occupation: 2.000 Energy: -0.51773 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. 7:
--------------------------
Occupation: 2.000 Energy: -0.50968 a.u.
( 1 O 1p-1: 0.46) ( 1 O 2p-1: 0.27) ( 2 O 1p-1: 0.46)
( 2 O 2p-1: 0.27)
Molecular Orbital No. 8:
--------------------------
Occupation: 2.000 Energy: -0.24089 a.u.
( 1 O 1p-1: 0.53) ( 1 O 2p-1: 0.41) ( 2 O 1p-1: -0.53)
( 2 O 2p-1: -0.41)
Molecular Orbital No. 9:
--------------------------
Occupation: 0.000 Energy: -0.17013 a.u.
( 1 O 1p+1: -0.52) ( 1 O 2p+1: -0.42) ( 2 O 1p+1: 0.52)
( 2 O 2p+1: 0.42)
Molecular Orbital No. 10:
--------------------------
Occupation: 0.000 Energy: 0.23373 a.u.
( 1 O 2s : -0.21) ( 1 O 3s : -1.05) ( 1 O 1p0 : -0.50)
( 1 O 2p0 : -0.99) ( 2 O 2s : 0.21) ( 2 O 3s : 1.05)
( 2 O 1p0 : -0.50) ( 2 O 2p0 : -0.99)
Molecular Orbital No. 11:
--------------------------
Occupation: 0.000 Energy: 0.79320 a.u.
( 1 O 3s : 1.12) ( 1 O 1p0 : -0.65) ( 1 O 2p0 : 1.46)
( 2 O 3s : -1.12) ( 2 O 1p0 : -0.65) ( 2 O 2p0 : 1.46)
Molecular Orbital No. 12:
--------------------------
Occupation: 0.000 Energy: 0.82293 a.u.
( 1 O 1p+1: -0.67) ( 1 O 2p+1: 0.63) ( 2 O 1p+1: -0.67)
( 2 O 2p+1: 0.63)
Molecular Orbital No. 13:
--------------------------
Occupation: 0.000 Energy: 0.84078 a.u.
( 1 O 1p-1: -0.67) ( 1 O 2p-1: 0.63) ( 2 O 1p-1: -0.67)
( 2 O 2p-1: 0.63)
print("Closed-shell singlet relative to triplet ground state (in kcal/mol)")
print(
f" restricted triplet state: {(singlet_rdft_drv.get_scf_energy() - triplet_rodft_drv.get_scf_energy()) * au2kcal : 14.8f}"
)
print(
f" unrestricted triplet state: {(singlet_rdft_drv.get_scf_energy() - triplet_udft_drv.get_scf_energy()) * au2kcal : 14.8f}"
)
Closed-shell singlet relative to triplet ground state (in kcal/mol)
restricted triplet state: 36.92128398
unrestricted triplet state: 39.29472211
These energies match reasonably well with the higher \(b^1\Sigma^+_g\) singlet state.
There are two different ways to reach the lower singlet using DFT. The first one is to compute the (negative) excitation energy using TD-DFT.
rsp_drv = vlx.TDAExciDriver()
rsp_drv.update_settings({"nstates": 1}, {"xcfun": "B3LYP"})
rsp_results = rsp_drv.compute(molecule, basis, singlet_rdft_results)
exc_energy = rsp_results["eigenvalues"][0]
print(f"Excitation energy (in kcal/mol): {exc_energy * au2kcal : 14.8}\n")
print("Open-shell singlet energy relative to triplet ground state (in kcal/mol)")
print(
f" restricted triplet state: {(singlet_rdft_drv.get_scf_energy() + exc_energy - triplet_rodft_drv.get_scf_energy()) * au2kcal : 14.8f}"
)
print(
f" unrestricted triplet state: {(singlet_rdft_drv.get_scf_energy() + exc_energy - triplet_udft_drv.get_scf_energy()) * au2kcal : 14.8f}"
)
Excitation energy (in kcal/mol): -11.732361
Open-shell singlet energy relative to triplet ground state (in kcal/mol)
restricted triplet state: 25.18892342
unrestricted triplet state: 27.56236155
These energies is in better agreement with the lowest singlet state. We have here reached the open-shell singlet state by means of a TD-DFT calculation based on a closed-shell reference state that is higher in energy.
Let us instead compute the lower singlet state in a direct manner using unrestricted DFT.
singlet_udft_drv = vlx.ScfUnrestrictedDriver()
singlet_udft_drv.xcfun = "B3LYP"
singlet_udft_result = singlet_udft_drv.compute(molecule, basis)
print(f"RDFT singlet: {singlet_rdft_drv.get_scf_energy() : 14.8f}")
print(f"UDFT singlet: {singlet_udft_drv.get_scf_energy() : 14.8f}")
RDFT singlet: -150.27130784
UDFT singlet: -150.27130784
Since the closed-shell state represent a stable point, the UDFT solution will become identical.
In the absence of magnetic fields, there is no impetus for the MOs of the \(\alpha\) and \(\beta\) spin-orbitals to differ and the SCF optimization gets stuck in a symmetric solution. In the high-spin case, on the other hand, this symmetry was broken by the fact that we had an excess of \(\alpha\)-spin electrons.
However, we can break the symmetry by introducing different starting occupations for the \(\alpha\) and \(\beta\) spin-orbitals and remain in this region during the SCF optimization with the maximum-overlap method.
# A list of which orbitals are occupied in alpha and beta (starting from the triplet orbitals)
a_occ = [0, 1, 2, 3, 4, 5, 6, 7]
b_occ = [0, 1, 2, 3, 4, 5, 6, 8]
broken_udft_drv = vlx.ScfUnrestrictedDriver()
broken_udft_drv.xcfun = "B3LYP"
broken_udft_drv.maximum_overlap(
molecule, basis, triplet_rodft_drv.mol_orbs, a_occ, b_occ
)
broken_udft_results = broken_udft_drv.compute(molecule, basis)
We can look at the natural orbitals and see that the state indeed corresponds to two singly occupied \(\pi\)-orbitals.
natural_orbs = broken_udft_drv.natural_orbitals(broken_udft_results)
natural_orbs.print_orbitals(molecule, basis)
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 4:
--------------------------
Occupation: 2.000 Energy: -1.01237 a.u.
( 1 O 2s : -0.32) ( 1 O 3s : -0.38) ( 2 O 2s : -0.32)
( 2 O 3s : -0.38)
Molecular Orbital No. 5:
--------------------------
Occupation: 2.000 Energy: -0.80124 a.u.
( 1 O 2s : -0.36) ( 1 O 3s : -0.46) ( 2 O 2s : 0.36)
( 2 O 3s : 0.46)
Molecular Orbital No. 6:
--------------------------
Occupation: 1.999 Energy: -0.51076 a.u.
( 1 O 1p+1: -0.34) ( 1 O 1p-1: 0.30) ( 1 O 2p+1: -0.20)
( 1 O 2p-1: 0.18) ( 2 O 1p+1: -0.34) ( 2 O 1p-1: 0.30)
( 2 O 2p+1: -0.20) ( 2 O 2p-1: 0.18)
Molecular Orbital No. 7:
--------------------------
Occupation: 1.999 Energy: -0.51076 a.u.
( 1 O 1p+1: 0.30) ( 1 O 1p-1: 0.34) ( 1 O 2p+1: 0.18)
( 1 O 2p-1: 0.20) ( 2 O 1p+1: 0.30) ( 2 O 1p-1: 0.34)
( 2 O 2p+1: 0.18) ( 2 O 2p-1: 0.20)
Molecular Orbital No. 8:
--------------------------
Occupation: 1.000 Energy: -0.20217 a.u.
( 1 O 1p+1: -0.52) ( 1 O 2p+1: -0.38) ( 2 O 1p+1: 0.52)
( 2 O 2p+1: 0.38)
Molecular Orbital No. 9:
--------------------------
Occupation: 1.000 Energy: -0.20217 a.u.
( 1 O 1p-1: 0.52) ( 1 O 2p-1: 0.38) ( 2 O 1p-1: -0.52)
( 2 O 2p-1: -0.38)
Molecular Orbital No. 10:
--------------------------
Occupation: 0.001 Energy: 0.91738 a.u.
( 1 O 1p+1: 0.65) ( 1 O 1p-1: 0.18) ( 1 O 2p+1: -0.58)
( 1 O 2p-1: -0.16) ( 2 O 1p+1: 0.65) ( 2 O 1p-1: 0.18)
( 2 O 2p+1: -0.58) ( 2 O 2p-1: -0.16)
Molecular Orbital No. 11:
--------------------------
Occupation: 0.001 Energy: 0.91738 a.u.
( 1 O 1p+1: -0.18) ( 1 O 1p-1: 0.65) ( 1 O 2p+1: 0.16)
( 1 O 2p-1: -0.58) ( 2 O 1p+1: -0.18) ( 2 O 1p-1: 0.65)
( 2 O 2p+1: 0.16) ( 2 O 2p-1: -0.58)
Molecular Orbital No. 12:
--------------------------
Occupation: 0.000 Energy: 2.65807 a.u.
( 1 O 1d+2: 0.34) ( 1 O 1d-2: -0.64) ( 2 O 1d+2: -0.34)
( 2 O 1d-2: 0.64)
Molecular Orbital No. 13:
--------------------------
Occupation: 0.000 Energy: 2.34841 a.u.
( 1 O 1d+2: -0.33) ( 1 O 1d-2: 0.61) ( 2 O 1d+2: -0.33)
( 2 O 1d-2: 0.61)
The energy of this singlet relative to the triplet ground state becomes:
print(
f"Broken symmetry singlet-triplet gap: {(broken_udft_drv.get_scf_energy() - triplet_udft_drv.get_scf_energy()) * au2kcal : 10.8}"
)
Broken symmetry singlet-triplet gap: 10.308167
With this trick, we have thus found a lower energy singlet, albeit by introducing spin contamination.
print(
f"<S^2> UDFT = {broken_udft_drv.compute_s2(molecule, broken_udft_results) : 8.6f}"
)
<S^2> UDFT = 1.003282
We also note that the energy gap is now well below that of the experimental value. The reason is that, while high spin open-shells are single-determinant in nature, open-shell singlets and low-spin triplets are not.
A minimum of two determinant is needed for physically correct wave functions
Our single determinant reference state is seen to correspond to a 50/50 mix of the singlet and triplet wave functions as reflected in the spin contamination of the UDFT state.
True singlets and triplets have spin expectation values of \(0\) and \(2 \hbar^2\), respectively, whereas we obtained a value close to \(\hbar^2\).
However, knowing this fact, we can actually correct the energy. Since we know the energy of the triplet and of the broken-symmetry solution, we can estimate the energy of the true singlet.
print(
f"Estimated singlet-triplet gap: {2 * (broken_udft_drv.get_scf_energy() - triplet_udft_drv.get_scf_energy()) * au2kcal : 12.8f}"
)
Estimated singlet-triplet gap: 20.61633427
This is now very close to the experimental energy of 22.5 kcal/mol.
Our technique to obtain the “true” singlet energy is a special case of the more general “weighted-averaged broken symmetry” (WABS) method, which can in principle correct any broken symmetry solution with knowledge of the high-spin energy and the expectation value of the total spin operator.
Multi-reference solution#
Arguably a multi-reference approach is simpler. For O\(_2\), we create an active space with two electrons in two orbitals.
space = mtp.OrbSpace(molecule, singlet_rdft_drv.mol_orbs)
space.cas(2, 2)
singlet_mcscf_drv = mtp.McscfDriver()
mcscf_results = singlet_mcscf_drv.compute(molecule, basis, space)
print(
f"Lowest singlet relative energy: {(singlet_mcscf_drv.get_energy(0) - triplet_mcscf_drv.get_energy()) * au2kcal : 14.8f}"
)
Lowest singlet relative energy: 29.69991698
The CASSCF automatically finds the open-shell to be the lowest state.
The CASSCF code in MultiPsi even allows to compute singlets and triplets simultaneously. We only need to deactivate the spin restriction that is applied by defaults to singlets and compute several states at the same time.
space = mtp.OrbSpace(molecule, singlet_rdft_drv.mol_orbs)
space.spin_restricted = False
space.cas(2, 2)
all_mcscf_drv = mtp.McscfDriver()
mcscf_results = all_mcscf_drv.compute(molecule, basis, space, n_states=4)
print(
f"First singlet relative energy : {(all_mcscf_drv.get_energy(1) - all_mcscf_drv.get_energy(0)) * au2kcal : 14.8f}"
)
print(
f"First singlet relative energy : {(all_mcscf_drv.get_energy(2) - all_mcscf_drv.get_energy(0)) * au2kcal : 14.8f}"
)
print(
f"Second singlet relative energy: {(all_mcscf_drv.get_energy(3) - all_mcscf_drv.get_energy(0)) * au2kcal : 14.8f}"
)
First singlet relative energy : 29.55599147
First singlet relative energy : 29.55599147
Second singlet relative energy: 59.11198294
In a single calculation, the CASSCF finds the lowest triplet state and the two components of the doubly degenerate \(a^1\Delta_g\) state followed by a higher lying singlet.
The energies are not a great match with experiment due to the fact that a CASSCF calculation with so small an active space does not include correlation. A much better result is found by expanding the active space to include all atomic 2p orbitals.
space.cas(8, 6)
large_mcscf_drv = mtp.McscfDriver()
mcscf_results = large_mcscf_drv.compute(molecule, basis, space, n_states=4)
print(
f"First singlet relative energy : {(large_mcscf_drv.get_energy(1) - large_mcscf_drv.get_energy(0)) * au2kcal :14.8f}"
)
print(
f"First singlet relative energy : {(large_mcscf_drv.get_energy(2) - large_mcscf_drv.get_energy(0)) * au2kcal : 14.8f}"
)
print(
f"Second singlet relative energy: {(large_mcscf_drv.get_energy(3) - large_mcscf_drv.get_energy(0)) * au2kcal : 14.8f}"
)
First singlet relative energy : 21.00271262
First singlet relative energy : 21.00271262
Second singlet relative energy: 39.02178823
If we examine the CASSCF wave functions, we clearly see the multi-determinant nature of the states.
for i, vec in enumerate(all_mcscf_drv.CIVecs):
print("Wavefunction for state", i + 1)
print(vec)
Wavefunction for state 1
Determinant coef. weight
ab 0.707 0.500
ba -0.707 0.500
Wavefunction for state 2
Determinant coef. weight
20 -0.584 0.341
ab -0.398 0.159
ba -0.398 0.159
02 0.584 0.341
Wavefunction for state 3
Determinant coef. weight
20 -0.398 0.159
ab 0.584 0.341
ba 0.584 0.341
02 0.398 0.159
Wavefunction for state 4
Determinant coef. weight
20 -0.707 0.500
02 -0.707 0.500
The determinant notation lists the occupation of the two orbitals in order, with “0” meaning “empty”, “a” and “b” meaning respectively singly occupied with an \(\alpha\) or \(\beta\) electron, and “2” meaning doubly occupied. Clearly in this case, every single state in the calculation is an equal-weight superposition of two determinants.