diff --git a/.pylintdict b/.pylintdict
index d5213ea68..d1ef95662 100644
--- a/.pylintdict
+++ b/.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
diff --git a/qiskit_nature/second_q/drivers/electronic_structure_driver.py b/qiskit_nature/second_q/drivers/electronic_structure_driver.py
index 287fe374f..b464ed277 100644
--- a/qiskit_nature/second_q/drivers/electronic_structure_driver.py
+++ b/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)
diff --git a/qiskit_nature/second_q/drivers/gaussiand/gaussiandriver.py b/qiskit_nature/second_q/drivers/gaussiand/gaussiandriver.py
index 15e3a2e22..e60e48c66 100644
--- a/qiskit_nature/second_q/drivers/gaussiand/gaussiandriver.py
+++ b/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):
diff --git a/qiskit_nature/second_q/drivers/psi4d/_template.txt b/qiskit_nature/second_q/drivers/psi4d/_template.txt
index cc7b06710..166f09817 100644
--- a/qiskit_nature/second_q/drivers/psi4d/_template.txt
+++ 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())
diff --git a/qiskit_nature/second_q/drivers/pyscfd/pyscfdriver.py b/qiskit_nature/second_q/drivers/pyscfd/pyscfdriver.py
index 4764c48ae..a38b89496 100644
--- a/qiskit_nature/second_q/drivers/pyscfd/pyscfdriver.py
+++ b/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
         )
diff --git a/qiskit_nature/second_q/formats/__init__.py b/qiskit_nature/second_q/formats/__init__.py
index 4fe38bb58..4a639f110 100644
--- a/qiskit_nature/second_q/formats/__init__.py
+++ b/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",
 ]
diff --git a/qiskit_nature/second_q/formats/qcschema/qc_wavefunction.py b/qiskit_nature/second_q/formats/qcschema/qc_wavefunction.py
index a700dd89c..912d79bb6 100644
--- a/qiskit_nature/second_q/formats/qcschema/qc_wavefunction.py
+++ b/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
diff --git a/qiskit_nature/second_q/formats/qcschema_translator.py b/qiskit_nature/second_q/formats/qcschema_translator.py
index e5312c3ef..45e6c4168 100644
--- a/qiskit_nature/second_q/formats/qcschema_translator.py
+++ b/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
diff --git a/qiskit_nature/second_q/properties/angular_momentum.py b/qiskit_nature/second_q/properties/angular_momentum.py
index 877042f3e..8738bfd66 100644
--- a/qiskit_nature/second_q/properties/angular_momentum.py
+++ b/qiskit_nature/second_q/properties/angular_momentum.py
@@ -16,6 +16,8 @@
 
 from typing import Mapping
 
+import numpy as np
+
 import qiskit_nature  # pylint: disable=unused-import
 from qiskit_nature.second_q.operators import FermionicOp
 
@@ -23,7 +25,7 @@
 
 
 class AngularMomentum:
-    """The AngularMomentum property.
+    r"""The AngularMomentum property.
 
     The operator constructed by this property is the $S^2$ operator which is computed as:
 
@@ -31,6 +33,12 @@ class AngularMomentum:
 
        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())
 
diff --git a/qiskit_nature/second_q/properties/s_operators.py b/qiskit_nature/second_q/properties/s_operators.py
index 47a3e5455..c9d69c6e9 100644
--- a/qiskit_nature/second_q/properties/s_operators.py
+++ b/qiskit_nature/second_q/properties/s_operators.py
@@ -10,66 +10,141 @@
 # copyright notice, and modified files need to carry a notice indicating
 # that they have been altered from the originals.
 
