Handle orbital overlap in spin operators (backport #1275)

.pylintdict

@@ -365,6 +365,7 @@ møller
 namelist
 nan
 nannichini
+nao
 natom
 nbasis
 nd
@@ -410,6 +411,7 @@ orthonormal
 otimes
 outpath
 overline
+ovlp
 pac
 param
 parameterized

qiskit_nature/second_q/drivers/electronic_structure_driver.py

@@ -53,6 +53,7 @@ class _QCSchemaData:
     eri_mo_bb: np.ndarray | None = None
     e_nuc: float | None = None
     e_ref: float | None = None
+    overlap: np.ndarray | None = None
     mo_coeff: np.ndarray | None = None
     mo_coeff_b: np.ndarray | None = None
     mo_energy: np.ndarray | None = None
@@ -242,6 +243,9 @@ def format_np_array(arr):
             return arr.ravel().tolist()
 
         wavefunction = QCWavefunction(basis=data.basis)
+        if data.overlap is not None:
+            wavefunction.overlap = "scf_overlap"
+            wavefunction.scf_overlap = format_np_array(data.overlap)
         if data.mo_coeff is not None:
             wavefunction.orbitals_a = "scf_orbitals_a"
             wavefunction.scf_orbitals_a = format_np_array(data.mo_coeff)

qiskit_nature/second_q/drivers/gaussiand/gaussiandriver.py

@@ -297,6 +297,7 @@ def _qcschema_from_matrix_file(mel: MatEl, *, include_dipole: bool = True) -> QC
         einsum_func, _ = get_einsum()
         data = _QCSchemaData()
 
+        data.overlap = GaussianDriver._get_matrix(mel, "OVERLAP")
         data.mo_coeff = GaussianDriver._get_matrix(mel, "ALPHA MO COEFFICIENTS")
         data.mo_coeff_b = GaussianDriver._get_matrix(mel, "BETA MO COEFFICIENTS")
         if np.array_equal(data.mo_coeff, data.mo_coeff_b):

qiskit_nature/second_q/drivers/psi4d/_template.txt

@@ -36,6 +36,7 @@ if _has_B:
 
 data.e_nuc = _q_mol.nuclear_repulsion_energy()
 data.e_ref = _q_hf_energy
+data.overlap = np.asarray(_q_mints.ao_overlap())
 data.mo_coeff = np.asarray(_q_hf_wavefn.Ca())
 data.mo_coeff_b = np.asarray(_q_hf_wavefn.Cb()) if _has_B else None
 data.mo_energy = np.asarray(_q_hf_wavefn.epsilon_a())

qiskit_nature/second_q/drivers/pyscfd/pyscfdriver.py

@@ -512,6 +512,7 @@ def to_qcschema(self, *, include_dipole: bool = True) -> QCSchema:
         einsum_func, _ = get_einsum()
         data = _QCSchemaData()
 
+        data.overlap = self._calc.get_ovlp()
         data.mo_coeff, data.mo_coeff_b = self._expand_mo_object(
             self._calc.mo_coeff, array_dimension=3
         )

qiskit_nature/second_q/formats/__init__.py

@@ -36,12 +36,17 @@
    watson_to_problem
    qcschema_to_problem
    get_ao_to_mo_from_qcschema
+   get_overlap_ab_from_qcschema
 
 """
 
 from .molecule_info import MoleculeInfo
 from .fcidump_translator import fcidump_to_problem
-from .qcschema_translator import qcschema_to_problem, get_ao_to_mo_from_qcschema
+from .qcschema_translator import (
+    qcschema_to_problem,
+    get_ao_to_mo_from_qcschema,
+    get_overlap_ab_from_qcschema,
+)
 from .watson_translator import watson_to_problem
 
 __all__ = [
@@ -50,4 +55,5 @@
     "watson_to_problem",
     "qcschema_to_problem",
     "get_ao_to_mo_from_qcschema",
+    "get_overlap_ab_from_qcschema",
 ]

qiskit_nature/second_q/formats/qcschema/qc_wavefunction.py

@@ -44,6 +44,8 @@ class QCWavefunction(_QCBase):
     Qiskit Nature.
     """
 
+    overlap: str | None = None
+    """The name of the overlap matrix in the AO basis."""
     orbitals_a: str | None = None
     """The name of the alpha-spin orbitals in the AO basis."""
     orbitals_b: str | None = None
@@ -101,6 +103,8 @@ class QCWavefunction(_QCBase):
     dipole_mo_z_b: str | None = None
     """The name of beta-spin z-axis dipole moment integrals in the MO basis."""
 
+    scf_overlap: Sequence[float] | None = None
+    """The SCF overlap matrix in the AO basis."""
     scf_orbitals_a: Sequence[float] | None = None
     """The SCF alpha-spin orbitals in the AO basis."""
     scf_orbitals_b: Sequence[float] | None = None

qiskit_nature/second_q/formats/qcschema_translator.py

@@ -14,6 +14,7 @@
 
 from __future__ import annotations
 
+import logging
 from typing import cast
 
 import numpy as np
@@ -38,6 +39,8 @@
 from .molecule_info import MoleculeInfo
 from .qcschema import QCSchema
 
+LOGGER = logging.getLogger(__name__)
+
 
 def qcschema_to_problem(
     qcschema: QCSchema,
@@ -69,10 +72,15 @@ def qcschema_to_problem(
     dipole_x: ElectronicIntegrals | None = None
     dipole_y: ElectronicIntegrals | None = None
     dipole_z: ElectronicIntegrals | None = None
+    overlap: np.ndarray | None = None
 
     if basis == ElectronicBasis.AO:
         hamiltonian = _get_ao_hamiltonian(qcschema)
 
+        if qcschema.wavefunction.scf_overlap is not None:
+            nao = hamiltonian.register_length
+            overlap = _reshape_2(qcschema.wavefunction.scf_overlap, nao, nao)
+
         if include_dipole:
             dipole_x = _get_ao_dipole(qcschema, "x")
             dipole_y = _get_ao_dipole(qcschema, "y")
@@ -81,6 +89,21 @@ def qcschema_to_problem(
     elif basis == ElectronicBasis.MO:
         hamiltonian = _get_mo_hamiltonian(qcschema)
 
+        if qcschema.wavefunction.scf_overlap is not None:
+            try:
+                overlap = get_overlap_ab_from_qcschema(qcschema)
+            except AttributeError:
+                if not hamiltonian.electronic_integrals.beta.is_empty() and not np.allclose(
+                    hamiltonian.electronic_integrals.alpha["+-"],
+                    hamiltonian.electronic_integrals.beta["+-"],
+                ):
+                    LOGGER.warning(
+                        "Your MO coefficients appear to differ for the alpha- and beta-spin "
+                        "orbitals but you did not provide any orbital overlap data. This may lead "
+                        "to faulty expectation values of observables which mix alpha- and beta-spin"
+                        " components (for example the AngularMomentum)."
+                    )
+
         if include_dipole:
             dipole_x = _get_mo_dipole(qcschema, "x")
             dipole_y = _get_mo_dipole(qcschema, "y")
@@ -110,7 +133,7 @@ def qcschema_to_problem(
     problem.num_particles = num_particles
     problem.num_spatial_orbitals = norb
     problem.reference_energy = qcschema.properties.return_energy
-    problem.properties.angular_momentum = AngularMomentum(norb)
+    problem.properties.angular_momentum = AngularMomentum(norb, overlap)
     problem.properties.magnetization = Magnetization(norb)
     problem.properties.particle_number = ParticleNumber(norb)
 
@@ -256,3 +279,36 @@ def get_ao_to_mo_from_qcschema(qcschema: QCSchema) -> BasisTransformer:
     )
 
     return transformer
+
+
+def get_overlap_ab_from_qcschema(qcschema: QCSchema) -> np.ndarray | None:
+    """Builds the alpha-beta spin orbital overlap matrix from a :class:`.QCSchema` instance.
+
+    Args:
+        qcschema: the :class:`.QCSchema` object from which to build the problem.
+
+    Raises:
+        AttributeError: when the overlap or MO coefficients are missing.
+
+    Returns:
+        The overlap matrix as a 2-dimensional numpy array or `None` if the overlap is equal to the
+        identity matrix.
+    """
+    if qcschema.wavefunction.scf_orbitals_a is None or qcschema.wavefunction.scf_overlap is None:
+        raise AttributeError
+
+    nmo = qcschema.properties.calcinfo_nmo
+    nao = len(qcschema.wavefunction.scf_orbitals_a) // nmo
+    coeff_a = _reshape_2(qcschema.wavefunction.scf_orbitals_a, nao, nmo)
+    coeff_b = None
+    if qcschema.wavefunction.scf_orbitals_b is not None:
+        coeff_b = _reshape_2(qcschema.wavefunction.scf_orbitals_b, nao, nmo)
+
+    if coeff_b is None:
+        # when coeff_b is None, that means it is equal to coeff_a, which in turn means that the
+        # overlap is equal to the identity
+        return None
+
+    overlap = _reshape_2(qcschema.wavefunction.scf_overlap, nao, nao)
+
+    return coeff_a.T @ overlap @ coeff_b

qiskit_nature/second_q/properties/angular_momentum.py

@@ -16,21 +16,29 @@
 
 from typing import Mapping
 
+import numpy as np
+
 import qiskit_nature  # pylint: disable=unused-import
 from qiskit_nature.second_q.operators import FermionicOp
 
 from .s_operators import s_minus_operator, s_plus_operator, s_z_operator
 
 
 class AngularMomentum:
-    """The AngularMomentum property.
+    r"""The AngularMomentum property.
 
     The operator constructed by this property is the $S^2$ operator which is computed as:
 
     .. math::
 
        S^2 = S^- S^+ + S^z (S^z + 1)
 
+    .. warning::
+
+       If you are working with a non-orthogonal basis, you _must_ provide the ``overlap`` attribute
+       in order to obtain the correct expectation value of this observable. Refer to the more
+       extensive documentation of the :mod:`.s_operators` module for more details.
+
     See also:
         - the $S^z$ operator: :func:`.s_z_operator`
         - the $S^+$ operator: :func:`.s_plus_operator`
@@ -41,14 +49,19 @@ class AngularMomentum:
 
     Attributes:
         num_spatial_orbitals (int): the number of spatial orbitals.
+        overlap (np.ndarray | None): the overlap-matrix between the $\alpha$- and $\beta$-spin
+            orbitals. When this is `None`, the overlap-matrix is assumed to be identity.
     """
 
-    def __init__(self, num_spatial_orbitals: int) -> None:
-        """
+    def __init__(self, num_spatial_orbitals: int, overlap: np.ndarray | None = None) -> None:
+        r"""
         Args:
             num_spatial_orbitals: the number of spatial orbitals in the system.
+            overlap: the overlap-matrix between the $\alpha$- and $\beta$-spin orbitals. When this
+                is `None`, the overlap-matrix is assumed to be identity.
         """
         self.num_spatial_orbitals = num_spatial_orbitals
+        self.overlap = overlap
 
     def second_q_ops(self) -> Mapping[str, FermionicOp]:
         """Returns the second quantized angular momentum operator.
@@ -57,8 +70,10 @@ def second_q_ops(self) -> Mapping[str, FermionicOp]:
             A mapping of strings to `FermionicOp` objects.
         """
         s_z = s_z_operator(self.num_spatial_orbitals)
-        s_p = s_plus_operator(self.num_spatial_orbitals)
-        s_m = s_minus_operator(self.num_spatial_orbitals)
+        overlap_ab = self.overlap
+        s_p = s_plus_operator(self.num_spatial_orbitals, overlap=overlap_ab)
+        overlap_ba = overlap_ab.T if overlap_ab is not None else None
+        s_m = s_minus_operator(self.num_spatial_orbitals, overlap=overlap_ba)
 
         op = s_m @ s_p + s_z @ (s_z + FermionicOp.one())