From 5e2f2b60a65222a21bacbb465c87b3992fd5d6d9 Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Sat, 7 Dec 2024 14:17:27 +0400
Subject: [PATCH 01/14] add `Heisenberg`

---
 src/qibo/hamiltonians/__init__.py |   2 +-
 src/qibo/hamiltonians/models.py   | 267 ++++++++++++++++++++----------
 2 files changed, 183 insertions(+), 86 deletions(-)

diff --git a/src/qibo/hamiltonians/__init__.py b/src/qibo/hamiltonians/__init__.py
index 65bb0891a4..9826e01833 100644
--- a/src/qibo/hamiltonians/__init__.py
+++ b/src/qibo/hamiltonians/__init__.py
@@ -1,2 +1,2 @@
 from qibo.hamiltonians.hamiltonians import *
-from qibo.hamiltonians.models import TFIM, XXZ, MaxCut, X, Y, Z
+from qibo.hamiltonians.models import Heisenberg, TFIM, XXZ, MaxCut, X, Y, Z
diff --git a/src/qibo/hamiltonians/models.py b/src/qibo/hamiltonians/models.py
index a268c26da2..30c22a491e 100644
--- a/src/qibo/hamiltonians/models.py
+++ b/src/qibo/hamiltonians/models.py
@@ -1,94 +1,15 @@
+from typing import Union, List, Tuple
+
 from functools import reduce
 
+import numpy as np
+
 from qibo.backends import matrices
 from qibo.config import raise_error
 from qibo.hamiltonians.hamiltonians import Hamiltonian, SymbolicHamiltonian
 from qibo.hamiltonians.terms import HamiltonianTerm
 
 
-def _multikron(matrix_list):
-    """Calculates Kronecker product of a list of matrices.
-
-    Args:
-        matrix_list (list): List of matrices as ``ndarray``.
-
-    Returns:
-        ndarray: Kronecker product of all matrices in ``matrix_list``.
-    """
-    import numpy as np
-
-    return reduce(np.kron, matrix_list)
-
-
-def _build_spin_model(nqubits, matrix, condition):
-    """Helper method for building nearest-neighbor spin model Hamiltonians."""
-    h = sum(
-        _multikron(matrix if condition(i, j) else matrices.I for j in range(nqubits))
-        for i in range(nqubits)
-    )
-    return h
-
-
-def XXZ(nqubits, delta=0.5, dense: bool = True, backend=None):
-    """Heisenberg :math:`\\mathrm{XXZ}` model with periodic boundary conditions.
-
-    .. math::
-        H = \\sum _{k=0}^N \\, \\left( X_{k} \\, X_{k + 1} + Y_{k} \\, Y_{k + 1}
-            + \\delta Z_{k} \\, Z_{k + 1} \\right) \\, .
-
-    Args:
-        nqubits (int): number of qubits.
-        delta (float, optional): coefficient for the :math:`Z` component.
-            Defaults to :math:`0.5`.
-        dense (bool, optional): If ``True``, creates the Hamiltonian as a
-            :class:`qibo.core.hamiltonians.Hamiltonian`, otherwise it creates
-            a :class:`qibo.core.hamiltonians.SymbolicHamiltonian`.
-            Defaults to ``True``.
-        backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be used
-            in the execution. If ``None``, it uses the current backend.
-            Defaults to ``None``.
-
-    Example:
-        .. testcode::
-
-            from qibo.hamiltonians import XXZ
-            h = XXZ(3) # initialized XXZ model with 3 qubits
-    """
-    if nqubits < 2:
-        raise_error(ValueError, "Number of qubits must be larger than one.")
-    if dense:
-        condition = lambda i, j: i in {j % nqubits, (j + 1) % nqubits}
-        hx = _build_spin_model(nqubits, matrices.X, condition)
-        hy = _build_spin_model(nqubits, matrices.Y, condition)
-        hz = _build_spin_model(nqubits, matrices.Z, condition)
-        matrix = hx + hy + delta * hz
-        return Hamiltonian(nqubits, matrix, backend=backend)
-
-    hx = _multikron([matrices.X, matrices.X])
-    hy = _multikron([matrices.Y, matrices.Y])
-    hz = _multikron([matrices.Z, matrices.Z])
-    matrix = hx + hy + delta * hz
-    terms = [HamiltonianTerm(matrix, i, i + 1) for i in range(nqubits - 1)]
-    terms.append(HamiltonianTerm(matrix, nqubits - 1, 0))
-    ham = SymbolicHamiltonian(backend=backend)
-    ham.terms = terms
-    return ham
-
-
-def _OneBodyPauli(nqubits, matrix, dense: bool = True, backend=None):
-    """Helper method for constracting non-interacting :math:`X`, :math:`Y`, and :math:`Z` Hamiltonians."""
-    if dense:
-        condition = lambda i, j: i == j % nqubits
-        ham = -_build_spin_model(nqubits, matrix, condition)
-        return Hamiltonian(nqubits, ham, backend=backend)
-
-    matrix = -matrix
-    terms = [HamiltonianTerm(matrix, i) for i in range(nqubits)]
-    ham = SymbolicHamiltonian(backend=backend)
-    ham.terms = terms
-    return ham
-
-
 def X(nqubits, dense: bool = True, backend=None):
     """Non-interacting Pauli-:math:`X` Hamiltonian.
 
@@ -194,7 +115,6 @@ def MaxCut(nqubits, dense: bool = True, backend=None):
             Defaults to ``None``.
     """
     import sympy as sp