-"""Generator functions for various spin-related operators."""
+r"""Generator functions for various spin-related operators.
+
+When dealing with non-orthonormal orbitals, you need to make sure that you include the `overlap`
+matrices when using the methods below. This ensures that the operators can resolve any spin
+contamination that may be present in your orbitals.
+
+The overlap matrices that you provide have to be computed in the same basis in which the spin
+operator is encoded. If you are working in the molecular orbital (MO) basis, the overlap can be
+easily constructed starting from the atomic orbital (AO) overlap matrix, which can be obtained from
+any standard quantum chemistry program (for example from the `get_oplp()` method in PySCF). This AO
+overlap matrix can be transformed to the MO basis using the AO-to-MO transformation matrix, $C$,
+according to the following equation:
+
+.. math::
+
+   s^{MO} = C^T s^{AO} C.
+
+For restricted spin orbitals (i.e. :math:`C_\alpha == C_\beta`), the equation above simplifies to
+the identity matrix (because the MOs will be orthonormal), in which case you can omit the `overlap`
+arguments below). Otherwise, you must include the correct overlap. For example, the overlap-matrix
+between the $`\alpha`$- and $`\beta`$-spin orbitals is:
+
+.. math::
+
+   s^{\alpha,\beta} = C_\alpha^T s^{AO} C_\beta.
+"""
 
 from __future__ import annotations
 
+import numpy as np
+
 from qiskit_nature.second_q.operators import FermionicOp
 
 
-def s_plus_operator(num_spatial_orbitals: int) -> FermionicOp:
+def s_plus_operator(num_spatial_orbitals: int, overlap: np.ndarray | None = None) -> FermionicOp:
     r"""Constructs the $S^+$ operator.
 
     The $S^+$ operator is defined as:
 
     .. math::
 
-       S^+ = \sum_{i=1}^{n} \hat{a}_{i}^{\dagger} \hat{a}_{i+n}
+       S^+ = \sum_{i,j} s_{ij}^{\alpha,\beta} \hat{a}_{i}^{\dagger} \hat{a}_{j},
+
+    where $s$ denotes the overlap-matrix between the $\alpha$- and $\beta$-spin orbitals.
 
-    Here, $n$ denotes the index of the *spatial* orbital. Since Qiskit Nature employs the blocked
+    Note that for orthonormal orbitals this overlap-matrix will become the identity matrix,
+    simplifying the operator above to become:
+
+    .. math::
+
+       S^+ = \sum_{i=1}^{n} \hat{a}_{i}^{\dagger} \hat{a}_{i+n},
+
+    where, $n$ denotes the index of the *spatial* orbital. Since Qiskit Nature employs the blocked
     spin-ordering convention, the creation operator above is applied to the :math:`\alpha`-spin
     orbital and the annihilation operator is applied to the corresponding :math:`\beta`-spin
     orbital.
 
     Args:
         num_spatial_orbitals: the size of the operator which to generate.
+        overlap: the overlap-matrix between the $\alpha$- and $\beta$-spin orbitals. When this is
+            `None`, the overlap-matrix is assumed to be identity, resulting in the second definition
+            above.
 
     Returns:
         The $S^+$ operator of the requested size.
     """
-    op = FermionicOp(
-        {f"+_{orb} -_{orb + num_spatial_orbitals}": 1.0 for orb in range(num_spatial_orbitals)}
-    )
-    return op
-
-
-def s_minus_operator(num_spatial_orbitals: int) -> FermionicOp:
+    if overlap is None:
+        op = FermionicOp(
+            {f"+_{orb} -_{orb + num_spatial_orbitals}": 1.0 for orb in range(num_spatial_orbitals)}
+        )
+    else:
+        op = FermionicOp(
+            {
+                f"+_{idx[0]} -_{idx[1] + num_spatial_orbitals}": overlap[idx]
+                for idx in np.ndindex(*overlap.shape)
+            }
+        )
+    return op.simplify()
+
+
+def s_minus_operator(num_spatial_orbitals: int, overlap: np.ndarray | None = None) -> FermionicOp:
     r"""Constructs the $S^-$ operator.
 
     The $S^-$ operator is defined as:
 
     .. math::
 
+       S^- = \sum_{i,j} s_{ij}^{\beta,\alpha} \hat{a}_{i}^{\dagger} \hat{a}_{j},
+
+    where $s$ denotes the overlap-matrix between the $\beta$- and $\alpha$-spin orbitals.
+
+    .. note::
+
+       The `overlap` input to this method is related to the input of the other methods
+       (:meth:`s_plus_operator`, :meth:`s_x_operator`, and :meth:`s_y_operator`) by its transpose,
+       because the following relation holds:
+
+       .. math::
+
+          s_{ij}^{\beta,\alpha} = \left(s_{ij}^{\alpha,\beta}\right)^T.
+
+    Note that for orthonormal orbitals this overlap-matrix will become the identity matrix,
+    simplifying the operator above to become:
+
        S^- = \sum_{i=1}^{n} \hat{a}_{i+n}^{\dagger} \hat{a}_{i}
 
-    Here, $n$ denotes the index of the *spatial* orbital. Since Qiskit Nature employs the blocked
+    where, $n$ denotes the index of the *spatial* orbital. Since Qiskit Nature employs the blocked
     spin-ordering convention, the creation operator above is applied to the :math:`\beta`-spin
     orbital and the annihilation operator is applied to the corresponding :math:`\alpha`-spin
     orbital.
 
     Args:
         num_spatial_orbitals: the size of the operator which to generate.
+        overlap: the overlap-matrix between the $\beta$- and $\alpha$-spin orbitals. When this is
+            `None`, the overlap-matrix is assumed to be identity, resulting in the second definition
+            above.
 
     Returns:
         The $S^-$ operator of the requested size.
     """
