add _fermion_operator_ for df hamiltonian (#308)

* add _fermion_operator_ for df hamiltonian

* increase test tolerance

python/ffsim/hamiltonians/double_factorized_hamiltonian.py

@@ -20,6 +20,8 @@
 from ffsim.contract.num_op_sum import num_op_sum_linop
 from ffsim.hamiltonians.molecular_hamiltonian import MolecularHamiltonian
 from ffsim.linalg import double_factorized
+from ffsim.operators import FermionOperator
+from ffsim.protocols import fermion_operator
 from ffsim.states import dim
 
 
@@ -270,6 +272,10 @@ def _diag_(self, norb: int, nelec: tuple[int, int]) -> np.ndarray:
         """Return the diagonal entries of the Hamiltonian."""
         return self.to_molecular_hamiltonian()._diag_(norb, nelec)
 
+    def _fermion_operator_(self) -> FermionOperator:
+        """Return a FermionOperator representing the object."""
+        return fermion_operator(self.to_molecular_hamiltonian())
+
 
 def _df_z_representation(
     diag_coulomb_mats: np.ndarray, orbital_rotations: np.ndarray

python/ffsim/random/random.py

@@ -526,7 +526,9 @@ def random_double_factorized_hamiltonian(
     diag_coulomb_mats = np.stack(
         [random_real_symmetric_matrix(norb, seed=rng) for _ in range(rank)]
     )
-    orbital_rotations = np.stack([random_unitary(norb, seed=rng) for _ in range(rank)])
+    orbital_rotations = np.stack(
+        [random_orthogonal(norb, seed=rng) for _ in range(rank)]
+    )
     constant = rng.standard_normal()
     return hamiltonians.DoubleFactorizedHamiltonian(
         one_body_tensor=one_body_tensor,

tests/python/hamiltonians/double_factorized_hamiltonian_test.py

@@ -113,3 +113,28 @@ def test_diag():
     np.testing.assert_allclose(
         ffsim.diag(hamiltonian, norb=norb, nelec=nelec), np.diag(hamiltonian_dense)
     )
+
+
+@pytest.mark.parametrize(
+    "norb, nelec",
+    [
+        (4, (2, 2)),
+        (4, (1, 2)),
+        (4, (0, 2)),
+        (4, (0, 0)),
+    ],
+)
+def test_fermion_operator(norb: int, nelec: tuple[int, int]):
+    """Test FermionOperator."""
+    rng = np.random.default_rng()
+
+    df_hamiltonian = ffsim.random.random_double_factorized_hamiltonian(norb, seed=rng)
+    vec = ffsim.random.random_state_vector(ffsim.dim(norb, nelec), seed=rng)
+
+    op = ffsim.fermion_operator(df_hamiltonian)
+    linop = ffsim.linear_operator(op, norb, nelec)
+    expected_linop = ffsim.linear_operator(df_hamiltonian, norb, nelec)
+
+    actual = linop @ vec
+    expected = expected_linop @ vec
+    np.testing.assert_allclose(actual, expected)

tests/python/trotter/double_factorized_test.py

@@ -26,7 +26,7 @@
         (3, (1, 1), 0.1, 10, 1, True, 3e-3),
         (4, (2, 1), 0.1, 1, 2, False, 3e-3),
         (4, (1, 2), 0.1, 3, 2, True, 3e-3),
-        (4, (2, 2), 0.1, 4, 1, False, 3e-3),
+        (4, (2, 2), 0.1, 4, 1, False, 4e-3),
         (5, (3, 2), 0.1, 5, 1, True, 3e-3),
     ],
 )