-    from numpy import ones
 
     Z = sp.symbols(f"Z:{nqubits}")
     V = sp.symbols(f"V:{nqubits**2}")
@@ -205,7 +125,7 @@ def MaxCut(nqubits, dense: bool = True, backend=None):
     )
     sham /= 2
 
-    v = ones(nqubits**2)
+    v = np.ones(nqubits**2)
     smap = {s: (i, matrices.Z) for i, s in enumerate(Z)}
     smap.update({s: (i, v[i]) for i, s in enumerate(V)})
 
@@ -213,3 +133,180 @@ def MaxCut(nqubits, dense: bool = True, backend=None):
     if dense:
         return ham.dense
     return ham
+
+
+def Heisenberg(
+    nqubits,
+    coupling_constants: Union[float, int, List[int], Tuple[int, ...]],
+    external_field_strength: Union[float, int],
+    dense: bool = True,
+    backend=None,
+):
+    """Heisenberg model on a :math:`1`-dimensional periodic lattice.
+
+    The general :math:`n`-qubit Hamiltonian is given by
+
+    .. math::
+        H = -\\sum_{k = 1}^{n} \\, \\left(
+            J_{x} \\, X_{k} \\, X_{k + 1}
+            + J_{y} \\, Y_{k} \\, Y_{k + 1}
+            + J_{z} \\, Z_{k} \\, Z_{k + 1} \\right)
+            - h \\, \\sum_{k = 1}^{n} \\left(X_{k} + Y_{k} + Z_{k}\\right) \\, ,
+
+    where :math:`\\{J_{x}, J_{y}, J_{z}\\}` are called the ``coupling constants``,
+    :math:`h` is called the ``external field strength``, and :math:`\\{X, Y, Z\\}`
+    are the usual Pauli operators.
+
+    Args:
+        nqubits (int): number of qubits.
+        coupling_constants (float or int or list or tuple): list or tuple with the
+            three coupling constants :math:`\\{J_{x}, J_{y}, J{z}\\}`.
+            If ``int`` or ``float``, then :math:`J_{x} = J_{y} = J_{z}`.
+        external_field_strength (float or int): external magnetic field strength :math:`h`.
+        dense (bool, optional): If ``True``, creates the Hamiltonian as a
+            :class:`qibo.core.hamiltonians.Hamiltonian`, otherwise it creates
+            a :class:`qibo.core.hamiltonians.SymbolicHamiltonian`.
+            Defaults to ``True``.
+        backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be used
+            in the execution. If ``None``, it uses the current backend.
+            Defaults to ``None``.
+
+    Returns:
+        :class:`qibo.hamiltonians.Hamiltonian` or :class:`qibo.hamiltonians.SymbolicHamiltonian`: 
+        Heisenberg Hamiltonian.
+    """
+    if isinstance(coupling_constants, (list, tuple)) and len(coupling_constants) != 3:
+        raise_error(
+            ValueError,
+            f"When `coupling_constants` is type `int` or `list`, it must have length 3.",
+        )
+    elif isinstance(coupling_constants, (float, int)):
+        coupling_constants = [coupling_constants] * 3
+
+    if dense:
+        condition = lambda i, j: i in {j % nqubits, (j + 1) % nqubits}
+        hx = _build_spin_model(nqubits, matrices.X, condition)
+        hy = _build_spin_model(nqubits, matrices.Y, condition)
+        hz = _build_spin_model(nqubits, matrices.Z, condition)
+        matrix = (
+            -coupling_constants[0] * hx
+            - coupling_constants[1] * hy
+            - coupling_constants[2] * hz
+        )
+
+        for pauli in [matrices.X, matrices.Y, matrices.Z]:
+            matrix += (
+                external_field_strength
+                * _OneBodyPauli(nqubits, pauli, dense, backend).matrix
+            )
+
+        return Hamiltonian(nqubits, matrix, backend=backend)
+
+    hx = _multikron([matrices.X, matrices.X])
+    hy = _multikron([matrices.Y, matrices.Y])
+    hz = _multikron([matrices.Z, matrices.Z])
+
+    matrix = (
+        -coupling_constants[0] * hx
+        - coupling_constants[1] * hy
+        - coupling_constants[2] * hz
+    )
+
+    terms = [HamiltonianTerm(matrix, i, i + 1) for i in range(nqubits - 1)]
+    terms.append(HamiltonianTerm(matrix, nqubits - 1, 0))
+
+    if external_field_strength != 0.0:
+        terms.extend(
+            [
+                -external_field_strength * HamiltonianTerm(pauli, qubit)
+                for qubit in range(nqubits)
+                for pauli in [matrices.X, matrices.Y, matrices.Z]
+            ]
+        )
+
+    ham = SymbolicHamiltonian(backend=backend)
+    ham.terms = terms
+
+    return ham
+
+
+def XXZ(nqubits, delta=0.5, dense: bool = True, backend=None):
+    """Heisenberg :math:`\\mathrm{XXZ}` model with periodic boundary conditions.
+
+    .. math::
+        H = \\sum _{k=0}^N \\, \\left( X_{k} \\, X_{k + 1} + Y_{k} \\, Y_{k + 1}
+            + \\delta Z_{k} \\, Z_{k + 1} \\right) \\, .
+
+    Args:
+        nqubits (int): number of qubits.
+        delta (float, optional): coefficient for the :math:`Z` component.
+            Defaults to :math:`0.5`.
+        dense (bool, optional): If ``True``, creates the Hamiltonian as a
+            :class:`qibo.core.hamiltonians.Hamiltonian`, otherwise it creates
+            a :class:`qibo.core.hamiltonians.SymbolicHamiltonian`.
+            Defaults to ``True``.
+        backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be used
+            in the execution. If ``None``, it uses the current backend.
+            Defaults to ``None``.
+
+    Example:
+        .. testcode::
+
+            from qibo.hamiltonians import XXZ
+            h = XXZ(3) # initialized XXZ model with 3 qubits
+    """
+    if nqubits < 2:
+        raise_error(ValueError, "Number of qubits must be larger than one.")
+    if dense:
+        condition = lambda i, j: i in {j % nqubits, (j + 1) % nqubits}
+        hx = _build_spin_model(nqubits, matrices.X, condition)
+        hy = _build_spin_model(nqubits, matrices.Y, condition)
+        hz = _build_spin_model(nqubits, matrices.Z, condition)
+        matrix = hx + hy + delta * hz
+        return Hamiltonian(nqubits, matrix, backend=backend)
+
+    hx = _multikron([matrices.X, matrices.X])
+    hy = _multikron([matrices.Y, matrices.Y])
+    hz = _multikron([matrices.Z, matrices.Z])
+    matrix = hx + hy + delta * hz
+    terms = [HamiltonianTerm(matrix, i, i + 1) for i in range(nqubits - 1)]
+    terms.append(HamiltonianTerm(matrix, nqubits - 1, 0))
+    ham = SymbolicHamiltonian(backend=backend)
+    ham.terms = terms
+    return ham
+
+
+def _multikron(matrix_list):
+    """Calculates Kronecker product of a list of matrices.
+
+    Args:
+        matrix_list (list): List of matrices as ``ndarray``.
+
+    Returns:
+        ndarray: Kronecker product of all matrices in ``matrix_list``.
+    """
+    return reduce(np.kron, matrix_list)
+
+
+def _build_spin_model(nqubits, matrix, condition):
+    """Helper method for building nearest-neighbor spin model Hamiltonians."""
+    h = sum(
+        _multikron(matrix if condition(i, j) else matrices.I for j in range(nqubits))
+        for i in range(nqubits)
+    )
+    return h
+
+
+def _OneBodyPauli(nqubits, matrix, dense: bool = True, backend=None):
+    """Helper method for constracting non-interacting
+    :math:`X`, :math:`Y`, and :math:`Z` Hamiltonians."""
+    if dense:
+        condition = lambda i, j: i == j % nqubits
+        ham = -_build_spin_model(nqubits, matrix, condition)
+        return Hamiltonian(nqubits, ham, backend=backend)
+
+    matrix = -matrix
+    terms = [HamiltonianTerm(matrix, i) for i in range(nqubits)]
+    ham = SymbolicHamiltonian(backend=backend)
+    ham.terms = terms
+    return ham

