support generating random complex df hamiltonian (#320)

python/ffsim/hamiltonians/double_factorized_hamiltonian.py

@@ -222,12 +222,12 @@ def to_molecular_hamiltonian(self):
         """Convert the DoubleFactorizedHamiltonian to a MolecularHamiltonian."""
         df_hamiltonian = self.to_number_representation()
         two_body_tensor = np.einsum(
-            "kpi,kqi,kij,krj,ksj->pqrs",
-            df_hamiltonian.orbital_rotations,
-            df_hamiltonian.orbital_rotations,
+            "kij,kpi,kqi,krj,ksj->pqrs",
             df_hamiltonian.diag_coulomb_mats,
             df_hamiltonian.orbital_rotations,
+            df_hamiltonian.orbital_rotations.conj(),
             df_hamiltonian.orbital_rotations,
+            df_hamiltonian.orbital_rotations.conj(),
         )
         one_body_tensor = df_hamiltonian.one_body_tensor + 0.5 * np.einsum(
             "prqr", two_body_tensor

python/ffsim/random/random.py

@@ -504,7 +504,12 @@ def random_diagonal_coulomb_hamiltonian(
 
 
 def random_double_factorized_hamiltonian(
-    norb: int, *, rank: int | None = None, z_representation: bool = False, seed=None
+    norb: int,
+    *,
+    rank: int | None = None,
+    z_representation: bool = False,
+    real: bool = False,
+    seed=None,
 ) -> hamiltonians.DoubleFactorizedHamiltonian:
     """Sample a random double-factorized Hamiltonian.
 
@@ -513,6 +518,7 @@ def random_double_factorized_hamiltonian(
         rank: The desired number of terms in the two-body part of the Hamiltonian.
             If not specified, it will be set to ``norb * (norb + 1) // 2``.
         z_representation: Whether to return a Hamiltonian in the "Z" representation.
+        real: Whether to sample a real-valued object rather than a complex-valued one.
         seed: A seed to initialize the pseudorandom number generator.
             Should be a valid input to ``np.random.default_rng``.
 
@@ -522,13 +528,19 @@ def random_double_factorized_hamiltonian(
     if rank is None:
         rank = norb * (norb + 1) // 2
     rng = np.random.default_rng(seed)
-    one_body_tensor = random_hermitian(norb, seed=rng)
+    if real:
+        one_body_tensor = random_real_symmetric_matrix(norb, seed=rng)
+        orbital_rotations = np.stack(
+            [random_orthogonal(norb, seed=rng) for _ in range(rank)]
+        )
+    else:
+        one_body_tensor = random_hermitian(norb, seed=rng)
+        orbital_rotations = np.stack(
+            [random_unitary(norb, seed=rng) for _ in range(rank)]
+        )
     diag_coulomb_mats = np.stack(
         [random_real_symmetric_matrix(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

@@ -12,6 +12,8 @@
 
 from __future__ import annotations
 
+import itertools
+
 import numpy as np
 import pytest
 
@@ -80,8 +82,8 @@ def test_z_representation_round_trip():
 
 
 @pytest.mark.parametrize("z_representation", [False, True])
-def test_to_molecular_hamiltonian(z_representation: bool):
-    """Test converting to molecular Hamiltonian"""
+def test_to_molecular_hamiltonian_roundtrip(z_representation: bool):
+    """Test converting to molecular Hamiltonian."""
     rng = np.random.default_rng(2787)
     norb = 5
     hamiltonian = ffsim.random.random_molecular_hamiltonian(norb, seed=rng, dtype=float)
@@ -99,6 +101,23 @@ def test_to_molecular_hamiltonian(z_representation: bool):
     np.testing.assert_allclose(hamiltonian_roundtrip.constant, hamiltonian.constant)
 
 
+@pytest.mark.parametrize("z_representation", [False, True])
+def test_to_molecular_hamiltonian_symmetry(z_representation: bool):
+    """Test molecular Hamiltonian symmetry."""
+    rng = np.random.default_rng(2787)
+    norb = 5
+    df_hamiltonian = ffsim.random.random_double_factorized_hamiltonian(
+        norb, z_representation=z_representation, seed=rng
+    )
+    mol_hamiltonian = df_hamiltonian.to_molecular_hamiltonian()
+    two_body_tensor = mol_hamiltonian.two_body_tensor
+    for i, j, k, ell in itertools.product(range(norb), repeat=4):
+        val = two_body_tensor[i, j, k, ell]
+        np.testing.assert_allclose(two_body_tensor[k, ell, i, j], val)
+        np.testing.assert_allclose(two_body_tensor[j, i, ell, k], val.conjugate())
+        np.testing.assert_allclose(two_body_tensor[ell, k, j, i], val.conjugate())
+
+
 def test_diag():
     """Test computing diagonal."""
     rng = np.random.default_rng(2222)
@@ -128,7 +147,9 @@ 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)
+    df_hamiltonian = ffsim.random.random_double_factorized_hamiltonian(
+        norb, real=True, seed=rng
+    )
     vec = ffsim.random.random_state_vector(ffsim.dim(norb, nelec), seed=rng)
 
     op = ffsim.fermion_operator(df_hamiltonian)

tests/python/trotter/double_factorized_test.py

@@ -27,7 +27,7 @@
         (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, 4e-3),
-        (5, (3, 2), 0.1, 5, 1, True, 3e-3),
+        (5, (3, 2), 0.1, 10, 1, True, 3e-3),
     ],
 )
 def test_random(