Support complex two-body tensor linear operator (#279)

* support linop for complex two-body tensor

* require pyscf >= 2.7

* restore single-factorized hamiltonian tests

single- or double- factorization does not support complex two-body operators

pyproject.toml

@@ -26,7 +26,7 @@ dependencies = [
     "numpy",
     "opt_einsum",
     "orjson",
-    "pyscf >= 2.4",
+    "pyscf >= 2.7",
     "qiskit >= 1.1",
     "scipy",
     "typing-extensions",

python/ffsim/hamiltonians/molecular_hamiltonian.py

@@ -18,7 +18,7 @@
 import pyscf.ao2mo
 import pyscf.tools
 from opt_einsum import contract
-from pyscf.fci.direct_nosym import absorb_h1e, contract_1e, contract_2e, make_hdiag
+from pyscf.fci.direct_nosym import absorb_h1e, contract_2e, make_hdiag
 from scipy.sparse.linalg import LinearOperator
 from typing_extensions import deprecated
 
@@ -122,34 +122,22 @@ def rotated(self, orbital_rotation: np.ndarray) -> MolecularHamiltonian:
 
     def _linear_operator_(self, norb: int, nelec: tuple[int, int]) -> LinearOperator:
         """Return a SciPy LinearOperator representing the object."""
-        if np.iscomplexobj(self.two_body_tensor):
-            raise NotImplementedError(
-                "This Hamiltonian has a complex-valued two-body tensor. "
-                "LinearOperator support for complex two-body tensors is not yet "
-                "implemented. See https://github.com/qiskit-community/ffsim/issues/81."
-            )
-
         n_alpha, n_beta = nelec
         linkstr_index_a = gen_linkstr_index(range(norb), n_alpha)
         linkstr_index_b = 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
+        )
 
-        if np.iscomplexobj(self.one_body_tensor):
-            return _linear_operator_complex(
-                one_body_tensor=self.one_body_tensor,
-                two_body_tensor=self.two_body_tensor,
-                constant=self.constant,
-                norb=norb,
-                nelec=nelec,
-                link_index=link_index,
-            )
-        return _linear_operator_real(
-            one_body_tensor=self.one_body_tensor,
-            two_body_tensor=self.two_body_tensor,
-            constant=self.constant,
-            norb=norb,
-            nelec=nelec,
-            link_index=link_index,
+        def matvec(vec: np.ndarray):
+            result = self.constant * vec.astype(complex, copy=False)
+            result += contract_2e(two_body, vec, norb, nelec, link_index=link_index)
+            return result
+
+        dim_ = dim(norb, nelec)
+        return LinearOperator(
+            shape=(dim_, dim_), matvec=matvec, rmatvec=matvec, dtype=complex
         )
 
     def _diag_(self, norb: int, nelec: tuple[int, int]) -> np.ndarray:
@@ -203,57 +191,3 @@ def _approx_eq_(self, other, rtol: float, atol: float) -> bool:
                 return False
             return True
         return NotImplemented
-
-
-def _linear_operator_complex(
-    one_body_tensor: np.ndarray,
-    two_body_tensor: np.ndarray,
-    constant: float,
-    norb: int,
-    nelec: tuple[int, int],
-    link_index: tuple[np.ndarray, np.ndarray],
-):
-    two_body = absorb_h1e(one_body_tensor.real, two_body_tensor, norb, nelec, 0.5)
-    dim_ = dim(norb, nelec)
-
-    def matvec(vec: np.ndarray):
-        result = constant * vec.astype(complex, copy=False)
-        result += 1j * contract_1e(
-            one_body_tensor.imag, vec.real, norb, nelec, link_index=link_index
-        )
-        result -= contract_1e(
-            one_body_tensor.imag, vec.imag, norb, nelec, link_index=link_index
-        )
-        result += contract_2e(two_body, vec.real, norb, nelec, link_index=link_index)
-        result += 1j * contract_2e(
-            two_body, vec.imag, norb, nelec, link_index=link_index
-        )
-        return result
-
-    return LinearOperator(
-        shape=(dim_, dim_), matvec=matvec, rmatvec=matvec, dtype=complex
-    )
-
-
-def _linear_operator_real(
-    one_body_tensor: np.ndarray,
-    two_body_tensor: np.ndarray,
-    constant: float,
-    norb: int,
-    nelec: tuple[int, int],
-    link_index: tuple[np.ndarray, np.ndarray],
-):
-    two_body = absorb_h1e(one_body_tensor, two_body_tensor, norb, nelec, 0.5)
-    dim_ = dim(norb, nelec)
-
-    def matvec(vec: np.ndarray):
-        result = constant * vec.astype(complex, copy=False)
-        result += contract_2e(two_body, vec.real, norb, nelec, link_index=link_index)
-        result += 1j * contract_2e(
-            two_body, vec.imag, norb, nelec, link_index=link_index
-        )
-        return result
-
-    return LinearOperator(
-        shape=(dim_, dim_), matvec=matvec, rmatvec=matvec, dtype=complex
-    )

tests/python/hamiltonians/molecular_hamiltonian_test.py

@@ -106,8 +106,7 @@ def test_fermion_operator(norb: int, nelec: tuple[int, int]):
     rng = np.random.default_rng()
 
     one_body_tensor = ffsim.random.random_hermitian(norb, seed=rng)
-    # TODO remove dtype=float after adding support for complex
-    two_body_tensor = ffsim.random.random_two_body_tensor(norb, seed=rng, dtype=float)
+    two_body_tensor = ffsim.random.random_two_body_tensor(norb, seed=rng)
     constant = rng.standard_normal()
     mol_hamiltonian = ffsim.MolecularHamiltonian(
         one_body_tensor, two_body_tensor, constant=constant
@@ -132,8 +131,7 @@ def test_rotated():
 
     # generate a random molecular Hamiltonian
     one_body_tensor = ffsim.random.random_hermitian(norb, seed=rng)
-    # TODO remove dtype=float after adding support for complex
-    two_body_tensor = ffsim.random.random_two_body_tensor(norb, seed=rng, dtype=float)
+    two_body_tensor = ffsim.random.random_two_body_tensor(norb, seed=rng)
     constant = rng.standard_normal()
     mol_hamiltonian = ffsim.MolecularHamiltonian(
         one_body_tensor, two_body_tensor, constant=constant

tests/python/hamiltonians/single_factorized_hamiltonian_test.py

@@ -33,7 +33,6 @@ def test_linear_operator(norb: int, nelec: tuple[int, int], cholesky: bool):
     """Test linear operator."""
     rng = np.random.default_rng(2474)
 
-    # TODO remove dtype=float once complex two-body is supported
     mol_hamiltonian = ffsim.random.random_molecular_hamiltonian(
         norb, seed=rng, dtype=float
     )
@@ -58,7 +57,6 @@ def test_reduced_matrix_product_states(norb: int, nelec: tuple[int, int]):
     """Test computing reduced matrix on product states."""
     rng = np.random.default_rng(7869)
 
-    # TODO remove dtype=float once complex two-body is supported
     mol_hamiltonian = ffsim.random.random_molecular_hamiltonian(
         norb, seed=rng, dtype=float
     )
@@ -104,7 +102,6 @@ def test_expectation_product_state_slater_determinant(
     """Test computing expectation value on Slater determinant product state."""
     rng = np.random.default_rng(3400)
 
-    # TODO remove dtype=float once complex two-body is supported
     mol_hamiltonian = ffsim.random.random_molecular_hamiltonian(
         norb, seed=rng, dtype=float
     )
@@ -137,7 +134,6 @@ def test_expectation_product_state(norb: int, nelec: tuple[int, int]):
     """Test computing expectation value on product state."""
     rng = np.random.default_rng(6775)
 
-    # TODO remove dtype=float once complex two-body is supported
     mol_hamiltonian = ffsim.random.random_molecular_hamiltonian(
         norb, seed=rng, dtype=float
     )