From 2a68633850d8a3d514934b8bf4dd4c1487950737 Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Sat, 7 Dec 2024 14:22:11 +0400
Subject: [PATCH 02/14] rewrite XXZ

---
 src/qibo/hamiltonians/models.py | 19 ++-----------------
 1 file changed, 2 insertions(+), 17 deletions(-)

diff --git a/src/qibo/hamiltonians/models.py b/src/qibo/hamiltonians/models.py
index 30c22a491e..8524591f85 100644
--- a/src/qibo/hamiltonians/models.py
+++ b/src/qibo/hamiltonians/models.py
@@ -172,7 +172,7 @@ def Heisenberg(
             Defaults to ``None``.
 
     Returns:
-        :class:`qibo.hamiltonians.Hamiltonian` or :class:`qibo.hamiltonians.SymbolicHamiltonian`: 
+        :class:`qibo.hamiltonians.Hamiltonian` or :class:`qibo.hamiltonians.SymbolicHamiltonian`:
         Heisenberg Hamiltonian.
     """
     if isinstance(coupling_constants, (list, tuple)) and len(coupling_constants) != 3:
@@ -257,23 +257,8 @@ def XXZ(nqubits, delta=0.5, dense: bool = True, backend=None):
     """
     if nqubits < 2:
         raise_error(ValueError, "Number of qubits must be larger than one.")
-    if dense:
-        condition = lambda i, j: i in {j % nqubits, (j + 1) % nqubits}
-        hx = _build_spin_model(nqubits, matrices.X, condition)
-        hy = _build_spin_model(nqubits, matrices.Y, condition)
-        hz = _build_spin_model(nqubits, matrices.Z, condition)
-        matrix = hx + hy + delta * hz
-        return Hamiltonian(nqubits, matrix, backend=backend)
 
-    hx = _multikron([matrices.X, matrices.X])
-    hy = _multikron([matrices.Y, matrices.Y])
-    hz = _multikron([matrices.Z, matrices.Z])
-    matrix = hx + hy + delta * hz
-    terms = [HamiltonianTerm(matrix, i, i + 1) for i in range(nqubits - 1)]
-    terms.append(HamiltonianTerm(matrix, nqubits - 1, 0))
-    ham = SymbolicHamiltonian(backend=backend)
-    ham.terms = terms
-    return ham
+    return Heisenberg(nqubits, [-1, -1, -delta], 0, dense=dense, backend=backend)
 
 
 def _multikron(matrix_list):

From 3e50d4d4c5df29f5ee98b42a63d900c3276998c5 Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Sat, 7 Dec 2024 14:24:56 +0400
Subject: [PATCH 03/14] reorder api doc

---
 doc/source/api-reference/qibo.rst | 23 ++++++++++++++++-------
 1 file changed, 16 insertions(+), 7 deletions(-)

diff --git a/doc/source/api-reference/qibo.rst b/doc/source/api-reference/qibo.rst
index 3ca9e40f37..0c4c714a12 100644
--- a/doc/source/api-reference/qibo.rst
+++ b/doc/source/api-reference/qibo.rst
@@ -1267,13 +1267,6 @@ In addition to the abstract Hamiltonian models, Qibo provides the following
 pre-coded Hamiltonians:
 
 
-Heisenberg XXZ
-^^^^^^^^^^^^^^
-
-.. autoclass:: qibo.hamiltonians.XXZ
-    :members:
-    :member-order: bysource
-
 Non-interacting Pauli-X
 ^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -1326,6 +1319,22 @@ Max Cut
     This is useful for systems that contain many qubits for which constructing
     the full matrix is intractable.
 
+
+Heisenberg model
+^^^^^^^^^^^^^^^^
+
+.. autoclass:: qibo.hamiltonians.Heisenberg
+    :members:
+    :member-order: bysource
+
+
+Heisenberg XXZ
+^^^^^^^^^^^^^^
+
+.. autoclass:: qibo.hamiltonians.XXZ
+    :members:
+    :member-order: bysource
+
 _______________________
 
 

From 84cd61b1a8f84df2a722843b7b533012d04ae68b Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Sat, 7 Dec 2024 10:28:16 +0000
Subject: [PATCH 04/14] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 src/qibo/hamiltonians/__init__.py | 2 +-
 src/qibo/hamiltonians/models.py   | 3 +--
 2 files changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/qibo/hamiltonians/__init__.py b/src/qibo/hamiltonians/__init__.py
index 9826e01833..ce21da1582 100644
--- a/src/qibo/hamiltonians/__init__.py
+++ b/src/qibo/hamiltonians/__init__.py
@@ -1,2 +1,2 @@
 from qibo.hamiltonians.hamiltonians import *
-from qibo.hamiltonians.models import Heisenberg, TFIM, XXZ, MaxCut, X, Y, Z
+from qibo.hamiltonians.models import TFIM, XXZ, Heisenberg, MaxCut, X, Y, Z
diff --git a/src/qibo/hamiltonians/models.py b/src/qibo/hamiltonians/models.py
index 8524591f85..8077eb0520 100644
--- a/src/qibo/hamiltonians/models.py
+++ b/src/qibo/hamiltonians/models.py
@@ -1,6 +1,5 @@
-from typing import Union, List, Tuple
-
 from functools import reduce
+from typing import List, Tuple, Union
 
 import numpy as np
 

From b914f0bbc3522dbd57bb6c96869f241d1109445d Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Sat, 7 Dec 2024 14:36:19 +0400
Subject: [PATCH 05/14] reorder tests

---
 tests/test_hamiltonians_models.py | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/tests/test_hamiltonians_models.py b/tests/test_hamiltonians_models.py
index 8539d0e116..be69065bd2 100644
--- a/tests/test_hamiltonians_models.py
+++ b/tests/test_hamiltonians_models.py
@@ -6,18 +6,18 @@
 from qibo import hamiltonians, matrices
 
 models_config = [
-    ("TFIM", {"nqubits": 3, "h": 0.0}, "tfim_N3h0.0.out"),
-    ("TFIM", {"nqubits": 3, "h": 0.5}, "tfim_N3h0.5.out"),
-    ("TFIM", {"nqubits": 3, "h": 1.0}, "tfim_N3h1.0.out"),
-    ("XXZ", {"nqubits": 3, "delta": 0.0}, "heisenberg_N3delta0.0.out"),
-    ("XXZ", {"nqubits": 3, "delta": 0.5}, "heisenberg_N3delta0.5.out"),
-    ("XXZ", {"nqubits": 3, "delta": 1.0}, "heisenberg_N3delta1.0.out"),
     ("X", {"nqubits": 3}, "x_N3.out"),
     ("Y", {"nqubits": 4}, "y_N4.out"),
     ("Z", {"nqubits": 5}, "z_N5.out"),
+    ("TFIM", {"nqubits": 3, "h": 0.0}, "tfim_N3h0.0.out"),
+    ("TFIM", {"nqubits": 3, "h": 0.5}, "tfim_N3h0.5.out"),
+    ("TFIM", {"nqubits": 3, "h": 1.0}, "tfim_N3h1.0.out"),
     ("MaxCut", {"nqubits": 3}, "maxcut_N3.out"),
     ("MaxCut", {"nqubits": 4}, "maxcut_N4.out"),
     ("MaxCut", {"nqubits": 5}, "maxcut_N5.out"),
+    ("XXZ", {"nqubits": 3, "delta": 0.0}, "heisenberg_N3delta0.0.out"),
+    ("XXZ", {"nqubits": 3, "delta": 0.5}, "heisenberg_N3delta0.5.out"),
+    ("XXZ", {"nqubits": 3, "delta": 1.0}, "heisenberg_N3delta1.0.out"),
 ]
 
 

From f484ed71650df83a9fa9283c643fc6fe476d421a Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Sat, 7 Dec 2024 14:57:13 +0400
Subject: [PATCH 06/14] fix backend bug

---
 src/qibo/hamiltonians/models.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/src/qibo/hamiltonians/models.py b/src/qibo/hamiltonians/models.py
index 8077eb0520..5c232b380e 100644
--- a/src/qibo/hamiltonians/models.py
+++ b/src/qibo/hamiltonians/models.py
@@ -3,7 +3,7 @@
 
 import numpy as np
 
-from qibo.backends import matrices
+from qibo.backends import matrices, _check_backend
 from qibo.config import raise_error
 from qibo.hamiltonians.hamiltonians import Hamiltonian, SymbolicHamiltonian
 from qibo.hamiltonians.terms import HamiltonianTerm
@@ -182,6 +182,8 @@ def Heisenberg(
     elif isinstance(coupling_constants, (float, int)):
         coupling_constants = [coupling_constants] * 3
 
+    backend = _check_backend(backend)
+
     if dense:
         condition = lambda i, j: i in {j % nqubits, (j + 1) % nqubits}
         hx = _build_spin_model(nqubits, matrices.X, condition)
@@ -192,9 +194,10 @@ def Heisenberg(
             - coupling_constants[1] * hy
             - coupling_constants[2] * hz
         )
+        matrix = backend.cast(matrix, dtype=matrix.dtype)
 
         for pauli in [matrices.X, matrices.Y, matrices.Z]:
-            matrix += (
+            matrix = matrix + (
                 external_field_strength
                 * _OneBodyPauli(nqubits, pauli, dense, backend).matrix
             )

From adb63b6544f010594739f7236b116c6cea273a45 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Sat, 7 Dec 2024 10:57:46 +0000
Subject: [PATCH 07/14] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 src/qibo/hamiltonians/models.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/qibo/hamiltonians/models.py b/src/qibo/hamiltonians/models.py
index 5c232b380e..303a422b34 100644
--- a/src/qibo/hamiltonians/models.py
+++ b/src/qibo/hamiltonians/models.py
@@ -3,7 +3,7 @@
 
 import numpy as np
 
-from qibo.backends import matrices, _check_backend
+from qibo.backends import _check_backend, matrices
 from qibo.config import raise_error
 from qibo.hamiltonians.hamiltonians import Hamiltonian, SymbolicHamiltonian
 from qibo.hamiltonians.terms import HamiltonianTerm

From c4508398af6b8c48a7cd8fdb4581ba53c6f0b3db Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Sun, 8 Dec 2024 11:30:38 +0400
Subject: [PATCH 08/14] api ref

---
 doc/source/api-reference/qibo.rst | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/doc/source/api-reference/qibo.rst b/doc/source/api-reference/qibo.rst
index 0c4c714a12..c8bf85a1e2 100644
--- a/doc/source/api-reference/qibo.rst
+++ b/doc/source/api-reference/qibo.rst
@@ -1328,6 +1328,14 @@ Heisenberg model
     :member-order: bysource
 
 
+Heisenberg XXX
+^^^^^^^^^^^^^^
+
+.. autoclass:: qibo.hamiltonians.XXX
+    :members:
+    :member-order: bysource
+
+
 Heisenberg XXZ
 ^^^^^^^^^^^^^^
 

From 5aa804cccd28d078247ad0884bd7b73924035e73 Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Sun, 8 Dec 2024 11:30:48 +0400
Subject: [PATCH 09/14] add __init__

---
 src/qibo/hamiltonians/__init__.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/qibo/hamiltonians/__init__.py b/src/qibo/hamiltonians/__init__.py
index ce21da1582..4cf40b7e01 100644
--- a/src/qibo/hamiltonians/__init__.py
+++ b/src/qibo/hamiltonians/__init__.py
@@ -1,2 +1,2 @@
-from qibo.hamiltonians.hamiltonians import *
-from qibo.hamiltonians.models import TFIM, XXZ, Heisenberg, MaxCut, X, Y, Z
+from qibo.hamiltonians.hamiltonians import Hamiltonian, SymbolicHamiltonian, TrotterHamiltonian
+from qibo.hamiltonians.models import TFIM, XXX, XXZ, Heisenberg, MaxCut, X, Y, Z

From 3655b466b8e77eef30fe7782d9f3829f1f337e2b Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Sun, 8 Dec 2024 11:30:53 +0400
Subject: [PATCH 10/14] add XXX

---
 src/qibo/hamiltonians/models.py | 153 ++++++++++++++++++++++++--------
 1 file changed, 116 insertions(+), 37 deletions(-)

diff --git a/src/qibo/hamiltonians/models.py b/src/qibo/hamiltonians/models.py
index 303a422b34..62f6c4955c 100644
--- a/src/qibo/hamiltonians/models.py
+++ b/src/qibo/hamiltonians/models.py
@@ -136,8 +136,8 @@ def MaxCut(nqubits, dense: bool = True, backend=None):
 
 def Heisenberg(
     nqubits,
-    coupling_constants: Union[float, int, List[int], Tuple[int, ...]],
-    external_field_strength: Union[float, int],
+    coupling_constants: Union[float, int, list, tuple],
+    external_field_strengths: Union[float, int],
     dense: bool = True,
     backend=None,
 ):
@@ -146,22 +146,28 @@ def Heisenberg(
     The general :math:`n`-qubit Hamiltonian is given by
 
     .. math::
-        H = -\\sum_{k = 1}^{n} \\, \\left(
+        \\begin{align}
+        H &= -\\sum_{k = 1}^{n} \\, \\left(
             J_{x} \\, X_{k} \\, X_{k + 1}
             + J_{y} \\, Y_{k} \\, Y_{k + 1}
-            + J_{z} \\, Z_{k} \\, Z_{k + 1} \\right)
-            - h \\, \\sum_{k = 1}^{n} \\left(X_{k} + Y_{k} + Z_{k}\\right) \\, ,
+            + J_{z} \\, Z_{k} \\, Z_{k + 1} \\right) \\\\
+        &\\quad\\,\\, - \\sum_{k = 1}^{n} \\left(
+            h_{x} \\, X_{k} + h_{y} \\, Y_{k} + h_{z} \\, Z_{k}
+            \\right) \\, ,
+        \\end{align}
 
-    where :math:`\\{J_{x}, J_{y}, J_{z}\\}` are called the ``coupling constants``,
-    :math:`h` is called the ``external field strength``, and :math:`\\{X, Y, Z\\}`
-    are the usual Pauli operators.
+    where :math:`\\{J_{x}, \\, J_{y}, \\, J_{z}\\}` are called the ``coupling constants``,
+    :math:`\\{h_{x}, \\, h_{y}, \\, h_{z}\\}` are called the ``external field strengths``, 
+    and :math:`\\{X, \\, Y, \\, Z\\}` are the usual Pauli operators.
 
     Args:
         nqubits (int): number of qubits.
-        coupling_constants (float or int or list or tuple): list or tuple with the
+        coupling_constants (list or tuple or float or int): list or tuple with the
             three coupling constants :math:`\\{J_{x}, J_{y}, J{z}\\}`.
             If ``int`` or ``float``, then :math:`J_{x} = J_{y} = J_{z}`.
-        external_field_strength (float or int): external magnetic field strength :math:`h`.
+        external_field_strength (list or tuple or float or int): list or tuple with the 
+            external magnetic field strengths :math:`\\{h_{x}, \\, h_{y}, \\, h_{z}\\}`.
+            If ``int`` or ``float``, then :math:`h_{x} = h_{y} = h_{z}`.
         dense (bool, optional): If ``True``, creates the Hamiltonian as a
             :class:`qibo.core.hamiltonians.Hamiltonian`, otherwise it creates
             a :class:`qibo.core.hamiltonians.SymbolicHamiltonian`.
@@ -177,28 +183,39 @@ def Heisenberg(
     if isinstance(coupling_constants, (list, tuple)) and len(coupling_constants) != 3:
         raise_error(
             ValueError,
-            f"When `coupling_constants` is type `int` or `list`, it must have length 3.",
+            f"When `coupling_constants` is type `int` or `list`, it must have length == 3.",
         )
-    elif isinstance(coupling_constants, (float, int)):
+
+    if (
+        isinstance(external_field_strengths, (list, tuple))
+        and len(external_field_strengths) != 3
+    ):
+        raise_error(
+            ValueError,
+            f"When `external_field_strengths` is type `int` or `list`, it must have length == 3.",
+        )
+
+    if isinstance(coupling_constants, (float, int)):
         coupling_constants = [coupling_constants] * 3
 
+    if isinstance(external_field_strengths, (float, int)):
+        external_field_strengths = [external_field_strengths] * 3
+
     backend = _check_backend(backend)
 
+    paulis = [matrices.X, matrices.Y, matrices.Z]
+
     if dense:
         condition = lambda i, j: i in {j % nqubits, (j + 1) % nqubits}
-        hx = _build_spin_model(nqubits, matrices.X, condition)
-        hy = _build_spin_model(nqubits, matrices.Y, condition)
-        hz = _build_spin_model(nqubits, matrices.Z, condition)
-        matrix = (
-            -coupling_constants[0] * hx
-            - coupling_constants[1] * hy
-            - coupling_constants[2] * hz
-        )
+        matrix = np.zeros((2**nqubits, 2**nqubits), dtype=complex)
         matrix = backend.cast(matrix, dtype=matrix.dtype)
-
-        for pauli in [matrices.X, matrices.Y, matrices.Z]:
-            matrix = matrix + (
-                external_field_strength
+        for ind, pauli in enumerate(paulis):
+            double_term = _build_spin_model(nqubits, pauli, condition)
+            double_term = backend.cast(double_term, dtype=double_term.dtype)
+            matrix = matrix - coupling_constants[ind] * double_term
+            matrix = (
+                matrix
+                + external_field_strengths[ind]
                 * _OneBodyPauli(nqubits, pauli, dense, backend).matrix
             )
 
@@ -217,14 +234,14 @@ def Heisenberg(
     terms = [HamiltonianTerm(matrix, i, i + 1) for i in range(nqubits - 1)]
     terms.append(HamiltonianTerm(matrix, nqubits - 1, 0))
 
-    if external_field_strength != 0.0:
-        terms.extend(
-            [
-                -external_field_strength * HamiltonianTerm(pauli, qubit)
-                for qubit in range(nqubits)
-                for pauli in [matrices.X, matrices.Y, matrices.Z]
-            ]
-        )
+    terms.extend(
+        [
+            -field_strength * HamiltonianTerm(pauli, qubit)
+            for qubit in range(nqubits)
+            for field_strength, pauli in zip(external_field_strengths, paulis)
+            if field_strength != 0.0
+        ]
+    )
 
     ham = SymbolicHamiltonian(backend=backend)
     ham.terms = terms
@@ -232,6 +249,64 @@ def Heisenberg(
     return ham
 
 
+def XXX(
+    nqubits,
+    coupling_constant: Union[float, int] = 1,
+    external_field_strengths: Union[float, int, list, tuple] = [0.5, 0, 0],
+    dense: bool = True,
+    backend=None,
+):
+    """Heisenberg :math:`\\mathrm{XXX}` model with periodic boundary conditions.
+
+    The :math:`n`-qubit :math:`\\mathrm{XXX}` Hamiltonian is given by
+
+    .. math::
+        H = -J \\, \\sum_{k = 1}^{n} \\, \\left(
+            X_{k} \\, X_{k + 1} + Y_{k} \\, Y_{k + 1} + Z_{k} \\, Z_{k + 1}
+            \\right)
+            - \\sum_{k = 1}^{n} \\left(h_{x} \\, X_{k} + h_{y} \\, Y_{k}
+            + h_{z} \\, Z_{k} \\right) \\, ,
+
+    where :math:`J` is the ``coupling_constant``, :math:`\\{h_{x}, \\, h_{y}, \\, h_{z}\\}`
+    are called the ``external field strengths``, and :math:`\\{X, \\, Y, \\, Z\\}` are the usual
+    Pauli operators.
+
+    Args:
+        nqubits (int): number of qubits.
+        coupling_constant (float or int, optional): coupling constant :math:`J`.
+            Defaults to :math:`1`.
+        external_field_strengths (float or int or list or tuple, optional): list or tuple with the
+            external magnetic field strengths :math:`\\{h_{x}, \\, h_{y}, \\, h_{z}\\}`.
+            If ``int`` or ``float``, then :math:`h_{x} = h_{y} = h_{z}`.
+            Defaults to :math:`[1/2, \\, 0, \\, 0]`.
+        dense (bool, optional): If ``True``, creates the Hamiltonian as a
+            :class:`qibo.core.hamiltonians.Hamiltonian`, otherwise it creates
+            a :class:`qibo.core.hamiltonians.SymbolicHamiltonian`.
+            Defaults to ``True``.
+        backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be used
+            in the execution. If ``None``, it uses the current backend.
+            Defaults to ``None``.
+
+    Returns:
+        :class:`qibo.hamiltonians.Hamiltonian` or :class:`qibo.hamiltonians.SymbolicHamiltonian`:
+        Heisenberg :math:`\\mathrm{XXX}` Hamiltonian.
+    """
+    if not isinstance(coupling_constant, (float, int)):
+        raise_error(
+            TypeError,
+            "`coupling_constant` must be type float or int, "
+            + f"but it is type {type(coupling_constant)}.",
+        )
+
+    return Heisenberg(
+        nqubits,
+        coupling_constants=[coupling_constant] * 3,
+        external_field_strengths=external_field_strengths,
+        dense=dense,
+        backend=backend,
+    )
+
+
 def XXZ(nqubits, delta=0.5, dense: bool = True, backend=None):
     """Heisenberg :math:`\\mathrm{XXZ}` model with periodic boundary conditions.
 
@@ -239,6 +314,12 @@ def XXZ(nqubits, delta=0.5, dense: bool = True, backend=None):
         H = \\sum _{k=0}^N \\, \\left( X_{k} \\, X_{k + 1} + Y_{k} \\, Y_{k + 1}
             + \\delta Z_{k} \\, Z_{k + 1} \\right) \\, .
 
+    Example:
+        .. testcode::
+
+            from qibo.hamiltonians import XXZ
+            h = XXZ(3) # initialized XXZ model with 3 qubits
+
     Args:
         nqubits (int): number of qubits.
         delta (float, optional): coefficient for the :math:`Z` component.
@@ -251,11 +332,9 @@ def XXZ(nqubits, delta=0.5, dense: bool = True, backend=None):
             in the execution. If ``None``, it uses the current backend.
             Defaults to ``None``.
 
-    Example:
-        .. testcode::
-
-            from qibo.hamiltonians import XXZ
-            h = XXZ(3) # initialized XXZ model with 3 qubits
+    Returns:
+        :class:`qibo.hamiltonians.Hamiltonian` or :class:`qibo.hamiltonians.SymbolicHamiltonian`:
+        Heisenberg :math:`\\mathrm{XXZ}` Hamiltonian.
     """
     if nqubits < 2:
         raise_error(ValueError, "Number of qubits must be larger than one.")

From 77600eff8fe8e991e5dfa770468548e820ff8fa7 Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Sun, 8 Dec 2024 11:30:58 +0400
Subject: [PATCH 11/14] add tests

---
 tests/test_hamiltonians_models.py | 23 +++++++++++++++++++++++
 1 file changed, 23 insertions(+)

diff --git a/tests/test_hamiltonians_models.py b/tests/test_hamiltonians_models.py
index be69065bd2..e20ebc9962 100644
--- a/tests/test_hamiltonians_models.py
+++ b/tests/test_hamiltonians_models.py
@@ -4,6 +4,7 @@
 import pytest
 
 from qibo import hamiltonians, matrices
+from qibo.hamiltonians.models import XXX, Heisenberg
 
 models_config = [
     ("X", {"nqubits": 3}, "x_N3.out"),
@@ -59,3 +60,25 @@ def test_maxcut(backend, nqubits, dense, calcterms):
 def test_missing_neighbour_qubit(backend, model):
     with pytest.raises(ValueError):
         H = getattr(hamiltonians, model)(nqubits=1, backend=backend)
+
+
+@pytest.mark.parametrize("dense", [True, False])
+def test_xxx(backend, dense):
+    nqubits = 2
+
+    with pytest.raises(ValueError):
+        test = XXX(
+            nqubits, external_field_strengths=[0, 1], dense=dense, backend=backend
+        )
+
+    with pytest.raises(TypeError):
+        test = XXX(nqubits, coupling_constant=[1], dense=dense, backend=backend)
+
+    with pytest.raises(ValueError):
+        test = Heisenberg(
+            nqubits,
+            coupling_constants=[0, 1],
+            external_field_strengths=1,
+            dense=dense,
+            backend=backend,
+        )

From 7cb95e2ddae6ea4123e7ccf75ef5d7bacc7ca54c Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Sun, 8 Dec 2024 07:31:31 +0000
Subject: [PATCH 12/14] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 src/qibo/hamiltonians/__init__.py | 6 +++++-
 src/qibo/hamiltonians/models.py   | 6 +++---
 2 files changed, 8 insertions(+), 4 deletions(-)

diff --git a/src/qibo/hamiltonians/__init__.py b/src/qibo/hamiltonians/__init__.py
index 4cf40b7e01..35c2673e32 100644
--- a/src/qibo/hamiltonians/__init__.py
+++ b/src/qibo/hamiltonians/__init__.py
@@ -1,2 +1,6 @@
-from qibo.hamiltonians.hamiltonians import Hamiltonian, SymbolicHamiltonian, TrotterHamiltonian
+from qibo.hamiltonians.hamiltonians import (
+    Hamiltonian,
+    SymbolicHamiltonian,
+    TrotterHamiltonian,
+)
 from qibo.hamiltonians.models import TFIM, XXX, XXZ, Heisenberg, MaxCut, X, Y, Z
diff --git a/src/qibo/hamiltonians/models.py b/src/qibo/hamiltonians/models.py
index 62f6c4955c..27a9e3c6a6 100644
--- a/src/qibo/hamiltonians/models.py
+++ b/src/qibo/hamiltonians/models.py
@@ -1,5 +1,5 @@
 from functools import reduce
-from typing import List, Tuple, Union
+from typing import Union
 
 import numpy as np
 
@@ -157,7 +157,7 @@ def Heisenberg(
         \\end{align}
 
     where :math:`\\{J_{x}, \\, J_{y}, \\, J_{z}\\}` are called the ``coupling constants``,
-    :math:`\\{h_{x}, \\, h_{y}, \\, h_{z}\\}` are called the ``external field strengths``, 
+    :math:`\\{h_{x}, \\, h_{y}, \\, h_{z}\\}` are called the ``external field strengths``,
     and :math:`\\{X, \\, Y, \\, Z\\}` are the usual Pauli operators.
 
     Args:
@@ -165,7 +165,7 @@ def Heisenberg(
         coupling_constants (list or tuple or float or int): list or tuple with the
             three coupling constants :math:`\\{J_{x}, J_{y}, J{z}\\}`.
             If ``int`` or ``float``, then :math:`J_{x} = J_{y} = J_{z}`.
-        external_field_strength (list or tuple or float or int): list or tuple with the 
+        external_field_strength (list or tuple or float or int): list or tuple with the
             external magnetic field strengths :math:`\\{h_{x}, \\, h_{y}, \\, h_{z}\\}`.
             If ``int`` or ``float``, then :math:`h_{x} = h_{y} = h_{z}`.
         dense (bool, optional): If ``True``, creates the Hamiltonian as a

From a2d631574cf97cde8cd2736285e57955ae68dfd9 Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Sun, 8 Dec 2024 12:06:09 +0400
Subject: [PATCH 13/14] fix coverage

---
 src/qibo/hamiltonians/models.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/qibo/hamiltonians/models.py b/src/qibo/hamiltonians/models.py
index 27a9e3c6a6..dac1c34b7f 100644
--- a/src/qibo/hamiltonians/models.py
+++ b/src/qibo/hamiltonians/models.py
@@ -300,7 +300,7 @@ def XXX(
 
     return Heisenberg(
         nqubits,
-        coupling_constants=[coupling_constant] * 3,
+        coupling_constants=coupling_constant,
         external_field_strengths=external_field_strengths,
         dense=dense,
         backend=backend,

From 707a106fce7d922fb555b491756fabece6840ec9 Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Sun, 8 Dec 2024 13:33:24 +0400
Subject: [PATCH 14/14] fix coverage

---
 src/qibo/hamiltonians/models.py   | 6 +++---
 tests/test_hamiltonians_models.py | 6 +++++-
 2 files changed, 8 insertions(+), 4 deletions(-)

diff --git a/src/qibo/hamiltonians/models.py b/src/qibo/hamiltonians/models.py
index dac1c34b7f..61fa918021 100644
--- a/src/qibo/hamiltonians/models.py
+++ b/src/qibo/hamiltonians/models.py
@@ -186,6 +186,9 @@ def Heisenberg(
             f"When `coupling_constants` is type `int` or `list`, it must have length == 3.",
         )
 
+    if isinstance(coupling_constants, (float, int)):
+        coupling_constants = [coupling_constants] * 3
+
     if (
         isinstance(external_field_strengths, (list, tuple))
         and len(external_field_strengths) != 3
@@ -195,9 +198,6 @@ def Heisenberg(
             f"When `external_field_strengths` is type `int` or `list`, it must have length == 3.",
         )
 
-    if isinstance(coupling_constants, (float, int)):
-        coupling_constants = [coupling_constants] * 3
-
     if isinstance(external_field_strengths, (float, int)):
         external_field_strengths = [external_field_strengths] * 3
 
diff --git a/tests/test_hamiltonians_models.py b/tests/test_hamiltonians_models.py
index e20ebc9962..5d2a2017d2 100644
--- a/tests/test_hamiltonians_models.py
+++ b/tests/test_hamiltonians_models.py
@@ -68,7 +68,11 @@ def test_xxx(backend, dense):
 
     with pytest.raises(ValueError):
         test = XXX(
-            nqubits, external_field_strengths=[0, 1], dense=dense, backend=backend
+            nqubits,
+            coupling_constant=1,
+            external_field_strengths=[0, 1],
+            dense=dense,
+            backend=backend,
         )
 
     with pytest.raises(TypeError):