{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to QuTiP\n",
    "Contact: Nathan Shammah nathan.shammah@gmail.com\n",
    "\n",
    "You can find this notebook at https://github.com/nathanshammah/interactive-notebooks\n",
    "\n",
    "**Run this notebook live** at https://mybinder.org/v2/gh/nathanshammah/interactive-notebooks/binder\n",
    "\n",
    "This notebook has been developed for the interactive lecture on [QuTiP](qutip.org) for the [iTHEMS](https://ithems.riken.jp/en/events/quantum-physics-and-science-with-open-source-software-qutip-the-quantum-toolbox-in-python) Quantum Computing and Information (QCoIn) working group."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from qutip import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qutip.about()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantum Objects\n",
    "\n",
    "### Operators and kets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's define a quantum object. We can start from a simple \n",
    "# bosonic destruction operator\n",
    "n_max = 2\n",
    "a = destroy(n_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# in QuTiP, all number (Fock) states are defined on the \"basis\". \n",
    "psi0 = basis(n_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "psi0 # it's a ket!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a*psi0 # destroing the ground state gives the vacuum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.dag()*psi0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "psi0 = basis(n_max)\n",
    "for n in range(1,n_max):\n",
    "    print(\"n = \",n)\n",
    "    psi = a.dag()*psi0/np.sqrt(n)\n",
    "    print(psi)\n",
    "    psi0 = psi\n",
    "psi0 = basis(n_max)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Density matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho0 = psi0*psi0.dag()\n",
    "rho0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho1 = fock(n_max,1)*fock(n_max,1).dag()\n",
    "rho1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mixed = (rho0 + rho1)/2\n",
    "mixed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's define the Schrodinger cat state\n",
    "cat = fock(n_max,0)+fock(n_max,1)\n",
    "cat = cat \n",
    "cat_rho = ket2dm(cat)/2\n",
    "cat_rho"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# trace Tr[rho]\n",
    "print(cat_rho.tr())\n",
    "print(mixed.tr())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# purity Tr[rho*rho]\n",
    "print((cat_rho*cat_rho).tr())\n",
    "print((mixed*mixed).tr())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hamiltonian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can now build a Hamiltonian, \n",
    "# that of the Quantum Harmonic Oscillator \n",
    "# exploring its features with QuTiP\n",
    "n_max = 6\n",
    "a = destroy(n_max)\n",
    "\n",
    "H = a.dag()*a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E_n , Psi_n = H.eigenstates()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E_n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Psi_n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = .6\n",
    "H_driven = a.dag()*a + g*(a + a.dag())\n",
    "En_driven , Psin_driven = H_driven.eigenstates()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(E_n,\"o-\",markersize = 10, label=\"QHO\")\n",
    "plt.plot(En_driven,\"s--\",markersize = 10, label=\"Driven QHO\")\n",
    "\n",
    "plt.ylabel(\"energy\")\n",
    "plt.xlabel(\"eigenvalue\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Open Quantum System: Lindblad Master Equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{eqnarray}\n",
    "\\dot{\\rho}&=&-i\\lbrack H, \\rho\\rbrack +\\kappa\\left(a\\rho a^\\dagger - a^\\dagger a\\rho- \\rho a^\\dagger a\\right)\\\\\n",
    "\\partial_t{\\rho}&=&-i\\lbrack H, \\rho\\rbrack +\\kappa\\mathcal{D}_{[a]}(\\rho)\n",
    "\\end{eqnarray}\n",
    "\n",
    "* Non-unitary dynamics: Non-Hermitian matrix,  $\\partial_t{\\rho}=\\mathcal{L}\\rho$\n",
    "* Contractive map \n",
    "* Linear map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho0 = ket2dm(basis(n_max,n_max-1))\n",
    "t = np.linspace(0,30,300)\n",
    "kappa = 0.1\n",
    "results = mesolve(H,rho0,t,[kappa*a])\n",
    "rho_t = results.states\n",
    "rho_t[0]\n",
    "\n",
    "rho_t[-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us consider a driven, lossy photonic cavity "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = mesolve(H_driven,rho0,t,c_ops=[kappa*a], e_ops=[a,a.dag()*a])\n",
    "at = results.expect[0]\n",
    "nt = results.expect[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(t,nt,\"-\",markersize = 10, label=\"$\\\\langle a^\\dagger a \\\\rangle$\")\n",
    "plt.ylabel(\"photons in cavity\")\n",
    "plt.xlabel(\"time\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cavity QED: Qubit in a cavity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho0qubit = ket2dm(basis(2,1))\n",
    "rho0_tot = tensor(rho0,rho0qubit)\n",
    "Om = 0.5\n",
    "Hint_tot = Om*tensor(a+a.dag(),sigmax())\n",
    "H0_tot = tensor(H,qeye(2))+tensor(qeye(n_max),sigmaz())\n",
    "H_tot=H0_tot+Hint_tot\n",
    "t = np.linspace(0,50,300)\n",
    "kappa = 0.1\n",
    "a_tot = tensor(a,qeye(2))\n",
    "sz_tot = tensor(qeye(n_max),sigmaz())\n",
    "sm_tot = tensor(qeye(n_max),sigmam())\n",
    "results_tot = mesolve(H_tot,rho0_tot,t,\n",
    "                      c_ops=[kappa*a_tot, kappa*sm_tot], \n",
    "                      e_ops= [a_tot.dag()*a_tot,sz_tot])\n",
    "nphot_tot_t = results_tot.expect[0]\n",
    "sz_tot_t = results_tot.expect[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(t,sz_tot_t,\"-\",markersize = 10, \n",
    "         label=\"$\\\\langle \\sigma_z \\\\rangle (t)$\")\n",
    "plt.plot(t,nphot_tot_t,\"--\",markersize = 10, \n",
    "         label=\"$\\\\langle a^\\dagger a \\\\rangle (t)$\")\n",
    "plt.ylabel(\"System excitation\")\n",
    "plt.xlabel(\"time\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bloch Sphere\n",
    "\n",
    "We can visualize a state on the Bloch sphere. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = Bloch()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b2 = Bloch()\n",
    "vec = [0,1,0]\n",
    "b2.add_vectors(vec)\n",
    "b2.render()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = (basis(2,0)+(1+0j)*basis(2,1)).unit() \n",
    "y = (basis(2,0)+(0+1j)*basis(2,1)).unit()\n",
    "z = (basis(2,0)+(0+0j)*basis(2,1)).unit() \n",
    "b3 = Bloch()      \n",
    "b3.add_states([x,y,z])\n",
    "b3.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wigner Function\n",
    "\n",
    "You can find a Gallery of Wigner functions at qutip.org/tutorials and at R.J. Johansson's  http://github.com/jrjohansson/qutip-lectures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_wigner_2d_3d(psi):\n",
    "    #fig, axes = plt.subplots(1, 2, subplot_kw={'projection': '3d'}, figsize=(12, 6))\n",
    "    fig = plt.figure(figsize=(17, 8))\n",
    "    \n",
    "    ax = fig.add_subplot(1, 2, 1)\n",
    "    plot_wigner(psi, fig=fig, ax=ax, alpha_max=6);\n",
    "\n",
    "    ax = fig.add_subplot(1, 2, 2, projection='3d')\n",
    "    plot_wigner(psi, fig=fig, ax=ax, projection='3d', alpha_max=6);\n",
    "    \n",
    "    plt.close(fig)\n",
    "    return fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Coherent state $|\\psi\\rangle=|\\alpha\\rangle$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 30\n",
    "psi = coherent(N, 2.0)\n",
    "plot_wigner_2d_3d(psi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Schrödinger Cat state $|\\psi\\rangle=\\frac{1}{\\sqrt{2}}(|\\alpha\\rangle+|-\\alpha\\rangle)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 30\n",
    "psi = (coherent(N, -2.0) + coherent(N, 2.0)) / np.sqrt(2)\n",
    "plot_wigner_2d_3d(psi)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Squeezed state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "psi = squeeze(N, 0.5) * basis(N, 0)\n",
    "display(plot_wigner_2d_3d(psi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "psi = (coherent(N, -2.0) + coherent(N, -1j) + coherent(N, 1j) + coherent(N, 2.0)).unit()\n",
    "plot_wigner_2d_3d(psi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NN = 8\n",
    "\n",
    "fig, axes = plt.subplots(NN, 1, figsize=(5, 5 * NN), sharex=True, sharey=True) \n",
    "for n in range(NN):\n",
    "    psi = sum([coherent(N, 2*np.exp(2j * np.pi * m / (n + 2))) for m in range(n + 2)]).unit()\n",
    "    plot_wigner(psi, fig=fig, ax=axes[n])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantum Circuits with QuTiP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = QubitCircuit(3, reverse_states=False)\n",
    "q.add_gate(\"TOFFOLI\", controls=[0, 2], targets=[1])\n",
    "q.add_gate(\"CNOT\", targets=[1], controls=[0])\n",
    "display(\"q\",q.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q.add_gate(\"TOFFOLI\", controls=[0, 2], targets=[1])\n",
    "display(\"q  updated\",q.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q2 = QubitCircuit(3, reverse_states=False)\n",
    "q2.add_gate(\"TOFFOLI\", controls=[0, 2], targets=[1])\n",
    "q2.add_gate(\"TOFFOLI\", controls=[0, 2], targets=[1])\n",
    "display(\"q2\",q2.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U = gate_sequence_product(q.propagators())\n",
    "U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qutip.about()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "[1] J. Robert Johansson, Paul D. Nation, and Franco Nori, QuTiP: An open-source Python framework for the dynamics of open quantum systems, Comp. Phys. Comm. 183, 1760 (2012); QuTiP 2: A Python framework for the dynamics of open quantum systems, Comp. Phys. Comm. . 184, 1234 (2013).\n",
    "\n",
    "[2] Nathan Shammah and Shahnawaz Ahmed, The rise of open source in quantum physics research, Nature’s Physics Blog “On your Wavelength”, January 9th, 2019, http://blogs.nature.com/onyourwavelength/2019/01/09/the-rise-of-open-source-in-quantum-physics-research/\n",
    "\n",
    "[3] Nathan Shammah, Neill Lambert, Franco Nori, and Simone De Liberato, Superradiance with local phase-breaking effects, Phys. Rev. A 96, 023863 (2017); Nathan Shammah, Shahnawaz Ahmed, Neill Lambert, Simone De Liberato, and Franco Nori, Open quantum systems with local and collective incoherent processes: Efficient numerical simulation using permutational invariance, Phys. Rev. A 98, 063815 (2018). \n",
    "\n",
    "[4] Nathan Shammah, A Guide to Building Your Open-Source Science Library, https://github.com/nathanshammah/opensource/blob/master/README.md"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernel_info": {
   "name": "python3"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "nteract": {
   "version": "0.14.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}