-    op = FermionicOp(
-        {f"+_{orb + num_spatial_orbitals} -_{orb}": 1.0 for orb in range(num_spatial_orbitals)}
-    )
-    return op
-
-
-def s_x_operator(num_spatial_orbitals: int) -> FermionicOp:
+    if overlap is None:
+        op = FermionicOp(
+            {f"+_{orb + num_spatial_orbitals} -_{orb}": 1.0 for orb in range(num_spatial_orbitals)}
+        )
+    else:
+        op = FermionicOp(
+            {
+                f"+_{idx[0] + num_spatial_orbitals} -_{idx[1]}": overlap[idx]
+                for idx in np.ndindex(*overlap.shape)
+            }
+        )
+    return op.simplify()
+
+
+def s_x_operator(num_spatial_orbitals: int, overlap: np.ndarray | None = None) -> FermionicOp:
     r"""Constructs the $S^x$ operator.
 
     The $S^x$ operator is defined as:
@@ -80,20 +155,28 @@ def s_x_operator(num_spatial_orbitals: int) -> FermionicOp:
 
     Args:
         num_spatial_orbitals: the size of the operator which to generate.
+        overlap: the overlap-matrix between the $\alpha$- and $\beta$-spin orbitals. When this is
+            `None`, the overlap-matrix is assumed to be identity.
 
     Returns:
         The $S^x$ operator of the requested size.
     """
-    op = FermionicOp(
-        {
-            f"+_{orb} -_{(orb + num_spatial_orbitals) % (2*num_spatial_orbitals)}": 0.5
-            for orb in range(2 * num_spatial_orbitals)
-        }
-    )
+    if overlap is None:
+        op = FermionicOp(
+            {
+                f"+_{orb} -_{(orb + num_spatial_orbitals) % (2*num_spatial_orbitals)}": 0.5
+                for orb in range(2 * num_spatial_orbitals)
+            }
+        )
+    else:
+        op = 0.5 * (
+            s_plus_operator(num_spatial_orbitals, overlap)
+            + s_minus_operator(num_spatial_orbitals, overlap.T)
+        )
     return op
 
 
-def s_y_operator(num_spatial_orbitals: int) -> FermionicOp:
+def s_y_operator(num_spatial_orbitals: int, overlap: np.ndarray | None = None) -> FermionicOp:
     r"""Constructs the $S^y$ operator.
 
     The $S^y$ operator is defined as:
@@ -104,17 +187,25 @@ def s_y_operator(num_spatial_orbitals: int) -> FermionicOp:
 
     Args:
         num_spatial_orbitals: the size of the operator which to generate.
+        overlap: the overlap-matrix between the $\alpha$- and $\beta$-spin orbitals. When this is
+            `None`, the overlap-matrix is assumed to be identity.
 
     Returns:
         The $S^y$ operator of the requested size.
     """
