Skip to content

Commit

Permalink
fix edge case in 1D Fermi-Hubbard model
Browse files Browse the repository at this point in the history
  • Loading branch information
bartandrews committed Apr 24, 2024
1 parent ed51b27 commit 7a4463b
Show file tree
Hide file tree
Showing 2 changed files with 77 additions and 4 deletions.
15 changes: 11 additions & 4 deletions python/ffsim/operators/fermi_hubbard.py
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down
66 changes: 66 additions & 0 deletions tests/python/operators/fermi_hubbard_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down

0 comments on commit 7a4463b

Please sign in to comment.