From b2fb884029b1adc4bb301d52bdf722ad9a322f43 Mon Sep 17 00:00:00 2001 From: Bart Andrews <38833761+bartandrews@users.noreply.github.com> Date: Wed, 24 Apr 2024 11:48:48 -0700 Subject: [PATCH] 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 | 12 ++-- tests/python/operators/fermi_hubbard_test.py | 70 ++++++++++++++++++++ 2 files changed, 77 insertions(+), 5 deletions(-) diff --git a/python/ffsim/operators/fermi_hubbard.py b/python/ffsim/operators/fermi_hubbard.py index 5f5eee6eb..d185973e2 100644 --- a/python/ffsim/operators/fermi_hubbard.py +++ b/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)) diff --git a/tests/python/operators/fermi_hubbard_test.py b/tests/python/operators/fermi_hubbard_test.py index 36e04a297..0b9c39b1e 100644 --- a/tests/python/operators/fermi_hubbard_test.py +++ b/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