feat: improve the implementation of the AngularMomentum operator (qis…

…kit-community#1239)

The AngularMomentum operator is actually just the $S^2$ operator.
However, in the code prior to this commit this was really not easy to
follow.

Instead, this commit changes the implementation to rely on the
sub-operators which the $S^2$ operator consists of. More specifically,
it adds generator functions for the $S^+$ and $S^-$ ladder operators as
well as the $S^x$, $S^y$, and $S^z$ operators. All of these may also be
useful to users who want to evaluate these separately.

Unittests are added to assert the relations of all of these operators.

docs/apidocs/qiskit_nature.second_q.properties.s_operators.rst

@@ -0,0 +1,28 @@
+
+
+s_operators
+=============================================
+
+.. automodule:: qiskit_nature.second_q.properties.s_operators
+
+
+
+   .. rubric:: Functions
+
+   .. autosummary::
+
+      s_minus_operator
+      s_plus_operator
+      s_x_operator
+      s_y_operator
+      s_z_operator
+
+
+
+
+
+
+
+
+
+

qiskit_nature/second_q/properties/__init__.py

@@ -36,6 +36,14 @@
    Magnetization
    ParticleNumber
 
+The following submodule provides multiple functions which can be used to generator more modular
+spin-operator components.
+
+.. autosummary::
+   :toctree:
+
+   s_operators
+
 Vibrational Properties
 ----------------------
 

qiskit_nature/second_q/properties/angular_momentum.py

@@ -16,18 +16,26 @@
 
 from typing import Mapping
 
-import itertools
-
-import numpy as np
-
 import qiskit_nature  # pylint: disable=unused-import
-from qiskit_nature.second_q.operators import FermionicOp, PolynomialTensor
-from qiskit_nature.second_q.operators.tensor_ordering import IndexType, to_physicist_ordering
+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.
 
+    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)
+
+    See also:
+        - the $S^z$ operator: :func:`.s_z_operator`
+        - the $S^+$ operator: :func:`.s_plus_operator`
+        - the $S^-$ operator: :func:`.s_minus_operator`
+
     The following attributes can be set via the initializer but can also be read and updated once
     the ``AngularMomentum`` object has been constructed.
 
