From 7a4463ba7487cb2392fd2a2b9e344b7c581106c7 Mon Sep 17 00:00:00 2001 From: bartandrews Date: Wed, 24 Apr 2024 08:31:00 -0700 Subject: [PATCH 1/5] fix edge case in 1D Fermi-Hubbard model --- python/ffsim/operators/fermi_hubbard.py | 15 +++-- tests/python/operators/fermi_hubbard_test.py | 66 ++++++++++++++++++++ 2 files changed, 77 insertions(+), 4 deletions(-) diff --git a/python/ffsim/operators/fermi_hubbard.py b/python/ffsim/operators/fermi_hubbard.py index 5f5eee6eb..d3f23b0bc 100644 --- a/python/ffsim/operators/fermi_hubbard.py +++ b/python/ffsim/operators/fermi_hubbard.py @@ -58,11 +58,18 @@ def fermi_hubbard_1d( """ coeffs: dict[tuple[tuple[bool, bool, int], ...], complex] = {} + # initialize tunneling keys 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))] = 0 + coeffs[(cre_b(p), des_b((p + 1) % norb))] = 0 + coeffs[(cre_a((p + 1) % norb), des_a(p))] = 0 + coeffs[(cre_b((p + 1) % norb), des_b(p))] = 0 + + 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 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..133af05f4 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,12 @@ 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 = fermi_hubbard_1d(norb=2, tunneling=1, interaction=0, periodic=True) + ham_periodic = ffsim.linear_operator(op_periodic, 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 +205,12 @@ 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 = fermi_hubbard_1d(norb=2, tunneling=1, interaction=2, periodic=True) + ham_periodic = ffsim.linear_operator(op_periodic, 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 +230,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 = fermi_hubbard_1d( + norb=2, tunneling=1, interaction=2, chemical_potential=3, periodic=True + ) + ham_periodic = ffsim.linear_operator(op_periodic, 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 +269,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 = 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, 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 From ed5993c0594f2d2c815ef10efb400ab6e134fd48 Mon Sep 17 00:00:00 2001 From: bartandrews Date: Wed, 24 Apr 2024 08:46:35 -0700 Subject: [PATCH 2/5] add populate keys label --- python/ffsim/operators/fermi_hubbard.py | 1 + 1 file changed, 1 insertion(+) diff --git a/python/ffsim/operators/fermi_hubbard.py b/python/ffsim/operators/fermi_hubbard.py index d3f23b0bc..00a107dc9 100644 --- a/python/ffsim/operators/fermi_hubbard.py +++ b/python/ffsim/operators/fermi_hubbard.py @@ -65,6 +65,7 @@ def fermi_hubbard_1d( coeffs[(cre_a((p + 1) % norb), des_a(p))] = 0 coeffs[(cre_b((p + 1) % norb), des_b(p))] = 0 + # populate keys 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 From d3b66cfa674e845f392826732444558090af8c33 Mon Sep 17 00:00:00 2001 From: bartandrews Date: Wed, 24 Apr 2024 09:21:27 -0700 Subject: [PATCH 3/5] rename to op_periodic_edge --- tests/python/operators/fermi_hubbard_test.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/tests/python/operators/fermi_hubbard_test.py b/tests/python/operators/fermi_hubbard_test.py index 133af05f4..0b9c39b1e 100644 --- a/tests/python/operators/fermi_hubbard_test.py +++ b/tests/python/operators/fermi_hubbard_test.py @@ -183,8 +183,10 @@ def test_non_interacting_fermi_hubbard_1d_eigenvalue(): np.testing.assert_allclose(eigs_periodic[0], -4.000000000000) # periodic boundary conditions (edge case) - op_periodic = fermi_hubbard_1d(norb=2, tunneling=1, interaction=0, periodic=True) - ham_periodic = ffsim.linear_operator(op_periodic, norb=2, nelec=(1, 1)) + 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) @@ -206,8 +208,10 @@ def test_fermi_hubbard_1d_with_interaction_eigenvalue(): np.testing.assert_allclose(eigs_periodic[0], -2.828427124746) # periodic boundary conditions (edge case) - op_periodic = fermi_hubbard_1d(norb=2, tunneling=1, interaction=2, periodic=True) - ham_periodic = ffsim.linear_operator(op_periodic, norb=2, nelec=(1, 1)) + 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) @@ -231,10 +235,10 @@ def test_fermi_hubbard_1d_with_chemical_potential_eigenvalue(): np.testing.assert_allclose(eigs_periodic[0], -14.828427124746) # periodic boundary conditions (edge case) - op_periodic = fermi_hubbard_1d( + op_periodic_edge = fermi_hubbard_1d( norb=2, tunneling=1, interaction=2, chemical_potential=3, periodic=True ) - ham_periodic = ffsim.linear_operator(op_periodic, norb=2, nelec=(1, 1)) + 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) @@ -270,7 +274,7 @@ def test_fermi_hubbard_1d_with_nearest_neighbor_interaction_eigenvalue(): np.testing.assert_allclose(eigs_periodic[0], -8.781962448006) # periodic boundary conditions (edge case) - op_periodic = fermi_hubbard_1d( + op_periodic_edge = fermi_hubbard_1d( norb=2, tunneling=1, interaction=2, @@ -278,7 +282,7 @@ def test_fermi_hubbard_1d_with_nearest_neighbor_interaction_eigenvalue(): nearest_neighbor_interaction=4, periodic=True, ) - ham_periodic = ffsim.linear_operator(op_periodic, norb=2, nelec=(1, 1)) + 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) From 3d54e2e1e23f8b46fc5d423274418d33c2e40671 Mon Sep 17 00:00:00 2001 From: bartandrews Date: Wed, 24 Apr 2024 09:57:21 -0700 Subject: [PATCH 4/5] add defaultdict --- python/ffsim/operators/fermi_hubbard.py | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/python/ffsim/operators/fermi_hubbard.py b/python/ffsim/operators/fermi_hubbard.py index 00a107dc9..c55476b49 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,16 +58,8 @@ def fermi_hubbard_1d( .. _The Hubbard Model: https://doi.org/10.1146/annurev-conmatphys-031620-102024 """ - coeffs: dict[tuple[tuple[bool, bool, int], ...], complex] = {} - - # initialize tunneling keys - for p in range(norb - 1 + periodic): - coeffs[(cre_a(p), des_a((p + 1) % norb))] = 0 - coeffs[(cre_b(p), des_b((p + 1) % norb))] = 0 - coeffs[(cre_a((p + 1) % norb), des_a(p))] = 0 - coeffs[(cre_b((p + 1) % norb), des_b(p))] = 0 + coeffs: dict[tuple[tuple[bool, bool, int], ...], complex] = defaultdict(lambda: 0) - # populate keys 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 From e4c8c995248bd5f458547b06eea26e2f305638a0 Mon Sep 17 00:00:00 2001 From: bartandrews Date: Wed, 24 Apr 2024 11:03:54 -0700 Subject: [PATCH 5/5] add defaultdict(float) --- python/ffsim/operators/fermi_hubbard.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/ffsim/operators/fermi_hubbard.py b/python/ffsim/operators/fermi_hubbard.py index c55476b49..d185973e2 100644 --- a/python/ffsim/operators/fermi_hubbard.py +++ b/python/ffsim/operators/fermi_hubbard.py @@ -58,7 +58,7 @@ def fermi_hubbard_1d( .. _The Hubbard Model: https://doi.org/10.1146/annurev-conmatphys-031620-102024 """ - coeffs: dict[tuple[tuple[bool, bool, int], ...], complex] = defaultdict(lambda: 0) + 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