diff --git a/CHANGELOG.md b/CHANGELOG.md index 4b119536..b0c7fab5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,7 +11,8 @@ This project adheres to [Semantic Versioning], with the exception that minor rel ### Added -- added deterministic ensemble evolution with optional autocorrelator and two-time correlator outputs, including periodic-wrap two-site observable support on `(L-1, 0)` ([#409]) ([**@Gauthameshwar**])\ +- added Fermionic and Jordan-Wigner MPO encodings of 1D Fermi-Hubbard model ([#220]) ([**@thilomueller**]) +- added deterministic ensemble evolution with optional autocorrelator and two-time correlator outputs, including periodic-wrap two-site observable support on `(L-1, 0)` ([#409]) ([**@Gauthameshwar**]) ### Changed @@ -100,7 +101,7 @@ _πŸ“š Refer to the [GitHub Release Notes](https://github.com/munich-quantum-tool -[Unreleased]: https://github.com/munich-quantum-toolkit/yaqs/compare/v0.4.0...HEAD +[Unreleased]: https://github.com/munich-quantum-toolkit/yaqs/compare/v0.5.0...HEAD [0.5.0]: https://github.com/munich-quantum-toolkit/yaqs/compare/v0.5.0 [0.4.0]: https://github.com/munich-quantum-toolkit/yaqs/releases/tag/v0.4.0 [0.3.3]: https://github.com/munich-quantum-toolkit/yaqs/releases/tag/v0.3.3 @@ -109,6 +110,7 @@ _πŸ“š Refer to the [GitHub Release Notes](https://github.com/munich-quantum-tool +[#220]: https://github.com/munich-quantum-toolkit/yaqs/pull/220 [#420]: https://github.com/munich-quantum-toolkit/yaqs/pull/420 [#409]: https://github.com/munich-quantum-toolkit/yaqs/pull/409 [#344]: https://github.com/munich-quantum-toolkit/yaqs/pull/344 diff --git a/docs/examples/fermi_hubbard_mpo.md b/docs/examples/fermi_hubbard_mpo.md new file mode 100644 index 00000000..8b8f424f --- /dev/null +++ b/docs/examples/fermi_hubbard_mpo.md @@ -0,0 +1,61 @@ +--- +file_format: mystnb +kernelspec: + name: python3 +mystnb: + number_source_lines: true + execution_timeout: 300 +--- + +# 1D Fermi-Hubbard MPO + +This example shows how to build a 1D Fermi-Hubbard Hamiltonian as an MPO using {class}`~mqt.yaqs.core.data_structures.networks.MPO.fermi_hubbard_1d`. + +YAQS supports two representations: + +- **Fermionic sites** (default): one site with local dimension 4 per physical lattice site. + Ladder operators act on a composite ↑/↓ basis per site; this is not a Jordan–Wigner qubit chain + across sites, but matches the standard tensor-product embedding of site Fock spaces. +- **Jordan-Wigner Pauli chain** (`jordan_wigner=True`): qubits in the order 1↑, 1↓, 2↑, 2↓, … with + local dimension 2 and full JW signs between spin orbitals. Use this mode for Pauli-string / + qubit simulators. + +The Hamiltonian (open boundaries, no chemical potential) is + +$$ +H = -t \sum_{i,\sigma} \left(c^\dagger_{i,\sigma} c_{i+1,\sigma} + \mathrm{h.c.}\right) ++ U \sum_i n_{i,\uparrow} n_{i,\downarrow}. +$$ + +## Fermionic MPO + +```{code-cell} ipython3 +from mqt.yaqs.core.data_structures.networks import MPO + +num_sites = 4 +t = 1.0 +u = 0.5 + +h_mpo = MPO.fermi_hubbard_1d(num_sites, t=t, u=u) +print(f"sites={h_mpo.length}, local dim={h_mpo.physical_dimension}, matrix shape={h_mpo.to_matrix().shape}") +``` + +The single-site basis is $|0\rangle, |\!\downarrow\rangle, |\!\uparrow\rangle, |\!\uparrow\downarrow\rangle$ (NumPy `kron` ordering for $|\!\uparrow\rangle \otimes |\!\downarrow\rangle$). + +## Jordan-Wigner MPO + +For the same model on $L$ physical sites, pass `length=2 * L` spin orbitals: + +```{code-cell} ipython3 +num_orbitals = 2 * num_sites + +h_jw = MPO.fermi_hubbard_1d(num_orbitals, t=t, u=u, jordan_wigner=True) +print(f"orbitals={h_jw.length}, local dim={h_jw.physical_dimension}, matrix shape={h_jw.to_matrix().shape}") +``` + +## Relation to the Trotter circuit helper + +{func}`~mqt.yaqs.core.libraries.circuit_library.create_1d_fermi_hubbard_circuit` builds a **digital** Trotter circuit on separate ↑ and ↓ registers and can include a chemical potential $\mu$. +The MPO factories above target the **analog** Hamiltonian without $\mu$ and use either fermionic operators or an interleaved JW layout. + +For digital simulation of the circuit model, use the circuit API; for tensor-network evolution of the Hubbard Hamiltonian, use `MPO.fermi_hubbard_1d`. diff --git a/docs/index.md b/docs/index.md index af439f55..caf9031c 100644 --- a/docs/index.md +++ b/docs/index.md @@ -39,6 +39,7 @@ examples/ensemble_evolution examples/solver_comparison examples/scheduled_jumps examples/transmon_emulation +examples/fermi_hubbard_mpo examples/process_tomography examples/strong_circuit_simulation examples/sample_observable_digital_tjm diff --git a/src/mqt/yaqs/core/data_structures/networks.py b/src/mqt/yaqs/core/data_structures/networks.py index cc4cbeec..91321f02 100644 --- a/src/mqt/yaqs/core/data_structures/networks.py +++ b/src/mqt/yaqs/core/data_structures/networks.py @@ -1309,6 +1309,7 @@ class MPO: - ``MPO.ising(...)`` / ``MPO.heisenberg(...)``: qubit Pauli Hamiltonians. - ``MPO.hamiltonian(...)``: generic one-/two-body Pauli interactions. + - ``MPO.fermi_hubbard_1d(...)``: 1D Fermi-Hubbard (fermionic or Jordan-Wigner Pauli). - ``MPO.coupled_transmon(...)``: alternating qubit/resonator chain MPO. - ``from_pauli_sum(...)``: in-place build from a sum of Pauli-string terms. - ``identity(...)``, ``custom(...)``, ``finite_state_machine(...)``: in-place builders. @@ -1508,6 +1509,146 @@ def heisenberg( n_sweeps=n_sweeps, ) + @classmethod + def fermi_hubbard_1d( + cls, + length: int, + t: float, + u: float, + *, + jordan_wigner: bool = False, + ) -> MPO: + r"""Construct a 1D Fermi-Hubbard Hamiltonian MPO. + + Without ``jordan_wigner``, builds the standard fermionic MPO on sites with + local dimension 4. The single-site basis is + :math:`|0\\rangle, |\\!\\downarrow\\rangle, |\\!\\uparrow\\rangle, |\\!\\uparrow\\downarrow\\rangle` + (NumPy ``kron`` ordering for :math:`|\\!\\uparrow\\rangle \\otimes |\\!\\downarrow\\rangle`). + The Hamiltonian is + :math:`H = -t \\sum_{i,\\sigma} (c^\\dagger_{i,\\sigma} c_{i+1,\\sigma} + \\mathrm{h.c.}) + + U \\sum_i n_{i,\\uparrow} n_{i,\\downarrow}`. + + With ``jordan_wigner=True``, builds the Jordan-Wigner Pauli-string MPO on an + interleaved spin chain 1↑, 1↓, 2↑, 2↓, ... (local dimension 2): + + .. math:: + + U n_{i,\\uparrow} n_{i,\\downarrow} + = \\frac{U}{4} \\left(I - Z_{i,\\uparrow} - Z_{i,\\downarrow} + + Z_{i,\\uparrow} Z_{i,\\downarrow}\\right) + + H = \\sum_i \\frac{U}{4} \\left(I - Z_{i,\\uparrow} - Z_{i,\\downarrow} + + Z_{i,\\uparrow} Z_{i,\\downarrow}\\right) + - \\frac{t}{2} \\sum_i \\left( X_{\\uparrow,i} Z_{\\downarrow,i} X_{\\uparrow,i+1} + + Y_{\\uparrow,i} Z_{\\downarrow,i} Y_{\\uparrow,i+1} \\right) + - \\frac{t}{2} \\sum_i \\left( X_{\\downarrow,i} Z_{\\uparrow,i+1} X_{\\downarrow,i+1} + + Y_{\\downarrow,i} Z_{\\uparrow,i+1} Y_{\\downarrow,i+1} \\right) + + Without ``jordan_wigner``, the MPO uses fermionic ladder operators on composite + dimension-4 sites (hard-core constraint per site). Inter-site algebra matches + that embedding; use ``jordan_wigner=True`` for a Pauli-chain representation + with full Jordan-Wigner signs between spin orbitals. + + In JW mode ``length`` is the number of **spin orbitals** and must be even and + at least 2. + + Args: + length: Chain length. Number of fermionic sites if ``jordan_wigner`` is + False; number of spin orbitals (even) if True. + t: Hopping strength. + u: On-site interaction strength. + jordan_wigner: If True, use the JW-transformed Pauli MPO; otherwise use + the fermionic operator MPO. + + Returns: + An MPO representing the 1D Fermi-Hubbard Hamiltonian. + + Raises: + ValueError: If ``length`` is invalid for the chosen representation. + """ + if jordan_wigner: + if length % 2 != 0 or length < 2: + msg = "length must be an even integer β‰₯ 2 (ordering: 1↑,1↓,2↑,2↓,...)." + raise ValueError(msg) + return cls._fermi_hubbard_1d_jordan_wigner(length=length, t=t, u=u) + return cls._fermi_hubbard_1d_fermionic(length=length, t=t, u=u) + + @classmethod + def _fermi_hubbard_1d_fermionic(cls, length: int, t: float, u: float) -> MPO: + if length <= 0: + msg = "length must be positive." + raise ValueError(msg) + + physical_dimension = 4 + identity = np.eye(physical_dimension, dtype=complex) + zero = np.zeros_like(identity, dtype=complex) + c = np.array([[0, 1], [0, 0]], dtype=complex) + c_dag = np.array([[0, 0], [1, 0]], dtype=complex) + c_up = np.kron(c, np.eye(2, dtype=complex)) + c_down = np.kron(np.eye(2, dtype=complex), c) + c_up_dag = np.kron(c_dag, np.eye(2, dtype=complex)) + c_down_dag = np.kron(np.eye(2, dtype=complex), c_dag) + n_up = np.kron(c_dag @ c, np.eye(2, dtype=complex)) + n_down = np.kron(np.eye(2, dtype=complex), c_dag @ c) + onsite = u * n_up @ n_down + + # Bond layout matches ``bose_hubbard``: channels + # 0=identity, 1=c↑†, 2=c↓†, 3=c↑, 4=c↓, 5=accumulator. + tensor = np.empty((6, 6, physical_dimension, physical_dimension), dtype=object) + tensor[:, :] = [[zero for _ in range(6)] for _ in range(6)] + tensor[0, 0] = identity + tensor[0, 1] = c_up_dag + tensor[0, 2] = c_down_dag + tensor[0, 3] = c_up + tensor[0, 4] = c_down + tensor[0, 5] = onsite + tensor[1, 5] = -t * c_up + tensor[2, 5] = -t * c_down + tensor[3, 5] = -t * c_up_dag + tensor[4, 5] = -t * c_down_dag + tensor[5, 5] = identity + + tensors = [np.transpose(tensor.copy(), (2, 3, 0, 1)).astype(np.complex128) for _ in range(length)] + tensors[0] = tensors[0][:, :, 0:1, :] + if length == 1: + tensors[0] = tensors[0][:, :, :, 5:6] + else: + tensors[-1] = tensors[-1][:, :, :, 5:6] + + mpo = cls() + mpo.tensors = tensors + mpo.length = length + mpo.physical_dimension = physical_dimension + assert mpo.check_if_valid_mpo(), "MPO initialized wrong" + return mpo + + @classmethod + def _fermi_hubbard_1d_jordan_wigner(cls, length: int, t: float, u: float) -> MPO: + num_sites = length // 2 + terms: list[tuple[complex | float, str]] = [] + for site in range(num_sites): + up, down = 2 * site, 2 * site + 1 + terms.extend([ + (u / 4, ""), + (-u / 4, f"Z{up}"), + (-u / 4, f"Z{down}"), + (u / 4, f"Z{up} Z{down}"), + ]) + for site in range(num_sites - 1): + up, down = 2 * site, 2 * site + 1 + up_next = 2 * (site + 1) + down_next = 2 * (site + 1) + 1 + terms.extend([ + (-t / 2, f"X{up} Z{down} X{up_next}"), + (-t / 2, f"Y{up} Z{down} Y{up_next}"), + (-t / 2, f"X{down} Z{up_next} X{down_next}"), + (-t / 2, f"Y{down} Z{up_next} Y{down_next}"), + ]) + + mpo = cls() + mpo.from_pauli_sum(terms=terms, length=length, n_sweeps=0) + return mpo + @classmethod def coupled_transmon( cls, @@ -1692,12 +1833,11 @@ def bose_hubbard( # build the full tensor list tensors = [np.transpose(tensor.copy(), (2, 3, 0, 1)).astype(np.complex128) for _ in range(length)] - - # Left boundary: take only row 0 - tensors[0] = np.transpose(tensor.copy(), (2, 3, 0, 1))[:, :, 0:1, :].astype(np.complex128) - - # Right boundary: take only col 3 - tensors[-1] = np.transpose(tensor.copy(), (2, 3, 0, 1))[:, :, :, 3:4].astype(np.complex128) + tensors[0] = tensors[0][:, :, 0:1, :] + if length == 1: + tensors[0] = tensors[0][:, :, :, 3:4] + else: + tensors[-1] = tensors[-1][:, :, :, 3:4] mpo = cls() mpo.tensors = tensors diff --git a/tests/core/data_structures/test_networks.py b/tests/core/data_structures/test_networks.py index d247b07c..da532333 100644 --- a/tests/core/data_structures/test_networks.py +++ b/tests/core/data_structures/test_networks.py @@ -19,6 +19,7 @@ from __future__ import annotations import copy +import re from typing import TYPE_CHECKING import numpy as np @@ -272,6 +273,117 @@ def embed(op_list: list[np.ndarray]) -> np.ndarray: return H +def _embed_local_ops(length: int, local_dim: int, site_ops: list[np.ndarray]) -> np.ndarray: + """Embed a list of single-site operators into the full Hilbert space. + + Returns: + np.ndarray: The embedded operator on the full chain Hilbert space. + """ + identity = np.eye(local_dim, dtype=complex) + op_list = [identity] * length + for site, op in enumerate(site_ops): + op_list[site] = op + out = np.array([[1.0]], dtype=complex) + for op in op_list: + out = np.kron(out, op) + return out + + +def _fermi_hubbard_1d_fermionic_dense(length: int, t: float, u: float) -> np.ndarray: + r"""Dense 1D Fermi-Hubbard Hamiltonian matching ``MPO.fermi_hubbard_1d``. + + Uses open boundaries and + :math:`H = -t \\sum_{i,\\sigma} (c^\\dagger_{i,\\sigma} c_{i+1,\\sigma} + \\mathrm{h.c.}) + + U \\sum_i n_{i,\\uparrow} n_{i,\\downarrow}` on local dimension-4 sites + (basis :math:`|0\\rangle, |\\!\\downarrow\\rangle, |\\!\\uparrow\\rangle, |\\!\\uparrow\\downarrow\\rangle`). + + Returns: + np.ndarray: Dense Hamiltonian matrix of shape ``(4**length, 4**length)``. + """ + local_dim = 4 + c = np.array([[0, 1], [0, 0]], dtype=complex) + c_dag = np.array([[0, 0], [1, 0]], dtype=complex) + identity2 = np.eye(2, dtype=complex) + c_up = np.kron(c, identity2) + c_down = np.kron(identity2, c) + c_up_dag = np.kron(c_dag, identity2) + c_down_dag = np.kron(identity2, c_dag) + n_up = c_up_dag @ c_up + n_down = c_down_dag @ c_down + identity4 = np.eye(local_dim, dtype=complex) + onsite = u * n_up @ n_down + + dim = local_dim**length + h = np.zeros((dim, dim), dtype=complex) + + for site in range(length): + site_ops = [identity4] * length + site_ops[site] = onsite + h += _embed_local_ops(length, local_dim, site_ops) + + for site in range(length - 1): + for c_right, c_left in ((c_up, c_up_dag), (c_down, c_down_dag)): + hop_ops = [identity4] * length + hop_ops[site] = c_left + hop_ops[site + 1] = c_right + h += -t * _embed_local_ops(length, local_dim, hop_ops) + + hop_ops = [identity4] * length + hop_ops[site] = c_right + hop_ops[site + 1] = c_left + h += -t * _embed_local_ops(length, local_dim, hop_ops) + + return h + + +def _fermi_hubbard_1d_jordan_wigner_dense(num_orbitals: int, t: float, u: float) -> np.ndarray: + """Dense JW-transformed Fermi-Hubbard on an interleaved spin chain. + + ``num_orbitals`` must be even. Orbitals are ordered 1↑, 1↓, 2↑, 2↓, ... + and the Hamiltonian matches the docstring of ``MPO.fermi_hubbard_1d(..., jordan_wigner=True)``. + + Returns: + np.ndarray: Dense Hamiltonian matrix of shape ``(2**num_orbitals, 2**num_orbitals)``. + + Raises: + ValueError: If ``num_orbitals`` is odd. + """ + if num_orbitals % 2 != 0: + msg = "num_orbitals must be even." + raise ValueError(msg) + + num_sites = num_orbitals // 2 + dim = 2**num_orbitals + h = np.zeros((dim, dim), dtype=complex) + + def term(ops: list[tuple[int, np.ndarray]]) -> np.ndarray: + local = [_I2] * num_orbitals + for index, op in ops: + local[index] = op + return _embed_local_ops(num_orbitals, 2, local) + + for site in range(num_sites): + up = 2 * site + down = 2 * site + 1 + h += (u / 4) * term([(up, _I2), (down, _I2)]) + h += -(u / 4) * term([(up, _Z2)]) + h += -(u / 4) * term([(down, _Z2)]) + h += (u / 4) * term([(up, _Z2), (down, _Z2)]) + + for site in range(num_sites - 1): + up = 2 * site + down = 2 * site + 1 + up_next = 2 * (site + 1) + h += -(t / 2) * term([(up, _X2), (down, _Z2), (up_next, _X2)]) + h += -(t / 2) * term([(up, _Y2), (down, _Z2), (up_next, _Y2)]) + + down_next = 2 * (site + 1) + 1 + h += -(t / 2) * term([(down, _X2), (up_next, _Z2), (down_next, _X2)]) + h += -(t / 2) * term([(down, _Y2), (up_next, _Z2), (down_next, _Y2)]) + + return h + + def untranspose_block(mpo_tensor: NDArray[np.complex128]) -> NDArray[np.complex128]: """Reverse the transposition of an MPO tensor. @@ -409,6 +521,69 @@ def test_bose_hubbard_correct_operator() -> None: np.testing.assert_allclose(H_mpo, H_dense, atol=1e-8) +def test_fermi_hubbard_1d_correct_operator() -> None: + """Verify the fermionic 1D Fermi-Hubbard MPO matches the dense Hamiltonian.""" + length = 3 + u, t = 0.5, 1.0 + + mpo = MPO.fermi_hubbard_1d(length, t, u) + + assert mpo.length == length + assert mpo.physical_dimension == 4 + assert len(mpo.tensors) == length + assert all(tensor.shape[2] <= 6 and tensor.shape[3] <= 6 for tensor in mpo.tensors) + + h_dense = _fermi_hubbard_1d_fermionic_dense(length, t, u) + np.testing.assert_allclose(mpo.to_matrix(), h_dense, atol=1e-10) + + +def test_fermi_hubbard_1d_jordan_wigner_correct_operator() -> None: + """Verify the JW 1D Fermi-Hubbard MPO matches the dense Pauli Hamiltonian.""" + num_orbitals = 4 + u, t = 0.5, 1.0 + + mpo = MPO.fermi_hubbard_1d(num_orbitals, t, u, jordan_wigner=True) + + assert mpo.length == num_orbitals + assert mpo.physical_dimension == 2 + assert len(mpo.tensors) == num_orbitals + assert all(tensor.shape[2] <= 16 and tensor.shape[3] <= 16 for tensor in mpo.tensors) + + h_dense = _fermi_hubbard_1d_jordan_wigner_dense(num_orbitals, t, u) + np.testing.assert_allclose(mpo.to_matrix(), h_dense, atol=1e-10) + + with pytest.raises(ValueError, match=re.escape("length must be an even integer β‰₯ 2 (ordering: 1↑,1↓,2↑,2↓,...).")): + MPO.fermi_hubbard_1d(length=5, t=t, u=u, jordan_wigner=True) + + +def test_fermi_hubbard_1d_length_one() -> None: + """Verify a single fermionic site MPO matches the dense reference.""" + length = 1 + u, t = 0.5, 1.0 + + mpo = MPO.fermi_hubbard_1d(length, t, u) + + assert mpo.length == length + assert mpo.physical_dimension == 4 + assert mpo.tensors[0].shape == (4, 4, 1, 1) + + h_dense = _fermi_hubbard_1d_fermionic_dense(length, t, u) + np.testing.assert_allclose(mpo.to_matrix(), h_dense, atol=1e-10) + + +def test_fermi_hubbard_1d_cross_representation() -> None: + """Onsite terms agree between fermionic and JW MPOs under the site basis map. + + Hopping terms differ between representations (composite fermionic sites vs JW + qubit chain), so this test uses ``t=0`` to compare only the interaction part. + """ + u = 0.5 + for length in (1, 2, 3): + h_ferm = MPO.fermi_hubbard_1d(length, t=0.0, u=u).to_matrix() + h_jw = MPO.fermi_hubbard_1d(2 * length, t=0.0, u=u, jordan_wigner=True).to_matrix() + np.testing.assert_allclose(h_ferm, h_jw, atol=1e-10) + + def test_identity() -> None: """Test that identity initializes an identity MPO correctly.