1273 unfortunately has not been fixed properly.

This issue resurfaced while I was looking into some more extended tests. Specifically, using the overlap matrix which I obtained for a (6,4) active space of O2:

ovlpab = np.asarray(
    [
        [0.987256, -0.001123, 0.00006, -0.0],
        [-0.001123, -0.987256, -0.0, -0.00006],
        [0.000019, 0.000055, -0.3195, -0.931662],
        [0.000056, -0.000019, -0.931662, 0.3195],
    ],
)

Basically, comparing the computed S^2 values for different occupations obtained from PySCF's FCI module and Qiskit Nature still produce vastly different results.

How can we reproduce the issue?

import numpy as np

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 JordanWignerMapper
from qiskit_nature.second_q.properties import AngularMomentum, ElectronicDensity

norb = 4
nelec = (4, 2)
ovlpab = np.asarray(
    [
        [0.987256, -0.001123, 0.00006, -0.0],
        [-0.001123, -0.987256, -0.0, -0.00006],
        [0.000019, 0.000055, -0.3195, -0.931662],
        [0.000056, -0.000019, -0.931662, 0.3195],
    ],
)

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

ovlpaa = np.eye(ovlpab.shape[0])
ovlpbb = np.eye(ovlpab.shape[0])
ovlpba = ovlpab.T

# manual computation just like done inside of `pyscf.fci.spin_op`
ssz = (
    np.einsum("ijkl,ij,kl->", dm2aa, ovlpaa, ovlpaa)
    - np.einsum("ijkl,ij,kl->", dm2ab, ovlpaa, ovlpbb)
    + np.einsum("ijkl,ij,kl->", dm2bb, ovlpbb, ovlpbb)
    - np.einsum("ijkl,ij,kl->", dm2ab, ovlpaa, ovlpbb)
    + np.einsum("ji,ij->", dm1a, ovlpaa)
    + np.einsum("ji,ij->", dm1b, ovlpbb)
) * 0.25

dm2abba = -dm2ab.transpose(0, 3, 2, 1)  # alpha^+ beta^+ alpha beta
dm2baab = -dm2ab.transpose(2, 1, 0, 3)  # beta^+ alpha^+ beta alpha
ssxy = (
    np.einsum("ijkl,ij,kl->", dm2abba, ovlpab, ovlpba)
    + np.einsum("ijkl,ij,kl->", dm2baab, ovlpba, ovlpab)
    + np.einsum("ji,ij->", dm1a, ovlpaa)
    + np.einsum("ji,ij->", dm1b, ovlpbb)
) * 0.5
spin_square = ssxy + ssz

s = np.sqrt(spin_square + 0.25) - 0.5
mult = s * 2 + 1
print("FCI", spin_square, mult)

ang_mom = AngularMomentum(norb, ovlpab).second_q_ops()

mapper = JordanWignerMapper()
qubit_op = mapper.map(ang_mom)

hf_state = HartreeFock(norb, nelec, mapper)

result = estimate_observables(Estimator(), hf_state, qubit_op)
new_ss = result["AngularMomentum"][0]
new_s = np.sqrt(new_ss + 0.25) - 0.5
new_mult = new_s * 2 + 1
print("Qiskit", new_ss, new_mult)

This results in:

FCI 2.050648651787 3.0335778557914086
Qiskit 2.0 3.0

What should happen?

Qiskit should produce the same result as the PySCF implementation.

Any suggestions?

Basically, there are multiple problems:

the formula $S^2 = S^- S^+ + S^z (S^z + 1)$ assumes that there are more alpha than beta spin particles. One needs to take the symmetric version to allow either case: $S^2 = 0.5 * (S^- S^+ + S^+ S^-) + S^z S^z$
on top of that, the logic that was added in #1275 is actually flawed. Through mostly trial and error, I found a patch last night which I admittedly did not understand completely but it made all my test cases pass :eyes: Essentially, I ended up overwriting the coefficients of the one-body terms which result from normal-ordering the $S^+ S^-$ (and v.v.) product to all be 1 (instead of the coefficient that resulted from the overlap squared).