-    op = FermionicOp(
-        {
-            f"+_{orb} -_{(orb + num_spatial_orbitals) % (2*num_spatial_orbitals)}": 0.5j
-            * (-1.0) ** (orb < num_spatial_orbitals)
-            for orb in range(2 * num_spatial_orbitals)
-        }
-    )
+    if overlap is None:
+        op = FermionicOp(
+            {
+                f"+_{orb} -_{(orb + num_spatial_orbitals) % (2*num_spatial_orbitals)}": 0.5j
+                * (-1.0) ** (orb < num_spatial_orbitals)
+                for orb in range(2 * num_spatial_orbitals)
+            }
+        )
+    else:
+        op = -0.5j * (
+            s_plus_operator(num_spatial_orbitals, overlap)
+            - s_minus_operator(num_spatial_orbitals, overlap.T)
+        )
     return op
 
 
@@ -127,13 +218,19 @@ def s_z_operator(num_spatial_orbitals: int) -> FermionicOp:
 
        S^z = \frac{1}{2} \sum_{i=1}^{n} \left(
         \hat{a}_{i}^{\dagger}\hat{a}_{i} - \hat{a}_{i+n}^{\dagger}\hat{a}_{i+n}
-       \right)
+       \right),
 
-    Here, $n$ denotes the index of the *spatial* orbital. Since Qiskit Nature employs the blocked
+    where, $n$ denotes the index of the *spatial* orbital. Since Qiskit Nature employs the blocked
     spin-ordering convention, this means that the above corresponds to evaluating the number
     operator (:math:`\hat{a}^{\dagger}\hat{a}`) once on the :math:`\alpha`-spin orbital and once on
     the :math:`\beta`-spin orbital and taking their difference.
 
+    .. note::
+
+       Contrary to the other methods in this module, this one does not require the inclusion of an
+       overlap-matrix for non-orthonormal orbitals, because it does not mix the $\alpha$- and
+       $\beta$-spin contributions.
+
     Args:
         num_spatial_orbitals: the size of the operator which to generate.
 
diff --git a/qiskit_nature/second_q/transformers/active_space_transformer.py b/qiskit_nature/second_q/transformers/active_space_transformer.py
index a9cde774a..1b9f659ca 100644
--- a/qiskit_nature/second_q/transformers/active_space_transformer.py
+++ b/qiskit_nature/second_q/transformers/active_space_transformer.py
@@ -22,7 +22,7 @@
 
 from qiskit_nature import QiskitNatureError
 from qiskit_nature.second_q.hamiltonians import ElectronicEnergy, Hamiltonian
