Repository: qutip/qutip-qip (language: python)

Issue:
Compile the documentation with qutip-v5
Fix #243 
And compile all the documentation and doctest using qutip-v5


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/black.yml b/.github/workflows/black.yml
index 21a42176..ab6aaa90 100644
--- a/.github/workflows/black.yml
+++ b/.github/workflows/black.yml
@@ -6,7 +6,7 @@ jobs:
   lint:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - uses: psf/black@stable
         with:
           options: "--check --diff"
diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index 4f47f623..0af13222 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -46,9 +46,9 @@ jobs:
       OVERRIDE_VERSION: ${{ github.event.inputs.override_version }}
 
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
 
-      - uses: actions/setup-python@v4
+      - uses: actions/setup-python@v5
         name: Install Python
         with:
           # For the sdist we should be as conservative as possible with our
@@ -68,7 +68,7 @@ jobs:
           # dependencies are specified by our setup code.
           python -m build --sdist .
 
-      - uses: actions/upload-artifact@v3
+      - uses: actions/upload-artifact@v4
         with:
           name: sdist
           path: dist/*.tar.gz
@@ -83,9 +83,9 @@ jobs:
       OVERRIDE_VERSION: ${{ github.event.inputs.override_version }}
 
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
 
-      - uses: actions/setup-python@v4
+      - uses: actions/setup-python@v5
         name: Install Python
         with:
           # This is about the build environment, not the released wheel version.
@@ -105,7 +105,7 @@ jobs:
           # dependencies are specified by our setup code.
           python -m build --wheel --outdir wheelhouse .
 
-      - uses: actions/upload-artifact@v3
+      - uses: actions/upload-artifact@v4
         with:
           name: wheels
           path: ./wheelhouse/*.whl
@@ -124,9 +124,9 @@ jobs:
 
     steps:
       - name: Download build artifacts to local runner
-        uses: actions/download-artifact@v3
+        uses: actions/download-artifact@v4
 
-      - uses: actions/setup-python@v4
+      - uses: actions/setup-python@v5
         name: Install Python
         with:
           python-version: '3.11'
diff --git a/.github/workflows/build_documentation.yml b/.github/workflows/build_documentation.yml
index 63da4050..88b56c5d 100644
--- a/.github/workflows/build_documentation.yml
+++ b/.github/workflows/build_documentation.yml
@@ -9,12 +9,12 @@ jobs:
     runs-on: ubuntu-latest
 
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
 
-      - uses: actions/setup-python@v4
+      - uses: actions/setup-python@v5
         name: Install Python
         with:
-          python-version: '3.9'
+          python-version: '3.11'
           
       - name: Install Pandoc
         run: |
@@ -27,7 +27,6 @@ jobs:
 
       - name: Install qutip-qip from GitHub
         run: |
-          python -mpip install qutip==4.7.3
           python -mpip install -e .
           # Install in editable mode so it doesn't matter if we import from
           # inside the installation directory, otherwise we can get some errors
@@ -45,7 +44,7 @@ jobs:
           #   -T : display a full traceback if a Python exception occurs
 
       - name: Upload built files
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v4
         with:
           name: qutip_qip_html_docs
           path: doc/_build/html/*
diff --git a/VERSION b/VERSION
index 9e11b32f..d15723fb 100644
--- a/VERSION
+++ b/VERSION
@@ -1,1 +1,1 @@
-0.3.1
+0.3.2
diff --git a/doc/pulse-paper/dj_algorithm.py b/doc/pulse-paper/dj_algorithm.py
index d3eebce5..d9fc93b4 100644
--- a/doc/pulse-paper/dj_algorithm.py
+++ b/doc/pulse-paper/dj_algorithm.py
@@ -136,8 +136,8 @@
 first_two_qubits = result1.states[-1].ptrace([0,1])
 _, ax4 = hinton(
     first_two_qubits, ax=ax4,
-    xlabels=['$|00\\rangle$', '$|01\\rangle$', '$|10\\rangle$', '$|11\\rangle$'],
-    ylabels=['$|00\\rangle$', '$|01\\rangle$', '$|10\\rangle$', '$|11\\rangle$'])
+    x_basis=['$|00\\rangle$', '$|01\\rangle$', '$|10\\rangle$', '$|11\\rangle$'],
+    y_basis=['$|00\\rangle$', '$|01\\rangle$', '$|10\\rangle$', '$|11\\rangle$'])
 fig4.savefig("hinton.pdf")
 fig4.show()
 
diff --git a/doc/pulse-paper/qft.py b/doc/pulse-paper/qft.py
index 81620afd..7a5bca64 100644
--- a/doc/pulse-paper/qft.py
+++ b/doc/pulse-paper/qft.py
@@ -23,7 +23,7 @@
 # Pulse-level simulation
 processor = LinearSpinChain(num_qubits)
 processor.load_circuit(qc)
-options = Options(max_step=5000)
+options = Options(max_step=5000, rtol=1.e-8)
 state2 = processor.run_state(
     basis([2]*num_qubits, [0]*num_qubits),
     options=options
diff --git a/doc/requirements.txt b/doc/requirements.txt
index e81a103c..4094f966 100644
--- a/doc/requirements.txt
+++ b/doc/requirements.txt
@@ -1,15 +1,16 @@
-numpy==1.22.4
-scipy==1.10.0
-qutip==4.7.3
-cython==0.29.30
-sphinx==5.0.2
-sphinx_rtd_theme==1.0.0
-readthedocs-sphinx-search==0.1.2
-numpydoc==1.4.0
-matplotlib==3.5.2
-docutils==0.17.1
-sphinxcontrib-bibtex==2.4.2
+numpy==1.26.4
+scipy==1.13.1
+qutip==5.0.3.post1
+qutip-qtrl==0.1.3
+cython==3.0.11
+sphinx==7.4.7
+sphinx_rtd_theme==2.0.0
+readthedocs-sphinx-search==0.3.2
+numpydoc==1.8.0
+matplotlib==3.9.2
+docutils==0.20.1
+sphinxcontrib-bibtex==2.6.2
 pyqir-generator==0.6.2
 pyqir-parser==0.6.2
-qiskit==0.46.1
-qiskit-aer==0.14.0.1
+qiskit==0.46.2
+qiskit-aer==0.14.2
\ No newline at end of file
diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst
index 500900d5..255b9c3a 100644
--- a/doc/source/changelog.rst
+++ b/doc/source/changelog.rst
@@ -3,6 +3,21 @@ Changelog
 *********
 
 
+Version 0.3.2 (Sept 3, 2024)
+++++++++++++++++++++++++++++
+
+Improvements
+------------
+- Add RZX gate to the SCQubits model. Fix bugs on the ZX strength. (`#245 <https://github.com/qutip/qutip-qip/pull/245>`_)
+- Compile the documentation using qutip v5 (`#246 <https://github.com/qutip/qutip-qip/pull/246>`_)
+
+Bug Fixes
+---------
+- Add phi parameter to expand_operator call. (`#237 <https://github.com/qutip/qutip-qip/pull/237>`_ by Neill Lambert)
+- Fix the warning on MacOS. (`#240 <https://github.com/qutip/qutip-qip/pull/240>`_)
+- Fix a bug in :class:`~.QubitCircuit.add_1q_gate`. (`#250 <https://github.com/qutip/qutip-qip/pull/250>`_ issue raised by hatlabcz Chao Zhou)
+
+
 Version 0.3.1 (April 21, 2024)
 ++++++++++++++++++++++++++++++
 
diff --git a/doc/source/qip-basics.rst b/doc/source/qip-basics.rst
index 2d94a15b..2668cceb 100644
--- a/doc/source/qip-basics.rst
+++ b/doc/source/qip-basics.rst
@@ -75,22 +75,22 @@ full dimension of the circuit:
 
 .. testoutput::
 
-  [Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True
+  [Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=CSR, isherm=True
   Qobj data =
   [[1. 0. 0. 0.]
    [0. 0. 1. 0.]
    [0. 1. 0. 0.]
-   [0. 0. 0. 1.]], Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True
+   [0. 0. 0. 1.]], Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=CSR, isherm=True
   Qobj data =
   [[1. 0. 0. 0.]
    [0. 1. 0. 0.]
    [0. 0. 0. 1.]
-   [0. 0. 1. 0.]], Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True
+   [0. 0. 1. 0.]], Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=CSR, isherm=True
   Qobj data =
   [[0. 0. 1. 0.]
    [0. 0. 0. 1.]
    [1. 0. 0. 0.]
-   [0. 1. 0. 0.]], Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True
+   [0. 1. 0. 0.]], Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=CSR, isherm=True
   Qobj data =
   [[1. 0. 0. 0.]
    [0. 0. 1. 0.]
@@ -110,25 +110,25 @@ can be achieved with the argument ``expand=False`` specified to the
 
 .. testoutput::
 
-  [Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True
-  Qobj data =
-  [[1. 0. 0. 0.]
-   [0. 0. 1. 0.]
-   [0. 1. 0. 0.]
-   [0. 0. 0. 1.]], Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True
-  Qobj data =
-  [[1. 0. 0. 0.]
-   [0. 1. 0. 0.]
-   [0. 0. 0. 1.]
-   [0. 0. 1. 0.]], Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True
-  Qobj data =
-  [[0. 1.]
-   [1. 0.]], Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True
-  Qobj data =
-  [[1. 0. 0. 0.]
-   [0. 0. 1. 0.]
-   [0. 1. 0. 0.]
-   [0. 0. 0. 1.]]]
+    [Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=Dense, isherm=True
+    Qobj data =
+    [[1. 0. 0. 0.]
+     [0. 0. 1. 0.]
+     [0. 1. 0. 0.]
+     [0. 0. 0. 1.]], Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=Dense, isherm=True
+    Qobj data =
+    [[1. 0. 0. 0.]
+     [0. 1. 0. 0.]
+     [0. 0. 0. 1.]
+     [0. 0. 1. 0.]], Quantum object: dims=[[2], [2]], shape=(2, 2), type='oper', dtype=Dense, isherm=True
+    Qobj data =
+    [[0. 1.]
+     [1. 0.]], Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=Dense, isherm=True
+    Qobj data =
+    [[1. 0. 0. 0.]
+     [0. 0. 1. 0.]
+     [0. 1. 0. 0.]
+     [0. 0. 0. 1.]]]
 
 .. _quantum_gates:
 
@@ -195,7 +195,7 @@ gate function returning a :class:`qutip.Qobj` and save it in the attribute ``use
 
       def user_gate1(arg_value):
            # controlled rotation X
-           mat = np.zeros((4, 4), dtype=np.complex)
+           mat = np.zeros((4, 4), dtype=np.complex128)
            mat[0, 0] = mat[1, 1] = 1.
            mat[2:4, 2:4] = rx(arg_value).full()
            return Qobj(mat, dims=[[2, 2], [2, 2]])
@@ -229,7 +229,7 @@ gate function returning a :class:`qutip.Qobj` and save it in the attribute ``use
 .. testoutput::
   :options: +NORMALIZE_WHITESPACE
 
-  Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = False
+  Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=CSR, isherm=False
   Qobj data =
     [[1.     +0.j      0.     +0.j      0.     +0.j      0.     +0.j     ] 
      [0.     +0.j      1.     +0.j      0.     +0.j      0.     +0.j     ] 
@@ -246,7 +246,7 @@ gate function returning a :class:`qutip.Qobj` and save it in the attribute ``use
   :options: +NORMALIZE_WHITESPACE
 
 
-  Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = False
+  Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=CSR, isherm=False
   Qobj data =
   [[1.     +0.j      0.     +0.j      0.     +0.j      0.     +0.j     ]
    [0.     +0.j      0.70711+0.j      0.     +0.j      0.     -0.70711j]
@@ -263,12 +263,12 @@ gate function returning a :class:`qutip.Qobj` and save it in the attribute ``use
 .. testoutput::
   :options: +NORMALIZE_WHITESPACE
 
-  Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = False
-  Qobj data =
-  [[1.+0.j 0.+0.j 0.+0.j 0.+0.j]
-   [0.+0.j 0.+1.j 0.+0.j 0.+0.j]
-   [0.+0.j 0.+0.j 1.+0.j 0.+0.j]
-   [0.+0.j 0.+0.j 0.+0.j 0.+1.j]]
+    Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=CSR, isherm=False
+    Qobj data =
+    [[1.+0.j 0.+0.j 0.+0.j 0.+0.j]
+     [0.+0.j 0.+1.j 0.+0.j 0.+0.j]
+     [0.+0.j 0.+0.j 1.+0.j 0.+0.j]
+     [0.+0.j 0.+0.j 0.+0.j 0.+1.j]]
 
 .. _quantum_circuit_plots:
 
diff --git a/doc/source/qip-processor.rst b/doc/source/qip-processor.rst
index 8c374568..eb1ca3cc 100644
--- a/doc/source/qip-processor.rst
+++ b/doc/source/qip-processor.rst
@@ -350,15 +350,16 @@ The following example demonstrates a biased Gaussian noise on the pulse amplitud
     # Plot the ideal pulse
     fig1, axis1 = processor.plot_pulses(
         title="Original control amplitude", figsize=(5,3),
-        use_control_latex=False)
+        use_control_latex=False, dpi=200)
 
     # Plot the noisy pulse
     qobjevo, _ = processor.get_qobjevo(noisy=True)
-    noisy_coeff = qobjevo.to_list()[1][1] + qobjevo.to_list()[2][1]
     fig2, axis2 = processor.plot_pulses(
         title="Noisy control amplitude", figsize=(5,3),
-        use_control_latex=False)
-    axis2[0].step(qobjevo.tlist, noisy_coeff)
+        use_control_latex=False, dpi=200)
+    tlist = np.linspace(0., processor.get_full_tlist()[-1], 300)
+    noisy_coeff = [qobjevo(t)[0, 0] for t in tlist]
+    axis2[0].step(tlist, noisy_coeff)
 
 Lindblad noise
 --------------
diff --git a/doc/source/qip-simulator.rst b/doc/source/qip-simulator.rst
index 553652a3..b766c673 100644
--- a/doc/source/qip-simulator.rst
+++ b/doc/source/qip-simulator.rst
@@ -52,7 +52,7 @@ method.
 .. testoutput::
   :options: +NORMALIZE_WHITESPACE
 
-  Quantum object: dims = [[2, 2, 2], [1, 1, 1]], shape = (8, 1), type = ket
+  Quantum object: dims=[[2, 2, 2], [1, 1, 1]], shape=(8, 1), type='ket', dtype=Dense
   Qobj data =
   [[0.     ]
    [0.57735]
@@ -95,7 +95,7 @@ outputs, we can use the :meth:`.QubitCircuit.run_statistics` function:
   :options: +NORMALIZE_WHITESPACE
 
   State:
-  Quantum object: dims = [[2, 2, 2], [1, 1, 1]], shape = (8, 1), type = ket
+  Quantum object: dims=[[2, 2, 2], [1, 1, 1]], shape=(8, 1), type='ket', dtype=Dense
   Qobj data =
   [[0.]
    [1.]
@@ -107,7 +107,7 @@ outputs, we can use the :meth:`.QubitCircuit.run_statistics` function:
    [0.]]
   with probability 0.33333
   State:
-  Quantum object: dims = [[2, 2, 2], [1, 1, 1]], shape = (8, 1), type = ket
+  Quantum object: dims=[[2, 2, 2], [1, 1, 1]], shape=(8, 1), type='ket', dtype=Dense
   Qobj data =
   [[0.]
    [0.]
@@ -119,7 +119,7 @@ outputs, we can use the :meth:`.QubitCircuit.run_statistics` function:
    [0.]]
   with probability 0.33333
   State:
-  Quantum object: dims = [[2, 2, 2], [1, 1, 1]], shape = (8, 1), type = ket
+  Quantum object: dims=[[2, 2, 2], [1, 1, 1]], shape=(8, 1), type='ket', dtype=Dense
   Qobj data =
   [[0.]
    [0.]
@@ -176,7 +176,7 @@ The :class:`.CircuitSimulator` class also enables stepping through the circuit:
 .. testoutput::
   :options: +NORMALIZE_WHITESPACE
 
-  Quantum object: dims = [[2, 2, 2], [1, 1, 1]], shape = (8, 1), type = ket
+  Quantum object: dims=[[2, 2, 2], [1, 1, 1]], shape=(8, 1), type='ket', dtype=Dense
   Qobj data =
   [[0.57735]
    [0.     ]
@@ -208,19 +208,19 @@ precomputes the product of the unitaries (in between the measurements):
 .. testoutput::
   :options: +NORMALIZE_WHITESPACE
 
-  [Quantum object: dims = [[2, 2, 2], [2, 2, 2]], shape = (8, 8), type = oper, isherm = False
-    Qobj data =
-    [[ 0.       0.57735  0.      -0.57735  0.       0.40825  0.      -0.40825]     
-     [ 0.57735  0.      -0.57735  0.       0.40825  0.      -0.40825  0.     ]     
-     [ 0.57735  0.       0.57735  0.       0.40825  0.       0.40825  0.     ]     
-     [ 0.       0.57735  0.       0.57735  0.       0.40825  0.       0.40825]     
-     [ 0.57735  0.       0.       0.      -0.8165   0.       0.       0.     ]     
-     [ 0.       0.57735  0.       0.       0.      -0.8165   0.       0.     ]     
-     [ 0.       0.       0.57735  0.       0.       0.      -0.8165   0.     ]     
-     [ 0.       0.       0.       0.57735  0.       0.       0.      -0.8165 ]],
-       Measurement(M0, target=[0], classical_store=0),
-       Measurement(M1, target=[1], classical_store=1),
-       Measurement(M2, target=[2], classical_store=2)]
+  [Quantum object: dims=[[2, 2, 2], [2, 2, 2]], shape=(8, 8), type='oper', dtype=CSR, isherm=False
+  Qobj data =
+  [[ 0.       0.57735  0.      -0.57735  0.       0.40825  0.      -0.40825]
+    [ 0.57735  0.      -0.57735  0.       0.40825  0.      -0.40825  0.     ]
+    [ 0.57735  0.       0.57735  0.       0.40825  0.       0.40825  0.     ]
+    [ 0.       0.57735  0.       0.57735  0.       0.40825  0.       0.40825]
+    [ 0.57735  0.       0.       0.      -0.8165   0.       0.       0.     ]
+    [ 0.       0.57735  0.       0.       0.      -0.8165   0.       0.     ]
+    [ 0.       0.       0.57735  0.       0.       0.      -0.8165   0.     ]
+    [ 0.       0.       0.       0.57735  0.       0.       0.      -0.8165 ]],
+    Measurement(M0, target=[0], classical_store=0),
+    Measurement(M1, target=[1], classical_store=1),
+    Measurement(M2, target=[2], classical_store=2)]
 
 
 Here, ``sim.ops`` stores all the circuit operations that are going to be applied during
@@ -264,7 +264,7 @@ just by measurement on the first qubit:
 .. testoutput::
   :options: +NORMALIZE_WHITESPACE
 
-    Quantum object: dims = [[2, 2, 2], [2, 2, 2]], shape = (8, 8), type = oper, isherm = True
+    Quantum object: dims=[[2, 2, 2], [2, 2, 2]], shape=(8, 8), type='oper', dtype=Dense, isherm=True
     Qobj data =
     [[0.      0.      0.      0.      0.      0.      0.      0.     ]
      [0.      0.33333 0.      0.      0.      0.      0.      0.     ]
diff --git a/src/qutip_qip/circuit/_decompose.py b/src/qutip_qip/circuit/_decompose.py
index 324d1b31..2eae07e2 100644
--- a/src/qutip_qip/circuit/_decompose.py
+++ b/src/qutip_qip/circuit/_decompose.py
@@ -59,6 +59,9 @@ def _gate_SNOT(gate, temp_resolved):
     )
 
 
+_gate_H = _gate_SNOT
+
+
 def _gate_PHASEGATE(gate, temp_resolved):
     temp_resolved.append(
         Gate(
diff --git a/src/qutip_qip/circuit/circuit.py b/src/qutip_qip/circuit/circuit.py
index d74ad665..eb6197a6 100644
--- a/src/qutip_qip/circuit/circuit.py
+++ b/src/qutip_qip/circuit/circuit.py
@@ -7,6 +7,7 @@
 import inspect
 import os
 from functools import partialmethod
+from typing import Optional, Union, Tuple, List, Dict, Any
 
 import numpy as np
 from copy import deepcopy
@@ -67,6 +68,7 @@ class QubitCircuit:
 
     Examples
     --------
+    >>> from qutip_qip.circuit import QubitCircuit
     >>> def user_gate():
     ...     mat = np.array([[1.,   0],
     ...                     [0., 1.j]])
@@ -271,11 +273,7 @@ def add_1q_gate(
         start=0,
         end=None,
         qubits=None,
-        arg_value=None,
-        arg_label=None,
-        classical_controls=None,
-        control_value=None,
-        classical_control_value=None,
+        **kwargs,
     ):
         """
         Adds a single qubit gate with specified parameters on a variable
@@ -285,45 +283,21 @@ def add_1q_gate(
         Parameters
         ----------
         name : string
-            Gate name.
+            Gate name or the :class:`~.operations.Gate` object.
         start : int
             Starting location of qubits.
         end : int
             Last qubit for the gate.
         qubits : list
             Specific qubits for applying gates.
-        arg_value : float
-            Argument value(phi).
-        arg_label : string
-            Label for gate representation.
+        kwargs : dict
+            Keyword arguments for the gate, except for `targets`.
+            See :class:`~.QubitCircuit.add_gate`.
         """
-        if qubits is not None:
-            for _, i in enumerate(qubits):
-                gate = GATE_CLASS_MAP[name](
-                    targets=qubits[i],
-                    controls=None,
-                    arg_value=arg_value,
-                    arg_label=arg_label,
-                    classical_controls=classical_controls,
-                    control_value=control_value,
-                    classical_control_value=classical_control_value,
-                )
-                self.gates.append(gate)
-
-        else:
-            if end is None:
-                end = self.N - 1
-            for i in range(start, end + 1):
-                gate = GATE_CLASS_MAP[name](
-                    targets=i,
-                    controls=None,
-                    arg_value=arg_value,
-                    arg_label=arg_label,
-                    classical_controls=classical_controls,
-                    control_value=control_value,
-                    classical_control_value=classical_control_value,
-                )
-                self.gates.append(gate)
+        if qubits is None:
+            qubits = range(start, end + 1)
+        for q in qubits:
+            self.add_gate(name, targets=q, **kwargs)
 
     def add_circuit(self, qc, start=0, overwrite_user_gates=False):
         """
@@ -595,9 +569,6 @@ def resolve_gates(self, basis=["CNOT", "RX", "RY", "RZ"]):
                     basis_1q.append(gate)
                 else:
                     pass
-                    raise NotImplementedError(
-                        "%s is not a valid basis gate" % gate
-                    )
             if len(basis_1q) == 1:
                 raise ValueError("Not sufficient single-qubit gates in basis")
             if len(basis_1q) == 0:
@@ -620,9 +591,12 @@ def resolve_gates(self, basis=["CNOT", "RX", "RY", "RZ"]):
                 )
             try:
                 _resolve_to_universal(gate, temp_resolved, basis_1q, basis_2q)
-            except AttributeError:
-                exception = f"Gate {gate.name} cannot be resolved."
-                raise NotImplementedError(exception)
+            except KeyError:
+                if gate.name in basis:
+                    temp_resolved.append(gate)
+                else:
+                    exception = f"Gate {gate.name} cannot be resolved."
+                    raise NotImplementedError(exception)
 
         match = False
         for basis_unit in ["CSIGN", "ISWAP", "SQRTSWAP", "SQRTISWAP"]:
diff --git a/src/qutip_qip/compiler/circuitqedcompiler.py b/src/qutip_qip/compiler/circuitqedcompiler.py
index 52c0f1d7..12434a0f 100644
--- a/src/qutip_qip/compiler/circuitqedcompiler.py
+++ b/src/qutip_qip/compiler/circuitqedcompiler.py
@@ -91,6 +91,7 @@ def __init__(self, num_qubits, params):
                 "RY": self.ry_compiler,
                 "RX": self.rx_compiler,
                 "CNOT": self.cnot_compiler,
+                "RZX": self.rzx_compiler,
             }
         )
         self.args = {  # Default configuration
@@ -130,10 +131,22 @@ def _rotation_compiler(self, gate, op_label, param_label, args):
             maximum=self.params[param_label][targets[0]],
             area=gate.arg_value / 2.0 / np.pi,
         )
+        f = 2 * np.pi * self.params["wq"][targets[0]]
         if args["DRAG"]:
             pulse_info = self._drag_pulse(op_label, coeff, tlist, targets[0])
+        elif op_label == "sx":
+            pulse_info = [
+                ("sx" + str(targets[0]), coeff),
+                # Add zero here just to make it easier to add the driving frequency later.
+                ("sy" + str(targets[0]), np.zeros(len(coeff))),
+            ]
+        elif op_label == "sy":
+            pulse_info = [
+                ("sx" + str(targets[0]), np.zeros(len(coeff))),
+                ("sy" + str(targets[0]), coeff),
+            ]
         else:
-            pulse_info = [(op_label + str(targets[0]), coeff)]
+            raise RuntimeError("Unknown label.")
         return [Instruction(gate, tlist, pulse_info)]
 
     def _drag_pulse(self, op_label, coeff, tlist, target):
@@ -198,6 +211,41 @@ def rx_compiler(self, gate, args):
         """
         return self._rotation_compiler(gate, "sx", "omega_single", args)
 
+    def rzx_compiler(self, gate, args):
+        """
+        Cross-Resonance RZX rotation, building block for the CNOT gate.
+
+        Parameters
+        ----------
+        gate : :obj:`~.operations.Gate`:
+            The quantum gate to be compiled.
+        args : dict
+            The compilation configuration defined in the attributes
+            :obj:`.GateCompiler.args` or given as a parameter in
+            :obj:`.GateCompiler.compile`.
+
+        Returns
+        -------
+        A list of :obj:`.Instruction`, including the compiled pulse
+        information for this gate.
+        """
+        result = []
+        q1, q2 = gate.targets
+        if q1 < q2:
+            zx_coeff = self.params["zx_coeff"][2 * q1]
+        else:
+            zx_coeff = self.params["zx_coeff"][2 * q1 - 1]
+        area = 0.5
+        coeff, tlist = self.generate_pulse_shape(
+            args["shape"], args["num_samples"], maximum=zx_coeff, area=area
+        )
+        area_rescale_factor = np.sqrt(np.abs(gate.arg_value) / (np.pi / 2))
+        tlist *= area_rescale_factor
+        coeff *= area_rescale_factor
+        pulse_info = [("zx" + str(q1) + str(q2), coeff)]
+        result += [Instruction(gate, tlist, pulse_info)]
+        return result
+
     def cnot_compiler(self, gate, args):
         """
         Compiler for CNOT gate using the cross resonance iteraction.
@@ -226,13 +274,8 @@ def cnot_compiler(self, gate, args):
         gate1 = Gate("RX", q2, arg_value=-np.pi / 2)
         result += self.gate_compiler[gate1.name](gate1, args)
 
-        zx_coeff = self.params["zx_coeff"][q1]
-        area = 1 / 2
-        coeff, tlist = self.generate_pulse_shape(
-            args["shape"], args["num_samples"], maximum=zx_coeff, area=area
-        )
-        pulse_info = [("zx" + str(q1) + str(q2), coeff)]
-        result += [Instruction(gate, tlist, pulse_info)]
+        gate2 = Gate("RZX", targets=[q1, q2], arg_value=np.pi / 2)
+        result += self.rzx_compiler(gate2, args)
 
         gate3 = Gate("RX", q1, arg_value=-np.pi / 2)
         result += self.gate_compiler[gate3.name](gate3, args)
diff --git a/src/qutip_qip/decompose/decompose_single_qubit_gate.py b/src/qutip_qip/decompose/decompose_single_qubit_gate.py
index 7c0d5058..a00ca5ab 100644
--- a/src/qutip_qip/decompose/decompose_single_qubit_gate.py
+++ b/src/qutip_qip/decompose/decompose_single_qubit_gate.py
@@ -1,3 +1,4 @@
+import warnings
 import numpy as np
 import cmath
 
@@ -23,7 +24,14 @@ def _angles_for_ZYZ(input_gate):
         The gate matrix that's supposed to be decomposed should be a Qobj.
     """
     input_array = input_gate.full()
-    normalization_constant = np.sqrt(np.linalg.det(input_array))
+
+    # On MacOS it raises a warning about dvived by 0 with numpy 1.26.4.
+    # It may be related to some openblas issue.
+    # https://github.com/numpy/numpy/issues/25214.
+    # Remove this later if the issue is fixed.
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        normalization_constant = np.sqrt(np.linalg.det(input_array))
     global_phase_angle = -cmath.phase(1 / normalization_constant)
     input_array = input_array * (1 / normalization_constant)
 
diff --git a/src/qutip_qip/device/circuitqed.py b/src/qutip_qip/device/circuitqed.py
index 2e8bc1fd..8f6a5828 100644
--- a/src/qutip_qip/device/circuitqed.py
+++ b/src/qutip_qip/device/circuitqed.py
@@ -73,7 +73,7 @@ def __init__(self, num_qubits, dims=None, zz_crosstalk=False, **params):
             **params,
         )
         super(SCQubits, self).__init__(model=model)
-        self.native_gates = ["RX", "RY", "CNOT"]
+        self.native_gates = ["RX", "RY", "CNOT", "RZX"]
         self._default_compiler = SCQubitsCompiler
         self.pulse_mode = "continuous"
 
@@ -318,13 +318,12 @@ def _compute_params(self):
                 )
             )
             zx_coeff.append(tmp)
-        for i in range(num_qubits - 1, 0, -1):
             tmp = (
-                J[i - 1]
-                * omega_cr[i]
+                J[i]
+                * omega_cr[i + 1]
                 * (
-                    1 / (wq[i] - wq[i - 1] + alpha[i])
-                    - 1 / (wq[i] - wq[i - 1])
+                    1 / (wq[i + 1] - wq[i] + alpha[i + 1])
+                    - 1 / (wq[i + 1] - wq[i])
                 )
             )
             zx_coeff.append(tmp)
diff --git a/src/qutip_qip/device/optpulseprocessor.py b/src/qutip_qip/device/optpulseprocessor.py
index b454fb9f..dbdead4e 100644
--- a/src/qutip_qip/device/optpulseprocessor.py
+++ b/src/qutip_qip/device/optpulseprocessor.py
@@ -5,7 +5,6 @@
 import numpy as np
 
 from qutip import Qobj, identity, tensor, mesolve
-import qutip.control.pulseoptim as cpo
 from ..circuit import QubitCircuit
 from .processor import Processor
 from ..operations import gate_sequence_product, expand_operator
@@ -195,6 +194,8 @@ def load_circuit(
                 ),
             )
 
+            import qutip.control.pulseoptim as cpo
+
             result = cpo.optimize_pulse_unitary(
                 full_drift_ham, full_ctrls_hams, U_0, U_targ, **kwargs
             )
diff --git a/src/qutip_qip/device/processor.py b/src/qutip_qip/device/processor.py
index cc6e06fd..bf7b4f1c 100644
--- a/src/qutip_qip/device/processor.py
+++ b/src/qutip_qip/device/processor.py
@@ -272,8 +272,8 @@ def add_control(
         >>> processor.get_control_labels()
         ['sx']
         >>> processor.get_control("sx") # doctest: +NORMALIZE_WHITESPACE
-        (Quantum object: dims = [[2], [2]], shape = (2, 2),
-        type = oper, isherm = True
+        (Quantum object: dims=[[2], [2]], shape=(2, 2),
+        type='oper', dtype=CSR, isherm=True
         Qobj data =
         [[0. 1.]
         [1. 0.]], [0])
@@ -310,8 +310,8 @@ def get_control(self, label):
         >>> processor.get_control_labels()
         ['sx0', 'sz0']
         >>> processor.get_control('sz0') # doctest: +NORMALIZE_WHITESPACE
-        (Quantum object: dims = [[2], [2]], shape = (2, 2),
-        type = oper, isherm = True
+        (Quantum object: dims=[[2], [2]], shape=(2, 2),
+        type='oper', dtype=CSR, isherm=True
         Qobj data =
         [[ 6.28319  0.     ]
          [ 0.      -6.28319]], 0)
@@ -997,10 +997,11 @@ def get_qobjevo(self, args=None, noisy=False):
                 qu = pulse.get_ideal_qobjevo(dims=self.dims)
             qu_list.append(qu)
 
-        final_qu = _merge_qobjevo(qu_list)
         if is_qutip5:
+            final_qu = sum(qu_list)
             final_qu.arguments(args)
         else:
+            final_qu = _merge_qobjevo(qu_list)
             final_qu.args.update(args)
 
         # bring all c_ops to the same tlist, won't need it in QuTiP 5
diff --git a/src/qutip_qip/noise.py b/src/qutip_qip/noise.py
index 94e61c5b..cb8635ce 100644
--- a/src/qutip_qip/noise.py
+++ b/src/qutip_qip/noise.py
@@ -470,7 +470,7 @@ def __init__(self, dt, rand_gen, indices=None, **kwargs):
         self.rand_gen = rand_gen
         self.kwargs = kwargs
         if "size" in kwargs:
-            raise ValueError("size is preditermined inside the noise object.")
+            raise ValueError("size is predetermined inside the noise object.")
         self.dt = dt
         self.indices = indices
 
@@ -510,7 +510,7 @@ def get_noisy_pulses(self, dims=None, pulses=None, systematic_noise=None):
         # create new tlist and random coeff
         num_rand = int(np.floor((t_max - t_min) / self.dt)) + 1
         tlist = np.arange(0, self.dt * num_rand, self.dt)[:num_rand] + t_min
-        # [:num_rand] for round of error like 0.2*6=1.2000000000002
+        # [:num_rand] for round off error like 0.2*6=1.2000000000002
 
         for i in indices:
             pulse = pulses[i]
diff --git a/src/qutip_qip/operations/gateclass.py b/src/qutip_qip/operations/gateclass.py
index a125bad8..b7817db7 100644
--- a/src/qutip_qip/operations/gateclass.py
+++ b/src/qutip_qip/operations/gateclass.py
@@ -49,6 +49,13 @@
     expand_operator,
 )
 
+from packaging.version import parse as parse_version
+
+if parse_version(qutip.__version__) >= parse_version("5.dev"):
+    is_qutip5 = True
+else:
+    is_qutip5 = False
+
 __all__ = [
     "Gate",
     "GATE_CLASS_MAP",
@@ -86,6 +93,7 @@
     "CS",
     "CT",
     "CPHASE",
+    "RZX",
 ]
 
 """
@@ -422,7 +430,7 @@ class X(SingleQubitGate):
     --------
     >>> from qutip_qip.operations import X
     >>> X(0).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True
+    Quantum object: dims=[[2], [2]], shape=(2, 2), type='oper', dtype=Dense, isherm=True
     Qobj data =
     [[0. 1.]
      [1. 0.]]
@@ -433,7 +441,7 @@ def __init__(self, targets, **kwargs):
         self.latex_str = r"X"
 
     def get_compact_qobj(self):
-        return qutip.sigmax()
+        return qutip.sigmax() if not is_qutip5 else qutip.sigmax(dtype="dense")
 
 
 class Y(SingleQubitGate):
@@ -444,7 +452,7 @@ class Y(SingleQubitGate):
     --------
     >>> from qutip_qip.operations import Y
     >>> Y(0).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True
+    Quantum object: dims=[[2], [2]], shape=(2, 2), type='oper', dtype=Dense, isherm=True
     Qobj data =
     [[0.+0.j 0.-1.j]
      [0.+1.j 0.+0.j]]
@@ -455,7 +463,7 @@ def __init__(self, targets, **kwargs):
         self.latex_str = r"Y"
 
     def get_compact_qobj(self):
-        return qutip.sigmay()
+        return qutip.sigmay() if not is_qutip5 else qutip.sigmay(dtype="dense")
 
 
 class Z(SingleQubitGate):
@@ -466,7 +474,7 @@ class Z(SingleQubitGate):
     --------
     >>> from qutip_qip.operations import Z
     >>> Z(0).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True
+    Quantum object: dims=[[2], [2]], shape=(2, 2), type='oper', dtype=Dense, isherm=True
     Qobj data =
     [[ 1.  0.]
      [ 0. -1.]]
@@ -477,7 +485,7 @@ def __init__(self, targets, **kwargs):
         self.latex_str = r"Z"
 
     def get_compact_qobj(self):
-        return qutip.sigmaz()
+        return qutip.sigmaz() if not is_qutip5 else qutip.sigmaz(dtype="dense")
 
 
 class RX(SingleQubitGate):
@@ -488,7 +496,7 @@ class RX(SingleQubitGate):
     --------
     >>> from qutip_qip.operations import RX
     >>> RX(0, 3.14159/2).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = False
+    Quantum object: dims=[[2], [2]], shape=(2, 2), type='oper', dtype=Dense, isherm=False
     Qobj data =
     [[0.70711+0.j      0.     -0.70711j]
      [0.     -0.70711j 0.70711+0.j     ]]
@@ -510,7 +518,7 @@ class RY(SingleQubitGate):
     --------
     >>> from qutip_qip.operations import RY
     >>> RY(0, 3.14159/2).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = False
+    Quantum object: dims=[[2], [2]], shape=(2, 2), type='oper', dtype=Dense, isherm=False
     Qobj data =
     [[ 0.70711 -0.70711]
      [ 0.70711  0.70711]]
@@ -532,7 +540,7 @@ class RZ(SingleQubitGate):
     --------
     >>> from qutip_qip.operations import RZ
     >>> RZ(0, 3.14159/2).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = False
+    Quantum object: dims=[[2], [2]], shape=(2, 2), type='oper', dtype=Dense, isherm=False
     Qobj data =
     [[0.70711-0.70711j 0.     +0.j     ]
      [0.     +0.j      0.70711+0.70711j]]
@@ -554,7 +562,7 @@ class H(SingleQubitGate):
     --------
     >>> from qutip_qip.operations import H
     >>> H(0).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True
+    Quantum object: dims=[[2], [2]], shape=(2, 2), type='oper', dtype=Dense, isherm=True
     Qobj data =
     [[ 0.70711  0.70711]
      [ 0.70711 -0.70711]]
@@ -580,7 +588,7 @@ class SQRTNOT(SingleQubitGate):
     --------
     >>> from qutip_qip.operations import SQRTNOT
     >>> SQRTNOT(0).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = False
+    Quantum object: dims=[[2], [2]], shape=(2, 2), type='oper', dtype=Dense, isherm=False
     Qobj data =
     [[0.5+0.5j 0.5-0.5j]
      [0.5-0.5j 0.5+0.5j]]
@@ -602,7 +610,7 @@ class S(SingleQubitGate):
     --------
     >>> from qutip_qip.operations import S
     >>> S(0).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = False
+    Quantum object: dims=[[2], [2]], shape=(2, 2), type='oper', dtype=Dense, isherm=False
     Qobj data =
     [[1.+0.j 0.+0.j]
      [0.+0.j 0.+1.j]]
@@ -624,7 +632,7 @@ class T(SingleQubitGate):
     --------
     >>> from qutip_qip.operations import T
     >>> T(0).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = False
+    Quantum object: dims=[[2], [2]], shape=(2, 2), type='oper', dtype=Dense, isherm=False
     Qobj data =
     [[1.     +0.j      0.     +0.j     ]
      [0.     +0.j      0.70711+0.70711j]]
@@ -653,7 +661,7 @@ class R(SingleQubitGate):
     --------
     >>> from qutip_qip.operations import R
     >>> R(0, (np.pi/2, np.pi/2)).get_compact_qobj().tidyup() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = False
+    Quantum object: dims=[[2], [2]], shape=(2, 2), type='oper', dtype=Dense, isherm=False
     Qobj data =
     [[ 0.70711 -0.70711]
      [ 0.70711  0.70711]]
@@ -678,7 +686,7 @@ class QASMU(SingleQubitGate):
     --------
     >>> from qutip_qip.operations import QASMU
     >>> QASMU(0, (np.pi/2, np.pi, np.pi/2)).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = False
+    Quantum object: dims=[[2], [2]], shape=(2, 2), type='oper', dtype=Dense, isherm=False
     Qobj data =
     [[-0.5-0.5j -0.5+0.5j]
      [ 0.5+0.5j -0.5+0.5j]]
@@ -711,7 +719,7 @@ class SWAP(TwoQubitGate):
     --------
     >>> from qutip_qip.operations import SWAP
     >>> SWAP([0, 1]).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True
+    Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=Dense, isherm=True
     Qobj data =
     [[1. 0. 0. 0.]
      [0. 0. 1. 0.]
@@ -735,7 +743,7 @@ class ISWAP(TwoQubitGate):
     --------
     >>> from qutip_qip.operations import ISWAP
     >>> ISWAP([0, 1]).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = False
+    Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=Dense, isherm=False
     Qobj data =
     [[1.+0.j 0.+0.j 0.+0.j 0.+0.j]
      [0.+0.j 0.+0.j 0.+1.j 0.+0.j]
@@ -759,7 +767,7 @@ class SQRTSWAP(TwoQubitGate):
     --------
     >>> from qutip_qip.operations import SQRTSWAP
     >>> SQRTSWAP([0, 1]).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = False
+    Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=Dense, isherm=False
     Qobj data =
     [[1. +0.j  0. +0.j  0. +0.j  0. +0.j ]
      [0. +0.j  0.5+0.5j 0.5-0.5j 0. +0.j ]
@@ -783,7 +791,7 @@ class SQRTISWAP(TwoQubitGate):
     --------
     >>> from qutip_qip.operations import SQRTISWAP
     >>> SQRTISWAP([0, 1]).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = False
+    Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=Dense, isherm=False
     Qobj data =
     [[1.     +0.j      0.     +0.j      0.     +0.j      0.     +0.j     ]
      [0.     +0.j      0.70711+0.j      0.     +0.70711j 0.     +0.j     ]
@@ -816,7 +824,7 @@ class BERKELEY(TwoQubitGate):
     --------
     >>> from qutip_qip.operations import BERKELEY
     >>> BERKELEY([0, 1]).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = False
+    Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=Dense, isherm=False
     Qobj data =
     [[0.92388+0.j      0.     +0.j      0.     +0.j      0.     +0.38268j]
      [0.     +0.j      0.38268+0.j      0.     +0.92388j 0.     +0.j     ]
@@ -849,7 +857,7 @@ class SWAPALPHA(TwoQubitGate):
     --------
     >>> from qutip_qip.operations import SWAPALPHA
     >>> SWAPALPHA([0, 1], 0.5).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = False
+    Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=Dense, isherm=False
     Qobj data =
     [[1. +0.j  0. +0.j  0. +0.j  0. +0.j ]
      [0. +0.j  0.5+0.5j 0.5-0.5j 0. +0.j ]
@@ -882,7 +890,7 @@ class MS(TwoQubitGate):
     --------
     >>> from qutip_qip.operations import MS
     >>> MS([0, 1], (np.pi/2, 0)).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = False
+    Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=Dense, isherm=False
     Qobj data =
     [[0.70711+0.j      0.     +0.j      0.     +0.j      0.     -0.70711j]
      [0.     +0.j      0.70711+0.j      0.     -0.70711j 0.     +0.j     ]
@@ -906,7 +914,7 @@ class TOFFOLI(Gate):
     --------
     >>> from qutip_qip.operations import TOFFOLI
     >>> TOFFOLI([0, 1, 2]).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2, 2, 2], [2, 2, 2]], shape = (8, 8), type = oper, isherm = True
+    Quantum object: dims=[[2, 2, 2], [2, 2, 2]], shape=(8, 8), type='oper', dtype=Dense, isherm=True
     Qobj data =
     [[1. 0. 0. 0. 0. 0. 0. 0.]
      [0. 1. 0. 0. 0. 0. 0. 0.]
@@ -934,7 +942,7 @@ class FREDKIN(Gate):
     --------
     >>> from qutip_qip.operations import FREDKIN
     >>> FREDKIN([0, 1, 2]).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2, 2, 2], [2, 2, 2]], shape = (8, 8), type = oper, isherm = True
+    Quantum object: dims=[[2, 2, 2], [2, 2, 2]], shape=(8, 8), type='oper', dtype=Dense, isherm=True
     Qobj data =
     [[1. 0. 0. 0. 0. 0. 0. 0.]
      [0. 1. 0. 0. 0. 0. 0. 0.]
@@ -1025,7 +1033,7 @@ class CNOT(_OneControlledGate):
     --------
     >>> from qutip_qip.operations import CNOT
     >>> CNOT(0, 1).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True
+    Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=Dense, isherm=True
     Qobj data =
     [[1. 0. 0. 0.]
      [0. 1. 0. 0.]
@@ -1055,7 +1063,7 @@ class CZ(_OneControlledGate):
     --------
     >>> from qutip_qip.operations import CZ
     >>> CSIGN(0, 1).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True
+    Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=Dense, isherm=True
     Qobj data =
     [[ 1.  0.  0.  0.]
      [ 0.  1.  0.  0.]
@@ -1084,7 +1092,7 @@ class CSIGN(_OneControlledGate):
     --------
     >>> from qutip_qip.operations import CSIGN
     >>> CSIGN(0, 1).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True
+    Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=Dense, isherm=True
     Qobj data =
     [[ 1.  0.  0.  0.]
      [ 0.  1.  0.  0.]
@@ -1122,7 +1130,7 @@ class CPHASE(_OneControlledGate):
     --------
     >>> from qutip_qip.operations import CPHASE
     >>> CPHASE(0, 1, np.pi/2).get_compact_qobj() # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = False
+    Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=Dense, isherm=False
     Qobj data =
     [[1.+0.j 0.+0.j 0.+0.j 0.+0.j]
      [0.+0.j 1.+0.j 0.+0.j 0.+0.j]
@@ -1143,7 +1151,56 @@ def __init__(
         )
 
     def get_compact_qobj(self):
-        return cphase(self.arg_value)
+        return cphase(self.arg_value).tidyup()
+
+
+class RZX(TwoQubitGate):
+    r"""
+    RZX gate.
+
+    .. math::
+
+        \begin{pmatrix}
+        \cos{\theta/2} & -i\sin{\theta/2} & 0 & 0 \\
+        -i\sin{\theta/2} & \cos{\theta/2} & 0 & 0 \\
+        0 & 0 & \cos{\theta/2} & i\sin{\theta/2} \\
+        0 & 0 & i\sin{\theta/2} & \cos{\theta/2} \\
+        \end{pmatrix}
+
+    Examples
+    --------
+    >>> from qutip_qip.operations import RZX
+    >>> RZX([0, 1], np.pi).get_compact_qobj().tidyup() # doctest: +NORMALIZE_WHITESPACE
+    Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=Dense, isherm=False
+    Qobj data =
+    [[0.+0.j 0.-1.j 0.+0.j 0.+0.j]
+    [0.-1.j 0.+0.j 0.+0.j 0.+0.j]
+    [0.+0.j 0.+0.j 0.+0.j 0.+1.j]
+    [0.+0.j 0.+0.j 0.+1.j 0.+0.j]]
+    """
+
+    def __init__(self, targets, arg_value, **kwargs):
+        self.target_gate = RZ
+        super().__init__(
+            targets=targets,
+            arg_value=arg_value,
+            target_gate=self.target_gate,
+            **kwargs,
+        )
+
+    def get_compact_qobj(self):
+        theta = self.arg_value
+        return Qobj(
+            np.array(
+                [
+                    [np.cos(theta / 2), -1.0j * np.sin(theta / 2), 0.0, 0.0],
+                    [-1.0j * np.sin(theta / 2), np.cos(theta / 2), 0.0, 0.0],
+                    [0.0, 0.0, np.cos(theta / 2), 1.0j * np.sin(theta / 2)],
+                    [0.0, 0.0, 1.0j * np.sin(theta / 2), np.cos(theta / 2)],
+                ]
+            ),
+            dims=[[2, 2], [2, 2]],
+        )
 
 
 CRY = partial(_OneControlledGate, target_gate=RY)
@@ -1198,4 +1255,5 @@ def get_compact_qobj(self):
     "CS": CS,
     "CT": CT,
     "CPHASE": CPHASE,
+    "RZX": RZX,
 }
diff --git a/src/qutip_qip/operations/gates.py b/src/qutip_qip/operations/gates.py
index 7a294f2f..5d078713 100644
--- a/src/qutip_qip/operations/gates.py
+++ b/src/qutip_qip/operations/gates.py
@@ -326,8 +326,8 @@ def snot(N=None, target=0):
     Examples
     --------
     >>> snot() # doctest: +SKIP
-    Quantum object: dims = [[2], [2]], \
-shape = [2, 2], type = oper, isHerm = True
+    Quantum object: dims=[[2], [2]], \
+shape = [2, 2], type='oper', dtype=Dense, isherm=True
     Qobj data =
     [[ 0.70710678+0.j  0.70710678+0.j]
      [ 0.70710678+0.j -0.70710678+0.j]]
@@ -356,8 +356,8 @@ def phasegate(theta, N=None, target=0):
     Examples
     --------
     >>> phasegate(pi/4) # doctest: +SKIP
-    Quantum object: dims = [[2], [2]], \
-shape = [2, 2], type = oper, isHerm = False
+    Quantum object: dims=[[2], [2]], \
+shape = [2, 2], type='oper', dtype=Dense, isherm=False
     Qobj data =
     [[ 1.00000000+0.j          0.00000000+0.j        ]
      [ 0.00000000+0.j          0.70710678+0.70710678j]]
@@ -502,8 +502,8 @@ def cnot(N=None, control=0, target=1):
     Examples
     --------
     >>> cnot() # doctest: +SKIP
-    Quantum object: dims = [[2, 2], [2, 2]], \
-shape = [4, 4], type = oper, isHerm = True
+    Quantum object: dims=[[2, 2], [2, 2]], \
+shape = [4, 4], type='oper', dtype=Dense, isherm=True
     Qobj data =
         [[ 1.+0.j  0.+0.j  0.+0.j  0.+0.j]
          [ 0.+0.j  1.+0.j  0.+0.j  0.+0.j]
@@ -535,8 +535,8 @@ def csign(N=None, control=0, target=1):
     Examples
     --------
     >>> csign() # doctest: +SKIP
-    Quantum object: dims = [[2, 2], [2, 2]], \
-shape = [4, 4], type = oper, isHerm = True
+    Quantum object: dims=[[2, 2], [2, 2]], \
+shape = [4, 4], type='oper', dtype=Dense, isherm=True
     Qobj data =
         [[ 1.+0.j  0.+0.j  0.+0.j  0.+0.j]
          [ 0.+0.j  1.+0.j  0.+0.j  0.+0.j]
@@ -568,8 +568,8 @@ def berkeley(N=None, targets=[0, 1]):
     Examples
     --------
     >>> berkeley() # doctest: +SKIP
-    Quantum object: dims = [[2, 2], [2, 2]], \
-shape = [4, 4], type = oper, isHerm = True
+    Quantum object: dims=[[2, 2], [2, 2]], \
+shape = [4, 4], type='oper', dtype=Dense, isherm=True
     Qobj data =
         [[ cos(pi/8).+0.j  0.+0.j           0.+0.j           0.+sin(pi/8).j]
          [ 0.+0.j          cos(3pi/8).+0.j  0.+sin(3pi/8).j  0.+0.j]
@@ -606,8 +606,8 @@ def swapalpha(alpha, N=None, targets=[0, 1]):
     Examples
     --------
     >>> swapalpha(alpha) # doctest: +SKIP
-    Quantum object: dims = [[2, 2], [2, 2]], \
-shape = [4, 4], type = oper, isHerm = True
+    Quantum object: dims=[[2, 2], [2, 2]], \
+shape = [4, 4], type='oper', dtype=Dense, isherm=True
     Qobj data =
     [[ 1.+0.j  0.+0.j                    0.+0.j                    0.+0.j]
      [ 0.+0.j  0.5*(1 + exp(j*pi*alpha)  0.5*(1 - exp(j*pi*alpha)  0.+0.j]
@@ -653,8 +653,8 @@ def swap(N=None, targets=[0, 1]):
     Examples
     --------
     >>> swap() # doctest: +SKIP
-    Quantum object: dims = [[2, 2], [2, 2]], \
-shape = [4, 4], type = oper, isHerm = True
+    Quantum object: dims=[[2, 2], [2, 2]], \
+shape = [4, 4], type='oper', dtype=Dense, isherm=True
     Qobj data =
     [[ 1.+0.j  0.+0.j  0.+0.j  0.+0.j]
      [ 0.+0.j  0.+0.j  1.+0.j  0.+0.j]
@@ -685,8 +685,8 @@ def iswap(N=None, targets=[0, 1]):
     Examples
     --------
     >>> iswap() # doctest: +SKIP
-    Quantum object: dims = [[2, 2], [2, 2]], \
-shape = [4, 4], type = oper, isHerm = False
+    Quantum object: dims=[[2, 2], [2, 2]], \
+shape = [4, 4], type='oper', dtype=Dense, isherm=False
     Qobj data =
     [[ 1.+0.j  0.+0.j  0.+0.j  0.+0.j]
      [ 0.+0.j  0.+0.j  0.+1.j  0.+0.j]
@@ -744,8 +744,8 @@ def sqrtiswap(N=None, targets=[0, 1]):
     Examples
     --------
     >>> sqrtiswap() # doctest: +SKIP
-    Quantum object: dims = [[2, 2], [2, 2]], \
-shape = [4, 4], type = oper, isHerm = False
+    Quantum object: dims=[[2, 2], [2, 2]], \
+shape = [4, 4], type='oper', dtype=Dense, isherm=False
     Qobj data =
     [[ 1.00000000+0.j   0.00000000+0.j   \
        0.00000000+0.j          0.00000000+0.j]
@@ -802,7 +802,7 @@ def molmer_sorensen(theta, phi=0.0, N=None, targets=[0, 1]):
     if N is not None:
         _deprecation_warnings_gate_expansion()
         return expand_operator(
-            molmer_sorensen(theta), dims=[2] * N, targets=targets
+            molmer_sorensen(theta, phi), dims=[2] * N, targets=targets
         )
 
     return Qobj(
@@ -842,8 +842,7 @@ def fredkin(N=None, control=0, targets=[1, 2]):
     Examples
     --------
     >>> fredkin() # doctest: +SKIP
-    Quantum object: dims = [[2, 2, 2], [2, 2, 2]], \
-shape = [8, 8], type = oper, isHerm = True
+    Quantum object: dims=[[2, 2, 2], [2, 2, 2]], shape = [8, 8], type='oper', dtype=Dense, isherm=True
     Qobj data =
         [[ 1.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j]
          [ 0.+0.j  1.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j]
@@ -889,8 +888,7 @@ def toffoli(N=None, controls=[0, 1], target=2):
     Examples
     --------
     >>> toffoli() # doctest: +SKIP
-    Quantum object: dims = [[2, 2, 2], [2, 2, 2]], \
-                    shape = [8, 8], type = oper, isHerm = True
+    Quantum object: dims=[[2, 2, 2], [2, 2, 2]], shape = [8, 8], type='oper', dtype=Dense, isherm=True
     Qobj data =
         [[ 1.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j]
          [ 0.+0.j  1.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j]
@@ -1018,8 +1016,8 @@ def globalphase(theta, N=1):
     Examples
     --------
     >>> phasegate(pi/4) # doctest: +SKIP
-    Quantum object: dims = [[2], [2]], \
-shape = [2, 2], type = oper, isHerm = False
+    Quantum object: dims=[[2], [2]], \
+shape = [2, 2], type='oper', dtype=Dense, isherm=False
     Qobj data =
     [[ 0.70710678+0.70710678j          0.00000000+0.j]
      [ 0.00000000+0.j          0.70710678+0.70710678j]]
@@ -1247,7 +1245,7 @@ def expand_operator(
     >>> from qutip_qip.operations import expand_operator, x_gate, cnot
     >>> import qutip
     >>> expand_operator(x_gate(), dims=[2,3], targets=[0]) # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2, 3], [2, 3]], shape = (6, 6), type = oper, isherm = True
+    Quantum object: dims=[[2, 3], [2, 3]], shape=(6, 6), type='oper', dtype=CSR, isherm=True
     Qobj data =
     [[0. 0. 0. 1. 0. 0.]
      [0. 0. 0. 0. 1. 0.]
@@ -1256,8 +1254,7 @@ def expand_operator(
      [0. 1. 0. 0. 0. 0.]
      [0. 0. 1. 0. 0. 0.]]
     >>> expand_operator(cnot(), dims=[2,2,2], targets=[1, 2]) # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2, 2, 2], [2, 2, 2]], shape = (8,
-    8), type = oper, isherm = True
+    Quantum object: dims=[[2, 2, 2], [2, 2, 2]], shape=(8, 8), type='oper', dtype=CSR, isherm=True
     Qobj data =
     [[1. 0. 0. 0. 0. 0. 0. 0.]
      [0. 1. 0. 0. 0. 0. 0. 0.]
@@ -1268,8 +1265,7 @@ def expand_operator(
      [0. 0. 0. 0. 0. 0. 0. 1.]
      [0. 0. 0. 0. 0. 0. 1. 0.]]
     >>> expand_operator(cnot(), dims=[2, 2, 2], targets=[2, 0]) # doctest: +NORMALIZE_WHITESPACE
-    Quantum object: dims = [[2, 2, 2], [2, 2, 2]], shape = (8,
-    8), type = oper, isherm = True
+    Quantum object: dims=[[2, 2, 2], [2, 2, 2]], shape=(8, 8), type='oper', dtype=CSR, isherm=True
     Qobj data =
     [[1. 0. 0. 0. 0. 0. 0. 0.]
      [0. 0. 0. 0. 0. 1. 0. 0.]
diff --git a/src/qutip_qip/pulse.py b/src/qutip_qip/pulse.py
index 675593cd..8497518c 100644
--- a/src/qutip_qip/pulse.py
+++ b/src/qutip_qip/pulse.py
@@ -417,15 +417,14 @@ def get_noisy_qobjevo(self, dims):
             noise.get_qobjevo(self.spline_kind, dims)
             for noise in self.coherent_noise
         ]
-        qu = _merge_qobjevo([ideal_qu] + noise_qu_list)
+        if parse_version(qutip.__version__) < parse_version("5.dev"):
+            qu = _merge_qobjevo([ideal_qu] + noise_qu_list)
+        else:
+            qu = sum(noise_qu_list, ideal_qu)
         c_ops = [
             noise.get_qobjevo(self.spline_kind, dims)
             for noise in self.lindblad_noise
         ]
-        full_tlist = self.get_full_tlist()
-        qu = _merge_qobjevo([qu], full_tlist)
-        for i, c_op in enumerate(c_ops):
-            c_ops[i] = _merge_qobjevo([c_op], full_tlist)
         return qu, c_ops
 
     def get_full_tlist(self, tol=1.0e-10):


Test patch:
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 3a375928..db96bc30 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -12,10 +12,9 @@ jobs:
     strategy:
       matrix:
         include:
-          - case-name: qutip@5.0.0
+          - case-name: qutip@5
             os: ubuntu-latest
-            qutip-version: '==5.0.0'
-            pyqir-version: ''
+            qutip-version: '==5.*'
             python-version: '3.11'
           - case-name: qutip@qutip-4.7.X
             os: ubuntu-latest
@@ -39,12 +38,6 @@ jobs:
             qiskit-version: '==0.46.*'
             pyqir-version: '==0.6.2'
             python-version: '3.9'
-          - case-name: qiskit+qir
-            os: macOS-latest
-            qutip-version: ''
-            qiskit-version: '==0.46.*'
-            pyqir-version: '==0.6.2'
-            python-version: '3.9'
           - case-name: qiskit+qir
             os: ubuntu-latest
             qutip-version: ''
@@ -53,9 +46,9 @@ jobs:
             python-version: '3.9'
 
     steps:
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - name: Set up Python
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
 
@@ -63,24 +56,30 @@ jobs:
       if: ${{ matrix.qutip-version != '' && ! startsWith( matrix.qutip-version, '@') }}
       run: python -m pip install 'qutip${{ matrix.qutip-version }}'
 
+    - name: Install PyQIR from PyPI
+      if: ${{ matrix.pyqir-version != '' }}
+      # We use each subpackage explicitly here; see https://github.com/qir-alliance/pyqir/issues/167.
+      run: python -m pip install 'pyqir-generator${{ matrix.pyqir-version }}' 'pyqir-parser${{ matrix.pyqir-version }}'
+
     - name: Install QuTiP from GitHub
       if: ${{ startsWith( matrix.qutip-version, '@') }}
       run: |
         python -m pip install 'git+https://github.com/qutip/qutip.git${{ matrix.qutip-version }}'
+        
+    - name: Install qutip-qtrl from GitHub
+      if: ${{ matrix.qutip-version == '@master' }}
+      run: |
+        python -m pip install 'git+https://github.com/qutip/qutip-qtrl.git${{ matrix.qutip-version }}'
 
     - name: Install Qiskit from PyPI
       if: ${{ matrix.qiskit-version != '' }}
       run: python -m pip install 'qiskit${{ matrix.qiskit-version }}' 'qiskit-aer==0.14.0.1'
 
-    - name: Install PyQIR from PyPI
-      if: ${{ matrix.pyqir-version != '' }}
-      # We use each subpackage explicitly here; see https://github.com/qir-alliance/pyqir/issues/167.
-      run: python -m pip install 'pyqir-generator${{ matrix.pyqir-version }}' 'pyqir-parser${{ matrix.pyqir-version }}'
-
     - name: Install qutip-qip
       # Installing in-place so that coveralls can locate the source code.
       run: |
-        pip install -e .[full]
+        pip install matplotlib pytest
+        pip install -e .
     - name: Test with pytest and generate coverage report
       run: |
         pip install pytest-cov coveralls
@@ -95,16 +94,15 @@ jobs:
     runs-on: ubuntu-latest
 
     steps:
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - name: Set up Python
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
-        python-version: 3.9
+        python-version: 3.11
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
         python -mpip install -r doc/requirements.txt
-        pip install qutip==4.7.3
         pip install .
     - name: Test code snippets in the documentation
       run: |
diff --git a/tests/test_circuit.py b/tests/test_circuit.py
index 98032bb2..019ba20a 100644
--- a/tests/test_circuit.py
+++ b/tests/test_circuit.py
@@ -156,6 +156,7 @@ def test_add_gate(self):
         qc.add_gate("SNOT", targets=[3])
         qc.add_gate(test_gate, index=[3])
         qc.add_1q_gate("RY", start=4, end=5, arg_value=1.570796)
+        qc.add_1q_gate("RX", qubits=[3], arg_value=-1.570796)
 
         # Test explicit gate addition
         assert qc.gates[0].name == "CNOT"
@@ -176,7 +177,12 @@ def test_add_gate(self):
         assert qc.gates[5].arg_value == 1.570796
         assert qc.gates[6].name == "RY"
         assert qc.gates[6].targets == [5]
-        assert qc.gates[5].arg_value == 1.570796
+        assert qc.gates[6].arg_value == 1.570796
+
+        # Test adding 1 qubit gate on qubits [3] 
+        assert qc.gates[7].name == "RX"
+        assert qc.gates[7].targets == [3]
+        assert qc.gates[7].arg_value == -1.570796
 
         dummy_gate1 = Gate("DUMMY1")
         inds = [1, 3, 4, 6]
@@ -197,6 +203,7 @@ def test_add_gate(self):
             'RY',       # 8
             'DUMMY1',   # 9
             'RY',       # 10
+            'RX',       # 11
         ]
         actual_gate_names = [gate.name for gate in qc.gates]
         assert actual_gate_names == expected_gate_names
@@ -220,6 +227,7 @@ def test_add_gate(self):
             'DUMMY1',   # 10
             'RY',       # 11
             'DUMMY2',   # 12
+            'RX',       # 13
         ]
         actual_gate_names = [gate.name for gate in qc.gates]
         assert actual_gate_names == expected_gate_names
diff --git a/tests/test_device.py b/tests/test_device.py
index 6170f2eb..96117185 100644
--- a/tests/test_device.py
+++ b/tests/test_device.py
@@ -8,7 +8,7 @@
 
 import qutip
 from qutip_qip.circuit import QubitCircuit
-from qutip_qip.operations import Gate, gate_sequence_product
+from qutip_qip.operations import Gate, gate_sequence_product, RZX
 from qutip_qip.device import (DispersiveCavityQED, LinearSpinChain,
                                 CircularSpinChain, SCQubits)
 
@@ -101,8 +101,21 @@ def test_analytical_evolution(num_qubits, gates, device_class, kwargs):
 @pytest.mark.filterwarnings("ignore:Not in the dispersive regime")
 @pytest.mark.parametrize(("num_qubits", "gates"), single_gate_tests)
 @pytest.mark.parametrize(("device_class", "kwargs"), device_lists_numeric)
-def test_numerical_evolution(
-    num_qubits, gates, device_class, kwargs):
+def test_numerical_evolution(num_qubits, gates, device_class, kwargs):
+    _test_numerical_evolution_helper(num_qubits, gates, device_class, kwargs)
+
+
+# Test for RZX gate, only available on SCQubits.
+_rzx = RZX([0, 1], arg_value=np.pi/2)
+@pytest.mark.parametrize(
+    ("num_qubits", "gates", "device_class", "kwargs"), 
+    [pytest.param(2, [_rzx], SCQubits, {}, id="RZX-SCQubits")]
+    )
+def test_numerical_evolution_zx(num_qubits, gates, device_class, kwargs):
+    _test_numerical_evolution_helper(num_qubits, gates, device_class, kwargs)
+
+
+def _test_numerical_evolution_helper(num_qubits, gates, device_class, kwargs):
     num_qubits = 2
     circuit = QubitCircuit(num_qubits)
     for gate in gates:
diff --git a/tests/test_gates.py b/tests/test_gates.py
index 49938105..7a2ac6e0 100644
--- a/tests/test_gates.py
+++ b/tests/test_gates.py
@@ -10,7 +10,7 @@
 from qutip_qip.operations import (
     X, Y, Z, RX, RY, RZ, H, SQRTNOT, S, T, QASMU, CNOT, CPHASE, ISWAP, SWAP,
     CZ, SQRTSWAP, SQRTISWAP, SWAPALPHA, SWAPALPHA, MS, TOFFOLI, FREDKIN,
-    BERKELEY, R, expand_operator)
+    BERKELEY, R, RZX, expand_operator)
 
 
 def _permutation_id(permutation):
@@ -382,6 +382,7 @@ def test_gates_class():
     circuit1.add_gate("TOFFOLI", [2, 0, 1])
     circuit1.add_gate("FREDKIN", [0, 1, 2])
     circuit1.add_gate("BERKELEY", [1, 0])
+    circuit1.add_gate("RZX", [1, 0], arg_value=1.)
     result1 = circuit1.run(init_state)
 
     circuit2 = QubitCircuit(3)
@@ -409,6 +410,7 @@ def test_gates_class():
     circuit2.add_gate(TOFFOLI([2, 0, 1]))
     circuit2.add_gate(FREDKIN([0, 1, 2]))
     circuit2.add_gate(BERKELEY([1, 0]))
+    circuit2.add_gate(RZX([1, 0], 1.))
     result2 = circuit2.run(init_state)
 
     assert pytest.approx(qutip.fidelity(result1, result2), 1.0e-6) == 1
diff --git a/tests/test_optpulseprocessor.py b/tests/test_optpulseprocessor.py
index b5fc9dd9..674c3c41 100644
--- a/tests/test_optpulseprocessor.py
+++ b/tests/test_optpulseprocessor.py
@@ -1,21 +1,30 @@
 import os
 
-from numpy.testing import (assert_, assert_allclose,
-                           assert_equal)
+from numpy.testing import assert_, assert_allclose, assert_equal
 import numpy as np
 import pytest
 
+pytest.importorskip("qutip_qtrl")
 from qutip_qip.device import OptPulseProcessor, SpinChainModel
 from qutip_qip.circuit import QubitCircuit
 from qutip_qip.qubits import qubit_states
 import qutip
-from qutip import (fidelity, Qobj, tensor, rand_ket, basis,  sigmaz,
-                    sigmax, sigmay, identity, destroy)
+from qutip import (
+    fidelity,
+    tensor,
+    rand_ket,
+    sigmaz,
+    sigmax,
+    sigmay,
+    identity,
+)
 from qutip_qip.operations import (
-    cnot, hadamard_transform, expand_operator, gate_sequence_product)
+    cnot,
+)
 import qutip
 from packaging.version import parse as parse_version
-if parse_version(qutip.__version__) < parse_version('5.dev'):
+
+if parse_version(qutip.__version__) < parse_version("5.dev"):
     from qutip import Options as SolverOptions
 else:
     from qutip import SolverOptions
@@ -38,7 +47,8 @@ def test_simple_hadamard(self):
         test = OptPulseProcessor(N, drift=H_d)
         test.add_control(H_c, targets=0)
         tlist, coeffs = test.load_circuit(
-            qc, num_tslots=num_tslots, evo_time=evo_time, verbose=True)
+            qc, num_tslots=num_tslots, evo_time=evo_time, verbose=True
+        )
 
         # test run_state
         rho0 = qubit_states(1, [0])
@@ -51,7 +61,7 @@ def test_multi_qubits(self):
         Test for multi-qubits system.
         """
         N = 3
-        H_d = tensor([sigmaz()]*3)
+        H_d = tensor([sigmaz()] * 3)
         H_c = []
 
         # test empty ctrls
@@ -59,28 +69,27 @@ def test_multi_qubits(self):
         evo_time = 10
         test = OptPulseProcessor(N)
         test.add_drift(H_d, [0, 1, 2])
-        test.add_control(tensor([sigmax(), sigmax()]),
-                          cyclic_permutation=True)
+        test.add_control(tensor([sigmax(), sigmax()]), cyclic_permutation=True)
         # test periodically adding ctrls
         sx = sigmax()
         iden = identity(2)
-        assert(len(test.get_control_labels()) == 3)
+        assert len(test.get_control_labels()) == 3
         test.add_control(sigmax(), cyclic_permutation=True)
         test.add_control(sigmay(), cyclic_permutation=True)
 
         # test pulse genration for cnot gate, with kwargs
         qc = [tensor([identity(2), cnot()])]
-        test.load_circuit(qc, num_tslots=num_tslots,
-                          evo_time=evo_time, min_fid_err=1.0e-6)
+        test.load_circuit(
+            qc, num_tslots=num_tslots, evo_time=evo_time, min_fid_err=1.0e-6
+        )
         rho0 = qubit_states(3, [1, 1, 1])
         rho1 = qubit_states(3, [1, 1, 0])
-        result = test.run_state(
-            rho0, options=SolverOptions(store_states=True))
-        assert_(fidelity(result.states[-1], rho1) > 1-1.0e-6)
+        result = test.run_state(rho0, options=SolverOptions(store_states=True))
+        assert_(fidelity(result.states[-1], rho1) > 1 - 1.0e-6)
 
     def test_multi_gates(self):
         N = 2
-        H_d = tensor([sigmaz()]*2)
+        H_d = tensor([sigmaz()] * 2)
         H_c = []
 
         test = OptPulseProcessor(N)
@@ -90,22 +99,23 @@ def test_multi_gates(self):
         test.add_control(tensor([sigmay(), sigmay()]))
 
         # qubits circuit with 3 gates
-        setting_args = {"SNOT": {"num_tslots": 10, "evo_time": 1},
-                        "SWAP": {"num_tslots": 30, "evo_time": 3},
-                        "CNOT": {"num_tslots": 30, "evo_time": 3}}
+        setting_args = {
+            "SNOT": {"num_tslots": 10, "evo_time": 1},
+            "SWAP": {"num_tslots": 30, "evo_time": 3},
+            "CNOT": {"num_tslots": 30, "evo_time": 3},
+        }
         qc = QubitCircuit(N)
         qc.add_gate("SNOT", 0)
         qc.add_gate("SWAP", targets=[0, 1])
-        qc.add_gate('CNOT', controls=1, targets=[0])
-        test.load_circuit(qc, setting_args=setting_args,
-                          merge_gates=False)
+        qc.add_gate("CNOT", controls=1, targets=[0])
+        test.load_circuit(qc, setting_args=setting_args, merge_gates=False)
 
         rho0 = rand_ket(4)  # use random generated ket state
         rho0.dims = [[2, 2], [1, 1]]
         U = qc.compute_unitary()
         rho1 = U * rho0
         result = test.run_state(rho0)
-        assert_(fidelity(result.states[-1], rho1) > 1-1.0e-6)
+        assert_(fidelity(result.states[-1], rho1) > 1 - 1.0e-6)
 
     def test_with_model(self):
         model = SpinChainModel(3, setup="linear")
@@ -117,7 +127,7 @@ def test_with_model(self):
             qc, merge_gates=True, num_tslots=10, evo_time=2.0
         )
 
-        if parse_version(qutip.__version__) < parse_version('5.dev'):
+        if parse_version(qutip.__version__) < parse_version("5.dev"):
             init_state = qutip.rand_ket(8, dims=[[2, 2, 2], [1, 1, 1]])
         else:
             init_state = qutip.rand_ket([2, 2, 2])