Today, I found out why (thanks also to @AlbertoBaiardi for the extensive discussions). The problem is, that $S^+$ and $S^-$ do NOT need to the include the overlap as coefficients! Rather, we need to adapt our anti-commutation relation and weight this with the overlap. This is explained quite nicely in sections 4.1 and 4.2 of this paper of which I will include the relevant screenshots below:

screenshot_1701338968 screenshot_1701338987 screenshot_1701339006

Long story short, to fix this issue we need to:

[x] use the symmetric formula for $S^2$
~revert the inclusion of the overlap as coefficients inside $S^+$ and $S^-$ (and $S^x$ and $S^y$ by extension)~ WRONG: see next comment
~integrate the non-unity anti-commutation relation directly into the FermionicOp~ also wrong

It turns out, that there are more things at play here...

The good news

The Qiskit calculation is mostly correct. In fact, all that is needed is to use a symmetric formula for $S^2$ which also allows the number of beta-spin particles to be larger than the alpha ones (rather than vice versa).

The "bad" news

The problem here is rather of conceptual nature. The problem is, that using unrestricted-spin orbitals for the active space, can result in the active alpha- and beta-spin orbitals to span different spaces. This then result in a non-unitary alpha-beta overlap matrix, which in turn actually breaks an assumption on which PySCF's implementation is based (which is where I obtained those numpy.einsum equations from). Fixing it is possible by simply changing the following two lines:

    + np.einsum("ji,ij->", dm1a, ovlpaa)
    + np.einsum("ji,ij->", dm1b, ovlpbb)

to become:

    + np.einsum("ji,ij->", dm1a, ovlpab @ ovlpba)
    + np.einsum("ji,ij->", dm1b, ovlpba @ ovlpab)

This properly accounts for the non-unitary nature of ovlpab. I tested this in a more extensive example in the snippet below.

A code snippet to compute the active and inactive $S^2$ expectation values in different active spaces of $O_2$.

