feat: additional ways to construct an ElectronicDensity

qiskit_nature/second_q/properties/electronic_density.py

@@ -21,6 +21,7 @@
 
 import qiskit_nature  # pylint: disable=unused-import
 from qiskit_nature.second_q.operators import ElectronicIntegrals, FermionicOp
+from qiskit_nature.utils import get_einsum
 
 
 class ElectronicDensity(ElectronicIntegrals):
@@ -33,6 +34,109 @@ class ElectronicDensity(ElectronicIntegrals):
     :meth:`~qiskit_nature.second_q.operators.ElectronicIntegrals.trace_spin`.
     """
 
+    @staticmethod
+    def _rdm2_from_rdm1s(
+        rdm1_a: np.ndarray, rdm1_b: np.ndarray, *, mixed_spin: bool = False
+    ) -> np.ndarray:
+        einsum_func, _ = get_einsum()
+        rdm2 = einsum_func("ij,kl->ijkl", rdm1_a, rdm1_b)
+        if not mixed_spin:
+            rdm2 -= einsum_func("ij,kl->ikjl", rdm1_a, rdm1_b)
+
+        return rdm2
+
+    @classmethod
+    def empty(cls, num_spatial_orbitals: int, *, include_rdm2: bool = True) -> ElectronicDensity:
+        """Initializes an empty (all-zero) ``ElectronicDensity``.
+
+        Args:
+            num_spatial_orbitals: the number of spatial orbitals.
+            include_rdm2: whether to include the 2-body RDMs.
+
+        Returns:
+            The resulting ``ElectronicDensity``.
+        """
+        rdm1 = np.zeros((num_spatial_orbitals, num_spatial_orbitals), dtype=float)
+        rdm2 = (
+            np.zeros(
+                (
+                    num_spatial_orbitals,
+                    num_spatial_orbitals,
+                    num_spatial_orbitals,
+                    num_spatial_orbitals,
+                ),
+                dtype=float,
+            )
+            if include_rdm2
+            else None
+        )
+
+        return ElectronicDensity.from_raw_integrals(
+            rdm1, rdm2, rdm1, rdm2, rdm2, auto_index_order=False
+        )
+
+    @classmethod
+    def identity(cls, num_spatial_orbitals: int, *, include_rdm2: bool = True) -> ElectronicDensity:
+        """Initializes an identity ``ElectronicDensity``.
+
+        Args:
+            num_spatial_orbitals: the number of spatial orbitals.
+            include_rdm2: whether to include the 2-body RDMs.
+
+        Returns:
+            The resulting ``ElectronicDensity``.
+        """
+        rdm1 = np.eye(num_spatial_orbitals, dtype=float)
+        rdm2: np.ndarray | None = None
+        rdm2_ba: np.ndarray | None = None
+        if include_rdm2:
+            rdm2 = ElectronicDensity._rdm2_from_rdm1s(rdm1, rdm1)
+            rdm2_ba = ElectronicDensity._rdm2_from_rdm1s(rdm1, rdm1, mixed_spin=True)
+
+        return ElectronicDensity.from_raw_integrals(
+            rdm1, rdm2, rdm1, rdm2, rdm2_ba, auto_index_order=False
+        )
+
+    @staticmethod
+    def from_particle_number(
+        num_spatial_orbitals: int,
+        num_particles: int | tuple[int, int],
+        *,
+        include_rdm2: bool = True,
+    ) -> ElectronicDensity:
+        """Initializes an ``ElectronicDensity`` from the provided number of particles.
+
+        Args:
+            num_spatial_orbitals: the number of spatial orbitals.
+            num_particles: the number of particles. If this is an integer it is interpreted as the
+                total number of particles. If it is a tuple of two integers, these are treated as
+                the number of alpha- and beta-spin particles, respectively.
+            include_rdm2: whether to include the 2-body RDMs.
+
+        Returns:
+            The resulting ``ElectronicDensity``.
+        """
+        if isinstance(num_particles, int):
+            num_beta = num_particles // 2
+            num_alpha = num_particles - num_beta
+        else:
+            num_alpha, num_beta = num_particles
+
+        rdm1_a = np.diag([1.0 if i < num_alpha else 0.0 for i in range(num_spatial_orbitals)])
+        rdm1_b = np.diag([1.0 if i < num_beta else 0.0 for i in range(num_spatial_orbitals)])
+
+        rdm2_aa: np.ndarray | None = None
+        rdm2_bb: np.ndarray | None = None
+        rdm2_ba: np.ndarray | None = None
+        if include_rdm2:
+            rdm2_aa = ElectronicDensity._rdm2_from_rdm1s(rdm1_a, rdm1_a)
+            rdm2_bb = ElectronicDensity._rdm2_from_rdm1s(rdm1_b, rdm1_b)
+            rdm2_ba = ElectronicDensity._rdm2_from_rdm1s(rdm1_b, rdm1_a, mixed_spin=True)
+
+        return ElectronicDensity.from_raw_integrals(
+            rdm1_a, rdm2_aa, rdm1_b, rdm2_bb, rdm2_ba, auto_index_order=False
+        )
+
     @staticmethod
     def from_orbital_occupation(
         alpha_occupation: Sequence[float],
@@ -60,13 +164,9 @@ def from_orbital_occupation(
         rdm2_bb: np.ndarray | None = None
         rdm2_ba: np.ndarray | None = None
         if include_rdm2:
-            rdm2_aa = np.einsum("ij,kl->ijkl", rdm1_a, rdm1_a) - np.einsum(
-                "ij,kl->iklj", rdm1_a, rdm1_a
-            )
-            rdm2_bb = np.einsum("ij,kl->ijkl", rdm1_b, rdm1_b) - np.einsum(
-                "ij,kl->iklj", rdm1_b, rdm1_b
-            )
-            rdm2_ba = np.einsum("ij,kl->ijkl", rdm1_a, rdm1_b)
+            rdm2_aa = ElectronicDensity._rdm2_from_rdm1s(rdm1_a, rdm1_a)
+            rdm2_bb = ElectronicDensity._rdm2_from_rdm1s(rdm1_b, rdm1_b)
+            rdm2_ba = ElectronicDensity._rdm2_from_rdm1s(rdm1_b, rdm1_a, mixed_spin=True)
 
         return ElectronicDensity.from_raw_integrals(
             rdm1_a, rdm2_aa, rdm1_b, rdm2_bb, rdm2_ba, auto_index_order=False

releasenotes/notes/feat-electronic-denisty-initializations-d1bd5131c5776418.yaml

@@ -0,0 +1,19 @@
+---
+features:
+  - |
+    Three news means of initializing an :class:`~qiskit_nature.second_q.properties.ElectronicDensity`
+    have been added.
+
+    .. code-block:: python
+
+      size = 4  # number of spatial orbitals
+
+      from qiskit_nature.second_q.properties import ElectronicDensity
+
+      empty = ElectronicDensity.empty(size, include_rdm2=True)
+
+      identity = ElectronicDensity.identity(size, include_rdm2=True)
+
+      num_particles = (2, 2)
+
+      two_and_two = ElectronicDensity.from_particle_number(size, num_particles)

test/second_q/properties/test_electronic_density.py

@@ -19,7 +19,7 @@
 from test import QiskitNatureTestCase
 
 import numpy as np
-from ddt import ddt, data
+from ddt import ddt, data, unpack
 
 import qiskit_nature.optionals as _optionals
 
@@ -39,6 +39,124 @@
 class TestElectronicDensity(QiskitNatureTestCase):
     """Test ElectronicDensity Property"""
 
+    @unpack
+    @data(
+        (1, True),
+        (1, False),
+        (2, True),
+        (2, False),
+        (3, True),
+        (3, False),
+    )
+    def test_empty(self, size: int, include_rdm2: bool):
+        """Test ElectronicDensity.empty."""
+        expected = ElectronicIntegrals.from_raw_integrals(
+            np.zeros((size, size)),
+            np.zeros((size, size, size, size)) if include_rdm2 else None,
+            np.zeros((size, size)),
+            np.zeros((size, size, size, size)) if include_rdm2 else None,
+            np.zeros((size, size, size, size)) if include_rdm2 else None,
+            auto_index_order=False,
+        )
+        density = ElectronicDensity.empty(size, include_rdm2=include_rdm2)
+        self.assertTrue(density.equiv(expected))
+
+    @unpack
+    @data(
+        (1, True),
+        (1, False),
+        (2, True),
+        (2, False),
+        (3, True),
+        (3, False),
+    )
+    def test_identity(self, size: int, include_rdm2: bool):
+        """Test ElectronicDensity.identity."""
+        rdm2: np.ndarray | None = None
+        rdm2_ba: np.ndarray | None = None
+        if include_rdm2:
+            rdm2 = np.zeros((size, size, size, size))
+            rdm2_ba = np.zeros((size, size, size, size))
+            if size == 1:
+                rdm2_ba[0, 0, 0, 0] = 1.0
+            elif size == 2:
+                rdm2[0, 0, 1, 1] = 1.0
+                rdm2[0, 1, 0, 1] = -1.0
+                rdm2[1, 0, 1, 0] = -1.0
+                rdm2[1, 1, 0, 0] = 1.0
+                rdm2_ba[0, 0, 0, 0] = 1.0
+                rdm2_ba[0, 0, 1, 1] = 1.0
+                rdm2_ba[1, 1, 0, 0] = 1.0
+                rdm2_ba[1, 1, 1, 1] = 1.0
+            elif size == 3:
+                rdm2[
+                    np.asarray([0, 0, 1, 1, 2, 2]),
+                    np.asarray([0, 0, 1, 1, 2, 2]),
+                    np.asarray([1, 2, 0, 2, 0, 1]),
+                    np.asarray([1, 2, 0, 2, 0, 1]),
+                ] = 1.0
+                rdm2[
+                    np.asarray([0, 0, 1, 1, 2, 2]),
+                    np.asarray([1, 2, 0, 2, 0, 1]),
+                    np.asarray([0, 0, 1, 1, 2, 2]),
+                    np.asarray([1, 2, 0, 2, 0, 1]),
+                ] = -1.0
+                rdm2_ba[
+                    np.asarray([0, 0, 0, 1, 1, 1, 2, 2, 2]),
+                    np.asarray([0, 0, 0, 1, 1, 1, 2, 2, 2]),
+                    np.asarray([0, 1, 2, 0, 1, 2, 0, 1, 2]),
+                    np.asarray([0, 1, 2, 0, 1, 2, 0, 1, 2]),
+                ] = 1.0
+        expected = ElectronicIntegrals.from_raw_integrals(
+            np.eye(size), rdm2, np.eye(size), rdm2, rdm2_ba, auto_index_order=False
+        )
+        density = ElectronicDensity.identity(size, include_rdm2=include_rdm2)
+        self.assertTrue(density.equiv(expected))
+
+    @unpack
+    @data(
+        (2, 2, True),
+        (2, 2, False),
+        (2, (1, 1), True),
+        (2, (1, 1), False),
+        (2, 3, True),
+        (2, 3, False),
+        (2, (2, 1), True),
+        (2, (2, 1), False),
+    )
+    def test_from_particle_number(
+        self, size: int, num_particles: int | tuple[int, int], include_rdm2: bool
+    ):
+        """Test ElectronicDensity.from_particle_number."""
+        if isinstance(num_particles, int):
+            num_beta = num_particles // 2
+            num_alpha = num_particles - num_beta
+        else:
+            num_alpha, num_beta = num_particles
+        rdm1_a = np.asarray([[1, 0], [0, 1 if num_alpha == 2 else 0]])
+        rdm1_b = np.asarray([[1, 0], [0, 0]])
+        rdm2_aa: np.ndarray | None = None
+        rdm2_bb: np.ndarray | None = None
+        rdm2_ba: np.ndarray | None = None
+        if include_rdm2:
+            rdm2_aa = np.zeros((size, size, size, size))
+            rdm2_bb = np.zeros((size, size, size, size))
+            rdm2_ba = np.zeros((size, size, size, size))
+            rdm2_ba[0, 0, 0, 0] = 1.0
+            if num_alpha == 2:
+                rdm2_aa[0, 0, 1, 1] = 1.0
+                rdm2_aa[0, 1, 0, 1] = -1.0
+                rdm2_aa[1, 0, 1, 0] = -1.0
+                rdm2_aa[1, 1, 0, 0] = 1.0
+                rdm2_ba[0, 0, 1, 1] = 1.0
+        expected = ElectronicIntegrals.from_raw_integrals(
+            rdm1_a, rdm2_aa, rdm1_b, rdm2_bb, rdm2_ba, auto_index_order=False
+        )
+        density = ElectronicDensity.from_particle_number(
+            size, num_particles, include_rdm2=include_rdm2
+        )
+        self.assertTrue(density.equiv(expected))
+
     @data(True, False)
     def test_from_orbital_occupation(self, include_rdm2: bool):
         """Test from_orbital_occupation."""