@@ -48,17 +56,11 @@ def second_q_ops(self) -> Mapping[str, FermionicOp]:
         Returns:
             A mapping of strings to `FermionicOp` objects.
         """
-        x_h1, x_h2 = _calc_s_x_squared_ints(self.num_spatial_orbitals)
-        y_h1, y_h2 = _calc_s_y_squared_ints(self.num_spatial_orbitals)
-        z_h1, z_h2 = _calc_s_z_squared_ints(self.num_spatial_orbitals)
-        h_1 = x_h1 + y_h1 + z_h1
-        h_2 = x_h2 + y_h2 + z_h2
-
-        tensor = PolynomialTensor(
-            {"+-": h_1, "++--": to_physicist_ordering(h_2, index_order=IndexType.CHEMIST)}
-        )
+        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)
 
-        op = FermionicOp.from_polynomial_tensor(tensor).simplify()
+        op = s_m @ s_p + s_z @ (s_z + FermionicOp.one())
 
         return {self.__class__.__name__: op}
 
@@ -86,144 +88,3 @@ def interpret(
                 result.total_angular_momentum.append(aux_op_eigenvalues[_key].real)
             else:
                 result.total_angular_momentum.append(None)
-
-
-def _calc_s_x_squared_ints(num_spatial_orbitals: int) -> tuple[np.ndarray, np.ndarray]:
-    return _calc_squared_ints(
-        num_spatial_orbitals, _modify_s_x_squared_ints_neq, _modify_s_x_squared_ints_eq
-    )
-
-
-def _calc_s_y_squared_ints(num_spatial_orbitals: int) -> tuple[np.ndarray, np.ndarray]:
-    return _calc_squared_ints(
-        num_spatial_orbitals, _modify_s_y_squared_ints_neq, _modify_s_y_squared_ints_eq
-    )
-
-
-def _calc_s_z_squared_ints(num_spatial_orbitals: int) -> tuple[np.ndarray, np.ndarray]:
-    return _calc_squared_ints(
-        num_spatial_orbitals, _modify_s_z_squared_ints_neq, _modify_s_z_squared_ints_eq
-    )
-
-
-def _calc_squared_ints(
-    num_spatial_orbitals: int, func_neq, func_eq
-) -> tuple[np.ndarray, np.ndarray]:
-    # calculates 1- and 2-body integrals for a given angular momentum axis (x or y or z,
-    # specified by func_neq and func_eq)
-    num_spin_orbitals = 2 * num_spatial_orbitals
-    h_1 = np.zeros((num_spin_orbitals, num_spin_orbitals))
-    h_2 = np.zeros((num_spin_orbitals, num_spin_orbitals, num_spin_orbitals, num_spin_orbitals))
-
-    # pylint: disable=invalid-name
-    for p, q in itertools.product(range(num_spatial_orbitals), repeat=2):
-        if p != q:
-            h_2 = func_neq(h_2, p, q, num_spatial_orbitals)
-        else:
-            h_2 = func_eq(h_2, p, num_spatial_orbitals)
-            h_1[p, p] += 1.0
-            h_1[p + num_spatial_orbitals, p + num_spatial_orbitals] += 1.0
-    h_1 *= 0.25
-    h_2 *= 0.25
-    return h_1, h_2
-
-
-def _modify_s_x_squared_ints_neq(
-    h_2: np.ndarray, p_ind: int, q_ind: int, num_spatial_orbitals: int
-) -> np.ndarray:
-    indices = [
-        (p_ind, p_ind + num_spatial_orbitals, q_ind, q_ind + num_spatial_orbitals),
-        (p_ind + num_spatial_orbitals, p_ind, q_ind, q_ind + num_spatial_orbitals),
-        (p_ind, p_ind + num_spatial_orbitals, q_ind + num_spatial_orbitals, q_ind),
-        (p_ind + num_spatial_orbitals, p_ind, q_ind + num_spatial_orbitals, q_ind),
-    ]
-    values = [1, 1, 1, 1]
-    # adds provided values to values of 2-body integrals (x axis of angular momentum) at given
-    # indices in case p not equal to q
-    return _add_values_to_s_squared_ints(h_2, indices, values)
-
-
-def _modify_s_x_squared_ints_eq(
-    h_2: np.ndarray, p_ind: int, num_spatial_orbitals: int
-) -> np.ndarray:
-    indices = [
-        (p_ind, p_ind + num_spatial_orbitals, p_ind, p_ind + num_spatial_orbitals),
-        (p_ind + num_spatial_orbitals, p_ind, p_ind + num_spatial_orbitals, p_ind),
-        (p_ind, p_ind, p_ind + num_spatial_orbitals, p_ind + num_spatial_orbitals),
-        (p_ind + num_spatial_orbitals, p_ind + num_spatial_orbitals, p_ind, p_ind),
-    ]
-    values = [-1, -1, -1, -1]
-    # adds provided values to values of 2-body integrals (x axis of angular momentum) at given
-    # indices in case p equal to q
-    return _add_values_to_s_squared_ints(h_2, indices, values)
-
-
-def _modify_s_y_squared_ints_neq(
-    h_2: np.ndarray, p_ind: int, q_ind: int, num_spatial_orbitals: int
-) -> np.ndarray:
-    indices = [
-        (p_ind, p_ind + num_spatial_orbitals, q_ind, q_ind + num_spatial_orbitals),
-        (p_ind + num_spatial_orbitals, p_ind, q_ind, q_ind + num_spatial_orbitals),
-        (p_ind, p_ind + num_spatial_orbitals, q_ind + num_spatial_orbitals, q_ind),
-        (p_ind + num_spatial_orbitals, p_ind, q_ind + num_spatial_orbitals, q_ind),
-    ]
-    values = [-1, 1, 1, -1]
-    # adds provided values to values of 2-body integrals (y axis of angular momentum) at given
-    # indices in case p not equal to q
-    return _add_values_to_s_squared_ints(h_2, indices, values)
-
-
-def _modify_s_y_squared_ints_eq(
-    h_2: np.ndarray, p_ind: int, num_spatial_orbitals: int
-) -> np.ndarray:
-    indices = [
-        (p_ind, p_ind + num_spatial_orbitals, p_ind, p_ind + num_spatial_orbitals),
-        (p_ind + num_spatial_orbitals, p_ind, p_ind + num_spatial_orbitals, p_ind),
-        (p_ind, p_ind, p_ind + num_spatial_orbitals, p_ind + num_spatial_orbitals),
-        (p_ind + num_spatial_orbitals, p_ind + num_spatial_orbitals, p_ind, p_ind),
-    ]
-    values = [1, 1, -1, -1]
-    # adds provided values to values of 2-body integrals (y axis of angular momentum) at given
-    # indices in case p equal to q
-    return _add_values_to_s_squared_ints(h_2, indices, values)
-
-
-def _modify_s_z_squared_ints_neq(
-    h_2: np.ndarray, p_ind: int, q_ind: int, num_spatial_orbitals: int
-) -> np.ndarray:
-    indices = [
-        (p_ind, p_ind, q_ind, q_ind),
-        (p_ind + num_spatial_orbitals, p_ind + num_spatial_orbitals, q_ind, q_ind),
-        (p_ind, p_ind, q_ind + num_spatial_orbitals, q_ind + num_spatial_orbitals),
-        (
-            p_ind + num_spatial_orbitals,
-            p_ind + num_spatial_orbitals,
-            q_ind + num_spatial_orbitals,
-            q_ind + num_spatial_orbitals,
-        ),
-    ]
-    values = [1, -1, -1, 1]
-    # adds provided values to values of 2-body integrals (z axis of angular momentum) at given
-    # indices in case p not equal to q
-    return _add_values_to_s_squared_ints(h_2, indices, values)
-
-
-def _modify_s_z_squared_ints_eq(
-    h_2: np.ndarray, p_ind: int, num_spatial_orbitals: int
-) -> np.ndarray:
-    indices = [
-        (p_ind, p_ind + num_spatial_orbitals, p_ind + num_spatial_orbitals, p_ind),
-        (p_ind + num_spatial_orbitals, p_ind, p_ind, p_ind + num_spatial_orbitals),
-    ]
-    values = [1, 1]
-    # adds provided values to values of 2-body integrals (z axis of angular momentum) at given
-    # indices in case p equal to q
-    return _add_values_to_s_squared_ints(h_2, indices, values)
-
-
-def _add_values_to_s_squared_ints(
-    h_2: np.ndarray, indices: list[tuple[int, int, int, int]], values: list[int]
-) -> np.ndarray:
-    for index, value in zip(indices, values):
-        h_2[index] += value
-    return h_2

qiskit_nature/second_q/properties/s_operators.py

@@ -0,0 +1,149 @@
+# This code is part of a Qiskit project.
+#
+# (C) Copyright IBM 2023.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# 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."""
+
+from __future__ import annotations
+
+from qiskit_nature.second_q.operators import FermionicOp
+
+
+def s_plus_operator(num_spatial_orbitals: int) -> 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}
+
+    Here, $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.
+
+    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:
+    r"""Constructs the $S^-$ operator.
+
+    The $S^-$ operator is defined as:
+
+    .. math::
+
+       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
+    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.
+
+    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:
+    r"""Constructs the $S^x$ operator.
+
+    The $S^x$ operator is defined as:
+
+    .. math::
+
+       S^x = \frac{1}{2} \left(S^+ + S^-\right)
+
+    Args:
+        num_spatial_orbitals: the size of the operator which to generate.
+
+    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)
+        }
+    )
+    return op
+
+
+def s_y_operator(num_spatial_orbitals: int) -> FermionicOp:
+    r"""Constructs the $S^y$ operator.
+
+    The $S^y$ operator is defined as:
+
+    .. math::
+
+       S^y = -\frac{i}{2} \left(S^+ - S^-\right)
+
+    Args:
+        num_spatial_orbitals: the size of the operator which to generate.
+
+    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)
+        }
+    )
+    return op
+
+
+def s_z_operator(num_spatial_orbitals: int) -> FermionicOp:
+    r"""Constructs the $S^z$ operator.
+
+    The $S^z$ operator is defined as:
+
+    .. math::
+
+       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)
+
+    Here, $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.
+
+    Args:
+        num_spatial_orbitals: the size of the operator which to generate.
+
+    Returns:
+        The $S^z$ operator of the requested size.
+    """
+    op = FermionicOp(
+        {
+            f"+_{orb} -_{orb}": 0.5 * (-1.0) ** (orb >= num_spatial_orbitals)
+            for orb in range(2 * num_spatial_orbitals)
+        }
+    )
+    return op