remove hamiltonian_linop

docs/tutorials/02-orbital-rotation.ipynb

@@ -118,9 +118,7 @@
     "one_body_tensor = ffsim.random.random_hermitian(norb, seed=1234)\n",
     "\n",
     "# Convert the Hamiltonian to a LinearOperator\n",
-    "hamiltonian = ffsim.contract.hamiltonian_linop(\n",
-    "    one_body_tensor=one_body_tensor, norb=norb, nelec=nelec\n",
-    ")\n",
+    "hamiltonian = ffsim.contract.one_body_linop(one_body_tensor, norb=norb, nelec=nelec)\n",
     "\n",
     "# Get the ground state of the Hamiltonian\n",
     "eigs, vecs = scipy.sparse.linalg.eigsh(hamiltonian, which=\"LA\", k=1)\n",
@@ -206,7 +204,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.11"
+   "version": "3.10.13"
   },
   "orig_nbformat": 4
  },

python/ffsim/contract/__init__.py

@@ -11,5 +11,5 @@
 """Functions for contracting tensors and constructing linear operators."""
 
 from ffsim.contract.diag_coulomb import contract_diag_coulomb, diag_coulomb_linop
-from ffsim.contract.hamiltonian import hamiltonian_linop, hamiltonian_trace
 from ffsim.contract.num_op_sum import contract_num_op_sum, num_op_sum_linop
+from ffsim.contract.one_body import contract_one_body, one_body_linop

python/ffsim/contract/one_body.py

@@ -0,0 +1,90 @@
+# (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.
+
+"""Contract one-body operator."""
+
+from __future__ import annotations
+
+import numpy as np
+import scipy.sparse.linalg
+from pyscf.fci import cistring
+from pyscf.fci.direct_nosym import contract_1e
+
+from ffsim.states import dim
+
+
+def contract_one_body(
+    vec: np.ndarray,
+    mat: np.ndarray,
+    norb: int,
+    nelec: tuple[int, int],
+    *,
+    linkstr_index_a: np.ndarray,
+    linkstr_index_b: np.ndarray,
+) -> np.ndarray:
+    link_index = (linkstr_index_a, linkstr_index_b)
+    result = contract_1e(mat.real, vec.real, norb, nelec, link_index=link_index).astype(
+        complex
+    )
+    result += 1j * contract_1e(mat.imag, vec.real, norb, nelec, link_index=link_index)
+    result += 1j * contract_1e(mat.real, vec.imag, norb, nelec, link_index=link_index)
+    result -= contract_1e(mat.imag, vec.imag, norb, nelec, link_index=link_index)
+    return result
+
+
+def one_body_linop(
+    mat: np.ndarray, norb: int, nelec: tuple[int, int]
+) -> scipy.sparse.linalg.LinearOperator:
+    r"""Convert a one-body tensor to a linear operator.
+
+    A one-body tensor has the form
+
+    .. math::
+
+        \sum_{ij} M_{ij} a^\dagger{i} a_j
+
+    where :math:`M` is a complex-valued matrix.
+
+    Args:
+        mat: The one-body tensor.
+        norb: The number of spatial orbitals.
+        nelec: The number of alpha and beta electrons.
+
+    Returns:
+        A LinearOperator that implements the action of the one-body tensor.
+    """
+    dim_ = dim(norb, nelec)
+    n_alpha, n_beta = nelec
+    linkstr_index_a = cistring.gen_linkstr_index(range(norb), n_alpha)
+    linkstr_index_b = cistring.gen_linkstr_index(range(norb), n_beta)
+
+    def matvec(vec: np.ndarray):
+        return contract_one_body(
+            vec,
+            mat,
+            norb=norb,
+            nelec=nelec,
+            linkstr_index_a=linkstr_index_a,
+            linkstr_index_b=linkstr_index_b,
+        )
+
+    def rmatvec(vec: np.ndarray):
+        return contract_one_body(
+            vec,
+            mat.T.conj(),
+            norb=norb,
+            nelec=nelec,
+            linkstr_index_a=linkstr_index_a,
+            linkstr_index_b=linkstr_index_b,
+        )
+
+    return scipy.sparse.linalg.LinearOperator(
+        shape=(dim_, dim_), matvec=matvec, rmatvec=rmatvec, dtype=complex
+    )

python/ffsim/hamiltonians/molecular_hamiltonian.py

@@ -13,9 +13,13 @@
 import dataclasses
 
 import numpy as np
+import scipy.sparse.linalg
+from pyscf.fci import cistring
+from pyscf.fci.direct_nosym import absorb_h1e, make_hdiag
+from pyscf.fci.fci_slow import contract_2e
 from scipy.sparse.linalg import LinearOperator
 
-from ffsim.contract.hamiltonian import hamiltonian_linop, hamiltonian_trace
+from ffsim.states import dim
 
 
 @dataclasses.dataclass
@@ -51,20 +55,26 @@ def norb(self):
 
     def _linear_operator_(self, norb: int, nelec: tuple[int, int]) -> LinearOperator:
         """Return a SciPy LinearOperator representing the object."""
