Add from_fermion_operator method to MoelcularHamiltonian class

python/ffsim/hamiltonians/diagonal_coulomb_hamiltonian.py

@@ -122,7 +122,24 @@ def _fermion_operator_(self) -> FermionOperator:
 
     @staticmethod
     def from_fermion_operator(op: FermionOperator) -> DiagonalCoulombHamiltonian:
-        """Convert a FermionOperator to a DiagonalCoulombHamiltonian."""
+        r"""Initialize a DiagonalCoulombHamiltonian from a FermionOperator.
+
+        The input operator must contain only terms of the following form:
+
+        - A real-valued constant
+        - :math:`a^\dagger_{\sigma, p} a_{\sigma, q}`
+        - :math:`n_{\sigma, p} n_{\tau, q}`
+
+        Any other terms will cause an error to be raised. No attempt will be made to
+        normal-order terms.
+
+        Args:
+            op: The FermionOperator from which to initialize the
+                DiagonalCoulombHamiltonian.
+
+        Returns:
+            The DiagonalCoulombHamiltonian represented by the input FermionOperator.
+        """
 
         # extract norb
         norb = 1 + max(orb for term in op for _, _, orb in term)
@@ -150,8 +167,8 @@ def from_fermion_operator(op: FermionOperator) -> DiagonalCoulombHamiltonian:
                 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}."
+                        f"DiagonalCoulombHamiltonian. The quadratic term {term} is not "
+                        r"of the form a^\dagger_{\sigma, p} a_{\sigma, q}."
                     )
             elif len(term) == 4:
                 # two-body term
@@ -171,8 +188,8 @@ def from_fermion_operator(op: FermionOperator) -> DiagonalCoulombHamiltonian:
                 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}."
+                        f"DiagonalCoulombHamiltonian. The quartic term {term} is not "
+                        r"of the form n_{\sigma, p} n_{\tau, q}."
                     )
             else:
                 raise ValueError(

python/ffsim/hamiltonians/molecular_hamiltonian.py

@@ -169,6 +169,83 @@ def _fermion_operator_(self) -> FermionOperator:
             )
         return op
 
+    @staticmethod
+    def from_fermion_operator(op: FermionOperator) -> MolecularHamiltonian:
+        r"""Initialize a MolecularHamiltonian from a FermionOperator.
+
+        The input operator must contain only terms of the following form:
+
+        - A real-valued constant
+        - :math:`a^\dagger_{\sigma, p} a_{\sigma, q}`
+        - :math:`a^\dagger_{\sigma,p}a^\dagger_{\tau,r}a_{\tau,s}a_{\sigma,q}`
+
+        Any other terms will cause an error to be raised. No attempt will be made to
+        normal-order terms.
+
+        Args:
+            op: The FermionOperator from which to initialize the MolecularHamiltonian.
+
+        Returns:
+            The MolecularHamiltonian represented by the input FermionOperator.
+        """
+        # extract number of spatial orbitals
+        norb = 1 + max(orb for term in op for _, _, orb in term)
+
+        # initialize constant, one‑ and two‑body tensors
+        constant: float = 0.0
+        one_body_tensor = np.zeros((norb, norb), dtype=complex)
+        two_body_tensor = np.zeros((norb, norb, norb, norb), dtype=complex)
+
+        for term, coeff in op.items():
+            # constant term: empty tuple
+            if not term:
+                if coeff.imag:
+                    raise ValueError(
+                        f"Constant term must be real. Instead, got {coeff}."
+                    )
+                constant = coeff.real
+            # one‑body term: a†_σ,p a_σ,q  (σ = α or β)
+            elif len(term) == 2:
+                (_, _, p), (_, _, q) = term
+                valid_one_body = [(cre_a(p), des_a(q)), (cre_b(p), des_b(q))]
+                if term in valid_one_body:
+                    one_body_tensor[p, q] += 0.5 * coeff
+                else:
+                    raise ValueError(
+                        "FermionOperator cannot be converted to MolecularHamiltonian. "
+                        f"The quadratic term {term} is not of the required form "
+                        r"a^\dagger_{\sigma, p} a_{\sigma, q}."
+                    )
+            # two‑body term: a†_σ,p a†_τ,r a_τ,s a_σ,q
+            elif len(term) == 4:
+                (_, _, p), (_, _, r), (_, _, s), (_, _, q) = term
+                valid_two_body = [
+                    (cre_a(p), cre_a(r), des_a(s), des_a(q)),
+                    (cre_a(p), cre_b(r), des_b(s), des_a(q)),
+                    (cre_b(p), cre_a(r), des_a(s), des_b(q)),
+                    (cre_b(p), cre_b(r), des_b(s), des_b(q)),
+                ]
+                if term not in valid_two_body:
+                    raise ValueError(
+                        "FermionOperator cannot be converted to MolecularHamiltonian. "
+                        f"The quartic term {term} is not of the required form "
+                        r"a^\dagger_{\sigma,p}a^\dagger_{\tau,r}a_{\tau,s}a_{\sigma,q}."
+                    )
+                two_body_tensor[p, q, r, s] += 0.5 * coeff
+            # other terms
+            else:
+                raise ValueError(
+                    "FermionOperator cannot be converted to MolecularHamiltonian."
+                    f" The term {term} is neither a constant, one-body, nor two-body "
+                    "term."
+                )
+
+        return MolecularHamiltonian(
+            one_body_tensor=one_body_tensor,
+            two_body_tensor=two_body_tensor,
+            constant=constant,
+        )
+
     def _approx_eq_(self, other, rtol: float, atol: float) -> bool:
         if isinstance(other, MolecularHamiltonian):
             if not np.allclose(self.constant, other.constant, rtol=rtol, atol=atol):

