Fix edge case in fermi_hubbard_1d (#142)

* fix edge case in 1D Fermi-Hubbard model

* add populate keys label

* rename to op_periodic_edge

* add defaultdict

* add defaultdict(float)

python/ffsim/operators/fermi_hubbard.py

@@ -10,6 +10,8 @@
 
 """Fermi-Hubbard model Hamiltonian."""
 
+from collections import defaultdict
+
 from ffsim._lib import FermionOperator
 from ffsim.operators.fermion_action import cre_a, cre_b, des_a, des_b
 
@@ -56,13 +58,13 @@ def fermi_hubbard_1d(
 
     .. _The Hubbard Model: https://doi.org/10.1146/annurev-conmatphys-031620-102024
     """
-    coeffs: dict[tuple[tuple[bool, bool, int], ...], complex] = {}
+    coeffs: dict[tuple[tuple[bool, bool, int], ...], complex] = defaultdict(float)
 
     for p in range(norb - 1 + periodic):
-        coeffs[(cre_a(p), des_a((p + 1) % norb))] = -tunneling
-        coeffs[(cre_b(p), des_b((p + 1) % norb))] = -tunneling
-        coeffs[(cre_a((p + 1) % norb), des_a(p))] = -tunneling
-        coeffs[(cre_b((p + 1) % norb), des_b(p))] = -tunneling
+        coeffs[(cre_a(p), des_a((p + 1) % norb))] -= tunneling
+        coeffs[(cre_b(p), des_b((p + 1) % norb))] -= tunneling
+        coeffs[(cre_a((p + 1) % norb), des_a(p))] -= tunneling
+        coeffs[(cre_b((p + 1) % norb), des_b(p))] -= tunneling
         if nearest_neighbor_interaction:
             coeffs[
                 (cre_a(p), des_a(p), cre_a((p + 1) % norb), des_a((p + 1) % norb))

tests/python/operators/fermi_hubbard_test.py

@@ -132,6 +132,39 @@ def test_fermi_hubbard_1d():
         },
     )
 
+    # periodic boundary conditions (edge case)
+    op_periodic_edge = fermi_hubbard_1d(
+        norb=2,
+        tunneling=1,
+        interaction=2,
+        chemical_potential=3,
+        nearest_neighbor_interaction=4,
+        periodic=True,
+    )
+    np.testing.assert_equal(
+        dict(op_periodic_edge),
+        {
+            (cre_a(0), des_a(1)): -2,
+            (cre_b(0), des_b(1)): -2,
+            (cre_a(1), des_a(0)): -2,
+            (cre_b(1), des_b(0)): -2,
+            (cre_a(0), des_a(0), cre_b(0), des_b(0)): 2,
+            (cre_a(1), des_a(1), cre_b(1), des_b(1)): 2,
+            (cre_a(0), des_a(0)): -3,
+            (cre_b(0), des_b(0)): -3,
+            (cre_a(1), des_a(1)): -3,
+            (cre_b(1), des_b(1)): -3,
+            (cre_a(0), des_a(0), cre_a(1), des_a(1)): 4,
+            (cre_a(0), des_a(0), cre_b(1), des_b(1)): 4,
+            (cre_b(0), des_b(0), cre_a(1), des_a(1)): 4,
+            (cre_b(0), des_b(0), cre_b(1), des_b(1)): 4,
+            (cre_a(1), des_a(1), cre_a(0), des_a(0)): 4,
+            (cre_a(1), des_a(1), cre_b(0), des_b(0)): 4,
+            (cre_b(1), des_b(1), cre_a(0), des_a(0)): 4,
+            (cre_b(1), des_b(1), cre_b(0), des_b(0)): 4,
+        },
+    )
+
 
 def test_non_interacting_fermi_hubbard_1d_eigenvalue():
     """Test ground-state eigenvalue of the non-interacting one-dimensional Fermi-Hubbard
@@ -149,6 +182,14 @@ def test_non_interacting_fermi_hubbard_1d_eigenvalue():
     eigs_periodic, _ = scipy.sparse.linalg.eigsh(ham_periodic, which="SA", k=1)
     np.testing.assert_allclose(eigs_periodic[0], -4.000000000000)
 
+    # periodic boundary conditions (edge case)
+    op_periodic_edge = fermi_hubbard_1d(
+        norb=2, tunneling=1, interaction=0, periodic=True
+    )
+    ham_periodic = ffsim.linear_operator(op_periodic_edge, norb=2, nelec=(1, 1))
+    eigs_periodic, _ = scipy.sparse.linalg.eigsh(ham_periodic, which="SA", k=1)
+    np.testing.assert_allclose(eigs_periodic[0], -4.000000000000)
+
 
 def test_fermi_hubbard_1d_with_interaction_eigenvalue():
     """Test ground-state eigenvalue of the one-dimensional Fermi-Hubbard model
@@ -166,6 +207,14 @@ def test_fermi_hubbard_1d_with_interaction_eigenvalue():
     eigs_periodic, _ = scipy.sparse.linalg.eigsh(ham_periodic, which="SA", k=1)
     np.testing.assert_allclose(eigs_periodic[0], -2.828427124746)
 
+    # periodic boundary conditions (edge case)
+    op_periodic_edge = fermi_hubbard_1d(
+        norb=2, tunneling=1, interaction=2, periodic=True
+    )
+    ham_periodic = ffsim.linear_operator(op_periodic_edge, norb=2, nelec=(1, 1))
+    eigs_periodic, _ = scipy.sparse.linalg.eigsh(ham_periodic, which="SA", k=1)
+    np.testing.assert_allclose(eigs_periodic[0], -3.123105625618)
+
 
 def test_fermi_hubbard_1d_with_chemical_potential_eigenvalue():
     """Test ground-state eigenvalue of the one-dimensional Fermi-Hubbard model
@@ -185,6 +234,14 @@ def test_fermi_hubbard_1d_with_chemical_potential_eigenvalue():
     eigs_periodic, _ = scipy.sparse.linalg.eigsh(ham_periodic, which="SA", k=1)
     np.testing.assert_allclose(eigs_periodic[0], -14.828427124746)
 
+    # periodic boundary conditions (edge case)
+    op_periodic_edge = fermi_hubbard_1d(
+        norb=2, tunneling=1, interaction=2, chemical_potential=3, periodic=True
+    )
+    ham_periodic = ffsim.linear_operator(op_periodic_edge, norb=2, nelec=(1, 1))
+    eigs_periodic, _ = scipy.sparse.linalg.eigsh(ham_periodic, which="SA", k=1)
+    np.testing.assert_allclose(eigs_periodic[0], -9.123105625618)
+
 
 def test_fermi_hubbard_1d_with_nearest_neighbor_interaction_eigenvalue():
     """Test ground-state eigenvalue of the one-dimensional Fermi-Hubbard model
@@ -216,6 +273,19 @@ def test_fermi_hubbard_1d_with_nearest_neighbor_interaction_eigenvalue():
     eigs_periodic, _ = scipy.sparse.linalg.eigsh(ham_periodic, which="SA", k=1)
     np.testing.assert_allclose(eigs_periodic[0], -8.781962448006)
 
+    # periodic boundary conditions (edge case)
+    op_periodic_edge = fermi_hubbard_1d(
+        norb=2,
+        tunneling=1,
+        interaction=2,
+        chemical_potential=3,
+        nearest_neighbor_interaction=4,
+        periodic=True,
+    )
+    ham_periodic = ffsim.linear_operator(op_periodic_edge, norb=2, nelec=(1, 1))
+    eigs_periodic, _ = scipy.sparse.linalg.eigsh(ham_periodic, which="SA", k=1)
+    np.testing.assert_allclose(eigs_periodic[0], -6.000000000000)
+
 
 def test_fermi_hubbard_1d_with_unequal_filling_eigenvalue():
     """Test ground-state eigenvalue of the one-dimensional Fermi-Hubbard model