-        return hamiltonian_linop(
-            norb=norb,
-            nelec=nelec,
-            one_body_tensor=self.one_body_tensor,
-            two_body_tensor=self.two_body_tensor,
-            constant=self.constant,
+        n_alpha, n_beta = nelec
+        linkstr_index_a = cistring.gen_linkstr_index(range(norb), n_alpha)
+        linkstr_index_b = cistring.gen_linkstr_index(range(norb), n_beta)
+        link_index = (linkstr_index_a, linkstr_index_b)
+        two_body = absorb_h1e(
+            self.one_body_tensor, self.two_body_tensor, norb, nelec, 0.5
+        )
+        dim_ = dim(norb, nelec)
+
+        def matvec(vec: np.ndarray):
+            return self.constant * vec + contract_2e(
+                two_body, vec, norb, nelec, link_index=link_index
+            )
+
+        return scipy.sparse.linalg.LinearOperator(
+            shape=(dim_, dim_), matvec=matvec, rmatvec=matvec, dtype=complex
         )
 
     def _trace_(self, norb: int, nelec: tuple[int, int]) -> float:
         """Return the trace of the object."""
-        return hamiltonian_trace(
-            norb=norb,
-            nelec=nelec,
-            one_body_tensor=self.one_body_tensor,
-            two_body_tensor=self.two_body_tensor,
-            constant=self.constant,
+        return self.constant * dim(norb, nelec) + np.sum(
+            make_hdiag(self.one_body_tensor, self.two_body_tensor, norb, nelec)
         )

tests/fermion_operator_test.py

@@ -507,8 +507,8 @@ def test_linear_operator():
     one_body_tensor = np.zeros((norb, norb))
     one_body_tensor[1, 2] = 0.5
     one_body_tensor[2, 1] = -0.5
-    expected_linop = ffsim.contract.hamiltonian_linop(
-        one_body_tensor=one_body_tensor, norb=norb, nelec=nelec
+    expected_linop = ffsim.contract.one_body_linop(
+        one_body_tensor, norb=norb, nelec=nelec
     )
 
     vec = ffsim.random.random_statevector(ffsim.dim(norb, nelec), seed=rng)

tests/gates/basic_gates_test.py

@@ -96,9 +96,7 @@ def test_apply_givens_rotation(norb: int, nelec: tuple[int, int]):
             a, b = target_orbs
             generator[a, b] = theta
             generator[b, a] = -theta
-            linop = ffsim.contract.hamiltonian_linop(
-                one_body_tensor=generator, norb=norb, nelec=nelec
-            )
+            linop = ffsim.contract.one_body_linop(generator, norb=norb, nelec=nelec)
             expected = scipy.sparse.linalg.expm_multiply(
                 linop, vec, traceA=np.sum(np.abs(generator))
             )
@@ -158,9 +156,7 @@ def test_apply_tunneling_interaction(norb: int, nelec: tuple[int, int]):
             generator = np.zeros((norb, norb))
             generator[i, j] = theta
             generator[j, i] = theta
-            linop = ffsim.contract.hamiltonian_linop(
-                one_body_tensor=generator, norb=norb, nelec=nelec
-            )
+            linop = ffsim.contract.one_body_linop(generator, norb=norb, nelec=nelec)
             expected = scipy.sparse.linalg.expm_multiply(
                 1j * linop, vec, traceA=np.sum(np.abs(generator))
             )
@@ -217,9 +213,7 @@ def test_apply_num_interaction(norb: int, nelec: tuple[int, int]):
         )
         generator = np.zeros((norb, norb))
         generator[target_orb, target_orb] = theta
-        linop = ffsim.contract.hamiltonian_linop(
-            one_body_tensor=generator, norb=norb, nelec=nelec
-        )
+        linop = ffsim.contract.one_body_linop(generator, norb=norb, nelec=nelec)
         expected = scipy.sparse.linalg.expm_multiply(
             1j * linop, vec, traceA=np.sum(np.abs(generator))
         )

tests/gates/diag_coulomb_test.py

@@ -50,8 +50,8 @@ def test_apply_diag_coulomb_evolution(z_representation: bool):
         op = ffsim.contract.diag_coulomb_linop(
             mat, norb=norb, nelec=nelec, z_representation=z_representation
         )
-        orbital_op = ffsim.contract.hamiltonian_linop(
-            one_body_tensor=scipy.linalg.logm(orbital_rotation), norb=norb, nelec=nelec
+        orbital_op = ffsim.contract.one_body_linop(
+            scipy.linalg.logm(orbital_rotation), norb=norb, nelec=nelec
         )
         expected = scipy.sparse.linalg.expm_multiply(
             -orbital_op, vec, traceA=np.sum(np.abs(orbital_rotation))
@@ -102,8 +102,8 @@ def test_apply_diag_coulomb_evolution_alpha_beta(z_representation: bool):
             mat_alpha_beta=mat_alpha_beta,
             z_representation=z_representation,
         )
-        orbital_op = ffsim.contract.hamiltonian_linop(
-            one_body_tensor=scipy.linalg.logm(orbital_rotation), norb=norb, nelec=nelec
+        orbital_op = ffsim.contract.one_body_linop(
+            scipy.linalg.logm(orbital_rotation), norb=norb, nelec=nelec
         )
         expected = scipy.sparse.linalg.expm_multiply(
             -orbital_op, vec, traceA=np.sum(np.abs(orbital_rotation))

tests/gates/num_op_sum_test.py

@@ -78,9 +78,7 @@ def test_apply_quadratic_hamiltonian_evolution():
         result = ffsim.apply_num_op_sum_evolution(
             vec, eigs, time, norb, nelec, orbital_rotation=vecs
         )
-        op = ffsim.contract.hamiltonian_linop(
-            one_body_tensor=mat, norb=norb, nelec=nelec
-        )
+        op = ffsim.contract.one_body_linop(mat, norb=norb, nelec=nelec)
         expected = scipy.sparse.linalg.expm_multiply(
             -1j * time * op, vec, traceA=np.sum(np.abs(mat))
         )