The AngularMomentum operator is unable to account for spin contamination

[mrossinek]

Environment

• Qiskit Nature version: main @ b0af762
• Python version: Any
• Operating system: Any

What is happening?

Evaluating the AngularMomentum operator (a.k.a. S^2 ) does not work as intended. As it is implemented right now, this operator has the baked-in assumption, that the orbitals form an orthonormal basis. While this is true for restricted-spin calculations, this does not necessarily hold for unrestricted-spin calculations. This reflects itself in the form of spin contamination as can be seen below.

How can we reproduce the issue?

import numpy as np
from pyscf import fci, gto, scf
from qiskit.primitives import Estimator
from qiskit_algorithms.observables_evaluator import estimate_observables
from qiskit_nature.second_q.circuit.library import HartreeFock
from qiskit_nature.second_q.mappers import ParityMapper
from qiskit_nature.second_q.operators import FermionicOp
from qiskit_nature.second_q.properties import AngularMomentum, ElectronicDensity
from qiskit_nature.second_q.properties.s_operators import s_z_operator


# we set up some non-singlet spin molecule
mol = gto.M(
    atom = 'Be  0.0      0.0      0.66242; \
            Be  0.0      0.0     -0.66242;',
    basis="sto-3g",
    spin = 2,
)

# we find its unrestricted-spin HF solution
hf = scf.UHF(mol)
hf.run()
# converged SCF energy = -28.538086239077  <S^2> = 2.1286098  2S+1 = 3.0845484
# in the output above, you can see the spin contamination because S^2 is not identical to 2.0

norb = mol.nao
nelec = hf.nelec
mo_coeff = hf.mo_coeff
ovlp = hf.get_ovlp()

density = ElectronicDensity.from_particle_number(norb, nelec)
dm1a = density.alpha["+-"]
dm1b = density.beta["+-"]
dm2aa = density.alpha["++--"]
dm2bb = density.beta["++--"]
dm2ab = density.alpha_beta["++--"]

# we can reproduce the correct S^2 eigenvalue above in the MO basis with PySCF like so:
ss, ms = fci.spin_op.spin_square_general(dm1a, dm1b, dm2aa, dm2ab, dm2bb, mo_coeff, ovlp)
print(ss, ms)
# 2.1286097507302992 3.0845484277153434

# Qiskit Nature's AngularMomentum operator will give a wrong result as we will see below
ops = AngularMomentum(norb).second_q_ops()

# to find the correct result we need to take the non-unity overlap into account:
ovlpab = mo_coeff[0].T @ ovlp @ mo_coeff[1]
ovlpba = mo_coeff[1].T @ ovlp @ mo_coeff[0]

# these overlaps now form the new coefficients for S^+ and S^-
s_p = FermionicOp(
    {f"+_{idx[0]} -_{idx[1] + norb}": ovlpab[idx] for idx in np.ndindex(*ovlpab.shape)}
)
s_m = FermionicOp(
    {f"+_{idx[0] + norb} -_{idx[1]}": ovlpba[idx] for idx in np.ndindex(*ovlpba.shape)}
)

# which we can then use to build the correct S^2 operator
s_z = s_z_operator(norb)
op = s_m @ s_p + s_z @ (s_z + FermionicOp.one())

ops["CorrectAngularMomentum"] = op.simplify()

mapper = ParityMapper(nelec)
qops = mapper.map(ops)

ansatz = HartreeFock(norb, nelec, mapper)

# evaluating these observables on a pure HF initial state, should give exactly the same values as produced by PySCF above
result = estimate_observables(Estimator(), ansatz, qops)
print(result)
# {'AngularMomentum': (1.9999999999999991, {}), 'CorrectAngularMomentum': (2.12860975073022, {})}

# for completeness, we can reproduce the current faulty Qiskit Nature result with PySCF, by simply injecting identity MO coefficients:
ss, ms = fci.spin_op.spin_square_general(dm1a, dm1b, dm2aa, dm2ab, dm2bb, np.eye(norb))
print(ss, ms)
# 2.0 3.0

What should happen?

We should be able to resolve the correct spin contamination on the pure HF ground state, just like PySCF can reproduce it. To do so, we need to take the non-trivial overlap between the alpha- and beta-spin orbitals into account when constructing the S^+ and S^- operators which go into the AngularMomentum operator.

Any suggestions?

The code snippet above already contains the solution. The API of the following methods/classes will need to be adapted slightly to make this a complete fix:

• AngularMomentum
• s_plus_operator
• s_minus_operator
• s_x_operator (can be built from S^+ and S^-)
• s_y_operator (can be built from S^+ and S^-)