From 1c4264ce17db15d8240f397e2298d33ede7f2aa0 Mon Sep 17 00:00:00 2001 From: Max Rossmannek Date: Fri, 1 Dec 2023 16:46:23 +0100 Subject: [PATCH 1/7] fix: use the symmetric formula for S^2 Prior to this fix, the `AngularMomentum` operator was working on the assumption that the number of alpha-spin particles would always exceed the number of beta-spin particles (for non-singlet systems). However, this is not guaranteed to be the case. This commit fixes this problem by using a symmetric formula for S^2 which is the average of the two cases (alpha- vs. beta-spin particles exceeding the other). --- qiskit_nature/second_q/properties/angular_momentum.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/qiskit_nature/second_q/properties/angular_momentum.py b/qiskit_nature/second_q/properties/angular_momentum.py index 8738bfd66..dd6806c71 100644 --- a/qiskit_nature/second_q/properties/angular_momentum.py +++ b/qiskit_nature/second_q/properties/angular_momentum.py @@ -31,7 +31,7 @@ class AngularMomentum: .. math:: - S^2 = S^- S^+ + S^z (S^z + 1) + S^2 = (S^+ S^- + S^- S^+) / 2 + S^z S^z .. warning:: @@ -75,7 +75,8 @@ def second_q_ops(self) -> Mapping[str, FermionicOp]: 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()) + spm_smp = (s_p @ s_m + s_m @ s_p).normal_order() + op = 0.5 * spm_smp + s_z @ s_z return {self.__class__.__name__: op} From c2d426adba3e24b3da37401d666a0c254ddd49ad Mon Sep 17 00:00:00 2001 From: Max Rossmannek Date: Fri, 1 Dec 2023 16:48:41 +0100 Subject: [PATCH 2/7] feat: log a warning when the alpha-beta overlap is non-unitary When dealing with active spaces obtained from unrestricted-spin orbitals, the alpha-beta overlap matrix which is used to construct the S^2 operator can become non-unitary (for example, when the active set of alpha- and beta-spin orbitals do not span the same space). This can result in an value measured on the active space that is largely different from (e.g.) the UHF value. More importantly, when this is the case, the difference between these two values is "hidden" in the inactive+inactive and inactive+active interactions. Especially when spin contamination is present, this can lead to vastly unexpected results and may indicate a poor choice of active space. --- .../second_q/properties/angular_momentum.py | 40 ++++++++++++++++++- 1 file changed, 38 insertions(+), 2 deletions(-) diff --git a/qiskit_nature/second_q/properties/angular_momentum.py b/qiskit_nature/second_q/properties/angular_momentum.py index dd6806c71..aa4e96c60 100644 --- a/qiskit_nature/second_q/properties/angular_momentum.py +++ b/qiskit_nature/second_q/properties/angular_momentum.py @@ -14,6 +14,7 @@ from __future__ import annotations +import logging from typing import Mapping import numpy as np @@ -23,6 +24,8 @@ from .s_operators import s_minus_operator, s_plus_operator, s_z_operator +LOGGER = logging.getLogger(__name__) + class AngularMomentum: r"""The AngularMomentum property. @@ -49,8 +52,6 @@ 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, overlap: np.ndarray | None = None) -> None: @@ -61,8 +62,43 @@ def __init__(self, num_spatial_orbitals: int, overlap: np.ndarray | None = None) is `None`, the overlap-matrix is assumed to be identity. """ self.num_spatial_orbitals = num_spatial_orbitals + self._overlap: np.ndarray | None = None self.overlap = overlap + @property + def overlap(self) -> np.ndarray | None: + r"""The overlap-matrix between the $\alpha$- and $\beta$-spin orbitals. + + When this is ``None``, the overlap-matrix is assumed to be identity. + """ + return self._overlap + + @overlap.setter + def overlap(self, overlap: np.ndarray | None) -> None: + self._overlap = overlap + + if overlap is not None: + norb = self.num_spatial_orbitals + delta = np.eye(2 * norb) + delta[:norb,:norb] -= overlap.T @ overlap + delta[norb:,norb:] -= overlap @ overlap.T + summed = np.einsum("ij->", np.abs(delta)) + if not np.isclose(summed, 0.0): + LOGGER.warning( + "The provided alpha-beta overlap matrix is NOT unitary! This can happen when " + "the alpha- and beta-spin orbitals do not span the same space. To provide an " + "example of what this means, consider an active space chosen from unrestricted-" + "spin orbitals. Computing within this active space may not result in the " + "same value as obtained on the single-reference starting point. More " + "importantly, this implies that the inactive subspace will account for the " + "difference between these two values, possibly resulting in significant " + "spin contamination in both subspaces. You should verify whether this is " + "intentional/acceptable or whether your choice of active space can be improved." + " As a reference, here is the summed-absolute deviation of `S^T @ S` from the " + f"identity: {summed}" + ) + + def second_q_ops(self) -> Mapping[str, FermionicOp]: """Returns the second quantized angular momentum operator. From abc6e0eda6ce3717846e8bf5e1a101fe9e96f22e Mon Sep 17 00:00:00 2001 From: Max Rossmannek Date: Fri, 1 Dec 2023 17:05:24 +0100 Subject: [PATCH 3/7] Add reno --- .../notes/fix-angular-momentum-2-cebabde075b61a4e.yaml | 10 ++++++++++ 1 file changed, 10 insertions(+) create mode 100644 releasenotes/notes/fix-angular-momentum-2-cebabde075b61a4e.yaml diff --git a/releasenotes/notes/fix-angular-momentum-2-cebabde075b61a4e.yaml b/releasenotes/notes/fix-angular-momentum-2-cebabde075b61a4e.yaml new file mode 100644 index 000000000..d8efec56d --- /dev/null +++ b/releasenotes/notes/fix-angular-momentum-2-cebabde075b61a4e.yaml @@ -0,0 +1,10 @@ +--- +features: + - | + The :meth:`.AngularMomentum.overlap` property will now warn the user, when + this matrix is non-unitary. +fixes: + - | + Fixes the :class:`.AngularMomentum` operator further to also support cases + where the number of beta-spin particles exceeds the number of alpha-spin + particles. From 89f15aea7b2847da1da0934d285783e9d0878182 Mon Sep 17 00:00:00 2001 From: Max Rossmannek Date: Fri, 1 Dec 2023 17:15:16 +0100 Subject: [PATCH 4/7] test: add a regression test --- .../properties/test_angular_momentum.py | 27 +++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/test/second_q/properties/test_angular_momentum.py b/test/second_q/properties/test_angular_momentum.py index 32dc7576c..f6ca0baf1 100644 --- a/test/second_q/properties/test_angular_momentum.py +++ b/test/second_q/properties/test_angular_momentum.py @@ -104,6 +104,33 @@ def test_with_overlap(self): result = estimate_observables(Estimator(), hf_state, qubit_op) self.assertAlmostEqual(result["AngularMomentum"][0], 0.29663167846210015) + def test_with_non_unitary_overlap(self): + """Tests the result with a non-unitary overlap. + + This is a regression test against + https://github.com/qiskit-community/qiskit-nature/issues/1291. + """ + norb = 4 + nelec = (4, 2) + ovlpab = np.asarray( + [ + [0.987256, -0.001123, 0.00006, -0.0], + [-0.001123, -0.987256, -0.0, -0.00006], + [0.000019, 0.000055, -0.3195, -0.931662], + [0.000056, -0.000019, -0.931662, 0.3195], + ], + ) + + 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], 1.9700743392855005) + if __name__ == "__main__": unittest.main() From 93b8bb6377f6ca4562e82ec3a867bc8ee03b4428 Mon Sep 17 00:00:00 2001 From: Max Rossmannek Date: Fri, 1 Dec 2023 17:15:29 +0100 Subject: [PATCH 5/7] fix: update warning tolerance --- qiskit_nature/second_q/properties/angular_momentum.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qiskit_nature/second_q/properties/angular_momentum.py b/qiskit_nature/second_q/properties/angular_momentum.py index aa4e96c60..f86574a57 100644 --- a/qiskit_nature/second_q/properties/angular_momentum.py +++ b/qiskit_nature/second_q/properties/angular_momentum.py @@ -83,7 +83,7 @@ def overlap(self, overlap: np.ndarray | None) -> None: delta[:norb,:norb] -= overlap.T @ overlap delta[norb:,norb:] -= overlap @ overlap.T summed = np.einsum("ij->", np.abs(delta)) - if not np.isclose(summed, 0.0): + if not np.isclose(summed, 0.0, atol=1e-6): LOGGER.warning( "The provided alpha-beta overlap matrix is NOT unitary! This can happen when " "the alpha- and beta-spin orbitals do not span the same space. To provide an " From 212533d4565d07870ea91988926426d5956c21ff Mon Sep 17 00:00:00 2001 From: Max Rossmannek Date: Fri, 1 Dec 2023 17:43:34 +0100 Subject: [PATCH 6/7] Linting and formatting --- qiskit_nature/second_q/properties/angular_momentum.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/qiskit_nature/second_q/properties/angular_momentum.py b/qiskit_nature/second_q/properties/angular_momentum.py index f86574a57..686361e08 100644 --- a/qiskit_nature/second_q/properties/angular_momentum.py +++ b/qiskit_nature/second_q/properties/angular_momentum.py @@ -80,8 +80,8 @@ def overlap(self, overlap: np.ndarray | None) -> None: if overlap is not None: norb = self.num_spatial_orbitals delta = np.eye(2 * norb) - delta[:norb,:norb] -= overlap.T @ overlap - delta[norb:,norb:] -= overlap @ overlap.T + delta[:norb, :norb] -= overlap.T @ overlap + delta[norb:, norb:] -= overlap @ overlap.T summed = np.einsum("ij->", np.abs(delta)) if not np.isclose(summed, 0.0, atol=1e-6): LOGGER.warning( @@ -95,10 +95,10 @@ def overlap(self, overlap: np.ndarray | None) -> None: "spin contamination in both subspaces. You should verify whether this is " "intentional/acceptable or whether your choice of active space can be improved." " As a reference, here is the summed-absolute deviation of `S^T @ S` from the " - f"identity: {summed}" + "identity: %s", + str(summed), ) - def second_q_ops(self) -> Mapping[str, FermionicOp]: """Returns the second quantized angular momentum operator. From 531bb81f769fb5126c913ed9f3acf0995ca1918e Mon Sep 17 00:00:00 2001 From: Max Rossmannek Date: Fri, 1 Dec 2023 17:50:15 +0100 Subject: [PATCH 7/7] docs: fix verbatim in RST --- qiskit_nature/second_q/properties/angular_momentum.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qiskit_nature/second_q/properties/angular_momentum.py b/qiskit_nature/second_q/properties/angular_momentum.py index 686361e08..db5cf8c57 100644 --- a/qiskit_nature/second_q/properties/angular_momentum.py +++ b/qiskit_nature/second_q/properties/angular_momentum.py @@ -59,7 +59,7 @@ def __init__(self, num_spatial_orbitals: int, overlap: np.ndarray | None = None) 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. + is ``None``, the overlap-matrix is assumed to be identity. """ self.num_spatial_orbitals = num_spatial_orbitals self._overlap: np.ndarray | None = None