-from qiskit_nature.second_q.operators import ElectronicIntegrals
+from qiskit_nature.second_q.operators import ElectronicIntegrals, Tensor
 from qiskit_nature.second_q.problems import BaseProblem, ElectronicBasis, ElectronicStructureProblem
 from qiskit_nature.second_q.properties import (
     AngularMomentum,
@@ -279,7 +279,26 @@ def _transform_electronic_structure_problem(
                 new_problem.properties.electronic_density = ElectronicDensity(
                     transformed.alpha, transformed.beta, transformed.beta_alpha
                 )
-            elif isinstance(prop, (AngularMomentum, Magnetization, ParticleNumber)):
+            elif isinstance(prop, AngularMomentum):
+                if prop.overlap is None:
+                    # only the size needs to be changed
+                    new_problem.properties.add(prop.__class__(new_problem.num_spatial_orbitals))
+                    continue
+
+                if isinstance(self._active_basis.coefficients, ElectronicIntegrals):
+                    coeff_alpha = self._active_basis.coefficients.alpha["+-"]
+                    coeff_beta: Tensor
+                    if self._active_basis.coefficients.beta.is_empty():
+                        coeff_beta = coeff_alpha
+                    else:
+                        coeff_beta = self._active_basis.coefficients.beta["+-"]
+
+                    new_problem.properties.angular_momentum = AngularMomentum(
+                        new_problem.num_spatial_orbitals,
+                        coeff_alpha.transpose() @ prop.overlap @ coeff_beta,
+                    )
+
+            elif isinstance(prop, (Magnetization, ParticleNumber)):
                 new_problem.properties.add(prop.__class__(new_problem.num_spatial_orbitals))
             else:
                 LOGGER.warning("Encountered an unsupported property of type '%s'.", type(prop))
diff --git a/qiskit_nature/second_q/transformers/basis_transformer.py b/qiskit_nature/second_q/transformers/basis_transformer.py
index 2ab335be0..26c1df469 100644
--- a/qiskit_nature/second_q/transformers/basis_transformer.py
+++ b/qiskit_nature/second_q/transformers/basis_transformer.py
@@ -21,7 +21,7 @@
 
 from qiskit_nature.exceptions import QiskitNatureError
 from qiskit_nature.second_q.hamiltonians import ElectronicEnergy, Hamiltonian
-from qiskit_nature.second_q.operators import ElectronicIntegrals, PolynomialTensor
+from qiskit_nature.second_q.operators import ElectronicIntegrals, PolynomialTensor, Tensor
 from qiskit_nature.second_q.problems import BaseProblem, ElectronicBasis, ElectronicStructureProblem
 from qiskit_nature.second_q.properties import (
     AngularMomentum,
@@ -152,7 +152,26 @@ def _transform_electronic_structure_problem(
                 new_problem.properties.electronic_density = self.transform_electronic_integrals(
                     prop
                 )
-            elif isinstance(prop, (AngularMomentum, Magnetization, ParticleNumber)):
+            elif isinstance(prop, AngularMomentum):
+                if prop.overlap is None:
+                    # nothing needs to be changed
+                    new_problem.properties.add(prop)
+                    continue
+
+                if isinstance(self.coefficients, ElectronicIntegrals):
+                    coeff_alpha = self.coefficients.alpha["+-"]
+                    coeff_beta: Tensor
+                    if self.coefficients.beta.is_empty():
+                        coeff_beta = coeff_alpha
+                    else:
+                        coeff_beta = self.coefficients.beta["+-"]
+
+                    new_problem.properties.angular_momentum = AngularMomentum(
+                        prop.num_spatial_orbitals,
+                        coeff_alpha.transpose() @ prop.overlap @ coeff_beta,
+                    )
+
+            elif isinstance(prop, (Magnetization, ParticleNumber)):
                 new_problem.properties.add(prop)
             else:
                 LOGGER.warning("Encountered an unsupported property of type '%s'.", type(prop))
diff --git a/releasenotes/notes/fix-angular-momentum-73eca0a5afb70679.yaml b/releasenotes/notes/fix-angular-momentum-73eca0a5afb70679.yaml
new file mode 100644
index 000000000..30d393b8d
--- /dev/null
+++ b/releasenotes/notes/fix-angular-momentum-73eca0a5afb70679.yaml
@@ -0,0 +1,22 @@
+---
+features:
+  - |
+    Added the :func:`.get_overlap_ab_from_qcschema` function which extracts the
+    alpha-beta spin orbital overlap matrix from a :class:`.QCSchema` instance.
+fixes:
+  - |
+    Fixes the following operators when dealing with non-orthonormal orbitals
+    (for example when using unrestricted spin orbitals):
+    - :class:`.AnglarMomentum`
+    - :func:`.s_plus_operator`
+    - :func:`.s_minus_operator`
+    - :func:`.s_x_operator`
+    - :func:`.s_y_operator`
+
+    To make the fix take effect, the new additional `overlap` argument needs to
+    be provided to all of these operators.
+
+    Prior to this fix, none of the operators above were able to resolve any spin
+    contamination and would yield misleadingly "clean" expectation values. See
+    `this issue <https://github.com/qiskit-community/qiskit-nature/issues/1273>`_
+    for a more complete discussion.
diff --git a/test/second_q/drivers/pyscfd/test_driver_methods_pyscf.py b/test/second_q/drivers/pyscfd/test_driver_methods_pyscf.py
index adc3d451c..7f4a521e5 100644
--- a/test/second_q/drivers/pyscfd/test_driver_methods_pyscf.py
+++ b/test/second_q/drivers/pyscfd/test_driver_methods_pyscf.py
@@ -15,7 +15,9 @@
 import unittest
 
 from test.second_q.drivers.test_driver_methods_gsc import TestDriverMethods
+from qiskit.quantum_info import Statevector
 from qiskit.test import slow_test
+from qiskit_nature.second_q.circuit.library import HartreeFock
 from qiskit_nature.units import DistanceUnit
 from qiskit_nature.second_q.drivers import PySCFDriver, MethodType
 from qiskit_nature.second_q.mappers import BravyiKitaevMapper, ParityMapper
@@ -235,6 +237,34 @@ def test_oh_uhf_bk(self):
         result = self._run_driver(driver, mapper=BravyiKitaevMapper())
         self._assert_energy_and_dipole(result, "oh")
 
+    @slow_test
+    def test_be2_spin_contaminaton(self):
+        """Test the spin contamination is preserved.
+
+        This is a regression test against
+        https://github.com/qiskit-community/qiskit-nature/issues/1273.
+        """
+        driver = PySCFDriver(
+            atom="Be 0.0 0.0 0.66242; Be 0.0 0.0 -0.66242",
+            unit=DistanceUnit.ANGSTROM,
+            charge=0,
+            spin=2,
+            basis="sto-3g",
+            method=MethodType.UHF,
+        )
+        problem = driver.run()
+        _, ops = problem.second_q_ops()
+
+        mapper = ParityMapper(problem.num_particles)
+        qop = mapper.map(ops["AngularMomentum"])
+
+        hf_state = Statevector(
+            HartreeFock(problem.num_spatial_orbitals, problem.num_particles, mapper)
+        )
+
+        result = hf_state.expectation_value(qop)
+        self.assertAlmostEqual(result, 2.1286097507302286)
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/test/second_q/properties/test_angular_momentum.py b/test/second_q/properties/test_angular_momentum.py
index 1629f6863..32dc7576c 100644
--- a/test/second_q/properties/test_angular_momentum.py
+++ b/test/second_q/properties/test_angular_momentum.py
@@ -16,6 +16,12 @@
 import json
 from test.second_q.properties.property_test import PropertyTest
 
+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 ParityMapper
 from qiskit_nature.second_q.properties import AngularMomentum
 from qiskit_nature.second_q.operators import FermionicOp
 
@@ -42,6 +48,62 @@ def test_second_q_ops(self):
 
         self.assertEqual(op.normal_order(), expected_op.normal_order())
 
+    def test_with_overlap(self):
+        """Test that the overlap is taken into account.
+
+        The idea of this test is simple. I found the MO coefficients of H2 @ 0.735 Angstrom using a
+        RHF calculation (``mo_coeff`` in the code below). Being restricted-spin coefficients for a
+        singlet-spin solution, these alone would give a total angular momentum of 0. However, we can
+        take the matrix and rotate it by some angle to form a spin contaminated pair of alpha- and
+        beta-spin coefficients. These should result in a non-zero angular momentum.
+
+        Below is the code to reproduce the RHF MO coefficients and overlap matrix with PySCF:
+
+        .. code:: python
+
+           from pyscf import gto, scf
+
+           mol = gto.M(atom="H 0.0 0.0 0.0; H 0.0 0.0 0.735;", basis="sto-3g")
+
+           hf = scf.RHF(mol)
+           hf.run()
+
+           norb = mol.nao
+           nelec = mol.nelec
+           mo_coeff = hf.mo_coeff
+           ovlp = hf.get_ovlp()
+
+        """
+        norb = 2
+        nelec = (1, 1)
+
+        mo_coeff = np.asarray([[0.54830202, 1.21832731], [0.54830202, -1.21832731]])
+        ovlp = np.asarray([[1.0, 0.66314574], [0.66314574, 1.0]])
+
+        # first, we ensure that our restricted-spin MO coefficients and overlap match
+        self.assertTrue(np.allclose(mo_coeff.T @ ovlp @ mo_coeff, np.eye(norb)))
+
+        theta = 33
+        rot = np.asarray(
+            [
+                [np.cos(np.deg2rad(theta)), -np.sin(np.deg2rad(theta))],
+                [np.sin(np.deg2rad(theta)), np.cos(np.deg2rad(theta))],
+            ],
+        )
+        mo_coeff_rot = mo_coeff @ rot
+
+        ovlpab = mo_coeff.T @ ovlp @ mo_coeff_rot
+
+        ang_mom = AngularMomentum(norb, ovlpab).second_q_ops()
+
+        mapper = ParityMapper(nelec)
+        qubit_op = mapper.map(ang_mom)
+
+        hf_state = HartreeFock(norb, nelec, mapper)
+
+        result = estimate_observables(Estimator(), hf_state, qubit_op)
+        self.assertAlmostEqual(result["AngularMomentum"][0], 0.29663167846210015)
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/test/second_q/properties/test_s_operators.py b/test/second_q/properties/test_s_operators.py
index 8b5514c0e..30936a01e 100644
--- a/test/second_q/properties/test_s_operators.py
+++ b/test/second_q/properties/test_s_operators.py
@@ -15,6 +15,12 @@
 import unittest
 from test import QiskitNatureTestCase
 
+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 ParityMapper
 from qiskit_nature.second_q.operators import FermionicOp
 from qiskit_nature.second_q.operators.commutators import commutator
 from qiskit_nature.second_q.properties import AngularMomentum
@@ -103,5 +109,66 @@ def test_commutator_s2z(self) -> None:
         self.assertEqual(commutator(s_2, s_z), FermionicOp.zero())
 
 
+class TestSOperatorsWithOverlap(QiskitNatureTestCase):
+    """Tests for the spin operator generator functions with non-identity overlap.
+
+    See also ``TestAngularMomentum.test_with_overlap``.
+    """
+
+    def setUp(self) -> None:
+        super().setUp()
+
+        self.norb = 2
+        self.nelec = (1, 1)
+
+        mo_coeff = np.asarray([[0.54830202, 1.21832731], [0.54830202, -1.21832731]])
+        ovlp = np.asarray([[1.0, 0.66314574], [0.66314574, 1.0]])
+
+        # first, we ensure that our restricted-spin MO coefficients and overlap match
+        self.assertTrue(np.allclose(mo_coeff.T @ ovlp @ mo_coeff, np.eye(self.norb)))
+
+        theta = 33
+        rot = np.asarray(
+            [
+                [np.cos(np.deg2rad(theta)), -np.sin(np.deg2rad(theta))],
+                [np.sin(np.deg2rad(theta)), np.cos(np.deg2rad(theta))],
+            ],
+        )
+        mo_coeff_rot = mo_coeff @ rot
+
+        self.ovlpab = mo_coeff.T @ ovlp @ mo_coeff_rot
+
+        self.mapper = ParityMapper(self.nelec)
+        self.hf_state = HartreeFock(self.norb, self.nelec, self.mapper)
+
+    def test_s_plus_operator(self) -> None:
+        """Tests the $S^+$ operator with non-identity overlap."""
+        s_p = s_plus_operator(self.norb, self.ovlpab)
+        qubit_op = self.mapper.map(s_p)
+        result = estimate_observables(Estimator(), self.hf_state, {"S+": qubit_op})
+        self.assertAlmostEqual(result["S+"][0], -0.2723195166091395 + 0.2723195166091395j)
+
+    def test_s_minus_operator(self) -> None:
+        """Tests the $S^-$ operator with non-identity overlap."""
+        s_m = s_minus_operator(self.norb, self.ovlpab.T)
+        qubit_op = self.mapper.map(s_m)
+        result = estimate_observables(Estimator(), self.hf_state, {"S-": qubit_op})
+        self.assertAlmostEqual(result["S-"][0], -0.2723195166091395 - 0.2723195166091395j)
+
+    def test_s_x_operator(self) -> None:
+        """Tests the $S^x$ operator with non-identity overlap."""
+        s_x = s_x_operator(self.norb, self.ovlpab)
+        qubit_op = self.mapper.map(s_x)
+        result = estimate_observables(Estimator(), self.hf_state, {"Sx": qubit_op})
+        self.assertAlmostEqual(result["Sx"][0], -0.27231951660913956)
+
+    def test_s_y_operator(self) -> None:
+        """Tests the $S^y$ operator with non-identity overlap."""
+        s_y = s_y_operator(self.norb, self.ovlpab)
+        qubit_op = self.mapper.map(s_y)
+        result = estimate_observables(Estimator(), self.hf_state, {"Sy": qubit_op})
+        self.assertAlmostEqual(result["Sy"][0], 0.27231951660913956)
+
+
 if __name__ == "__main__":
     unittest.main()