diff --git a/.pylintdict b/.pylintdict index 3536edbd8a..f3fabf3645 100644 --- a/.pylintdict +++ b/.pylintdict @@ -460,6 +460,7 @@ radians rangle rcccx rccx +rdm readme readthedocs realise diff --git a/qiskit_nature/second_q/properties/electronic_density.py b/qiskit_nature/second_q/properties/electronic_density.py index 54791f1d00..a38a179ffc 100644 --- a/qiskit_nature/second_q/properties/electronic_density.py +++ b/qiskit_nature/second_q/properties/electronic_density.py @@ -1,6 +1,6 @@ # This code is part of Qiskit. # -# (C) Copyright IBM 2022. +# (C) Copyright IBM 2022, 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 @@ -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 diff --git a/releasenotes/notes/feat-electronic-denisty-initializations-d1bd5131c5776418.yaml b/releasenotes/notes/feat-electronic-denisty-initializations-d1bd5131c5776418.yaml new file mode 100644 index 0000000000..07dcac02c3 --- /dev/null +++ b/releasenotes/notes/feat-electronic-denisty-initializations-d1bd5131c5776418.yaml @@ -0,0 +1,40 @@ +--- +features: + - | + Three new methods for creating instances + :class:`~qiskit_nature.second_q.properties.ElectronicDensity` have been + added: + + 1. constructing an empty (or all-zero) density of a given size: + + .. code-block:: python + + empty = ElectronicDensity.empty(num_spatial_orbitals=4) + + 2. constructing an idenity density, meaning that the 1-body matrices are + initialized with identity matrices + + .. code-block:: python + + identity = ElectronicDensity.identity(num_spatial_orbitals=4) + + 3. constructing from a provided number of particles. This is a shorter + variant of the already existing ``from_orbital_occupation`` method for the + most common use-case. + + .. code-block:: python + + num_spatial_orbitals = 4 + num_particles = (2, 2) + + two_and_two = ElectronicDensity.from_particle_number(num_spatial_orbitals, num_particles) + + # for example now the 1-body matrices will be: + # [[1, 0, 0, 0], + # [0, 1, 0, 0], + # [0, 0, 0, 0], + # [0, 0, 0, 0]] + + All of the methods above take the optional keyword-argument ``include_rdm2`` + which determines whether or not the 2-body matrices are computed based on the + constructed 1-body matrices. By default, this is set to ``True``. diff --git a/test/second_q/properties/test_electronic_density.py b/test/second_q/properties/test_electronic_density.py index ebace43672..ed7a6cac46 100644 --- a/test/second_q/properties/test_electronic_density.py +++ b/test/second_q/properties/test_electronic_density.py @@ -1,6 +1,6 @@ # This code is part of Qiskit. # -# (C) Copyright IBM 2022. +# (C) Copyright IBM 2022, 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 @@ -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."""