```python from copy import copy import numpy as np from pyscf import 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 JordanWignerMapper, ParityMapper from qiskit_nature.second_q.properties import AngularMomentum, ElectronicDensity from qiskit_nature.second_q.properties.s_operators import s_z_operator from qiskit_nature.second_q.operators import FermionicOp def diff_operators(total_op, active_op, padding): """Takes the difference of two operators by padding the latter to the size of the former.""" ops = {} for name in total_op.keys(): padded_op = FermionicOp( { " ".join( f"{term[0]}_{term[1] + padding * (1 + (term[1] // (active_op[name].num_spin_orbitals // 2)))}" for term in terms ): coeff for terms, coeff in active_op[name].terms() }, num_spin_orbitals=total_op[name].num_spin_orbitals, ) ops[name] = (total_op[name] - padded_op).simplify() return ops def compute_s2_qiskit(norb, nelec, ops_2nd): """Computes S^2 on the HF state using Qiskit (Nature).""" mapper = JordanWignerMapper() qubit_op = mapper.map(ops_2nd) hf_state = HartreeFock(norb, nelec, mapper) result = estimate_observables(Estimator(), hf_state, qubit_op) ss = result["AngularMomentum"][0] ssz = result["S_z^2"][0] return ss, ssz def compute_s2_numpy(hf_density, ovlpab, ovlpaa=None, ovlpbb=None): """Computes S^2 on the HF density and overlap using numpy.""" dm1a = hf_density.alpha["+-"] dm1b = hf_density.beta["+-"] dm2aa = hf_density.alpha["++--"] dm2bb = hf_density.beta["++--"] dm2ab = hf_density.alpha_beta["++--"] if ovlpaa is None: ovlpaa = np.eye(ovlpab.shape[0]) if ovlpbb is None: ovlpbb = np.eye(ovlpab.shape[1]) ovlpba = ovlpab.T ssz = ( np.einsum("ijkl,ij,kl->", dm2aa, ovlpaa, ovlpaa) - np.einsum("ijkl,ij,kl->", dm2ab, ovlpaa, ovlpbb) + np.einsum("ijkl,ij,kl->", dm2bb, ovlpbb, ovlpbb) - np.einsum("ijkl,ij,kl->", dm2ab, ovlpaa, ovlpbb) + np.einsum("ji,ij->", dm1a, ovlpaa) + np.einsum("ji,ij->", dm1b, ovlpbb) ) * 0.25 dm2abba = -dm2ab.transpose(0, 3, 2, 1) # alpha^+ beta^+ alpha beta dm2baab = -dm2ab.transpose(2, 1, 0, 3) # beta^+ alpha^+ beta alpha ssxy = ( np.einsum("ijkl,ij,kl->", dm2abba, ovlpab, ovlpba) + np.einsum("ijkl,ij,kl->", dm2baab, ovlpba, ovlpab) + np.einsum("ji,ij->", dm1a, ovlpab @ ovlpba) # ovlpaa) + np.einsum("ji,ij->", dm1b, ovlpba @ ovlpab) # ovlpbb) ) * 0.5 ss = ssxy + ssz return ss, ssz # we use the triplet O2 as a test case m = gto.M(atom="O 0 0 0; O 0 0 1.3", basis="sto3g", spin=2) # UHF has a little bit of spin contamination for this system: h = scf.UHF(m).run() # converged SCF energy = -147.620355111774 = 2.0033214 2S+1 = 3.0022135 ref_ss, _ = h.spin_square() def verify_results(spin_square, source, label): """Verifies the result against the reference values.""" if np.isclose(spin_square, ref_ss): print(f" The {label} computed by {source} is correct!") return True else: print(f" ERROR: The {label} computed by {source} is WRONG!") print(f" = {spin_square:.5f} Δ = {ref_ss - spin_square:.5f}") return False mo_a = h.mo_coeff[0] mo_b = h.mo_coeff[1] occ_a = h.mo_occ[0] occ_b = h.mo_occ[1] ovlp = h.get_ovlp() ovlpab_tot = mo_a.T @ ovlp @ mo_b nelec_tot = (int(sum(occ_a)), int(sum(occ_b))) norb_tot = ovlpab_tot.shape[0] print("\n===== References ======") # as a reference, we (re-)compute the total S^2 values tot_2nd_ops = AngularMomentum(norb_tot, ovlpab_tot).second_q_ops() tot_s_z = s_z_operator(norb_tot) tot_2nd_ops["S_z^2"] = tot_s_z @ tot_s_z tot_ss_qiskit, tot_ssz_qiskit = compute_s2_qiskit(norb_tot, nelec_tot, tot_2nd_ops) _ = verify_results(tot_ss_qiskit, "Qiskit", "total") density_tot = ElectronicDensity.from_particle_number(norb_tot, nelec_tot, include_rdm2=True) tot_ss_numpy, tot_ssz_numpy = compute_s2_numpy(density_tot, ovlpab_tot) _ = verify_results(tot_ss_numpy, "Numpy", "total") # we now look towards active spaces and gradually shrink it by removing more and more inner orbitals for ninact in range(10): print(f"\n--- {ninact} Inactive Orbitals ---", flush=True) # lists of inactive and active orbital indices inact = list(range(ninact)) act = list(range(ninact, norb_tot)) # active parameters norb_act = len(act) nelec_act = (int(sum(occ_a[act])), int(sum(occ_b[act]))) ovlpab_act = mo_a[:, act].T @ ovlp @ mo_b[:, act] # compute the active value with Qiskit active_2nd_ops = AngularMomentum(norb_act, ovlpab_act).second_q_ops() act_s_z = s_z_operator(norb_act) active_2nd_ops["S_z^2"] = act_s_z @ act_s_z act_ss_qiskit, act_ssz_qiskit = compute_s2_qiskit(norb_act, nelec_act, active_2nd_ops) # we can compute the missing based on the difference of the total and active operators diff_ops = diff_operators(tot_2nd_ops, active_2nd_ops, ninact) diff_ss_qiskit, diff_ssz_qiskit = compute_s2_qiskit(norb_tot, nelec_tot, diff_ops) # we can now sum the active and inactive to verify against the reference sum_ss_qiskit = act_ss_qiskit + diff_ss_qiskit _ = verify_results(sum_ss_qiskit, "Qiskit", "active + inactive") # now we compute the active value with Numpy density_act = ElectronicDensity.from_particle_number(norb_act, nelec_act, include_rdm2=True) act_ss_numpy, act_ssz_numpy = compute_s2_numpy(density_act, ovlpab_act) # to obtain the inactive contribution while also including the active+inactive terms, we # must remove the active+active sector from the overlap used by Numpy ovlpaa_diff = np.eye(norb_tot) ovlpaa_diff[ninact:, ninact:] = 0 ovlpbb_diff = np.eye(norb_tot) ovlpbb_diff[ninact:, ninact:] = 0 ovlpab_diff = copy(ovlpab_tot) ovlpab_diff[ninact:, ninact:] = 0 diff_ss_numpy, diff_ssz_numpy = compute_s2_numpy(density_tot, ovlpab_diff, ovlpaa_diff, ovlpbb_diff) sum_ss_numpy = act_ss_numpy + diff_ss_numpy correct = verify_results(sum_ss_numpy, "Numpy", "active + inactive") if not correct: print("\tThe error lies in the value:") print(f"\t total Qiskit = {tot_ssz_qiskit}") print(f"\t active Qiskit = {act_ssz_qiskit}") print(f"\t inactive Qiskit = {diff_ssz_qiskit}") print(f"\t total Numpy = {tot_ssz_numpy}") print(f"\t active Numpy = {act_ssz_numpy}") print(f"\t inactive Numpy = {diff_ssz_numpy}") print(f"\tThe number of active electrons is: {nelec_act}") nelec_inact = (nelec_tot[0] - nelec_act[0], nelec_tot[1] - nelec_act[1]) print(f"\tThe number of inactive electrons is: {nelec_inact}") ```

