Repository: qiskit-community/ffsim (python)

Issue:
Add DiagonalCoulombHamiltonian
This is a Hamiltonian with a diagonal Coulomb term. See Eq. (1) of https://arxiv.org/pdf/1711.04789.pdf. We can absorb the "number" term into the one body-part.

We can use the corresponding class in [OpenFermion](https://github.com/quantumlib/OpenFermion/blob/master/src/openfermion/ops/representations/diagonal_coulomb_hamiltonian.py) as a reference.

The implementation should be similar to [MolecularHamiltonian](https://github.com/qiskit-community/ffsim/blob/main/python/ffsim/hamiltonians/molecular_hamiltonian.py) or [DoubleFactorizedHamiltonian](https://github.com/qiskit-community/ffsim/blob/main/python/ffsim/hamiltonians/double_factorized_hamiltonian.py). Something like this:

```python
@dataclasses.dataclass(frozen=True)
class DiagonalCoulombHamiltonian:
    one_body_tensor: np.ndarray
    diag_coulomb_mat: np.ndarray
```

We can add a method
```python
@staticmethod
def from_fermion_operator(op: FermionOperator) -> DiagonalCoulombHamiltonian
```
to initialize from a compatible FermionOperator, like the Hubbard model. Passing a FermionOperator with incompatible terms should raise an error.

Pull request description:
Add diagonal Coulomb Hamiltonian class
Fixes #117 

Patch:
diff --git a/python/ffsim/__init__.py b/python/ffsim/__init__.py
index 3d3ccc1b..12f7139a 100644
--- a/python/ffsim/__init__.py
+++ b/python/ffsim/__init__.py
@@ -26,6 +26,7 @@ from ffsim.gates import (
     apply_tunneling_interaction,
 )
 from ffsim.hamiltonians import (
+    DiagonalCoulombHamiltonian,
     DoubleFactorizedHamiltonian,
     MolecularHamiltonian,
     SingleFactorizedHamiltonian,
@@ -94,6 +95,7 @@ from ffsim.variational import (
 )
 
 __all__ = [
+    "DiagonalCoulombHamiltonian",
     "DoubleFactorizedHamiltonian",
     "FermionAction",
     "FermionOperator",
diff --git a/python/ffsim/hamiltonians/__init__.py b/python/ffsim/hamiltonians/__init__.py
index b3184e32..dfa7e79b 100644
--- a/python/ffsim/hamiltonians/__init__.py
+++ b/python/ffsim/hamiltonians/__init__.py
@@ -10,11 +10,13 @@
 
 """Classes for representing Hamiltonians."""
 
+from ffsim.hamiltonians.diagonal_coulomb_hamiltonian import DiagonalCoulombHamiltonian
 from ffsim.hamiltonians.double_factorized_hamiltonian import DoubleFactorizedHamiltonian
 from ffsim.hamiltonians.molecular_hamiltonian import MolecularHamiltonian
 from ffsim.hamiltonians.single_factorized_hamiltonian import SingleFactorizedHamiltonian
 
 __all__ = [
+    "DiagonalCoulombHamiltonian",
     "DoubleFactorizedHamiltonian",
     "MolecularHamiltonian",
     "SingleFactorizedHamiltonian",
diff --git a/python/ffsim/hamiltonians/diagonal_coulomb_hamiltonian.py b/python/ffsim/hamiltonians/diagonal_coulomb_hamiltonian.py
new file mode 100644
index 00000000..75f1b4dd
--- /dev/null
+++ b/python/ffsim/hamiltonians/diagonal_coulomb_hamiltonian.py
@@ -0,0 +1,175 @@
+# (C) Copyright IBM 2024.
+#
+# 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.
+
+from __future__ import annotations
+
+import dataclasses
+import itertools
+
+import numpy as np
+import scipy.linalg
+from scipy.sparse.linalg import LinearOperator
+
+from ffsim._lib import FermionOperator
+from ffsim.contract.diag_coulomb import diag_coulomb_linop
+from ffsim.contract.num_op_sum import num_op_sum_linop
+from ffsim.operators.fermion_action import cre_a, cre_b, des_a, des_b
+from ffsim.states import dim
+
+
+@dataclasses.dataclass(frozen=True)
+class DiagonalCoulombHamiltonian:
+    r"""A diagonal Coulomb Hamiltonian.
+
+    A Hamiltonian of the form
+
+    .. math::
+
+        H = \sum_{\sigma, pq} h_{pq} a^\dagger_{\sigma, p} a_{\sigma, q}
+            + \frac12 \sum_{\sigma \tau, pq} V_{(\sigma \tau), pq} n_{\sigma, p}
+            n_{\tau, q} + \text{constant}.
+
+    where :math:`n_{\sigma, p} = a_{\sigma, p}^\dagger a_{\sigma, p}` is the number
+    operator on orbital :math:`p` with spin :math:`\sigma`.
+
+    Here :math:`h_{pq}` is called the one-body tensor and :math:`V_{(\sigma \tau), pq}`
+    are called the diagonal Coulomb matrices. The brackets indicate that
+    :math:`V_{(\sigma \tau)}` is a circulant matrix, which satisfies
+    :math:`V_{\alpha\alpha}=V_{\beta\beta}` and :math:`V_{\alpha\beta}=V_{\beta\alpha}`.
+
+    Attributes:
+        one_body_tensor (np.ndarray): The one-body tensor :math:`h`.
+        diag_coulomb_mats (np.ndarray): The diagonal Coulomb matrices
+            :math:`V_{(\sigma \tau)}`, given as a pair of Numpy arrays specifying
+            independent coefficients for alpha-alpha and alpha-beta interactions (in
+            that order).
+        constant (float): The constant.
+    """
+
+    one_body_tensor: np.ndarray
+    diag_coulomb_mats: np.ndarray
+    constant: float = 0.0
+
+    @property
+    def norb(self) -> int:
+        """The number of spatial orbitals."""
+        return self.one_body_tensor.shape[0]
+
+    def _linear_operator_(self, norb: int, nelec: tuple[int, int]) -> LinearOperator:
+        """Return a SciPy LinearOperator representing the object."""
+        dim_ = dim(norb, nelec)
+        eigs, vecs = scipy.linalg.eigh(self.one_body_tensor)
+        num_linop = num_op_sum_linop(eigs, norb, nelec, orbital_rotation=vecs)
+        dc_linop = diag_coulomb_linop(
+            (
+                self.diag_coulomb_mats[0],
+                self.diag_coulomb_mats[1],
+                self.diag_coulomb_mats[0],
+            ),
+            norb,
+            nelec,
+        )
+
+        def matvec(vec: np.ndarray):
+            result = self.constant * vec
+            result += num_linop @ vec
+            result += dc_linop @ vec
+            return result
+
+        return LinearOperator(
+            shape=(dim_, dim_), matvec=matvec, rmatvec=matvec, dtype=complex
+        )
+
+    def _fermion_operator_(self) -> FermionOperator:
+        """Return a FermionOperator representing the object."""
+        op = FermionOperator({(): self.constant})
+        for p, q in itertools.product(range(self.norb), repeat=2):
+            coeff = self.one_body_tensor[p, q]
+            op += FermionOperator(
+                {
+                    (cre_a(p), des_a(q)): coeff,
+                    (cre_b(p), des_b(q)): coeff,
+                }
+            )
+        for p, q in itertools.product(range(self.norb), repeat=2):
+            (coeff_same_spin, coeff_diff_spin) = 0.5 * self.diag_coulomb_mats[:, p, q]
+            op += FermionOperator(
+                {
+                    (cre_a(p), des_a(p), cre_a(q), des_a(q)): coeff_same_spin,
+                    (cre_a(p), des_a(p), cre_b(q), des_b(q)): coeff_diff_spin,
+                    (cre_b(p), des_b(p), cre_a(q), des_a(q)): coeff_diff_spin,
+                    (cre_b(p), des_b(p), cre_b(q), des_b(q)): coeff_same_spin,
+                }
+            )
+        return op
+
+    @staticmethod
+    def from_fermion_operator(op: FermionOperator) -> DiagonalCoulombHamiltonian:
+        """Convert a FermionOperator to a DiagonalCoulombHamiltonian."""
+
+        # extract norb
+        norb = 1 + max(orb for term in op for _, _, orb in term)
+
+        # initialize variables
+        constant: float = 0
+        one_body_tensor = np.zeros((norb, norb), dtype=complex)
+        diag_coulomb_mats = np.zeros((2, norb, norb), dtype=float)
+
+        # populate the tensors
+        for term, coeff in op.items():
+            if not term:  # constant term
+                constant = coeff.real
+            elif len(term) == 2:  # one-body term candidate
+                (_, _, p), (_, _, q) = term
+                valid_terms = [(cre_a(p), des_a(q)), (cre_b(p), des_b(q))]
+                if term in valid_terms:
+                    one_body_tensor[p, q] += 0.5 * coeff
+                else:
+                    raise ValueError(
+                        "FermionOperator cannot be converted to "
+                        f"DiagonalCoulombHamiltonian. The one-body term {term} is not "
+                        "of the form a^\\dagger_{\\sigma, p} a_{\\sigma, q}."
+                    )
+            elif len(term) == 4:  # two-body term candidate
+                (_, _, p), _, (_, _, q), _ = term
+                valid_terms_same_spin = [
+                    (cre_a(p), des_a(p), cre_a(q), des_a(q)),
+                    (cre_b(p), des_b(p), cre_b(q), des_b(q)),
+                ]
+                valid_terms_diff_spin = [
+                    (cre_a(p), des_a(p), cre_b(q), des_b(q)),
+                    (cre_b(p), des_b(p), cre_a(q), des_a(q)),
+                ]
+                if term in valid_terms_same_spin:
+                    diag_coulomb_mats[0][p, q] += coeff.real
+                elif term in valid_terms_diff_spin:
+                    diag_coulomb_mats[1][p, q] += coeff.real
+                else:
+                    raise ValueError(
+                        "FermionOperator cannot be converted to "
+                        f"DiagonalCoulombHamiltonian. The two-body term {term} is not "
+                        "of the form n_{\\sigma, p} n_{\\tau, q}."
+                    )
+            else:
+                raise ValueError(
+                    "FermionOperator cannot be converted to DiagonalCoulombHamiltonian."
+                    f" The term {term} is neither a constant, one-body, or two-body "
+                    "term."
+                )
+
+        # ensure diag_coulomb_mats symmetry
+        diag_coulomb_mats += diag_coulomb_mats.transpose(0, 2, 1)
+        diag_coulomb_mats *= 0.5
+
+        return DiagonalCoulombHamiltonian(
+            one_body_tensor=one_body_tensor,
+            diag_coulomb_mats=diag_coulomb_mats,
+            constant=constant,
+        )


Test patch:
diff --git a/tests/python/hamiltonians/diagonal_coulomb_hamiltonian_test.py b/tests/python/hamiltonians/diagonal_coulomb_hamiltonian_test.py
new file mode 100644
index 00000000..cd2a8d7c
--- /dev/null
+++ b/tests/python/hamiltonians/diagonal_coulomb_hamiltonian_test.py
@@ -0,0 +1,240 @@
+# (C) Copyright IBM 2024.
+#
+# 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.
+
+"""Tests for diagonal Coulomb Hamiltonian."""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+import ffsim
+from ffsim import fermi_hubbard_1d, fermi_hubbard_2d
+from ffsim.hamiltonians.diagonal_coulomb_hamiltonian import DiagonalCoulombHamiltonian
+
+TEST_ARGS = pytest.mark.parametrize(
+    "norb, nelec",
+    [
+        (4, (2, 2)),
+        (4, (1, 2)),
+        (4, (0, 2)),
+        (4, (0, 0)),
+    ],
+)
+
+
+@TEST_ARGS
+def test_linear_operator(norb: int, nelec: tuple[int, int]):
+    """Test linear_operator method."""
+    dim = ffsim.dim(norb, nelec)
+    rng = np.random.default_rng()
+
+    one_body_tensor = ffsim.random.random_hermitian(norb, seed=rng)
+    diag_coulomb_mat = ffsim.random.random_real_symmetric_matrix(
+        norb, seed=rng, dtype=float
+    )
+    diag_coulomb_mats = np.array([diag_coulomb_mat, diag_coulomb_mat])
+    constant = rng.standard_normal()
+
+    dc_hamiltonian = ffsim.DiagonalCoulombHamiltonian(
+        one_body_tensor, diag_coulomb_mats, constant
+    )
+    df_hamiltonian = ffsim.DoubleFactorizedHamiltonian(
+        one_body_tensor,
+        np.array([diag_coulomb_mat]),
+        np.array([np.eye(norb)]),
+        constant,
+        z_representation=False,
+    )
+
+    actual_linop = ffsim.linear_operator(dc_hamiltonian, norb, nelec)
+    expected_linop = ffsim.linear_operator(df_hamiltonian, norb, nelec)
+
+    vec = ffsim.random.random_statevector(dim, seed=rng)
+    actual = actual_linop @ vec
+    expected = expected_linop @ vec
+    np.testing.assert_allclose(actual, expected)
+
+
+@TEST_ARGS
+def test_fermion_operator(norb: int, nelec: tuple[int, int]):
+    """Test fermion_operator method."""
+    dim = ffsim.dim(norb, nelec)
+    rng = np.random.default_rng()
+
+    one_body_tensor = ffsim.random.random_hermitian(norb, seed=rng)
+    diag_coulomb_mats = np.empty((2, norb, norb), dtype=float)
+    for i in range(2):
+        diag_coulomb_mat = ffsim.random.random_real_symmetric_matrix(
+            norb, seed=rng, dtype=float
+        )
+        diag_coulomb_mats[i] = diag_coulomb_mat
+    constant = rng.standard_normal()
+
+    dc_hamiltonian = ffsim.DiagonalCoulombHamiltonian(
+        one_body_tensor, diag_coulomb_mats, constant=constant
+    )
+
+    op = ffsim.fermion_operator(dc_hamiltonian)
+    actual_linop = ffsim.linear_operator(op, norb, nelec)
+    expected_linop = ffsim.linear_operator(dc_hamiltonian, norb, nelec)
+
+    vec = ffsim.random.random_statevector(dim, seed=rng)
+    actual = actual_linop @ vec
+    expected = expected_linop @ vec
+    np.testing.assert_allclose(actual, expected)
+
+
+@TEST_ARGS
+def test_from_fermion_operator(norb: int, nelec: tuple[int, int]):
+    """Test from_fermion_operator method."""
+    dim = ffsim.dim(norb, nelec)
+    rng = np.random.default_rng()
+
+    one_body_tensor = ffsim.random.random_hermitian(norb, seed=rng)
+    diag_coulomb_mats = np.empty((2, norb, norb), dtype=float)
+    for i in range(2):
+        diag_coulomb_mat = ffsim.random.random_real_symmetric_matrix(
+            norb, seed=rng, dtype=float
+        )
+        diag_coulomb_mats[i] = diag_coulomb_mat
+    constant = rng.standard_normal()
+
+    dc_hamiltonian = ffsim.DiagonalCoulombHamiltonian(
+        one_body_tensor, diag_coulomb_mats, constant=constant
+    )
+
+    op = ffsim.fermion_operator(dc_hamiltonian)
+    dc_hamiltonian_from_op = DiagonalCoulombHamiltonian.from_fermion_operator(op)
+    actual_linop = ffsim.linear_operator(dc_hamiltonian_from_op, norb, nelec)
+    expected_linop = ffsim.linear_operator(op, norb, nelec)
+
+    vec = ffsim.random.random_statevector(dim, seed=rng)
+    actual = actual_linop @ vec
+    expected = expected_linop @ vec
+    np.testing.assert_allclose(actual, expected)
+
+
+def test_from_fermion_operator_failure():
+    """Test from_fermion_operator method failure."""
+    norb = 4
+    rng = np.random.default_rng()
+
+    one_body_tensor = ffsim.random.random_hermitian(norb, seed=rng)
+    two_body_tensor = ffsim.random.random_two_body_tensor(norb, seed=rng, dtype=float)
+    constant = rng.standard_normal()
+
+    mol_hamiltonian = ffsim.MolecularHamiltonian(
+        one_body_tensor, two_body_tensor, constant=constant
+    )
+
+    op = ffsim.fermion_operator(mol_hamiltonian)
+
+    with pytest.raises(
+        ValueError,
+        match="FermionOperator cannot be converted to DiagonalCoulombHamiltonian",
+    ):
+        DiagonalCoulombHamiltonian.from_fermion_operator(op)
+
+
+@TEST_ARGS
+def test_from_fermion_operator_fermi_hubbard_1d(norb: int, nelec: tuple[int, int]):
+    """Test from_fermion_operator method with the fermi_hubbard_1d model."""
+    dim = ffsim.dim(norb, nelec)
+    rng = np.random.default_rng()
+    vec = ffsim.random.random_statevector(dim, seed=rng)
+
+    # open boundary conditions
+    op = fermi_hubbard_1d(
+        norb=4,
+        tunneling=1,
+        interaction=2,
+        chemical_potential=3,
+        nearest_neighbor_interaction=4,
+    )
+
+    dc_hamiltonian_from_op = DiagonalCoulombHamiltonian.from_fermion_operator(op)
+    actual_linop = ffsim.linear_operator(dc_hamiltonian_from_op, norb, nelec)
+    expected_linop = ffsim.linear_operator(op, norb, nelec)
+
+    actual = actual_linop @ vec
+    expected = expected_linop @ vec
+    np.testing.assert_allclose(actual, expected)
+
+    # periodic boundary conditions
+    op_periodic = fermi_hubbard_1d(
+        norb=4,
+        tunneling=1,
+        interaction=2,
+        chemical_potential=3,
+        nearest_neighbor_interaction=4,
+        periodic=True,
+    )
+
+    dc_hamiltonian_from_op_periodic = DiagonalCoulombHamiltonian.from_fermion_operator(
+        op_periodic
+    )
+    actual_linop_periodic = ffsim.linear_operator(
+        dc_hamiltonian_from_op_periodic, norb, nelec
+    )
+    expected_linop_periodic = ffsim.linear_operator(op_periodic, norb, nelec)
+
+    actual_periodic = actual_linop_periodic @ vec
+    expected_periodic = expected_linop_periodic @ vec
+    np.testing.assert_allclose(actual_periodic, expected_periodic)
+
+
+@TEST_ARGS
+def test_from_fermion_operator_fermi_hubbard_2d(norb: int, nelec: tuple[int, int]):
+    """Test from_fermion_operator method with the fermi_hubbard_2d model."""
+    dim = ffsim.dim(norb, nelec)
+    rng = np.random.default_rng()
+    vec = ffsim.random.random_statevector(dim, seed=rng)
+
+    # open boundary conditions
+    op = fermi_hubbard_2d(
+        norb_x=2,
+        norb_y=2,
+        tunneling=1,
+        interaction=2,
+        chemical_potential=3,
+        nearest_neighbor_interaction=4,
+    )
+
+    dc_hamiltonian_from_op = DiagonalCoulombHamiltonian.from_fermion_operator(op)
+    actual_linop = ffsim.linear_operator(dc_hamiltonian_from_op, norb, nelec)
+    expected_linop = ffsim.linear_operator(op, norb, nelec)
+
+    actual = actual_linop @ vec
+    expected = expected_linop @ vec
+    np.testing.assert_allclose(actual, expected)
+
+    # periodic boundary conditions
+    op_periodic = fermi_hubbard_2d(
+        norb_x=2,
+        norb_y=2,
+        tunneling=1,
+        interaction=2,
+        chemical_potential=3,
+        nearest_neighbor_interaction=4,
+        periodic=True,
+    )
+
+    dc_hamiltonian_from_op_periodic = DiagonalCoulombHamiltonian.from_fermion_operator(
+        op_periodic
+    )
+    actual_linop_periodic = ffsim.linear_operator(
+        dc_hamiltonian_from_op_periodic, norb, nelec
+    )
+    expected_linop_periodic = ffsim.linear_operator(op_periodic, norb, nelec)
+
+    actual_periodic = actual_linop_periodic @ vec
+    expected_periodic = expected_linop_periodic @ vec
+    np.testing.assert_allclose(actual_periodic, expected_periodic)