tests/python/hamiltonians/diagonal_coulomb_hamiltonian_test.py

@@ -118,7 +118,7 @@ def test_from_fermion_operator_errors():
 
     mol_hamiltonian = ffsim.random.random_molecular_hamiltonian(norb, seed=rng)
     op = ffsim.fermion_operator(mol_hamiltonian)
-    with pytest.raises(ValueError, match="two-body"):
+    with pytest.raises(ValueError, match="quartic"):
         _ = ffsim.DiagonalCoulombHamiltonian.from_fermion_operator(op)
 
     dc_hamiltonian = ffsim.random.random_diagonal_coulomb_hamiltonian(norb, seed=rng)

tests/python/hamiltonians/molecular_hamiltonian_test.py

@@ -330,3 +330,29 @@ def test_molecular_hamiltonian_spinless_rotated(norb: int, nelec: int):
     original_expectation = np.vdot(vec, linop @ vec)
     rotated_expectation = np.vdot(rotated_vec, linop_rotated @ rotated_vec)
     np.testing.assert_allclose(original_expectation, rotated_expectation)
+
+
+@pytest.mark.parametrize("norb", range(1, 5))
+def test_from_fermion_operator_roundtrip(norb: int):
+    """Test converting fermion operator to molecular Hamiltonian."""
+    mol_ham = ffsim.random.random_molecular_hamiltonian(norb=norb, seed=RNG)
+    roundtripped = ffsim.MolecularHamiltonian.from_fermion_operator(
+        ffsim.fermion_operator(mol_ham)
+    )
+    assert ffsim.approx_eq(roundtripped, mol_ham, atol=0)
+
+
+@pytest.mark.parametrize("norb", range(1, 5))
+def test_from_fermion_operator_invalid(norb: int):
+    """Test converting fermion operator with invalid terms."""
+    op = ffsim.FermionOperator({(ffsim.cre_a(3), ffsim.cre_b(2)): 1.0})
+    with pytest.raises(ValueError, match="quadratic"):
+        _ = ffsim.MolecularHamiltonian.from_fermion_operator(op)
+    op = ffsim.FermionOperator(
+        {(ffsim.cre_a(3), ffsim.cre_b(2), ffsim.cre_a(3), ffsim.cre_b(2)): 1.0}
+    )
+    with pytest.raises(ValueError, match="quartic"):
+        _ = ffsim.MolecularHamiltonian.from_fermion_operator(op)
+    op = ffsim.FermionOperator({(ffsim.cre_a(3),): 1.0})
+    with pytest.raises(ValueError, match="term"):
+        _ = ffsim.MolecularHamiltonian.from_fermion_operator(op)