One final problem remains, which is for the case where both, the inactive and active sector include unpaired electrons. Here, the numpy.einsum code does NOT produce the correct expectation value of $S_z$. It is possible to get the correct inactive value from the difference of the total and active terms, but computing the inactive value directly I have still not figured out yet. However, this is not really a problem because such a choice of active space is not really relevant.

Output of the above snippet

``` converged SCF energy = -147.620355111774 = 2.0033214 2S+1 = 3.0022135 ===== References ====== The total computed by Qiskit is correct! The total computed by Numpy is correct! --- 0 Inactive Orbitals --- The active + inactive computed by Qiskit is correct! The active + inactive computed by Numpy is correct! --- 1 Inactive Orbitals --- The active + inactive computed by Qiskit is correct! The active + inactive computed by Numpy is correct! --- 2 Inactive Orbitals --- The active + inactive computed by Qiskit is correct! The active + inactive computed by Numpy is correct! --- 3 Inactive Orbitals --- The active + inactive computed by Qiskit is correct! The active + inactive computed by Numpy is correct! --- 4 Inactive Orbitals --- The active + inactive computed by Qiskit is correct! The active + inactive computed by Numpy is correct! --- 5 Inactive Orbitals --- The active + inactive computed by Qiskit is correct! The active + inactive computed by Numpy is correct! --- 6 Inactive Orbitals --- The active + inactive computed by Qiskit is correct! The active + inactive computed by Numpy is correct! --- 7 Inactive Orbitals --- The active + inactive computed by Qiskit is correct! The active + inactive computed by Numpy is correct! --- 8 Inactive Orbitals --- The active + inactive computed by Qiskit is correct! ERROR: The active + inactive computed by Numpy is WRONG! = 1.50332 Δ = 0.50000 The error lies in the value: total Qiskit = 1.0 active Qiskit = 0.25 inactive Qiskit = 0.75 total Numpy = 1.0 active Numpy = 0.25 inactive Numpy = 0.25 The number of active electrons is: (1, 0) The number of inactive electrons is: (8, 7) --- 9 Inactive Orbitals --- The active + inactive computed by Qiskit is correct! The active + inactive computed by Numpy is correct! ```

Finally, this means that the implementation of #1275 was largely correct. This is also good news, because integrating the overlap-matrix into the FermionicOp would not really work as I had originally suggested. A PR to make some final improvements to fix this issue will follow shortly.

qiskit-community / qiskit-nature

The orbital overlap is still not handled correctly by the `AngularMomentum` operator #1291

Environment

What is happening?