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

Issue:
Links main repo in docs html
Fixes #57 

Adds a link to the main repo on each page of the documentation. 
![image](https://user-images.githubusercontent.com/66048318/122238578-dbe32200-ce85-11eb-828f-9e5e515e2b27.png)

I think the two colors don't match exactly. I did use the same color `#6d6d6d` for the octocat figure but the [`fork_me` ribbon ](https://github.blog/2008-12-19-github-ribbons/)appears a bit darker.



Discussion / hints:

Fix patch:
diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
new file mode 100644
index 00000000..f07470db
--- /dev/null
+++ b/.github/workflows/build.yml
@@ -0,0 +1,148 @@
+name: Build wheels, optionally deploy to PyPI
+
+on:
+  workflow_dispatch:
+    inputs:
+      confirm_ref:
+        description: "Confirm chosen branch name to deploy to PyPI (optional):"
+        default: ""
+      override_version:
+        description: "Override version number (optional):"
+        default: ""
+
+
+jobs:
+  # The deploy_test job is part of the test of whether we should deploy to PyPI.
+  # The job will succeed if either the confirmation reference is empty or if the
+  # confirmation is the selected branch or tag name.  It will fail if it is
+  # nonempty and does not match.  All later jobs depend on this job, so that
+  # they will be immediately cancelled if the confirmation is bad.  The
+  # dependency is currently necessary (2021-03) because GitHub Actions does not
+  # have a simpler method of cancelling an entire workflow---the normal use-case
+  # expects to try and run as much as possible despite one or two failures.
+  deploy_test:
+    name: Verify PyPI deployment confirmation
+    runs-on: ubuntu-latest
+    env:
+      GITHUB_REF: ${{ github.ref }}
+      CONFIRM_REF: ${{ github.event.inputs.confirm_ref }}
+    steps:
+      - name: Compare confirmation to current reference
+        shell: bash
+        run: |
+          [[ -z $CONFIRM_REF || $GITHUB_REF =~ ^refs/(heads|tags)/$CONFIRM_REF$ ]]
+          if [[ -z $CONFIRM_REF ]]; then
+            echo "Build only.  Nothing will be uploaded to PyPI."
+          else
+            echo "Full build and deploy.  Wheels and source will be uploaded to PyPI."
+          fi
+
+
+  build_sdist:
+    name: Build sdist on Ubuntu
+    needs: deploy_test
+    runs-on: ubuntu-latest
+    env:
+      OVERRIDE_VERSION: ${{ github.event.inputs.override_version }}
+
+    steps:
+      - uses: actions/checkout@v2
+
+      - uses: actions/setup-python@v2
+        name: Install Python
+        with:
+          # For the sdist we should be as conservative as possible with our
+          # Python version.  This should be the lowest supported version.  This
+          # means that no unsupported syntax can sneak through.
+          python-version: '3.6'
+
+      - name: Install pip build
+        run: |
+          python -m pip install 'build'
+
+      - name: Build sdist tarball
+        shell: bash
+        run: |
+          if [[ ! -z "$OVERRIDE_VERSION" ]]; then echo "$OVERRIDE_VERSION" > VERSION; fi
+          # The build package is the reference PEP 517 package builder.  All
+          # dependencies are specified by our setup code.
+          python -m build --sdist .
+
+      - uses: actions/upload-artifact@v2
+        with:
+          name: sdist
+          path: dist/*.tar.gz
+          if-no-files-found: error
+
+
+  build_wheels:
+    name: Build wheels
+    needs: deploy_test
+    runs-on: ubuntu-latest
+    env:
+      OVERRIDE_VERSION: ${{ github.event.inputs.override_version }}
+
+    steps:
+      - uses: actions/checkout@v2
+
+      - uses: actions/setup-python@v2
+        name: Install Python
+        with:
+          # This is about the build environment, not the released wheel version.
+          python-version: '3.7'
+
+      - name: Install pip build
+        run: |
+          python -m pip install 'build'
+
+      - name: Build wheels
+        shell: bash
+        run: |
+          # If the version override was specified, then write it the VERSION
+          # file with it.
+          if [[ ! -z "$OVERRIDE_VERSION" ]]; then echo "$OVERRIDE_VERSION" > VERSION; fi
+          # The build package is the reference PEP 517 package builder.  All
+          # dependencies are specified by our setup code.
+          python -m build --wheel --outdir wheelhouse .
+
+      - uses: actions/upload-artifact@v2
+        with:
+          name: wheels
+          path: ./wheelhouse/*.whl
+          if-no-files-found: error
+
+
+  deploy:
+    name: "Deploy to PyPI if desired"
+    # The confirmation is tested explicitly in `deploy_test`, so we know it is
+    # either a missing confirmation (so we shouldn't run this job) or a valid
+    # confirmation.  We don't need to retest the value of the confirmation,
+    # beyond checking that one existed.
+    if: ${{ github.event.inputs.confirm_ref != '' }}
+    needs: [deploy_test, build_sdist, build_wheels]
+    runs-on: ubuntu-latest
+
+    steps:
+      - name: Download build artifacts to local runner
+        uses: actions/download-artifact@v2
+
+      - uses: actions/setup-python@v2
+        name: Install Python
+        with:
+          python-version: '3.7'
+
+      - name: Verify this is not a dev version
+        shell: bash
+        run: |
+          python -m pip install wheels/*.whl
+          python -c 'import qutip_qip; print(qutip_qip.__version__); assert "dev" not in qutip_qip.__version__; assert "+" not in qutip_qip.__version__'
+
+      - name: Upload sdist and wheels to PyPI
+        env:
+          TWINE_USERNAME: __token__
+          TWINE_PASSWORD: ${{ secrets.PYPI_API_TOKEN }}
+          TWINE_NON_INTERACTIVE: 1
+          TWINE_REPOSITORY: pypi
+        run: |
+          python -m pip install "twine"
+          python -m twine upload --verbose wheels/*.whl sdist/*.tar.gz
diff --git a/.github/workflows/build_documentation.yml b/.github/workflows/build_documentation.yml
new file mode 100644
index 00000000..1a427ef5
--- /dev/null
+++ b/.github/workflows/build_documentation.yml
@@ -0,0 +1,46 @@
+name: Build HTML documentation
+
+on:
+  [push, pull_request]
+
+jobs:
+  build:
+    name: Build documentation
+    runs-on: ubuntu-latest
+
+    steps:
+      - uses: actions/checkout@v2
+
+      - uses: actions/setup-python@v2
+        name: Install Python
+        with:
+          python-version: '3.8'
+
+      - name: Install documentation dependencies
+        run: |
+          python -mpip install -r doc/requirements.txt
+
+      - name: Install qutip-qip from GitHub
+        run: |
+          python -mpip install -e .[full]
+          # Install in editable mode so it doesn't matter if we import from
+          # inside the installation directory, otherwise we can get some errors
+          # because we're importing from the wrong location.
+          python -c 'import qutip_qip; print("QuTiP QIP Version:      %s" % qutip_qip.__version__)'
+          python -c 'import qutip; qutip.about()'
+
+      - name: Build documentation
+        working-directory: doc
+        run: |
+          make html SPHINXOPTS="-W --keep-going -T"
+          # Above flags are:
+          #   -W : turn warnings into errors
+          #   --keep-going : do not stop after the first error
+          #   -T : display a full traceback if a Python exception occurs
+
+      - name: Upload built files
+        uses: actions/upload-artifact@v2
+        with:
+          name: qutip_qip_html_docs
+          path: doc/_build/html/*
+          if-no-files-found: error
diff --git a/.gitignore b/.gitignore
index baab02aa..d775c834 100644
--- a/.gitignore
+++ b/.gitignore
@@ -71,7 +71,7 @@ instance/
 .scrapy
 
 # Sphinx documentation
-docs/_build/
+doc/_build/
 
 # PyBuilder
 target/
diff --git a/.readthedocs.yaml b/.readthedocs.yaml
index cbdf05c8..9324641d 100644
--- a/.readthedocs.yaml
+++ b/.readthedocs.yaml
@@ -2,14 +2,14 @@
 
 version: 2
 
-# Build from the docs/ directory with Sphinx
+# Build from the doc/ directory with Sphinx
 sphinx:
-  configuration: docs/source/conf.py
+  configuration: doc/source/conf.py
 
 # Explicitly set the version of Python and its requirements
 python:
   version: 3.8
   install:
-    - requirements: docs/requirements.txt
+    - requirements: doc/requirements.txt
     - method: pip
       path: .
diff --git a/README.md b/README.md
old mode 100755
new mode 100644
index 70318f5e..0acf0173
--- a/README.md
+++ b/README.md
@@ -2,6 +2,8 @@
 
 [![build](https://github.com/qutip/qutip-qip/workflows/Tests/badge.svg)](https://github.com/qutip/qutip-qip/actions)
 [![Documentation Status](https://readthedocs.org/projects/qutip-qip/badge/?version=latest)](https://qutip-qip.readthedocs.io/en/latest/)
+[![PyPI version](https://badge.fury.io/py/qutip-qip.svg)](https://badge.fury.io/py/qutip-qip)
+[![arXiv paper](https://img.shields.io/badge/arXiv-2105.09902-<COLOR>.svg)](https://arxiv.org/abs/2105.09902)
 
 The qutip-qip package used to be a module ``qutip.qip`` under [QuTiP (Quantum Toolbox in Python)](http://qutip.org/index.html).
 From QuTiP 5.0, the community has decided to decrease the size of the core QuTiP package by reducing the external dependencies, in order to simplify maintenance.
@@ -14,25 +16,42 @@ The package offers two different approaches for simulating quantum circuits, one
 
 If you would like to know the future development plan and ideas, have a look at the [discussion panel](https://github.com/qutip/qutip-qip/discussions) as well as the [qutip documentation for ideas](https://github.com/qutip/qutip-doc/tree/master/development/ideas).
 
-Installation
-------------
-To install the package, download to source code and run
+Quick start
+-----------
+To install the package, use
 ```
-pip install qutip_qip
+pip install qutip-qip
 ```
 
-If you want to edit the source code, please download the source code and run the following command under the folder with `setup.cfg`
+Migrating from ``qutip.qip``
+--------------------------
+As the introduction suggested, this package is based on a module in the [QuTiP](http://qutip.org/docs/latest/) package `qutip.qip`.
+If you were using the `qutip` package and now want to try out the new features included in this package, you can simply install this package and replace all the `qutip.qip` in your import statement with `qutip_qip`. Everything should work smoothly as usual.
+
+Documentation and tutorials
+-------------
+
+The documentation of `qutip-qip` updated to the latest development version is hosted at [qutip-qip.readthedocs.io/](https://qutip-qip.readthedocs.io/en/latest/).
+Tutorials related to using quantum gates and circuits in `qutip-qip` can be found [*here*](https://qutip.org/tutorials#quantum-information-processing) and those related to using noise simulators areavailable at [*this link*](https://qutip.org/tutorials#nisq). 
+
+Code examples used in the preprint [*Pulse-level noisy quantum circuits with QuTiP*](https://arxiv.org/abs/2105.09902), updated for the latest code version, are hosted in [this folder](https://github.com/qutip/qutip-qip/tree/master/doc/pulse-paper).
+
+Installation from source
+------------------------
+If you want to edit the source code, please download the source code and run the following command under the root `qutip-qip` folder,
 ```
+pip install --upgrade pip
 pip install -e .
 ```
+which makes sure that you are up to date with the latest `pip` version. Contribution guidelines are available [*here*](https://qutip-qip.readthedocs.io/en/latest/contribution-code.html). 
 
 To build and test the documentation, additional packages need to be installed:
 
 ```
-pip install matplotlib sphinx numpydoc sphinx_rtd_theme
+pip install pytest matplotlib sphinx numpydoc sphinx_rtd_theme
 ```
 
-Under the `docs` directory, use
+Under the `doc` directory, use
 ```
 make html
 ```
@@ -42,18 +61,26 @@ make doctest
 ```
 to test the code in the documentation.
 
-Documentation
--------------
-
-The documentation of `qutip-qip` updated to the latest development version is hosted at [qutip-qip.readthedocs.io/](https://qutip-qip.readthedocs.io/en/latest/).
-
 Testing
 ------------
 To test the installation from a download of the source code, run from the `qutip-qip` directory
 ```
-    pytest tests
+pytest tests
+```
+Citing `qutip-qip`
+------------
+If you use `qutip-qip` in your research, the  [*preprint*](https://arxiv.org/abs/2105.09902) can be cited as
+
+```bibtex
+@misc{li2021pulselevel,
+      title={Pulse-level noisy quantum circuits with QuTiP},
+      author={Boxi Li and Shahnawaz Ahmed and Sidhant Saraogi and Neill Lambert and Franco Nori and Alexander Pitchford and Nathan Shammah},
+      year={2021},
+      eprint={2105.09902},
+      archivePrefix={arXiv},
+      primaryClass={quant-ph}
+}
 ```
-
 Support
 -------
 This package is supported and maintained by the same developers group as QuTiP.
diff --git a/VERSION b/VERSION
new file mode 100644
index 00000000..17e51c38
--- /dev/null
+++ b/VERSION
@@ -0,0 +1,1 @@
+0.1.1
diff --git a/docs/Makefile b/doc/Makefile
similarity index 96%
rename from docs/Makefile
rename to doc/Makefile
index d0c3cbf1..92dd33a1 100644
--- a/docs/Makefile
+++ b/doc/Makefile
@@ -6,7 +6,7 @@
 SPHINXOPTS    ?=
 SPHINXBUILD   ?= sphinx-build
 SOURCEDIR     = source
-BUILDDIR      = build
+BUILDDIR      = _build
 
 # Put it first so that "make" without argument is like "make help".
 help:
diff --git a/docs/make.bat b/doc/make.bat
similarity index 97%
rename from docs/make.bat
rename to doc/make.bat
index 9534b018..5001b3f6 100644
--- a/docs/make.bat
+++ b/doc/make.bat
@@ -8,7 +8,7 @@ if "%SPHINXBUILD%" == "" (
 	set SPHINXBUILD=sphinx-build
 )
 set SOURCEDIR=source
-set BUILDDIR=build
+set BUILDDIR=_build
 
 if "%1" == "" goto help
 
diff --git a/doc/pulse-paper/README.md b/doc/pulse-paper/README.md
new file mode 100644
index 00000000..a66ff71d
--- /dev/null
+++ b/doc/pulse-paper/README.md
@@ -0,0 +1,15 @@
+This folder contains code examples used in the publication [*Pulse-level noisy quantum circuits with QuTiP*](https://arxiv.org/abs/2105.09902). To run the examples, please first install the software package qutip-qip
+```
+pip install qutip_qip[full] joblib
+```
+All examples are self-contained and running the code should reproduce the plots used in the paper.
+
+The following table summarizes the sections in the paper and the corresponding code examples:
+
+| Section | Code example |
+| ----------- | ----------- |
+| Appendix A | `pulse_simulation.py` |
+| Appendix B | `dj_algorithm.py` |
+| Fig.4 and Appendix C | `customize.py` |
+| Fig.5 in Section 4 | `decoherence.py` |
+| Appendix D | `deutsch_jozsa.qasm` and `deutsch_jozsa-qasm.py` |
diff --git a/doc/pulse-paper/customize.py b/doc/pulse-paper/customize.py
new file mode 100644
index 00000000..45a14a7a
--- /dev/null
+++ b/doc/pulse-paper/customize.py
@@ -0,0 +1,287 @@
+TEXTWIDTH = 7.1398920714
+LINEWIDTH = 3.48692403487
+
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+
+from scipy.optimize import curve_fit
+
+try:
+    global_setup(fontsize=10)
+except:
+    pass
+from joblib import Parallel, delayed  # for parallel simulations
+
+import numpy as np
+from qutip import sigmax, sigmay, sigmaz, basis, qeye, tensor, Qobj, fock_dm
+from qutip_qip.circuit import QubitCircuit, Gate
+from qutip_qip.device import ModelProcessor
+from qutip_qip.compiler import GateCompiler, Instruction
+from qutip import Options
+from qutip_qip.noise import Noise
+
+
+class MyProcessor(ModelProcessor):
+    """Custom processor built using ModelProcessor as the base class.
+
+    This custom processor will inherit all the methods of the base class 
+    such as setting up of the T1 and T2 decoherence rates in the simulations.
+
+    In addition, it is possible to write your own functions to add control
+    pulses.
+
+    Args:
+        num_qubits (int): Number of qubits in the processor.
+        t1, t2 (float or list): The T1 and T2 decoherence rates for the
+                                qubit. If it is a list, then it is assumed that
+                                each element of the list corresponds to the rates
+                                for each of the qubits.
+    """
+
+    def __init__(self, num_qubits, t1=None, t2=None):
+        super().__init__(num_qubits, t1=t1, t2=t2)
+        self.pulse_mode = "discrete"  # set the control pulse as discrete or continuous.
+        self.set_up_ops()  # set up the available Hamiltonians.
+        self.dims = [2] * num_qubits  # dimension of the quantum system.
+        self.num_qubits = num_qubits
+        self.native_gates = ["RX", "RY"]
+
+    def set_up_ops(self):
+        """Sets up the single qubit control operators for each qubit."""
+        for m in range(self.num_qubits):
+            self.add_control(2 * np.pi * sigmax() / 2, m, label="sx" + str(m))
+        for m in range(self.num_qubits):
+            self.add_control(2 * np.pi * sigmay() / 2, m, label="sy" + str(m))
+
+
+class MyCompiler(GateCompiler):
+    """Custom compiler for generating pulses from gates using the base class 
+    GateCompiler.
+
+    Args:
+        num_qubits (int): The number of qubits in the processor
+        params (dict): A dictionary of parameters for gate pulses such as
+                       the pulse amplitude.
+    """
+
+    def __init__(self, num_qubits, params):
+        super().__init__(num_qubits, params=params)
+        self.params = params
+        self.gate_compiler = {
+            "ROT": self.rotation_with_phase_compiler,
+            "RX": self.single_qubit_gate_compiler,
+            "RY": self.single_qubit_gate_compiler,
+        }
+
+    def generate_pulse(self, gate, tlist, coeff, phase=0.0):
+        """Generates the pulses.
+
+        Args:
+            gate (qutip_qip.circuit.Gate): A qutip Gate object.
+            tlist (array): A list of times for the evolution.
+            coeff (array): An array of coefficients for the gate pulses
+            phase (float): The value of the phase for the gate.
+
+        Returns:
+            Instruction (qutip_qip.compiler.instruction.Instruction): An instruction
+            to implement a gate containing the control pulses.                                               
+        """
+        pulse_info = [
+            ("sx" + str(gate.targets[0]), np.cos(phase) * coeff),
+            ("sy" + str(gate.targets[0]), np.sin(phase) * coeff),
+        ]
+        return [Instruction(gate, tlist, pulse_info)]
+
+    def single_qubit_gate_compiler(self, gate, args):
+        """Compiles single qubit gates to pulses.
+        
+        Args:
+            gate (qutip_qip.circuit.Gate): A qutip Gate object.
+        
+        Returns:
+            Instruction (qutip_qip.compiler.instruction.Instruction): An instruction
+            to implement a gate containing the control pulses.
+        """
+        # gate.arg_value is the rotation angle
+        tlist = np.abs(gate.arg_value) / self.params["pulse_amplitude"]
+        coeff = self.params["pulse_amplitude"] * np.sign(gate.arg_value)
+        if gate.name == "RX":
+            return self.generate_pulse(gate, tlist, coeff, phase=0.0)
+        elif gate.name == "RY":
+            return self.generate_pulse(gate, tlist, coeff, phase=np.pi / 2)
+
+    def rotation_with_phase_compiler(self, gate, args):
+        """Compiles gates with a phase term.
+
+        Args:
+            gate (qutip_qip.circuit.Gate): A qutip Gate object.
+        
+        Returns:
+            Instruction (qutip_qip.compiler.instruction.Instruction): An instruction
+            to implement a gate containing the control pulses.
+        """
+        # gate.arg_value is the pulse phase
+        tlist = self.params["duration"]
+        coeff = self.params["pulse_amplitude"]
+        return self.generate_pulse(gate, tlist, coeff, phase=gate.arg_value)
+
+
+# Define a circuit and run the simulation
+num_qubits = 1
+
+circuit = QubitCircuit(1)
+circuit.add_gate("RX", targets=0, arg_value=np.pi / 2)
+circuit.add_gate("Z", targets=0)
+
+myprocessor = MyProcessor(1)
+
+mycompiler = MyCompiler(num_qubits, {"pulse_amplitude": 0.02})
+
+myprocessor.load_circuit(circuit, compiler=mycompiler)
+result = myprocessor.run_state(basis(2, 0))
+
+fig, ax = myprocessor.plot_pulses(
+    figsize=(LINEWIDTH * 0.7, LINEWIDTH / 2 * 0.7), dpi=200
+)
+ax[-1].set_xlabel("Time")
+fig.tight_layout()
+fig.savefig("custom_compiler_pulse.pdf")
+fig.show()
+
+
+class ClassicalCrossTalk(Noise):
+    def __init__(self, ratio):
+        self.ratio = ratio
+
+    def get_noisy_dynamics(self, dims=None, pulses=None, systematic_noise=None):
+        """Adds noise to the control pulses.
+        
+        Args:
+            dims: Dimension of the system, e.g., [2,2,2,...] for qubits.
+            pulses: A list of Pulse objects, representing the compiled pulses.
+            systematic_noise: A Pulse object with no ideal control, used to represent
+            pulse-independent noise such as decoherence (not used in this example).
+        Returns:
+            pulses: The list of modified pulses according to the noise model.
+            systematic_noise: A Pulse object (not used in this example). 
+        """
+        for i, pulse in enumerate(pulses):
+            if "sx" or "sy" not in pulse.label:
+                pass  # filter out other pulses, e.g. drift
+            target = pulse.targets[0]
+            if target != 0:  # add pulse to the left neighbour
+                pulses[i].add_control_noise(
+                    self.ratio * pulse.qobj,
+                    targets=[target - 1],
+                    coeff=pulse.coeff,
+                    tlist=pulse.tlist,
+                )
+            if target != len(dims) - 1:  # add pulse to the right neighbour
+                pulses[i].add_control_noise(
+                    self.ratio * pulse.qobj,
+                    targets=[target + 1],
+                    coeff=pulse.coeff,
+                    tlist=pulse.tlist,
+                )
+        return pulses, systematic_noise
+
+
+def single_crosstalk_simulation(num_gates):
+    """ A single simulation, with num_gates representing the number of rotations.
+
+    Args:
+        num_gates (int): The random number of gates to add in the simulation.
+
+    Returns:
+        result (qutip.solver.Result): A qutip Result object obtained from any of the
+                                      solver methods such as mesolve.
+    """
+    num_qubits = 2  # Qubit-0 is the target qubit. Qubit-1 suffers from crosstalk.
+    myprocessor = MyProcessor(num_qubits)
+    # Add qubit frequency detuning 1.852MHz for the second qubit.
+    myprocessor.add_drift(2 * np.pi * (sigmaz() + 1) / 2 * 1.852, targets=1)
+    myprocessor.native_gates = None  # Remove the native gates
+    mycompiler = MyCompiler(num_qubits, {"pulse_amplitude": 0.02, "duration": 25})
+    myprocessor.add_noise(ClassicalCrossTalk(1.0))
+    # Define a randome circuit.
+    gates_set = [
+        Gate("ROT", 0, arg_value=0),
+        Gate("ROT", 0, arg_value=np.pi / 2),
+        Gate("ROT", 0, arg_value=np.pi),
+        Gate("ROT", 0, arg_value=np.pi / 2 * 3),
+    ]
+    circuit = QubitCircuit(num_qubits)
+    for ind in np.random.randint(0, 4, num_gates):
+        circuit.add_gate(gates_set[ind])
+    # Simulate the circuit.
+    myprocessor.load_circuit(circuit, compiler=mycompiler)
+    init_state = tensor(
+        [Qobj([[init_fid, 0], [0, 0.025]]), Qobj([[init_fid, 0], [0, 0.025]])]
+    )
+    options = Options(nsteps=10000)  # increase the maximal allowed steps
+    e_ops = [tensor([qeye(2), fock_dm(2)])]  # observable
+
+    # compute results of the run using a solver of choice with custom options
+    result = myprocessor.run_state(init_state, solver="mesolve",
+        options=options, e_ops=e_ops)
+    result = result.expect[0][-1]  # measured expectation value at the end
+    return result
+
+
+num_qubits = 2
+num_sample = 2
+# num_sample = 1600
+fidelity = []
+fidelity_error = []
+init_fid = 0.975
+num_gates_list = [250]
+# num_gates_list = [250, 500, 750, 1000, 1250, 1500]
+
+# The full simulation may take several hours
+# so we just choose num_sample=2 and num_gates=250 as a test
+for num_gates in num_gates_list:
+    expect = Parallel(n_jobs=1)(delayed(single_crosstalk_simulation)(num_gates) for i in range(num_sample))
+    fidelity.append(np.mean(expect))
+    fidelity_error.append(np.std(expect)/np.sqrt(num_sample))
+
+# Recorded result
+num_gates_list = [250, 500, 750, 1000, 1250, 1500]
+data_y = [
+    0.9577285560461476,
+    0.9384849070716464,
+    0.9230217713086177,
+    0.9062344084919285,
+    0.889009550855518,
+    0.8749290612064392,
+]
+data_y_error = [
+    0.000431399017208067,
+    0.0008622091914303468,
+    0.0012216267555118497,
+    0.001537120153687202,
+    0.0018528957172559654,
+    0.0020169257334183596,
+]
+
+
+def linear(x, a):
+    return (1 / 2 + np.exp(-2 * a * x) / 2) * 0.975
+
+
+pos, cov = curve_fit(linear, num_gates_list, data_y, p0=[0.001])
+
+xline = np.linspace(0, 1700, 200)
+yline = linear(xline, *pos)
+
+fig, ax = plt.subplots(figsize=(LINEWIDTH, 0.65 * LINEWIDTH), dpi=200)
+ax.errorbar(
+    num_gates_list, data_y, yerr=data_y_error, fmt=".", capsize=2, color="slategrey"
+)
+ax.plot(xline, yline, color="slategrey")
+ax.set_ylabel("Average fidelity")
+ax.set_xlabel(r"Number of $\pi$ rotations")
+ax.set_xlim((0, 1700))
+
+fig.tight_layout()
+fig.savefig("fig4_cross_talk.pdf")
+fig.show()
diff --git a/doc/pulse-paper/decoherence.py b/doc/pulse-paper/decoherence.py
new file mode 100644
index 00000000..6eb9c434
--- /dev/null
+++ b/doc/pulse-paper/decoherence.py
@@ -0,0 +1,51 @@
+TEXTWIDTH = 7.1398920714
+LINEWIDTH = 3.48692403487
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+try:
+    from quantum_plots import global_setup
+    global_setup(fontsize = 10)
+except:
+    pass
+
+import matplotlib.pyplot as plt
+
+import numpy as np
+from qutip import sigmaz, basis
+from qutip_qip.operations import hadamard_transform
+from qutip_qip.pulse import Pulse
+from qutip_qip.device import LinearSpinChain
+from qutip_qip.circuit import QubitCircuit
+
+# define a circuit with only idling
+t = 30.
+circuit = QubitCircuit(1)
+circuit.add_gate("IDLE", 0, arg_value=t)
+# define a processor with t2 relaxation
+f = 0.5
+t2 = 10 / f
+processor = LinearSpinChain(1, t2=t2)
+processor.add_drift(
+    2*np.pi*sigmaz()/2*f, targets=[0])
+processor.load_circuit(circuit)
+# Record the expectation value
+plus_state = \
+    (basis(2,1) + basis(2,0)).unit()
+result = processor.run_state(
+    init_state=plus_state,
+    tlist = np.linspace(0., t, 1000),
+    # observable
+    e_ops=[plus_state*plus_state.dag()])
+
+tlist = np.linspace(0., t, 1000)
+fig, ax = plt.subplots(figsize = (LINEWIDTH, LINEWIDTH*0.65), dpi=200)
+# detail about lenght of tlist needs to be fixed
+ax.plot(tlist[:-1], result.expect[0][:-1], '-', label="Simulation", color="slategray")
+ax.plot(tlist[:-1], np.exp(-1./t2 * tlist[:-1])*0.5 + 0.5, '--', label="Theory", color="slategray")
+ax.set_xlabel(r"Time [$\mu$s]")
+ax.set_ylabel("Ramsey signal")
+ax.legend()
+ax.grid()
+fig.tight_layout()
+fig.savefig("fig5_decoherence.pdf")
+fig.show()
\ No newline at end of file
diff --git a/doc/pulse-paper/deutsch-jozsa.qasm b/doc/pulse-paper/deutsch-jozsa.qasm
new file mode 100644
index 00000000..564b166f
--- /dev/null
+++ b/doc/pulse-paper/deutsch-jozsa.qasm
@@ -0,0 +1,10 @@
+OPENQASM 2.0;
+include "qelib1.inc";
+
+qreg q[3];
+x q[2];
+h q; 
+cx q[0], q[2];
+cx q[1], q[2];
+h q[0];
+h q[1];
\ No newline at end of file
diff --git a/doc/pulse-paper/deutsch_jozsa-qasm.py b/doc/pulse-paper/deutsch_jozsa-qasm.py
new file mode 100644
index 00000000..cc435498
--- /dev/null
+++ b/doc/pulse-paper/deutsch_jozsa-qasm.py
@@ -0,0 +1,7 @@
+from qutip_qip.qasm import read_qasm
+
+qc = read_qasm("deutsch-jozsa.qasm")
+
+from qutip_qip.qasm import save_qasm
+
+save_qasm(qc, "deutsch-jozsa-qutip.qasm")
diff --git a/doc/pulse-paper/dj_algorithm.py b/doc/pulse-paper/dj_algorithm.py
new file mode 100644
index 00000000..a7d3f7af
--- /dev/null
+++ b/doc/pulse-paper/dj_algorithm.py
@@ -0,0 +1,151 @@
+TEXTWIDTH = 5.93
+LINEWIDTH = 3.22
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+try:
+    from quantum_plots import global_setup
+    global_setup(fontsize = 10)
+except:
+    pass
+plt.rcParams.update({"text.usetex": False})
+
+num_qubits = 3
+def get_operators_labels():
+    """
+    Get the labels for each Hamiltonian.
+    It is used in the method method :meth:`.Processor.plot_pulses`.
+    It is a 2-d nested list, in the plot,
+    a different color will be used for each sublist.
+    """
+    return ([[r"$\Omega^z_%d$" % n for n in range(num_qubits)],
+            [r"$\Omega^x_%d$" % n for n in range(num_qubits)],
+            [r"$g_{%d}$" % (n) for n in range(num_qubits - 1)],
+            ])
+                
+import numpy as np
+import matplotlib.pyplot as plt
+np.random.seed(0)
+from qutip_qip.device import OptPulseProcessor, LinearSpinChain, SCQubits
+from qutip_qip.circuit import QubitCircuit
+from qutip import sigmaz, sigmax, identity, tensor, basis
+
+# Deutsch-Josza algorithm
+dj_circuit = QubitCircuit(3)
+dj_circuit.add_gate("X", targets=2)
+dj_circuit.add_gate("SNOT", targets=0)
+dj_circuit.add_gate("SNOT", targets=1)
+dj_circuit.add_gate("SNOT", targets=2)
+
+# Oracle function f(x)
+dj_circuit.add_gate("CNOT", controls=0, targets=2)
+dj_circuit.add_gate("CNOT", controls=1, targets=2)
+
+dj_circuit.add_gate("SNOT", targets=0)
+dj_circuit.add_gate("SNOT", targets=1)
+
+# Spin chain model
+spinchain_processor = LinearSpinChain(3, t2=30)  # T2 = 30
+spinchain_processor.load_circuit(dj_circuit)
+initial_state = basis([2, 2, 2], [0, 0, 0])  # 3 qubits in the 000 state
+t_record = np.linspace(0, 20, 300)
+result1 = spinchain_processor.run_state(initial_state, tlist=t_record)
+
+# Superconducting qubits
+scqubits_processor = SCQubits(num_qubits)
+scqubits_processor.load_circuit(dj_circuit)
+initial_state = basis([3, 3, 3], [0, 0, 0])  #  3-level
+result2 = scqubits_processor.run_state(initial_state)
+
+# Optimal control model
+setting_args = {"SNOT": {"num_tslots": 6, "evo_time": 2},
+                "X": {"num_tslots": 1, "evo_time": 0.5},
+                "CNOT": {"num_tslots": 12, "evo_time": 5}}
+opt_processor = OptPulseProcessor(N=3)
+opt_processor.add_control(sigmax(), cyclic_permutation=True)
+opt_processor.add_control(sigmaz(), cyclic_permutation=True)
+opt_processor.add_control(tensor([sigmax(), sigmax(), identity(2)]))
+opt_processor.add_control(tensor([identity(2), sigmax(), sigmax()]))
+opt_processor.load_circuit(  # Provide parameters for the algorithm
+    dj_circuit, setting_args=setting_args, merge_gates=False,
+    verbose=True, amp_ubound=5, amp_lbound=0)
+initial_state = basis([2, 2, 2], [0, 0, 0])
+result3 = opt_processor.run_state(initial_state)
+
+# For plotting
+opt_processor.get_operators_labels = get_operators_labels
+width = TEXTWIDTH/3
+fig, ax = opt_processor.plot_pulses(figsize=(width, width*3/2.9), dpi=200);
+
+ax[0].set_ylabel(r"$\Omega^x_{0}$")
+ax[1].set_ylabel(r"$\Omega^x_{1}$")
+ax[2].set_ylabel(r"$\Omega^x_{2}$")
+ax[3].set_ylabel(r"$\Omega^z_{0}$")
+ax[4].set_ylabel(r"$\Omega^z_{1}$")
+ax[5].set_ylabel(r"$\Omega^z_{2}$")
+ax[6].set_ylabel(r"$g_{0}$")
+ax[7].set_ylabel(r"$g_{1}$")
+# ax[7].set_xlabel("Time")
+
+fig.tight_layout()
+fig.savefig("optimal_control_pulse.pdf")
+fig.show()
+
+width = TEXTWIDTH/3
+fig2, ax2 = spinchain_processor.plot_pulses(figsize=(width, width*3/2.9), dpi=200);
+fig2.tight_layout()
+ax2[0].set_ylabel(r"$\Omega^x_{0}$")
+ax2[1].set_ylabel(r"$\Omega^x_{1}$")
+ax2[2].set_ylabel(r"$\Omega^x_{2}$")
+ax2[3].set_ylabel(r"$\Omega^z_{0}$")
+ax2[4].set_ylabel(r"$\Omega^z_{1}$")
+ax2[5].set_ylabel(r"$\Omega^z_{2}$")
+ax2[6].set_ylabel(r"$g_{0}$")
+ax2[7].set_ylabel(r"$g_{1}$")
+# ax2[7].set_xlabel("Time")
+fig2.tight_layout()
+fig2.savefig("spin_chain_pulse.pdf")
+fig2.show()
+
+width = TEXTWIDTH/3
+# fig3, ax3 = scqubits_processor.plot_pulses(figsize=(width, width*3/2.43), dpi=200);
+fig3, ax3 = scqubits_processor.plot_pulses(figsize=(width, width*3/2.4), dpi=200);
+ax3[0].set_ylabel(r"$\Omega^x_{0}$")
+ax3[1].set_ylabel(r"$\Omega^x_{1}$")
+ax3[2].set_ylabel(r"$\Omega^x_{2}$")
+ax3[3].set_ylabel(r"$\Omega^y_{0}$")
+ax3[4].set_ylabel(r"$\Omega^y_{1}$")
+ax3[5].set_ylabel(r"$\Omega^y_{2}$")
+ax3[6].set_ylabel(r"$\Omega^{\rm{cr}1}_{0}$")
+ax3[7].set_ylabel(r"$\Omega^{\rm{cr}2}_{0}$")
+ax3[8].set_ylabel(r"$\Omega^{\rm{cr}1}_{1}$")
+ax3[9].set_ylabel(r"$\Omega^{\rm{cr}2}_{1}$")
+# ax3[9].set_xlabel("Time")
+fig3.tight_layout()
+fig3.savefig("transmon_pulse.pdf")
+fig3.show()
+
+
+plt.rcParams.update({"text.usetex": False})
+
+# Plot hinton
+from qutip import hinton
+fig4, ax4 = plt.subplots(figsize=(LINEWIDTH*0.9, LINEWIDTH*0.7), dpi=200)
+first_two_qubits = result1.states[-1].ptrace([0,1])
+_, ax4 = hinton(first_two_qubits, ax=ax4)
+fig4.savefig("hinton.pdf")
+fig4.show()
+
+# Plot trajectory
+expect = []
+for state in result1.states:
+    tmp = state.ptrace([0,1])
+    tmp = basis([2,2], [0,0]).dag() * tmp * basis([2,2], [0,0])
+    expect.append(np.real(tmp[0, 0]))
+# fig5, ax5 = plt.subplots(figsize=(LINEWIDTH*0.7, LINEWIDTH*0.7*0.6), dpi=200)
+fig5, ax5 = plt.subplots(figsize=(LINEWIDTH, LINEWIDTH*0.7), dpi=200)
+ax5.plot(t_record, expect, color="slategrey")
+ax5.set_ylabel(r"Population of $|00\rangle$")
+ax5.set_xlabel(r"Time [$\mu$s]")
+fig5.tight_layout()
+fig5.savefig("population.pdf")
+fig5.show()
diff --git a/doc/pulse-paper/pulse_simulation.py b/doc/pulse-paper/pulse_simulation.py
new file mode 100644
index 00000000..dbc7d07c
--- /dev/null
+++ b/doc/pulse-paper/pulse_simulation.py
@@ -0,0 +1,93 @@
+TEXTWIDTH = 7.1398920714
+LINEWIDTH = 3.48692403487
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+try:
+    from quantum_plots import global_setup
+    global_setup(fontsize = 10)
+except:
+    pass
+
+# discrete pulse
+import numpy as np
+from qutip_qip.device import Processor
+from qutip import sigmaz, sigmax, basis
+from qutip_qip.pulse import Pulse
+
+processor = Processor(1)
+coeff = np.array([1.])
+tlist = np.array([0., np.pi/2])
+pulse = Pulse(
+    sigmax(), targets=0, tlist=tlist,
+    coeff=coeff, label="sigmax")
+processor.add_pulse(pulse)
+plot_time = np.linspace(0, np.pi/2, 20)
+solver_result = processor.run_state(
+    init_state=basis(2, 0),
+    tlist=plot_time)
+
+mpl.rcParams['axes.prop_cycle'] = mpl.cycler(color=["#1f77b4"]) 
+fig1, ax1 = processor.plot_pulses(figsize=(LINEWIDTH/2, LINEWIDTH/4))
+ax1[0].set_xlim(-0.5, np.pi/2 + 0.5)
+ax1[0].set_ylim(-0.1, 1.1)
+ax1[0].axhline(0)
+ax1[0].set_ylabel(None)
+fig1.tight_layout()
+fig1.savefig("discrete_pulse.pdf")
+fig1.show()
+
+# continuous pulse
+processor = Processor(1)
+processor.pulse_mode = "continuous"
+processor.add_control(sigmax())
+tlist = np.linspace(0., np.pi/2, 21)
+coeff = np.array(np.sin(2*tlist) * np.pi/2)
+processor.pulses[0].tlist = tlist
+processor.pulses[0].coeff = coeff
+solver_result2 = processor.run_state(
+    init_state=basis(2, 0),
+    tlist=plot_time)
+
+mpl.rcParams['axes.prop_cycle'] = mpl.cycler(color=["#ff7f0e"]) 
+fig2, ax2 = processor.plot_pulses(figsize=(LINEWIDTH/2, LINEWIDTH/4))
+ax2[0].set_xlim(-0.5, np.pi/2 + 0.5)
+ax2[0].set_ylim(-0.1, 1.7)
+ax2[0].axhline(0)
+ax2[0].set_ylabel(None)
+fig2.tight_layout()
+fig2.savefig("continuous_pulse.pdf")
+fig2.show()
+
+# bloch sphere
+from qutip import Bloch
+fig3 = plt.figure(figsize = (LINEWIDTH, LINEWIDTH))
+b = Bloch(fig3)
+b.point_color = ["#1f77b4"]
+
+theta = 2 * np.arctan(np.asarray([np.abs(s[1,0]/s[0,0]) for s in solver_result.states]))
+phi = np.angle(np.asarray([s[1,0]/s[0,0] for s in solver_result.states]))
+xp = np.sin(theta) * np.cos(phi)
+yp = np.sin(theta) * np.sin(phi)
+zp = np.cos(theta)
+pnts = [xp, yp, zp]
+b.add_points(pnts)
+b.show()
+fig3.savefig("bloch_sphere1.pdf")
+fig3.show()
+
+from qutip import Bloch
+fig4 = plt.figure(figsize = (LINEWIDTH, LINEWIDTH))
+b = Bloch(fig4)
+b.point_color = ["#ff7f0e"]
+
+theta = 2 * np.arctan(np.asarray([np.abs(s[1,0]/s[0,0]) for s in solver_result2.states]))
+phi = np.angle(np.asarray([s[1,0]/s[0,0] for s in solver_result2.states]))
+xp = np.sin(theta) * np.cos(phi)
+yp = np.sin(theta) * np.sin(phi)
+zp = np.cos(theta)
+pnts = [xp, yp, zp]
+b.add_points(pnts)
+b.show()
+fig4.savefig("bloch_sphere2.pdf")
+fig4.show()
+
diff --git a/doc/pulse-paper/swap.qasm b/doc/pulse-paper/swap.qasm
new file mode 100644
index 00000000..341d7736
--- /dev/null
+++ b/doc/pulse-paper/swap.qasm
@@ -0,0 +1,8 @@
+OPENQASM 2.0;
+include "qelib1.inc";
+
+qreg q[2];
+
+cx q[0], q[1];
+cx q[1], q[0];
+cx q[0], q[1];
diff --git a/docs/requirements.txt b/doc/requirements.txt
similarity index 91%
rename from docs/requirements.txt
rename to doc/requirements.txt
index 68d93787..6603c3c6 100644
--- a/docs/requirements.txt
+++ b/doc/requirements.txt
@@ -7,3 +7,4 @@ sphinx_rtd_theme==0.5.1
 readthedocs-sphinx-search==0.1.0rc3
 numpydoc==1.1.0
 matplotlib==3.3.4
+docutils==0.16
diff --git a/docs/source/_apidoc/qutip_qip.algorithms.rst b/doc/source/_apidoc/qutip_qip.algorithms.rst
similarity index 100%
rename from docs/source/_apidoc/qutip_qip.algorithms.rst
rename to doc/source/_apidoc/qutip_qip.algorithms.rst
diff --git a/docs/source/_apidoc/qutip_qip.compiler.rst b/doc/source/_apidoc/qutip_qip.compiler.rst
similarity index 100%
rename from docs/source/_apidoc/qutip_qip.compiler.rst
rename to doc/source/_apidoc/qutip_qip.compiler.rst
diff --git a/docs/source/_apidoc/qutip_qip.device.rst b/doc/source/_apidoc/qutip_qip.device.rst
similarity index 100%
rename from docs/source/_apidoc/qutip_qip.device.rst
rename to doc/source/_apidoc/qutip_qip.device.rst
diff --git a/docs/source/_apidoc/qutip_qip.operations.rst b/doc/source/_apidoc/qutip_qip.operations.rst
similarity index 100%
rename from docs/source/_apidoc/qutip_qip.operations.rst
rename to doc/source/_apidoc/qutip_qip.operations.rst
diff --git a/docs/source/_apidoc/qutip_qip.rst b/doc/source/_apidoc/qutip_qip.rst
similarity index 100%
rename from docs/source/_apidoc/qutip_qip.rst
rename to doc/source/_apidoc/qutip_qip.rst
diff --git a/docs/source/_apidoc/qutip_qip.transpiler.rst b/doc/source/_apidoc/qutip_qip.transpiler.rst
similarity index 100%
rename from docs/source/_apidoc/qutip_qip.transpiler.rst
rename to doc/source/_apidoc/qutip_qip.transpiler.rst
diff --git a/doc/source/_templates/breadcrumbs.html b/doc/source/_templates/breadcrumbs.html
new file mode 100644
index 00000000..804ad69e
--- /dev/null
+++ b/doc/source/_templates/breadcrumbs.html
@@ -0,0 +1,7 @@
+{%- extends "sphinx_rtd_theme/breadcrumbs.html" %}
+
+{% block breadcrumbs_aside %}
+{% if not meta or meta.get('github_url') != 'hide' %}
+{{ super() }}
+{% endif %}
+{% endblock %}
diff --git a/doc/source/_templates/layout.html b/doc/source/_templates/layout.html
new file mode 100644
index 00000000..e0fcd3cf
--- /dev/null
+++ b/doc/source/_templates/layout.html
@@ -0,0 +1,49 @@
+<!-- https://github.com/tholman/github-corners -->
+<!-- https://stackoverflow.com/a/53329712/10241324 -->
+{% extends '!layout.html' %}
+{% block document %}
+{{ super() }}
+<a href="https://github.com/qutip/qutip-qip" class="github-corner" aria-label="View source on GitHub">
+  <svg width="80" height="80" viewBox="0 0 250 250" style="fill:#fff; color:#6d6d6d; position: absolute; top: 0; border: 0; right: 0;" aria-hidden="true">
+    <path d="M0,0 L115,115 L130,115 L142,142 L250,250 L250,0 Z"></path>
+    <path d="M128.3,109.0 C113.8,99.7 119.0,89.6 119.0,89.6 C122.0,82.7 120.5,78.6 120.5,78.6 C119.2,72.0 123.4,76.3 123.4,76.3 C127.3,80.9 125.5,87.3 125.5,87.3 C122.9,97.6 130.6,101.9 134.4,103.2" fill="currentColor"
+      style="transform-origin: 130px 106px;" class="octo-arm"></path>
+    <path
+      d="M115.0,115.0 C114.9,115.1 118.7,116.5 119.8,115.4 L133.7,101.6 C136.9,99.2 139.9,98.4 142.2,98.6 C133.8,88.0 127.5,74.4 143.8,58.0 C148.5,53.4 154.0,51.2 159.7,51.0 C160.3,49.4 163.2,43.6 171.4,40.1 C171.4,40.1 176.1,42.5 178.8,56.2 C183.1,58.6 187.2,61.8 190.9,65.4 C194.5,69.0 197.7,73.2 200.1,77.6 C213.8,80.2 216.3,84.9 216.3,84.9 C212.7,93.1 206.9,96.0 205.4,96.6 C205.1,102.4 203.0,107.8 198.3,112.5 C181.9,128.9 168.3,122.5 157.7,114.1 C157.9,116.9 156.7,120.9 152.7,124.9 L141.0,136.5 C139.8,137.7 141.6,141.9 141.8,141.8 Z"
+      fill="currentColor" class="octo-body"></path>
+  </svg>
+</a>
+<style>
+  .github-corner:hover .octo-arm {
+    animation: octocat-wave 560ms ease-in-out
+  }
+
+  @keyframes octocat-wave {
+
+    0%,
+    100% {
+      transform: rotate(0)
+    }
+
+    20%,
+    60% {
+      transform: rotate(-25deg)
+    }
+
+    40%,
+    80% {
+      transform: rotate(10deg)
+    }
+  }
+
+  @media (max-width:500px) {
+    .github-corner:hover .octo-arm {
+      animation: none
+    }
+
+    .github-corner .octo-arm {
+      animation: octocat-wave 560ms ease-in-out
+    }
+  }
+</style>
+{% endblock %}
diff --git a/docs/source/conf.py b/doc/source/conf.py
similarity index 86%
rename from docs/source/conf.py
rename to doc/source/conf.py
index e3ec4186..91b180b8 100644
--- a/docs/source/conf.py
+++ b/doc/source/conf.py
@@ -11,6 +11,7 @@
 # documentation root, use os.path.abspath to make it absolute, like shown here.
 #
 import os
+import pathlib
 import sys
 sys.path.insert(0, os.path.abspath('.'))
 
@@ -21,8 +22,24 @@
 copyright = '2021, QuTiP Community'
 author = 'QuTiP Community'
 
-# The full version, including alpha/beta/rc tags
-release = '0.1.0'
+
+def qutip_qip_version():
+    """ Retrieve the qutip-qip version from ``../../VERSION``.
+    """
+    src_folder_root = pathlib.Path(__file__).absolute().parent.parent.parent
+    version = src_folder_root.joinpath(
+        "VERSION"
+    ).read_text().strip()
+    return version
+
+
+# The version info for the project you're documenting, acts as replacement for
+# |version| and |release|, also used in various other places throughout the
+# built documents.
+# The full version, including alpha/beta/rc tags.
+release = qutip_qip_version()
+# The short X.Y version.
+version = ".".join(release.split(".")[:2])
 
 
 # -- General configuration ---------------------------------------------------
@@ -61,7 +78,7 @@
 # Add any paths that contain custom static files (such as style sheets) here,
 # relative to this directory. They are copied after the builtin static files,
 # so a file named "default.css" will overwrite the builtin "default.css".
-html_static_path = ['_static']
+html_static_path = []
 
 # -- Options for LaTeX output ---------------------------------------------
 
@@ -145,4 +162,4 @@
 
 numpydoc_show_class_members = False
 napoleon_numpy_docstring = True
-napoleon_use_admonition_for_notes = True
\ No newline at end of file
+napoleon_use_admonition_for_notes = True
diff --git a/doc/source/contribution-code.rst b/doc/source/contribution-code.rst
new file mode 100644
index 00000000..ed266c98
--- /dev/null
+++ b/doc/source/contribution-code.rst
@@ -0,0 +1,99 @@
+.. _contribute_code:
+
+*******************************
+Contributing to the source code
+*******************************
+
+Build up an development environment
+===================================
+
+Please follow the instruction on the `QuTiP contribution guide <https://qutip.org/docs/latest/development/contributing.html#building>`_ to
+build a conda environment.
+You don't need to build ``qutip`` in the editable mode unless you also want to contribute to `qutip`.
+Instead, you need to install ``qutip-qip`` by downloading the source code and run
+
+.. code-block:: bash
+
+    pip install -e .
+
+Docstrings for the code
+=======================
+
+Each class and function should be accompanied with a docstring
+explaining the functionality, including input parameters and returned values.
+The docstring should follow
+`NumPy Style Python Docstrings <https://www.sphinx-doc.org/en/master/usage/extensions/example_numpy.html>`_.
+
+Checking Code Style and Format
+==============================
+
+In order to check if your code in ``some_file.py`` follows `PEP8 <https://www.python.org/dev/peps/pep-0008/>`_
+style guidelines, `Black <https://black.readthedocs.io/en/stable/the_black_code_style/current_style.html>`_
+has to be installed.
+
+.. code-block::
+
+  pip install black
+
+In the directory that contains ``some_file.py``, use
+
+.. code-block::
+
+  black some_file.py --check
+  black some_file.py --diff --color
+  black some_file.py
+
+Using ``--check`` will show if any of the file will be reformatted or not.
+
+  * `Code 0 <https://black.readthedocs.io/en/stable/usage_and_configuration/the_basics.html#the-basics>`_ means nothing will be reformatted.
+  * Code 1 means one or more files could be reformatted. More than one files could
+    be reformatted if ``black some_directory --check`` is used.
+
+Using ``--diff --color`` will show a `difference <https://black.readthedocs.io/en/stable/usage_and_configuration/the_basics.html#diffs>`_ of
+the changes that will be made by ``Black``. If you would prefer these changes to be made, use the last line of above code block.
+
+.. note::
+  We are currently in the process of checking format of existing code in ``qutip-qip``.
+  Running ``black existing_file.py`` will attempt to format existing code. We
+  advise you to create a separate issue for ``existing_file.py`` or skip re-formatting
+  ``existing_file.py`` in the same PR as your new contribution.
+
+  It is advised to keep your new contribution ``PEP8`` compliant.
+
+Checking tests locally
+=======================
+
+You can run tests and generate code coverage report locally. First make sure
+required packages have been installed.
+
+.. code-block::
+
+  pip install pytest pytest-cov
+
+``pytest`` is used to test files containing tests. If you would like to test all the
+files contained in a directory then specify the path to this directory. In order to run
+tests in ``test_something.py`` then specify the exact path to this file for ``pytest``
+or navigate to the file before running the tests.
+
+.. code-block::
+
+  pytest path_to_some_directory
+  pytest /path_to_test_something/test_something.py
+  ~/path_to_test_something$ pytest test_something.py
+
+A code coverage report in ``html`` format  can be generated locally for
+``qutip-qip`` using the code line given below. By default the coverage report
+is generated in a temporary directory ``htmlcov``. The report can be output
+in `other formats <https://pytest-cov.readthedocs.io/en/latest/reporting.html>`_
+besides ``html``.
+
+.. code-block::
+
+  pytest --cov-report html --cov=qutip_qip tests/
+
+If you would prefer to check the code coverage of one specific file, specify
+the location of this file. Same as above the report can be accessed in ``htmlcov``.
+
+.. code-block::
+
+  pytest --cov-report html --cov=qutip_qip tests/test_something.py
diff --git a/docs/source/contribution-docs.rst b/doc/source/contribution-docs.rst
similarity index 96%
rename from docs/source/contribution-docs.rst
rename to doc/source/contribution-docs.rst
index 2e4ff25e..1e30f380 100644
--- a/docs/source/contribution-docs.rst
+++ b/doc/source/contribution-docs.rst
@@ -32,13 +32,13 @@ To build and test the documentation, the following packages are required:
 
     sphinx numpydoc sphinx_rtd_theme doctest
 
-Under the ``docs`` directory, use
+Under the ``doc`` directory, use
 
 .. code-block:: bash
 
     make html
 
-to build the documentation in html format. The build is saved under the directory ``docs/build/html``
+to build the documentation in html format. The build is saved under the directory ``doc/_build/html``
 
 Use the command
 
@@ -68,8 +68,8 @@ Doctest
 
 
 The doctest directive enables tests on interactive code examples. The simplest way
-to do this is by specifying a prompt along with it's respective output:
-::
+to do this is by specifying a prompt along with it's respective output: ::
+
     .. doctest::
 
         >>> a = 2
@@ -92,16 +92,16 @@ after it. This directive is ideally used when there are a number of examples tha
 need to be checked in quick succession.
 
 A different way to specify code examples (and test them) is using the associated
-**.. testcode::** directive which is effectively a code block:
-::
+**.. testcode::** directive which is effectively a code block: ::
+
     .. testcode::
 
         a = 2
         print(a)
 
 followed by  it's results. The result can be specified with the
-**.. testoutput::** block:
-::
+**.. testoutput::** block: ::
+
     .. testoutput::
 
         2
@@ -130,8 +130,8 @@ A few notes on using the doctest extension:
 
 - By default, each **testcode** and **doctest** block is run in a fresh namespace.
   To share a common namespace, we can specify a common group across the blocks
-  (within a single **.rst** file). For example,
-  ::
+  (within a single **.rst** file). For example, ::
+
         .. doctest:: [group_name]
 
           >>> a = 2
@@ -233,7 +233,7 @@ One can use the following code to generate API documentation:
 
 .. code-block::
 
-    cd docs
+    cd doc
     sphinx-apidoc -f -o source/_apidoc/ ../src/qutip_qip -T
 
 It scan the source code of qutip_qip and save the automatically generated API documentation under the path ``source/_apidoc/``.
diff --git a/docs/source/figures/illustration.png b/doc/source/figures/illustration.png
similarity index 100%
rename from docs/source/figures/illustration.png
rename to doc/source/figures/illustration.png
diff --git a/docs/source/figures/logo.png b/doc/source/figures/logo.png
similarity index 100%
rename from docs/source/figures/logo.png
rename to doc/source/figures/logo.png
diff --git a/docs/source/figures/processor-noise.png b/doc/source/figures/processor-noise.png
similarity index 100%
rename from docs/source/figures/processor-noise.png
rename to doc/source/figures/processor-noise.png
diff --git a/docs/source/figures/processor-workflow.png b/doc/source/figures/processor-workflow.png
similarity index 100%
rename from docs/source/figures/processor-workflow.png
rename to doc/source/figures/processor-workflow.png
diff --git a/docs/source/figures/quantum_circuit_example.png b/doc/source/figures/quantum_circuit_example.png
similarity index 100%
rename from docs/source/figures/quantum_circuit_example.png
rename to doc/source/figures/quantum_circuit_example.png
diff --git a/docs/source/figures/quantum_circuit_w_state.png b/doc/source/figures/quantum_circuit_w_state.png
similarity index 100%
rename from docs/source/figures/quantum_circuit_w_state.png
rename to doc/source/figures/quantum_circuit_w_state.png
diff --git a/docs/source/figures/workflow.png b/doc/source/figures/workflow.png
similarity index 100%
rename from docs/source/figures/workflow.png
rename to doc/source/figures/workflow.png
diff --git a/docs/source/index.rst b/doc/source/index.rst
similarity index 91%
rename from docs/source/index.rst
rename to doc/source/index.rst
index 05c31dd5..02bab36e 100644
--- a/docs/source/index.rst
+++ b/doc/source/index.rst
@@ -23,6 +23,12 @@ qutip-qip: QuTiP quantum information processing
     qip-simulator.rst
     qip-processor.rst
 
+.. toctree::
+    :maxdepth: 2
+    :caption: Tutorials
+
+    tutorials.rst
+
 .. toctree::
     :maxdepth: 2
     :caption: Contribute to qutip-qip
diff --git a/doc/source/installation.rst b/doc/source/installation.rst
new file mode 100644
index 00000000..f068cd4b
--- /dev/null
+++ b/doc/source/installation.rst
@@ -0,0 +1,136 @@
+************
+Installation
+************
+
+.. _quickstart:
+
+Quick start
+===========
+To install the package ``qutip-qip`` from PyPI, use
+
+.. code-block:: bash
+
+    pip install qutip-qip
+
+Migrating from ``qutip.qip``
+============================
+As the :ref:`introduction` suggested, this package is based on a module in the `QuTiP <http://qutip.org/docs/latest/>`_ package ``qutip.qip``.
+If you were using the ``qutip`` package and now want to try out the new features included in this package, you can simply install this package and replace all the ``qutip.qip`` in your import statement with ``qutip_qip``. Everything should work smoothly as usual.
+
+.. _prerequisites:
+
+Prerequisites
+=============
+This package is built upon QuTiP, of which the installation guide can be found at on `QuTiP Installation <http://qutip.org/docs/latest/installation.html>`_.
+The only difference is that C++ compilers are not required here
+since there is no run-time compiling for qutip-qip.
+The minimal Python version supported is Python 3.6.
+
+
+In particular, following packages are necessary for running qutip-qip
+
+.. code-block:: bash
+
+    numpy scipy cython qutip
+
+The following to packages are used for plotting and testing:
+
+.. code-block:: bash
+
+    matplotlib pytest
+
+In addition
+
+.. code-block:: bash
+
+    sphinx numpydoc sphinx_rtd_theme
+
+are used to build and test the documentation.
+
+.. _circuit_plot_packages:
+
+Additional software for Plotting Circuits
+=========================================
+In order to plot circuits, the following non-Python packages are needed:
+
+LaTeX
+-----
+The circuit plotting function in QuTiP uses LaTeX and a few LaTeX packages including ``qcircuit``, ``pdfcrop`` and ``braket``.
+
+**For Linux and Mac** :
+You need to install a TeX distribution such as `TeX Live <https://www.tug.org/texlive/>`_. You can either install it through ``apt-get``/``brew`` or using their installer on the website.
+If you would prefer to avoid installing the full ``texlive`` distribution (which is a few Gigabyte large), `this link <https://tex.stackexchange.com/a/504566/203959>`_
+has some useful discussion on selectively installing part of the packages. In particular:
+
+* ``braket`` will be installed as a part of ``texlive-latex-extra``
+* ``qcircuit`` will be installed as a part of ``texlive-pictures``
+* ``pdfcrop`` and ``pdflatex`` are installed when a minimal ``texlive`` is installed.
+
+.. note::
+  Do not install ``pdfcrop`` via pip because the package's old dependencies will clash
+  with newer dependencies of other programs installed in the ``qutip-qip`` virtual
+  environment.
+
+**For Windows** :
+We recommend installing `MiKTeX <https://miktex.org/>`_, which will automatically install necessary packages like ``qcircuit`` for you when it is used. It will only take a few more minutes in your first attempt at plotting a circuit.
+In addition, you also need to install perl for ``pdfcrop``.
+
+ImageMagick and Ghostscript
+---------------------------
+In order to display the circuit in Jupyter notebook, we need to convert it to png
+format. To do that, you will need to install `Ghostscript <https://www.ghostscript.com/doc/current/Make.htm>`_
+and `ImageMagick <https://imagemagick.org/script/install-source.php>`_.
+The first is responsible for reading the pdf file while the second will convert it to png.
+
+
+.. note::
+    You `might need to make changes <https://stackoverflow.com/a/52863413/10241324>`_ to ``policy.xml`` if the following error occurs :
+
+    .. code-block:: text
+
+        RuntimeError: convert-im6.q16: not authorized `qcirc.pdf' @ error/constitute.c/ReadImage/412.
+        convert-im6.q16: no images defined `qcirc.png' @ error/convert.c/ConvertImageCommand/3258.
+
+
+The output of a circuit plot will be output in a `Jupyter notebook <https://jupyterlab.readthedocs.io/en/stable/getting_started/installation.html>`_.
+In your virtual environment, start `a Jupyter notebook server <https://jupyter.readthedocs.io/en/latest/running.html#starting-the-notebook-server>`_ and run your examples
+in a notebook. If you try to access a circuit plot in a terminal or IPython console,
+you will only be able to access the location of this image in your device's memory.
+
+
+
+pdf2svg
+-------
+To convert the circuit into svg format, you will need to install ``pdf2svg``.
+Please visit `their website <https://github.com/dawbarton/pdf2sv>`_ for installation guide.
+
+.. note::
+    If you want to check whether all dependencies are installed,
+    see if the following three commands work correctly:
+    ``pdflatex``, ``pdfcrop`` and ``magick anypdf.pdf anypdf.png``,
+    where ``anypdf.pdf`` is any pdf file you have.
+
+.. _installation:
+
+Install qutip-qip from source code
+==================================
+
+To install the package, download to source code from `GitHub website <https://github.com/qutip/qutip-qip>`_ and run
+
+.. code-block:: bash
+
+    pip install .
+
+under the directory containing the ``setup.cfg`` file.
+
+If you want to edit the code, use instead
+
+.. code-block:: bash
+
+    pip install -e .
+
+To test the installation from a download of the source code, run from the `qutip-qip` directory
+
+```
+pytest tests
+```
diff --git a/docs/source/introduction.rst b/doc/source/introduction.rst
similarity index 61%
rename from docs/source/introduction.rst
rename to doc/source/introduction.rst
index 972cbddc..ff6e4ecf 100644
--- a/docs/source/introduction.rst
+++ b/doc/source/introduction.rst
@@ -1,4 +1,4 @@
-.. introduction:
+.. _introduction:
 
 ************
 Introduction
@@ -15,3 +15,27 @@ They are still maintained by the QuTiP team but hosted under different repositor
 The qutip-qip package, QuTiP quantum information processing, aims at providing basic tools for quantum computing simulation both for simple quantum algorithm design and for experimental realization.
 Compared to other libraries for quantum information processing, qutip-qip puts additional emphasis on the physics layer and the interaction with the QuTiP package.
 It offers two different approaches for simulating quantum circuits, one with :class:`.QubitCircuit` calculating unitary evolution under quantum gates by matrix product, another called :class:`.Processor` using open system solvers in QuTiP to simulate noisy quantum device.
+
+Citing
+===========
+
+If you use ``qutip.qip`` in your research, cite the `preprint <https://arxiv.org/abs/2105.09902>`_
+as
+
+.. code-block:: text
+
+  Boxi Li, Shahnawaz Ahmed, Sidhant Saraogi, Neill Lambert, Franco Nori, Alexander Pitchford, & Nathan Shammah. (2021). Pulse-level noisy quantum circuits with QuTiP.
+
+
+The bibtex can be downloaded directly :download:`here<qutip_qip.bib>` or
+copy-pasted using :
+
+.. code-block:: text
+
+  @misc{li2021pulselevel,
+        title={Pulse-level noisy quantum circuits with QuTiP},
+        author={Boxi Li and Shahnawaz Ahmed and Sidhant Saraogi and Neill Lambert and Franco Nori and Alexander Pitchford and Nathan Shammah},
+        year={2021},
+        eprint={2105.09902},
+        archivePrefix={arXiv},
+        primaryClass={quant-ph}}
diff --git a/docs/source/qip-basics.rst b/doc/source/qip-basics.rst
similarity index 92%
rename from docs/source/qip-basics.rst
rename to doc/source/qip-basics.rst
index 947d728f..09dbd374 100644
--- a/docs/source/qip-basics.rst
+++ b/doc/source/qip-basics.rst
@@ -284,20 +284,8 @@ Plotting a Quantum Circuit
 A quantum circuit (described above) can directly be plotted using the QCircuit library (https://github.com/CQuIC/qcircuit).
 QCiruit is a quantum circuit drawing application and is implemented directly into QuTiP.
 
-The circuit image visualization requires LaTeX and ImageMagick for display.
-The module automatically generates the LaTeX code for plotting the circuit,
-produces the pdf and converts it to the png format. On Mac and Linux,
-ImageMagick can be easily installed with the command conda install imagemagick if you have conda installed.
-Otherwise, please follow the installation instructions on the ImageMagick documentation.
-
-On windows, you need to download and install ImageMagick installer.
-In addition, you also need perl (for ``pdfcrop``) and
-Ghostscript (additional dependency of ImageMagick for png conversion).
-
-If you want to check whether all dependencies are installed,
-see if the following three commands work correctly:
-``pdflatex``, ``pdfcrop`` and ``magick anypdf.pdf anypdf.png``,
-where ``anypdf.pdf`` is any pdf file you have.
+More information related to installing these packages is also available in the
+installation guide (:ref:`circuit_plot_packages`).
 
 An example code for plotting the example quantum circuit from above is given:
 
diff --git a/docs/source/qip-processor.rst b/doc/source/qip-processor.rst
similarity index 100%
rename from docs/source/qip-processor.rst
rename to doc/source/qip-processor.rst
diff --git a/docs/source/qip-simulator.rst b/doc/source/qip-simulator.rst
similarity index 100%
rename from docs/source/qip-simulator.rst
rename to doc/source/qip-simulator.rst
diff --git a/doc/source/qutip_qip.bib b/doc/source/qutip_qip.bib
new file mode 100644
index 00000000..f2439a86
--- /dev/null
+++ b/doc/source/qutip_qip.bib
@@ -0,0 +1,7 @@
+  @misc{li2021pulselevel,
+      title={Pulse-level noisy quantum circuits with QuTiP},
+      author={Boxi Li and Shahnawaz Ahmed and Sidhant Saraogi and Neill Lambert and Franco Nori and Alexander Pitchford and Nathan Shammah},
+      year={2021},
+      eprint={2105.09902},
+      archivePrefix={arXiv},
+      primaryClass={quant-ph}}
diff --git a/doc/source/tutorials.rst b/doc/source/tutorials.rst
new file mode 100644
index 00000000..358df2e0
--- /dev/null
+++ b/doc/source/tutorials.rst
@@ -0,0 +1,10 @@
+.. _tutorials:
+
+************
+Tutorials
+************
+
+Tutorials related to using quantum gates and circuits in ``qutip-qip`` can be
+found `here <https://qutip.org/tutorials#quantum-information-processing>`_ and
+those related to using noise simulators are available at this
+`link <https://qutip.org/tutorials#nisq>`_.
diff --git a/docs/source/w-state.qasm b/doc/source/w-state.qasm
similarity index 100%
rename from docs/source/w-state.qasm
rename to doc/source/w-state.qasm
diff --git a/docs/source/contribution-code.rst b/docs/source/contribution-code.rst
deleted file mode 100644
index 10025f82..00000000
--- a/docs/source/contribution-code.rst
+++ /dev/null
@@ -1,25 +0,0 @@
-.. _contribute_code:
-
-*******************************
-Contributing to the source code
-*******************************
-
-Build up an development environment
-===================================
-
-Please follow the instruction on the `QuTiP contribution guide <https://github.com/qutip/qutip-doc/blob/master/CONTRIBUTING.md>`_ to 
-build a conda environment.
-Then install qutip-qip by downloading the source code and run
-
-.. code-block:: bash
-
-    pip install -e .
-
-Docstrings for the code
-=======================
-
-Each class and function should be accompanied with a docstring
-explaining the functionality, including input parameters and returned values.
-The docstring should follow
-`NumPy Style Python Docstrings <https://www.sphinx-doc.org/en/master/usage/extensions/example_numpy.html>`_.
-
diff --git a/docs/source/installation.rst b/docs/source/installation.rst
deleted file mode 100644
index 8c944100..00000000
--- a/docs/source/installation.rst
+++ /dev/null
@@ -1,62 +0,0 @@
-************
-Installation
-************
-
-Quick start
-===========
-To install the package ``qutip-qip`` from PyPI, use
-
-.. code-block:: bash
-
-    pip install qutip_qip
-
-.. _prerequisites:
-
-Prerequisites
-=============
-As the name indicates, this package is built upon QuTiP, the installation guide can be found at on `QuTiP Installation <http://qutip.org/docs/latest/installation.html>`_.
-The only difference is that C++ compilers are not required here
-since there is no run-time compiling for qutip-qip.
-The minimal Python version supported is Python 3.6.
-
-
-In particular, following packages are necessary for running qutip-qip
-
-.. code-block:: bash
-
-    numpy scipy cython qutip
-
-The following to packages are used for plotting and testing:
-
-.. code-block:: bash
-
-    matplotlib pytest
-
-In addition
-
-.. code-block:: bash
-
-    sphinx numpydoc sphinx_rtd_theme
-
-are used to build and test the documentation.
-
-A few other packages such as LaTeX is used for circuit plotting, please refer to the main documentation section for detailed instruction.
-
-.. _installation:
-
-Install qutip-qip from source code
-==================================
-
-To install the package, download to source code from `GitHub website <https://github.com/qutip/qutip-qip>`_ and run
-
-.. code-block:: bash
-
-    pip install .
-
-under the directory containing the ``setup.cfg`` file.
-
-If you want to edit the code, use instead
-
-.. code-block:: bash
-
-    pip install -e .
diff --git a/pyproject.toml b/pyproject.toml
index f6c16894..0f7f28f1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,7 @@
 [build-system]
 requires = [
     "setuptools>=42",
+    "packaging",
     "wheel",
 ]
 build-backend = "setuptools.build_meta"
diff --git a/setup.cfg b/setup.cfg
index 06f6a09f..40be6f8d 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -21,7 +21,6 @@ classifiers =
 author = Alexander Pitchford, Paul D. Nation, Robert J. Johansson, Chris Granade, Arne Grimsmo, Nathan Shammah, Shahnawaz Ahmed, Neill Lambert, Eric Giguere, Boxi Li, Jake Lishman
 author_email = qutip-admin@googlegroups.com
 platforms = Linux, Mac OSX, Unix, Windows
-version = 0.1.0
 
 [options]
 package_dir=
@@ -32,6 +31,8 @@ install_requires =
     numpy>=1.16.6
     scipy>=1.0
     qutip>=4.6,<5
+setup_requires =
+    packaging
 
 [options.packages.find]
 where = src
diff --git a/setup.py b/setup.py
new file mode 100755
index 00000000..6e45f3eb
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,96 @@
+#!/usr/bin/env python
+
+import os
+import subprocess
+import sys
+
+# Required third-party imports, must be specified in pyproject.toml.
+import packaging.version
+from setuptools import setup
+
+
+def process_options():
+    """
+    Determine all runtime options, returning a dictionary of the results.  The
+    keys are:
+        'rootdir': str
+            The root directory of the setup.  Almost certainly the directory
+            that this setup.py file is contained in.
+        'release': bool
+            Is this a release build (True) or a local development build (False)
+    """
+    options = {}
+    options['rootdir'] = os.path.dirname(os.path.abspath(__file__))
+    options = _determine_version(options)
+    return options
+
+
+def _determine_version(options):
+    """
+    Adds the 'short_version', 'version' and 'release' options.
+
+    Read from the VERSION file to discover the version.  This should be a
+    single line file containing valid Python package public identifier (see PEP
+    440), for example
+      4.5.2rc2
+      5.0.0
+      5.1.1a1
+    We do that here rather than in setup.cfg so we can apply the local
+    versioning number as well.
+    """
+    version_filename = os.path.join(options['rootdir'], 'VERSION')
+    with open(version_filename, "r") as version_file:
+        version_string = version_file.read().strip()
+    version = packaging.version.parse(version_string)
+    if isinstance(version, packaging.version.LegacyVersion):
+        raise ValueError("invalid version: " + version_string)
+    options['short_version'] = str(version.public)
+    options['release'] = not version.is_devrelease
+    if not options['release']:
+        # Put the version string into canonical form, if it wasn't already.
+        version_string = str(version)
+        version_string += "+"
+        try:
+            git_out = subprocess.run(
+                ('git', 'rev-parse', '--verify', '--short=7', 'HEAD'),
+                check=True,
+                stdout=subprocess.PIPE, stderr=subprocess.PIPE,
+            )
+            git_hash = git_out.stdout.decode(sys.stdout.encoding).strip()
+            version_string += git_hash or "nogit"
+        # CalledProcessError is for if the git command fails for internal
+        # reasons (e.g. we're not in a git repository), OSError is for if
+        # something goes wrong when trying to run git (e.g. it's not installed,
+        # or a permission error).
+        except (subprocess.CalledProcessError, OSError):
+            version_string += "nogit"
+    options['version'] = version_string
+    return options
+
+
+def create_version_py_file(options):
+    """
+    Generate and write out the file qutip/version.py, which is used to produce
+    the '__version__' information for the module.  This function will overwrite
+    an existing file at that location.
+    """
+    filename = os.path.join(
+        options['rootdir'], 'src', 'qutip_qip', 'version.py')
+    content = "\n".join([
+        "# This file is automatically generated by qutip-qip's setup.py.",
+        f"short_version = '{options['short_version']}'",
+        f"version = '{options['version']}'",
+        f"release = {options['release']}",
+    ])
+    with open(filename, 'w') as file:
+        print(content, file=file)
+
+
+if __name__ == "__main__":
+    options = process_options()
+    create_version_py_file(options)
+    # Most of the kwargs to setup are defined in setup.cfg; the only ones we
+    # keep here are ones that we have done some compile-time processing on.
+    setup(
+        version=options['version'],
+    )
diff --git a/src/qutip_qip/__init__.py b/src/qutip_qip/__init__.py
index 09e4688a..69a06464 100644
--- a/src/qutip_qip/__init__.py
+++ b/src/qutip_qip/__init__.py
@@ -30,4 +30,5 @@
 #    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 #    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 ###############################################################################
-__version__ = "0.1.0"
+
+from .version import version as __version__  # noqa
diff --git a/src/qutip_qip/circuit.py b/src/qutip_qip/circuit.py
index c8182a69..e75cb591 100644
--- a/src/qutip_qip/circuit.py
+++ b/src/qutip_qip/circuit.py
@@ -34,22 +34,20 @@
 from itertools import product
 import numbers
 
-import warnings
 import inspect
 
 import numpy as np
 from copy import deepcopy
 
 from . import circuit_latex as _latex
-from .operations.gates import (rx, ry, rz, sqrtnot, snot, phasegate,
-                                        x_gate, y_gate, z_gate, cy_gate,
-                                        cz_gate, s_gate, t_gate, cs_gate,
-                                        qasmu_gate, ct_gate, cphase, cnot,
-                                        csign, berkeley, swapalpha, swap,
-                                        iswap, sqrtswap, sqrtiswap, fredkin,
-                                        toffoli, controlled_gate, globalphase,
-                                        expand_operator, gate_sequence_product)
-from qutip import tensor, basis, identity, ket2dm, qeye
+from .operations import (
+    rx, ry, rz, sqrtnot, snot, phasegate,
+    x_gate, y_gate, z_gate, cy_gate, cz_gate, s_gate, t_gate, cs_gate,
+    qasmu_gate, ct_gate, cphase, cnot, csign, berkeley, swapalpha, swap,
+    iswap, sqrtswap, sqrtiswap, fredkin, toffoli, controlled_gate, globalphase,
+    expand_operator, gate_sequence_product
+)
+from qutip import basis, ket2dm, qeye
 from qutip.qobj import Qobj
 from qutip.measurement import measurement_statistics
 
@@ -413,7 +411,6 @@ def __init__(self, N, input_states=None, output_states=None,
         self.N = N
         self.reverse_states = reverse_states
         self.gates = []
-        self.U_list = []
         self.dims = dims
         self.num_cbits = num_cbits
 
@@ -555,8 +552,6 @@ def add_gate(self, gate, targets=None, controls=None, arg_value=None,
 
         else:
             for position in index:
-                num_mes = (sum(isinstance(op, Measurement) for op
-                               in self.gates[:position]))
                 gate = Gate(name, targets=targets, controls=controls,
                             arg_value=arg_value, arg_label=arg_label,
                             classical_controls=classical_controls,
@@ -780,8 +775,10 @@ def _gate_PHASEGATE(self, gate, temp_resolved):
     def _gate_NOTIMPLEMENTED(self, gate, temp_resolved):
         raise NotImplementedError("Cannot be resolved in this basis")
 
-    _gate_PHASEGATE = _gate_BERKELEY = _gate_SWAPalpha = _gate_NOTIMPLEMENTED
-    _gate_SQRTSWAP = _gate_SQRTISWAP = _gate_NOTIMPLEMENTED
+    _gate_BERKELEY = _gate_NOTIMPLEMENTED
+    _gate_SWAPalpha = _gate_NOTIMPLEMENTED
+    _gate_SQRTSWAP = _gate_NOTIMPLEMENTED
+    _gate_SQRTISWAP = _gate_NOTIMPLEMENTED
 
     def _gate_CSIGN(self, gate, temp_resolved):
         half_pi = np.pi / 2
@@ -1449,8 +1446,20 @@ def adjacent_gates(self):
 
     def propagators(self, expand=True):
         """
-        Propagator matrix calculator for N qubits returning the individual
+        Propagator matrix calculator returning the individual
         steps as unitary matrices operating from left to right.
+
+        Parameters
+        ----------
+        expand : bool, optional
+            Whether to expand the unitary matrices for the individual
+            steps to the full Hilbert space for N qubits.
+            Defaults to ``True``.
+            If ``False``, the unitary matrices will not be expanded and the
+            list of unitaries will need to be combined with the list of
+            gates in order to determine which qubits the unitaries should
+            act on.
+
         Returns
         -------
         U_list : list
@@ -1458,101 +1467,101 @@ def propagators(self, expand=True):
         """
 
         if not expand:
-            return self.propagators_no_expand()
+            return self._propagators_no_expand()
 
-        self.U_list = []
+        U_list = []
 
         gates = filter(lambda x: isinstance(x, Gate), self.gates)
 
         for gate in gates:
             if gate.name == "RX":
-                self.U_list.append(rx(
+                U_list.append(rx(
                     gate.arg_value, self.N, gate.targets[0]))
             elif gate.name == "RY":
-                self.U_list.append(ry(
+                U_list.append(ry(
                     gate.arg_value, self.N, gate.targets[0]))
             elif gate.name == "RZ":
-                self.U_list.append(rz(
+                U_list.append(rz(
                     gate.arg_value, self.N, gate.targets[0]))
             elif gate.name == "X":
-                self.U_list.append(x_gate(self.N, gate.targets[0]))
+                U_list.append(x_gate(self.N, gate.targets[0]))
             elif gate.name == "Y":
-                self.U_list.append(y_gate(self.N, gate.targets[0]))
+                U_list.append(y_gate(self.N, gate.targets[0]))
             elif gate.name == "CY":
-                self.U_list.append(cy_gate(
+                U_list.append(cy_gate(
                     self.N, gate.controls[0], gate.targets[0]))
             elif gate.name == "Z":
-                self.U_list.append(z_gate(self.N, gate.targets[0]))
+                U_list.append(z_gate(self.N, gate.targets[0]))
             elif gate.name == "CZ":
-                self.U_list.append(cz_gate(
+                U_list.append(cz_gate(
                     self.N, gate.controls[0], gate.targets[0]))
             elif gate.name == "T":
-                self.U_list.append(t_gate(self.N, gate.targets[0]))
+                U_list.append(t_gate(self.N, gate.targets[0]))
             elif gate.name == "CT":
-                self.U_list.append(ct_gate(
+                U_list.append(ct_gate(
                     self.N, gate.controls[0], gate.targets[0]))
             elif gate.name == "S":
-                self.U_list.append(s_gate(self.N, gate.targets[0]))
+                U_list.append(s_gate(self.N, gate.targets[0]))
             elif gate.name == "CS":
-                self.U_list.append(cs_gate(
+                U_list.append(cs_gate(
                     self.N, gate.controls[0], gate.targets[0]))
             elif gate.name == "SQRTNOT":
-                self.U_list.append(sqrtnot(self.N, gate.targets[0]))
+                U_list.append(sqrtnot(self.N, gate.targets[0]))
             elif gate.name == "SNOT":
-                self.U_list.append(snot(self.N, gate.targets[0]))
+                U_list.append(snot(self.N, gate.targets[0]))
             elif gate.name == "PHASEGATE":
-                self.U_list.append(phasegate(gate.arg_value, self.N,
-                                             gate.targets[0]))
+                U_list.append(phasegate(gate.arg_value, self.N,
+                                        gate.targets[0]))
             elif gate.name == "QASMU":
-                self.U_list.append(qasmu_gate(gate.arg_value, self.N,
-                                              gate.targets[0]))
+                U_list.append(qasmu_gate(gate.arg_value, self.N,
+                                         gate.targets[0]))
             elif gate.name == "CRX":
-                self.U_list.append(controlled_gate(rx(gate.arg_value),
-                                                   N=self.N,
-                                                   control=gate.controls[0],
-                                                   target=gate.targets[0]))
+                U_list.append(controlled_gate(rx(gate.arg_value),
+                                              N=self.N,
+                                              control=gate.controls[0],
+                                              target=gate.targets[0]))
             elif gate.name == "CRY":
-                self.U_list.append(controlled_gate(ry(gate.arg_value),
-                                                   N=self.N,
-                                                   control=gate.controls[0],
-                                                   target=gate.targets[0]))
+                U_list.append(controlled_gate(ry(gate.arg_value),
+                                              N=self.N,
+                                              control=gate.controls[0],
+                                              target=gate.targets[0]))
             elif gate.name == "CRZ":
-                self.U_list.append(controlled_gate(rz(gate.arg_value),
-                                                   N=self.N,
-                                                   control=gate.controls[0],
-                                                   target=gate.targets[0]))
+                U_list.append(controlled_gate(rz(gate.arg_value),
+                                              N=self.N,
+                                              control=gate.controls[0],
+                                              target=gate.targets[0]))
             elif gate.name == "CPHASE":
-                self.U_list.append(cphase(gate.arg_value, self.N,
-                                          gate.controls[0], gate.targets[0]))
+                U_list.append(cphase(gate.arg_value, self.N,
+                                     gate.controls[0], gate.targets[0]))
             elif gate.name == "CNOT":
-                self.U_list.append(cnot(self.N,
-                                        gate.controls[0], gate.targets[0]))
+                U_list.append(cnot(self.N,
+                                   gate.controls[0], gate.targets[0]))
             elif gate.name == "CSIGN":
-                self.U_list.append(csign(self.N,
-                                         gate.controls[0], gate.targets[0]))
+                U_list.append(csign(self.N,
+                                    gate.controls[0], gate.targets[0]))
             elif gate.name == "BERKELEY":
-                self.U_list.append(berkeley(self.N, gate.targets))
+                U_list.append(berkeley(self.N, gate.targets))
             elif gate.name == "SWAPalpha":
-                self.U_list.append(swapalpha(gate.arg_value, self.N,
-                                             gate.targets))
+                U_list.append(swapalpha(gate.arg_value, self.N,
+                                        gate.targets))
             elif gate.name == "SWAP":
-                self.U_list.append(swap(self.N, gate.targets))
+                U_list.append(swap(self.N, gate.targets))
             elif gate.name == "ISWAP":
-                self.U_list.append(iswap(self.N, gate.targets))
+                U_list.append(iswap(self.N, gate.targets))
             elif gate.name == "SQRTSWAP":
-                self.U_list.append(sqrtswap(self.N, gate.targets))
+                U_list.append(sqrtswap(self.N, gate.targets))
             elif gate.name == "SQRTISWAP":
-                self.U_list.append(sqrtiswap(self.N, gate.targets))
+                U_list.append(sqrtiswap(self.N, gate.targets))
             elif gate.name == "FREDKIN":
-                self.U_list.append(fredkin(self.N, gate.controls[0],
-                                           gate.targets))
+                U_list.append(fredkin(self.N, gate.controls[0],
+                                      gate.targets))
             elif gate.name == "TOFFOLI":
-                self.U_list.append(toffoli(self.N, gate.controls,
-                                           gate.targets[0]))
+                U_list.append(toffoli(self.N, gate.controls,
+                                      gate.targets[0]))
             elif gate.name == "GLOBALPHASE":
-                self.U_list.append(globalphase(gate.arg_value, self.N))
+                U_list.append(globalphase(gate.arg_value, self.N))
             elif gate.name == "IDLE":
-                self.U_list.append(qeye(self.N * [2]))
+                U_list.append(qeye(self.N * [2]))
             elif gate.name in self.user_gates:
                 if gate.controls is not None:
                     raise ValueError("A user defined gate {} takes only  "
@@ -1572,15 +1581,15 @@ def propagators(self, expand=True):
                     oper = func_or_oper
                 else:
                     raise ValueError("gate is neither function nor operator")
-                self.U_list.append(expand_operator(
+                U_list.append(expand_operator(
                     oper, N=self.N, targets=gate.targets, dims=self.dims))
             else:
                 raise NotImplementedError(
                     "{} gate is an unknown gate.".format(gate.name))
 
-        return self.U_list
+        return U_list
 
-    def propagators_no_expand(self):
+    def _propagators_no_expand(self):
         """
         Propagator matrix calculator for N qubits returning the individual
         steps as unitary matrices operating from left to right.
@@ -1591,75 +1600,75 @@ def propagators_no_expand(self):
             Return list of unitary matrices for the qubit circuit.
 
         """
-        self.U_list = []
+        U_list = []
 
         gates = filter(lambda x: isinstance(x, Gate), self.gates)
 
         for gate in gates:
             if gate.name == "RX":
-                self.U_list.append(rx(gate.arg_value))
+                U_list.append(rx(gate.arg_value))
             elif gate.name == "RY":
-                self.U_list.append(ry(gate.arg_value))
+                U_list.append(ry(gate.arg_value))
             elif gate.name == "RZ":
-                self.U_list.append(rz(gate.arg_value))
+                U_list.append(rz(gate.arg_value))
             elif gate.name == "X":
-                self.U_list.append(x_gate())
+                U_list.append(x_gate())
             elif gate.name == "Y":
-                self.U_list.append(y_gate())
+                U_list.append(y_gate())
             elif gate.name == "CY":
-                self.U_list.append(cy_gate())
+                U_list.append(cy_gate())
             elif gate.name == "Z":
-                self.U_list.append(z_gate())
+                U_list.append(z_gate())
             elif gate.name == "CZ":
-                self.U_list.append(cz_gate())
+                U_list.append(cz_gate())
             elif gate.name == "T":
-                self.U_list.append(t_gate())
+                U_list.append(t_gate())
             elif gate.name == "CT":
-                self.U_list.append(ct_gate())
+                U_list.append(ct_gate())
             elif gate.name == "S":
-                self.U_list.append(s_gate())
+                U_list.append(s_gate())
             elif gate.name == "CS":
-                self.U_list.append(cs_gate())
+                U_list.append(cs_gate())
             elif gate.name == "SQRTNOT":
-                self.U_list.append(sqrtnot())
+                U_list.append(sqrtnot())
             elif gate.name == "SNOT":
-                self.U_list.append(snot())
+                U_list.append(snot())
             elif gate.name == "PHASEGATE":
-                self.U_list.append(phasegate(gate.arg_value))
+                U_list.append(phasegate(gate.arg_value))
             elif gate.name == "QASMU":
-                self.U_list.append(qasmu_gate(gate.arg_value))
+                U_list.append(qasmu_gate(gate.arg_value))
             elif gate.name == "CRX":
-                self.U_list.append(controlled_gate(rx(gate.arg_value)))
+                U_list.append(controlled_gate(rx(gate.arg_value)))
             elif gate.name == "CRY":
-                self.U_list.append(controlled_gate(ry(gate.arg_value)))
+                U_list.append(controlled_gate(ry(gate.arg_value)))
             elif gate.name == "CRZ":
-                self.U_list.append(controlled_gate(rz(gate.arg_value)))
+                U_list.append(controlled_gate(rz(gate.arg_value)))
             elif gate.name == "CPHASE":
-                self.U_list.append(cphase(gate.arg_value))
+                U_list.append(cphase(gate.arg_value))
             elif gate.name == "CNOT":
-                self.U_list.append(cnot())
+                U_list.append(cnot())
             elif gate.name == "CSIGN":
-                self.U_list.append(csign())
+                U_list.append(csign())
             elif gate.name == "BERKELEY":
-                self.U_list.append(berkeley())
+                U_list.append(berkeley())
             elif gate.name == "SWAPalpha":
-                self.U_list.append(swapalpha(gate.arg_value))
+                U_list.append(swapalpha(gate.arg_value))
             elif gate.name == "SWAP":
-                self.U_list.append(swap())
+                U_list.append(swap())
             elif gate.name == "ISWAP":
-                self.U_list.append(iswap())
+                U_list.append(iswap())
             elif gate.name == "SQRTSWAP":
-                self.U_list.append(sqrtswap())
+                U_list.append(sqrtswap())
             elif gate.name == "SQRTISWAP":
-                self.U_list.append(sqrtiswap())
+                U_list.append(sqrtiswap())
             elif gate.name == "FREDKIN":
-                self.U_list.append(fredkin())
+                U_list.append(fredkin())
             elif gate.name == "TOFFOLI":
-                self.U_list.append(toffoli())
+                U_list.append(toffoli())
             elif gate.name == "GLOBALPHASE":
-                self.U_list.append(globalphase(gate.arg_value, n))
+                U_list.append(globalphase(gate.arg_value, self.N))
             elif gate.name == "IDLE":
-                self.U_list.append(qeye(2))
+                U_list.append(qeye(2))
             elif gate.name in self.user_gates:
                 if gate.controls is not None:
                     raise ValueError("A user defined gate {} takes only  "
@@ -1679,12 +1688,12 @@ def propagators_no_expand(self):
                     oper = func_or_oper
                 else:
                     raise ValueError("gate is neither function nor operator")
-                self.U_list.append(oper)
+                U_list.append(oper)
             else:
                 raise NotImplementedError(
                     "{} gate is an unknown gate.".format(gate.name))
 
-        return self.U_list
+        return U_list
 
     def latex_code(self):
         rows = []
diff --git a/src/qutip_qip/circuit_latex.py b/src/qutip_qip/circuit_latex.py
index cc9e32b5..6f8405e8 100644
--- a/src/qutip_qip/circuit_latex.py
+++ b/src/qutip_qip/circuit_latex.py
@@ -41,7 +41,8 @@
 
 _latex_template = r"""
 \documentclass{standalone}
-%s
+\usepackage[braket]{qcircuit}
+\renewcommand{\qswap}{*=<0em>{\times}}
 \begin{document}
 \Qcircuit @C=1cm @R=1cm {
 %s}
@@ -227,9 +228,17 @@ def image_from_latex(code, file_type="png"):
             try:
                 os.chdir(temporary_dir)
                 with open(filename + ".tex", "w") as file:
-                    file.write(_latex_template % (_qcircuit_latex_min, code))
-                _run_command((_pdflatex, '-interaction', 'batchmode',
-                              filename))
+                    file.write(_latex_template % code)
+                try:
+                    _run_command((_pdflatex, '-interaction', 'batchmode',
+                                  filename))
+                except RuntimeError as e:
+                    message = (
+                        "pdflatex failed."
+                        " Perhaps you do not have it installed, or you are"
+                        " missing the LaTeX package 'qcircuit'."
+                    )
+                    raise RuntimeError(message) from e
                 _crop_pdf(filename + ".pdf")
                 if file_type in _MISSING_CONVERTERS:
                     dependency = _MISSING_CONVERTERS[file_type]
@@ -251,50 +260,3 @@ def image_from_latex(code, file_type="png"):
 else:
     def image_from_latex(*args, **kwargs):
         raise RuntimeError("Could not find system 'pdflatex'.")
-
-
-_qcircuit_latex_min = r"""
-% Q-circuit version 2
-% Copyright (C) 2004 Steve Flammia & Bryan Eastin
-% Last modified on: 9/16/2011
-% License: http://www.gnu.org/licenses/gpl-2.0.html
-% Original file: http://physics.unm.edu/CQuIC/Qcircuit/Qcircuit.tex
-% Modified for QuTiP on: 2/19/2021
-\usepackage{xy}
-\xyoption{matrix}
-\xyoption{frame}
-\xyoption{arrow}
-\xyoption{arc}
-\usepackage{ifpdf}
-\entrymodifiers={!C\entrybox}
-\newcommand{\bra}[1]{{\left\langle{#1}\right\vert}}
-\newcommand{\ket}[1]{{\left\vert{#1}\right\rangle}}
-\newcommand{\qw}[1][-1]{\ar @{-} [0,#1]}
-\newcommand{\qwx}[1][-1]{\ar @{-} [#1,0]}
-\newcommand{\cw}[1][-1]{\ar @{=} [0,#1]}
-\newcommand{\cwx}[1][-1]{\ar @{=} [#1,0]}
-\newcommand{\gate}[1]{*+<.6em>{#1} \POS ="i","i"+UR;"i"+UL **\dir{-};"i"+DL **\dir{-};"i"+DR **\dir{-};"i"+UR **\dir{-},"i" \qw}
-\newcommand{\meter}{*=<1.8em,1.4em>{\xy ="j","j"-<.778em,.322em>;{"j"+<.778em,-.322em> \ellipse ur,_{}},"j"-<0em,.4em>;p+<.5em,.9em> **\dir{-},"j"+<2.2em,2.2em>*{},"j"-<2.2em,2.2em>*{} \endxy} \POS ="i","i"+UR;"i"+UL **\dir{-};"i"+DL **\dir{-};"i"+DR **\dir{-};"i"+UR **\dir{-},"i" \qw}
-\newcommand{\measure}[1]{*+[F-:<.9em>]{#1} \qw}
-\newcommand{\measuretab}[1]{*{\xy*+<.6em>{#1}="e";"e"+UL;"e"+UR **\dir{-};"e"+DR **\dir{-};"e"+DL **\dir{-};"e"+LC-<.5em,0em> **\dir{-};"e"+UL **\dir{-} \endxy} \qw}
-\newcommand{\measureD}[1]{*{\xy*+=<0em,.1em>{#1}="e";"e"+UR+<0em,.25em>;"e"+UL+<-.5em,.25em> **\dir{-};"e"+DL+<-.5em,-.25em> **\dir{-};"e"+DR+<0em,-.25em> **\dir{-};{"e"+UR+<0em,.25em>\ellipse^{}};"e"+C:,+(0,1)*{} \endxy} \qw}
-\newcommand{\multimeasure}[2]{*+<1em,.9em>{\hphantom{#2}} \qw \POS[0,0].[#1,0];p !C *{#2},p \drop\frm<.9em>{-}}
-\newcommand{\multimeasureD}[2]{*+<1em,.9em>{\hphantom{#2}} \POS [0,0]="i",[0,0].[#1,0]="e",!C *{#2},"e"+UR-<.8em,0em>;"e"+UL **\dir{-};"e"+DL **\dir{-};"e"+DR+<-.8em,0em> **\dir{-};{"e"+DR+<0em,.8em>\ellipse^{}};"e"+UR+<0em,-.8em> **\dir{-};{"e"+UR-<.8em,0em>\ellipse^{}},"i" \qw}
-\newcommand{\control}{*!<0em,.025em>-=-<.2em>{\bullet}}
-\newcommand{\controlo}{*+<.01em>{\xy -<.095em>*\xycircle<.19em>{} \endxy}}
-\newcommand{\ctrl}[1]{\control \qwx[#1] \qw}
-\newcommand{\ctrlo}[1]{\controlo \qwx[#1] \qw}
-\newcommand{\targ}{*+<.02em,.02em>{\xy ="i","i"-<.39em,0em>;"i"+<.39em,0em> **\dir{-}, "i"-<0em,.39em>;"i"+<0em,.39em> **\dir{-},"i"*\xycircle<.4em>{} \endxy} \qw}
-\newcommand{\qswap}{*=<0em>{\times}}
-\newcommand{\multigate}[2]{*+<1em,.9em>{\hphantom{#2}} \POS [0,0]="i",[0,0].[#1,0]="e",!C *{#2},"e"+UR;"e"+UL **\dir{-};"e"+DL **\dir{-};"e"+DR **\dir{-};"e"+UR **\dir{-},"i" \qw}
-\newcommand{\ghost}[1]{*+<1em,.9em>{\hphantom{#1}} \qw}
-\newcommand{\push}[1]{*{#1}}
-\newcommand{\gategroup}[6]{\POS"#1,#2"."#3,#2"."#1,#4"."#3,#4"!C*+<#5>\frm{#6}}
-\newcommand{\rstick}[1]{*!L!<-.5em,0em>=<0em>{#1}}
-\newcommand{\lstick}[1]{*!R!<.5em,0em>=<0em>{#1}}
-\newcommand{\ustick}[1]{*!D!<0em,-.5em>=<0em>{#1}}
-\newcommand{\dstick}[1]{*!U!<0em,.5em>=<0em>{#1}}
-\newcommand{\Qcircuit}{\xymatrix @*=<0em>}
-\newcommand{\link}[2]{\ar @{-} [#1,#2]}
-\newcommand{\pureghost}[1]{*+<1em,.9em>{\hphantom{#1}}}
-"""
diff --git a/src/qutip_qip/device/processor.py b/src/qutip_qip/device/processor.py
index 521aea57..95e13c96 100644
--- a/src/qutip_qip/device/processor.py
+++ b/src/qutip_qip/device/processor.py
@@ -698,19 +698,25 @@ def run_analytically(self, init_state=None, qc=None):
 
         Returns
         -------
-        evo_result : :class:`qutip.Result`
-            An instance of the class
-            :class:`qutip.Result` will be returned.
+        U_list: list
+            A list of propagators obtained for the physical implementation.
         """
         if init_state is not None:
             U_list = [init_state]
         else:
             U_list = []
         tlist = self.get_full_tlist()
-        # TODO replace this by get_complete_coeff
         coeffs = self.get_full_coeffs()
+
+        # Compute drift Hamiltonians
+        H_drift = 0
+        for drift_ham in self.drift.drift_hamiltonians:
+            H_drift += drift_ham.get_qobj(self.dims)
+
+        # Compute control Hamiltonians
         for n in range(len(tlist)-1):
-            H = sum([coeffs[m, n] * self.ctrls[m]
+            H = H_drift + sum(
+                [coeffs[m, n] * self.ctrls[m]
                     for m in range(len(self.ctrls))])
             dt = tlist[n + 1] - tlist[n]
             U = (-1j * H * dt).expm()
diff --git a/src/qutip_qip/noise.py b/src/qutip_qip/noise.py
index 7d5ddd8d..26bebaed 100644
--- a/src/qutip_qip/noise.py
+++ b/src/qutip_qip/noise.py
@@ -1,4 +1,5 @@
 import numbers
+import warnings
 from collections.abc import Iterable
 from copy import deepcopy
 import numpy as np
@@ -70,13 +71,12 @@ class Noise(object):
     The noise object can be added to :class:`.device.Processor` and
     contributes to evolution.
     """
-    def __init__(self):
-        pass
-
-    def get_noisy_dynamics(self, dims, pulses, systematic_noise):
+    def get_noisy_pulses(self, dims=None, pulses=None, systematic_noise=None):
         """
         Return the input pulses list with noise added and
-        the pulse independent noise in a dummy Pulse object.
+        the pulse independent noise in a dummy :class:`.Pulse` object.
+        This is a template method, a method with the same name and signatures
+        needs to be defined in the subclasses.
 
         Parameters
         ----------
@@ -85,7 +85,7 @@ def get_noisy_dynamics(self, dims, pulses, systematic_noise):
             [2,2...,2] for qubits system.
 
         pulses : list of :class:`.Pulse`
-            The input pulses, on which the noise object is to be applied.
+            The input pulses. The noise will be added to pulses in this list.
 
         systematic_noise : :class:`.Pulse`
             The dummy pulse with no ideal control element.
@@ -96,18 +96,26 @@ def get_noisy_dynamics(self, dims, pulses, systematic_noise):
             Noisy pulses.
 
         systematic_noise : :class:`.Pulse`
-            The dummy pulse representing pulse independent noise.
+            The dummy pulse representing pulse-independent noise.
         """
+        get_noisy_dynamics = getattr(self, "get_noisy_dynamics", None)
+        if get_noisy_dynamics is not None:
+            warnings.warn(
+                "Using get_noisy_dynamics as the hook function for custom "
+                "noise will be deprecated, "
+                "please use get_noisy_pulses instead.",
+                PendingDeprecationWarning)
+            return self.get_noisy_dynamics(dims, pulses, systematic_noise)
         raise NotImplementedError(
             "Subclass error needs a method"
-            "`get_noisy_dynamics` to process the noise.")
+            "`get_noisy_pulses` to process the noise.")
 
     def _apply_noise(self, pulses=None, systematic_noise=None, dims=None):
         """
         For backward compatibility, in case the method has no return value
         or only return the pulse.
         """
-        result = self.get_noisy_dynamics(
+        result = self.get_noisy_pulses(
             pulses=pulses, systematic_noise=systematic_noise, dims=dims)
         if result is None:  # in-place change
             pass
@@ -118,7 +126,7 @@ def _apply_noise(self, pulses=None, systematic_noise=None, dims=None):
             pulses = result
         else:
             raise TypeError(
-                "Returned value of get_noisy_dynamics not understood.")
+                "Returned value of get_noisy_pulses not understood.")
         return pulses, systematic_noise
 
 
@@ -172,8 +180,29 @@ def __init__(self, c_ops, targets=None, coeff=None, tlist=None,
                     "thus cannot be applied to all qubits")
         self.all_qubits = all_qubits
 
-    def get_noisy_dynamics(
+    def get_noisy_pulses(
             self, dims=None, pulses=None, systematic_noise=None):
+        """
+        Return the input pulses list with noise added and
+        the pulse independent noise in a dummy :class:`.Pulse` object.
+
+        Parameters
+        ----------
+        dims: list, optional
+            The dimension of the components system, the default value is
+            [2, 2, ..., 2] for a system of qubits.
+        pulses : list of :class:`.Pulse`
+            The input pulses. The noise will be added to pulses in this list.
+        systematic_noise : :class:`.Pulse`
+            The dummy pulse with no ideal control element.
+
+        Returns
+        -------
+        noisy_pulses: list of :class:`.Pulse`
+            Noisy pulses.
+        systematic_noise : :class:`.Pulse`
+            The dummy pulse representing pulse-independent noise.
+        """
         if systematic_noise is None:
             systematic_noise = Pulse(None, None, label="system")
         N = len(dims)
@@ -250,8 +279,29 @@ def _T_to_list(self, T, N):
                 "either the length is not equal to the number of qubits, "
                 "or T is not a positive number.".format(T))
 
-    def get_noisy_dynamics(
+    def get_noisy_pulses(
             self, dims=None, pulses=None, systematic_noise=None):
+        """
+        Return the input pulses list with noise added and
+        the pulse independent noise in a dummy :class:`.Pulse` object.
+
+        Parameters
+        ----------
+        dims: list, optional
+            The dimension of the components system, the default value is
+            [2,2...,2] for qubits system.
+        pulses : list of :class:`.Pulse`
+            The input pulses. The noise will be added to pulses in this list.
+        systematic_noise : :class:`.Pulse`
+            The dummy pulse with no ideal control element.
+
+        Returns
+        -------
+        noisy_pulses: list of :class:`.Pulse`
+            Noisy pulses.
+        systematic_noise : :class:`.Pulse`
+            The dummy pulse representing pulse-independent noise.
+        """
         if systematic_noise is None:
             systematic_noise = Pulse(None, None, label="system")
         N = len(dims)
@@ -319,7 +369,7 @@ def __init__(self, coeff, tlist=None, indices=None):
         self.tlist = tlist
         self.indices = indices
 
-    def get_noisy_dynamics(
+    def get_noisy_pulses(
             self, dims=None, pulses=None, systematic_noise=None):
         if pulses is None:
             pulses = []
@@ -389,8 +439,29 @@ def __init__(
         self.dt = dt
         self.indices = indices
 
-    def get_noisy_dynamics(
+    def get_noisy_pulses(
             self, dims=None, pulses=None, systematic_noise=None):
+        """
+        Return the input pulses list with noise added and
+        the pulse independent noise in a dummy :class:`.Pulse` object.
+
+        Parameters
+        ----------
+        dims: list, optional
+            The dimension of the components system, the default value is
+            [2,2...,2] for qubits system.
+        pulses : list of :class:`.Pulse`
+            The input pulses. The noise will be added to pulses in this list.
+        systematic_noise : :class:`.Pulse`
+            The dummy pulse with no ideal control element.
+
+        Returns
+        -------
+        noisy_pulses: list of :class:`.Pulse`
+            Noisy pulses.
+        systematic_noise : :class:`.Pulse`
+            The dummy pulse representing pulse-independent noise.
+        """
         if pulses is None:
             pulses = []
         if self.indices is None:
@@ -432,8 +503,29 @@ def __init__(
             self, params):
         self.params = params
 
-    def get_noisy_dynamics(
+    def get_noisy_pulses(
             self, dims=None, pulses=None, systematic_noise=None):
+        """
+        Return the input pulses list with noise added and
+        the pulse independent noise in a dummy :class:`.Pulse` object.
+
+        Parameters
+        ----------
+        dims: list, optional
+            The dimension of the components system, the default value is
+            [2,2...,2] for qubits system.
+        pulses : list of :class:`.Pulse`
+            The input pulses. The noise will be added to pulses in this list.
+        systematic_noise : :class:`.Pulse`
+            The dummy pulse with no ideal control element.
+
+        Returns
+        -------
+        noisy_pulses: list of :class:`.Pulse`
+            Noisy pulses.
+        systematic_noise : :class:`.Pulse`
+            The dummy pulse representing pulse-independent noise.
+        """
         J = self.params["J"]
         wr_dr = self.params["wr_dressed"]
         wr = self.params["wr"]
diff --git a/src/qutip_qip/operations/__init__.py b/src/qutip_qip/operations/__init__.py
index 74c71aca..8c7f8456 100644
--- a/src/qutip_qip/operations/__init__.py
+++ b/src/qutip_qip/operations/__init__.py
@@ -30,13 +30,4 @@
 #    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 #    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 ###############################################################################
-from .gates import (
-    rx, ry, rz, x_gate, y_gate, z_gate, s_gate, t_gate,
-    cy_gate, cz_gate, cs_gate, ct_gate, sqrtnot, snot,
-    phasegate, qrot, cphase, cnot,
-    csign, berkeley, swapalpha, swap, iswap, sqrtswap,
-    sqrtiswap, fredkin, molmer_sorensen,
-    toffoli, rotation, controlled_gate,
-    globalphase, hadamard_transform, gate_sequence_product,
-    gate_expand_1toN, gate_expand_2toN, gate_expand_3toN,
-    qubit_clifford_group, expand_operator)
+from .gates import *


Test patch:
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index ad151563..7623a72c 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -4,7 +4,6 @@ on:
   push:
     branches: [ master ]
   pull_request:
-    branches: [ master ]
 
 jobs:
   test:
@@ -32,7 +31,7 @@ jobs:
 
   doctest:
     runs-on: ubuntu-latest
-        
+
     steps:
     - uses: actions/checkout@v2
     - name: Set up Python
@@ -43,9 +42,14 @@ jobs:
       run: |
         python -m pip install --upgrade pip
         python -m pip install numpy scipy cython sphinx matplotlib numpydoc sphinx_rtd_theme
-        pip install git+https://github.com/qutip/qutip.git
+        pip install qutip
         pip install .
-    - name: Test with pytest
+    - name: Test code snippets in the documentation
       run: |
-        cd docs
+        cd doc
         make doctest
+    - name: Test code examples for the pulse paper
+      run: |
+        python -m pip install joblib pytest pytest-custom_exit_code
+        cd doc/pulse-paper
+        pytest *.py --suppress-no-test-exit-code
diff --git a/tests/test_qubitcircuit.py b/tests/test_circuit.py
similarity index 94%
rename from tests/test_qubitcircuit.py
rename to tests/test_circuit.py
index c5fcb05b..304731b2 100644
--- a/tests/test_qubitcircuit.py
+++ b/tests/test_circuit.py
@@ -45,6 +45,7 @@
                    rand_dm, bell_state, ket2dm)
 from qutip_qip.qasm import read_qasm
 from qutip_qip.operations.gates import gate_sequence_product
+import qutip as qp
 
 
 def _op_dist(A, B):
@@ -160,7 +161,6 @@ def testFREDKINdecompose(self):
         U2 = gates.gate_sequence_product(qc2.propagators())
         assert _op_dist(U1, U2) < 1e-12
 
-
     def testadjacentgates(self):
         """
         Adjacent Gates: compare unitary matrix for ISWAP and product of
@@ -396,6 +396,20 @@ def test_exceptions_controlled(self, gate):
                       targets=[1])
         pytest.raises(ValueError, qc.add_gate, gate)
 
+    def test_globalphase_gate_propagators(self):
+        qc = QubitCircuit(2)
+        qc.add_gate("GLOBALPHASE", arg_value=np.pi / 2)
+
+        [gate] = qc.gates
+        assert gate.name == "GLOBALPHASE"
+        assert gate.arg_value == np.pi / 2
+
+        [U_expanded] = qc.propagators()
+        assert U_expanded == 1j * qp.qeye([2, 2])
+
+        [U_unexpanded] = qc.propagators(expand=False)
+        assert U_unexpanded == 1j * qp.qeye([2, 2])
+
     def test_single_qubit_gates(self):
         """
         Text single qubit gates are added correctly
@@ -641,6 +655,19 @@ def test_wstate(self):
                 np.testing.assert_allclose(probs_initial, probs_final)
                 assert sum(result_cbits[i]) == 1
 
-
-if __name__ == "__main__":
-    run_module_suite()
+    def test_latex_code_teleportation_circuit(self):
+        qc = _teleportation_circuit()
+        latex = qc.latex_code()
+        assert latex == "\n".join([
+            r" & \lstick{c1} &  \qw  &  \qw  &  \qw  &  \qw"
+            r"  &  \qw \cwx[4]  &  \qw  &  \qw  &  \ctrl{2}  & \qw \\ ",
+            r" & \lstick{c0} &  \qw  &  \qw  &  \qw  &  \qw"
+            r"  &  \qw  &  \qw \cwx[2]  &  \ctrl{1}  &  \qw  & \qw \\ ",
+            r" & \lstick{\ket{0}} &  \qw  &  \targ  &  \qw  &  \qw"
+            r"  &  \qw  &  \qw  &  \gate{X}  &  \gate{Z}  & \qw \\ ",
+            r" & \lstick{\ket{0}} &  \gate{{\rm H}}  &  \ctrl{-1}  &"
+            r"  \targ  &  \qw  &  \qw  &  \meter &  \qw  &  \qw  & \qw \\ ",
+            r" & \lstick{\ket{q0}} &  \qw  &  \qw  &  \ctrl{-1}  &"
+            r"  \gate{{\rm H}}  &  \meter &  \qw  &  \qw  &  \qw  & \qw \\ ",
+            "",
+        ])
diff --git a/tests/test_noise.py b/tests/test_noise.py
index 35ffe12d..1236e289 100644
--- a/tests/test_noise.py
+++ b/tests/test_noise.py
@@ -1,5 +1,6 @@
 from numpy.testing import assert_, run_module_suite, assert_allclose
 import numpy as np
+import pytest
 
 from qutip import (
     tensor, qeye, sigmaz, sigmax, sigmay, destroy, identity, QobjEvo,
@@ -13,16 +14,6 @@
 from qutip_qip.circuit import QubitCircuit
 
 
-class DriftNoise(Noise):
-    def __init__(self, op):
-        self.qobj = op
-
-    def get_noisy_dynamics(self, dims, pulses, systematic_noise):
-        systematic_noise.add_coherent_noise(self.qobj, 0, coeff=True)
-        # test the backward compatibility
-        return pulses
-
-
 class TestNoise:
     def test_decoherence_noise(self):
         """
@@ -35,15 +26,15 @@ def test_decoherence_noise(self):
         decnoise = DecoherenceNoise(
             sigmaz(), coeff=coeff, tlist=tlist, targets=[1])
         dims = [2] * 2
-        pulses, systematic_noise = decnoise.get_noisy_dynamics(dims=dims)
+        pulses, systematic_noise = decnoise.get_noisy_pulses(dims=dims)
         noisy_qu, c_ops = systematic_noise.get_noisy_qobjevo(dims=dims)
         assert_allclose(c_ops[0].ops[0].qobj, tensor(qeye(2), sigmaz()))
         assert_allclose(c_ops[0].ops[0].coeff, coeff)
         assert_allclose(c_ops[0].tlist, tlist)
 
-        # Time-indenpendent and all qubits
+        # Time-independent and all qubits
         decnoise = DecoherenceNoise(sigmax(), all_qubits=True)
-        pulses, systematic_noise = decnoise.get_noisy_dynamics(dims=dims)
+        pulses, systematic_noise = decnoise.get_noisy_pulses(dims=dims)
         noisy_qu, c_ops = systematic_noise.get_noisy_qobjevo(dims=dims)
         c_ops = [qu.cte for qu in c_ops]
         assert_(tensor([qeye(2), sigmax()]) in c_ops)
@@ -52,7 +43,7 @@ def test_decoherence_noise(self):
         # Time-denpendent and all qubits
         decnoise = DecoherenceNoise(
             sigmax(), all_qubits=True, coeff=coeff*2, tlist=tlist)
-        pulses, systematic_noise = decnoise.get_noisy_dynamics(dims=dims)
+        pulses, systematic_noise = decnoise.get_noisy_pulses(dims=dims)
         noisy_qu, c_ops = systematic_noise.get_noisy_qobjevo(dims=dims)
         assert_allclose(c_ops[0].ops[0].qobj, tensor(sigmax(), qeye(2)))
         assert_allclose(c_ops[0].ops[0].coeff, coeff*2)
@@ -68,7 +59,7 @@ def test_relaxation_noise(self):
         dims = [2] * 3
         relnoise = RelaxationNoise(t1=[1., 1., 1.], t2=None)
         systematic_noise = Pulse(None, None, label="system")
-        pulses, systematic_noise = relnoise.get_noisy_dynamics(
+        pulses, systematic_noise = relnoise.get_noisy_pulses(
             dims=dims, systematic_noise=systematic_noise)
         noisy_qu, c_ops = systematic_noise.get_noisy_qobjevo(dims=dims)
         assert_(len(c_ops) == 3)
@@ -78,7 +69,7 @@ def test_relaxation_noise(self):
         dims = [2] * 2
         relnoise = RelaxationNoise(t1=None, t2=None)
         systematic_noise = Pulse(None, None, label="system")
-        pulses, systematic_noise = relnoise.get_noisy_dynamics(
+        pulses, systematic_noise = relnoise.get_noisy_pulses(
             dims=dims, systematic_noise=systematic_noise)
         noisy_qu, c_ops = systematic_noise.get_noisy_qobjevo(dims=dims)
         assert_(len(c_ops) == 0)
@@ -86,14 +77,14 @@ def test_relaxation_noise(self):
         # only t2
         relnoise = RelaxationNoise(t1=None, t2=[0.2, 0.7])
         systematic_noise = Pulse(None, None, label="system")
-        pulses, systematic_noise = relnoise.get_noisy_dynamics(
+        pulses, systematic_noise = relnoise.get_noisy_pulses(
             dims=dims, systematic_noise=systematic_noise)
         noisy_qu, c_ops = systematic_noise.get_noisy_qobjevo(dims=dims)
         assert_(len(c_ops) == 2)
 
         # t1+t2 and systematic_noise = None
         relnoise = RelaxationNoise(t1=[1., 1.], t2=[0.5, 0.5])
-        pulses, systematic_noise = relnoise.get_noisy_dynamics(dims=dims)
+        pulses, systematic_noise = relnoise.get_noisy_pulses(dims=dims)
         noisy_qu, c_ops = systematic_noise.get_noisy_qobjevo(dims=dims)
         assert_(len(c_ops) == 4)
 
@@ -108,7 +99,7 @@ def test_control_amplitude_noise(self):
         pulses = [Pulse(sigmaz(), 0, tlist, coeff)]
         connoise = ControlAmpNoise(coeff=coeff, tlist=tlist)
         noisy_pulses, systematic_noise = \
-            connoise.get_noisy_dynamics(pulses=pulses)
+            connoise.get_noisy_pulses(pulses=pulses)
         assert_allclose(pulses[0].coherent_noise[0].qobj, sigmaz())
         assert_allclose(noisy_pulses[0].coherent_noise[0].coeff, coeff)
 
@@ -129,7 +120,7 @@ def test_random_noise(self):
         gaussnoise = RandomNoise(
             dt=0.1, rand_gen=np.random.normal, loc=mean, scale=std)
         noisy_pulses, systematic_noise = \
-            gaussnoise.get_noisy_dynamics(pulses=pulses)
+            gaussnoise.get_noisy_pulses(pulses=pulses)
         assert_allclose(noisy_pulses[2].qobj, sigmay())
         assert_allclose(noisy_pulses[1].coherent_noise[0].qobj, sigmax())
         assert_allclose(
@@ -143,7 +134,7 @@ def test_random_noise(self):
         gaussnoise = RandomNoise(lam=0.1, dt=0.2, rand_gen=np.random.poisson)
         assert_(gaussnoise.rand_gen is np.random.poisson)
         noisy_pulses, systematic_noise = \
-            gaussnoise.get_noisy_dynamics(pulses=pulses)
+            gaussnoise.get_noisy_pulses(pulses=pulses)
         assert_allclose(
             noisy_pulses[0].coherent_noise[0].tlist,
             np.linspace(1, 6, int(5/0.2) + 1))
@@ -154,19 +145,6 @@ def test_random_noise(self):
             noisy_pulses[2].coherent_noise[0].tlist,
             np.linspace(1, 6, int(5/0.2) + 1))
 
-    def test_user_defined_noise(self):
-        """
-        Test for the user-defined noise object
-        """
-        dr_noise = DriftNoise(sigmax())
-        proc = Processor(1)
-        proc.add_noise(dr_noise)
-        tlist = np.array([0, np.pi/2.])
-        proc.add_pulse(Pulse(identity(2), 0, tlist, False))
-        result = proc.run_state(init_state=basis(2, 0))
-        assert_allclose(
-            fidelity(result.states[-1], basis(2, 1)), 1, rtol=1.0e-6)
-
     def test_zz_cross_talk(self):
         circuit = QubitCircuit(2)
         circuit.add_gate("X", 0)
@@ -177,3 +155,71 @@ def test_zz_cross_talk(self):
         for pulse in pulses:
             if not isinstance(pulse, Drift) and pulse.label=="systematic_noise":
                 assert(len(pulse.coherent_noise) == 1)
+
+
+class DriftNoise1(Noise):
+    """Standard defintion."""
+    def __init__(self, op):
+        self.qobj = op
+
+    def get_noisy_pulses(self, dims, pulses, systematic_noise):
+        systematic_noise.add_coherent_noise(self.qobj, 0, coeff=True)
+        return pulses, systematic_noise
+
+
+class DriftNoise2(Noise):
+    """Only pulses is returned, no system noise."""
+    def __init__(self, op):
+        self.qobj = op
+
+    def get_noisy_pulses(self, dims, pulses, systematic_noise):
+        systematic_noise.add_coherent_noise(self.qobj, 0, coeff=True)
+        return pulses
+
+
+class DriftNoiseOld1(Noise):
+    """Using get_noisy_dynamics as the hook function."""
+    def __init__(self, op):
+        self.qobj = op
+
+    def get_noisy_dynamics(self, dims, pulses, systematic_noise):
+        systematic_noise.add_coherent_noise(self.qobj, 0, coeff=True)
+        return pulses, systematic_noise
+
+
+class DriftNoiseError1(Noise):
+    """Missing hook function, check the error raised."""
+    def __init__(self, op):
+        self.qobj = op
+
+
+@pytest.mark.parametrize(
+        "noise_class, mode",
+        [(DriftNoise1, "pass"),
+        (DriftNoise2, "pass"),
+        (DriftNoiseOld1, "warning"),
+        (DriftNoiseError1, "error"),
+        ]
+    )
+def test_user_defined_noise(noise_class, mode):
+    """
+    Test for the user-defined noise object
+    """
+    dr_noise = noise_class(sigmax())
+    proc = Processor(1)
+    proc.add_noise(dr_noise)
+    tlist = np.array([0, np.pi/2.])
+    proc.add_pulse(Pulse(identity(2), 0, tlist, False))
+
+    if mode == "warning":
+        with pytest.warns(PendingDeprecationWarning):
+            result = proc.run_state(init_state=basis(2, 0))
+    elif mode == "error":
+        with pytest.raises(NotImplementedError):
+            result = proc.run_state(init_state=basis(2, 0))
+        return
+    else:
+        result = proc.run_state(init_state=basis(2, 0))
+
+    assert_allclose(
+        fidelity(result.states[-1], basis(2, 1)), 1, rtol=1.0e-6)
diff --git a/tests/test_processor.py b/tests/test_processor.py
index 63dccf3e..847e8fe5 100644
--- a/tests/test_processor.py
+++ b/tests/test_processor.py
@@ -324,11 +324,19 @@ def testDrift(self):
         Test for the drift Hamiltonian
         """
         processor = Processor(N=1)
-        processor.add_drift(sigmaz(), 0)
-        tlist = np.array([0., 1., 2.])
-        processor.add_pulse(Pulse(identity(2), 0, tlist, False))
+        processor.add_drift(sigmax() / 2, 0)
+        tlist = np.array([0., np.pi, 2*np.pi, 3*np.pi])
+        processor.add_pulse(Pulse(None, None, tlist, False))
         ideal_qobjevo, _ = processor.get_qobjevo(noisy=True)
-        assert_equal(ideal_qobjevo.cte, sigmaz())
+        assert_equal(ideal_qobjevo.cte, sigmax() / 2)
+
+        init_state = basis(2)
+        propagators = processor.run_analytically()
+        analytical_result = init_state
+        for unitary in propagators:
+            analytical_result = unitary * analytical_result
+        fid = fidelity(sigmax() * init_state, analytical_result)
+        assert((1 - fid) < 1.0e-6)
 
     def testChooseSolver(self):
         # setup and fidelity without noise