SFQ_single_and_two_qubit_gate_JK.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: TkAgg\n"
     ]
    }
   ],
   "source": [
    "%matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.constants import pi,h\n",
    "from qutip import *\n",
    "import functools\n",
    "import operator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Two-level single qubit gate with SFQ pulse train"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using a SFQ pulse train in H1(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check a Gaussian pulse train\n",
    "\n",
    "# parameters\n",
    "w_01 = 5 * 2 *pi # in GHz*2*pi \n",
    "tau = 2e-3 # in ns \n",
    "tc = 5*tau # in ns \n",
    "del_theta = 2*pi/360*2 # in rad\n",
    "n = 90\n",
    "\n",
    "t1 = np.linspace(0, n * 2*pi/w_01, 100001)\n",
    "y1 = functools.reduce(operator.add, [ np.exp(-(t1-tc -i*2*pi/w_01)**2/2/tau**2) * (np.heaviside(t1-i*2*pi/w_01,0)- np.heaviside(t1-2*tc-i*2*pi/w_01,0)) for i in range(n)])\n",
    "\n",
    "plt.plot(t1, y1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "def H1_coeff(t, args):\n",
    "    out = functools.reduce(operator.add, \n",
    "        [np.exp(-(t-args['tc']-i*2*pi/args['w_01'])**2/2/args['tau']**2) \n",
    "         * (np.heaviside(t-i*2*pi/args['w_01'],0) + np.heaviside(t+2*args['tc']-i*2*pi/args['w_01'],0)) \n",
    "         for i in range(args['n'])])\n",
    "    \n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tlist_pulseTrain = np.linspace(0, n * 2*pi/w_01, 10001)\n",
    "plt.plot(tlist_pulseTrain, H1_coeff(tlist_pulseTrain, args))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10.0%. Run time:   3.21s. Est. time left: 00:00:00:28\n",
      "20.0%. Run time:   7.03s. Est. time left: 00:00:00:28\n",
      "30.0%. Run time:  10.46s. Est. time left: 00:00:00:24\n",
      "40.0%. Run time:  14.00s. Est. time left: 00:00:00:20\n",
      "50.0%. Run time:  17.13s. Est. time left: 00:00:00:17\n",
      "60.0%. Run time:  20.07s. Est. time left: 00:00:00:13\n",
      "70.0%. Run time:  23.48s. Est. time left: 00:00:00:10\n",
      "80.0%. Run time:  26.76s. Est. time left: 00:00:00:06\n",
      "90.0%. Run time:  30.02s. Est. time left: 00:00:00:03\n",
      "Total run time:  33.24s\n"
     ]
    }
   ],
   "source": [
    "# parameters\n",
    "w_01 = 5 * 2 *pi # in GHz*2*pi \n",
    "tau = 2e-3 # in ns \n",
    "tc = 5*tau # in ns \n",
    "del_theta = 2*pi/360*4 # in rad\n",
    "n = 45\n",
    "\n",
    "# solve Lindblad master equation\n",
    "H0 = w_01/2*(qeye(2) - sigmaz()) \n",
    "H1 = del_theta/2/np.sqrt(2*pi)/tau*sigmay()\n",
    "H = [H0, [H1, H1_coeff]]\n",
    "\n",
    "psi0 = basis(2,0)\n",
    "tlist_pulseTrain = np.linspace(0, n * 2*pi/w_01, 5001)\n",
    "args = {'tau':tau, 'tc': tc, 'w_01': w_01,'n':n}\n",
    "\n",
    "output = mesolve(H, psi0, tlist_pulseTrain, [], [sigmax(), sigmay(),sigmaz()], args, progress_bar=True)  \n",
    "exp_x = output.expect[0]\n",
    "exp_y = output.expect[1]\n",
    "exp_z = output.expect[2]\n",
    "\n",
    "# plot\n",
    "fig,ax = plt.subplots(1,1, figsize=(10,6))\n",
    "ax.plot(tlist_pulseTrain, exp_z)\n",
    "ax.set_xlabel('Time (ns)', fontsize=18)\n",
    "ax.set_ylabel(r'$<\\sigma_z>$', fontsize=18)\n",
    "ax.tick_params(axis='x', labelsize=18)\n",
    "ax.tick_params(axis='y', labelsize=18)\n",
    "ax.grid('on')\n",
    "\n",
    "# # plot in Bloch sphere\n",
    "bloch = Bloch()\n",
    "bloch.add_points([exp_x, exp_y, exp_z], meth='l')\n",
    "bloch.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "bloch.add_points?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using evolution operator : faster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def H1_coeff_1(t, args):\n",
    "    \"\"\" Time-dependant coefficient of H1 Hamiltonian \"\"\"\n",
    "    return np.exp(-t**2/2/args['tau']**2)\n",
    "\n",
    "def U_f(t):\n",
    "    \"\"\" free time evolution operator \"\"\"\n",
    "    return (1j*w_01*t/2*sigmaz()).expm() \n",
    "\n",
    "def U_single_multiplier(U_single, n):\n",
    "    U = U_single\n",
    "    counter = 1\n",
    "    while True:\n",
    "        if counter > n: return\n",
    "        yield U \n",
    "        U *= U_single\n",
    "        counter += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters\n",
    "w_01 = 5 * 2 *pi # in GHz*2*pi \n",
    "tau = 4e-3 # in ns \n",
    "tc = 5*tau # in ns \n",
    "del_theta = 2*pi/360*2 # in rad\n",
    "n = 180\n",
    "\n",
    "# Compute evolution operator over one unit of SFQ pulse train,\n",
    "# i.e., one SFQ + free evolotuion\n",
    "H0 = w_01/2*(qeye(2) - sigmaz()) \n",
    "H1 = del_theta/2/np.sqrt(2*pi)/tau*sigmay() \n",
    "H = [H0, [H1, H1_coeff_1]]\n",
    "\n",
    "tlist = np.linspace(-tc, tc, 101)\n",
    "args = {'tau':tau, 'tc':tc, 'w_01':w_01}\n",
    "\n",
    "U_single_Gauss = propagator(H, tlist, [], args=args) # time evolution during a Gaussian pulse\n",
    "\n",
    "U_single = U_f(2*pi/w_01-tc) * U_single_Gauss[-1] * U_f(-tc) # Eq.26 in PRApplied paper.\n",
    "\n",
    "U_G_list = [i for i in U_single_multiplier(U_single, n)]\n",
    "\n",
    "# U_G_list = [U_single]\n",
    "# for i in range(n-1):\n",
    "#     U_G_list.append(U_single*U_G_list[-1])\n",
    "\n",
    "sigmaz_expect = expect(sigmaz(), U_G_list*basis(2,0)) # expectation value of sigma_z\n",
    "sigmay_expect = expect(sigmay(), U_G_list*basis(2,0))\n",
    "sigmax_expect = expect(sigmax(), U_G_list*basis(2,0))\n",
    "\n",
    "# plot\n",
    "tlist = [2*pi/w_01*i for i in range(n)]\n",
    "fig,ax = plt.subplots(1,1, figsize=(10,6))\n",
    "ax.plot(tlist, sigmaz_expect, label='<z>')\n",
    "ax.plot(tlist, sigmay_expect, label='<y>')\n",
    "ax.plot(tlist, sigmax_expect, label='<x>')\n",
    "ax.set_xlabel('Time (ns)', fontsize=18)\n",
    "ax.set_ylabel(r'$<\\sigma_z>$', fontsize=18)\n",
    "ax.tick_params(axis='x', labelsize=18)\n",
    "ax.tick_params(axis='y', labelsize=18)\n",
    "ax.grid('on')\n",
    "ax.legend(fontsize=18)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "U_single_Gauss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fidelity vs Number of SFQ pulses \n",
    "As in Fig.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters\n",
    "w_01 = 5 * 2 *pi # in GHz*2*pi \n",
    "tau = 4e-3 # in ns \n",
    "tc = 5*tau # in ns \n",
    "theta = pi/2\n",
    "\n",
    "H0 = w_01/2*(qeye(2) - sigmaz()) \n",
    "U_id = (-1j*theta*sigmay()/2).expm()\n",
    "\n",
    "N_list = np.arange(1, 300, 1)\n",
    "fid_list = []\n",
    "UG_list = []\n",
    "for n in N_list:\n",
    "    del_theta = theta / float(n)\n",
    "    H1 = del_theta/2/np.sqrt(2*pi)/tau*sigmay() \n",
    "    H = [H0, [H1, H1_coeff_1]]\n",
    "    \n",
    "    tlist = np.linspace(-tc, tc, 101)\n",
    "    args = {'tau':tau, 'tc':tc, 'w_01':w_01}\n",
    "\n",
    "    U_single_Gauss = propagator(H, tlist, [], args=args) # time evolution during a Gaussian pulse\n",
    "\n",
    "    U_single = U_f(2*pi/w_01-tc) * U_single_Gauss[-1] * U_f(-tc) # Eq.26 in PRApplied paper.\n",
    "\n",
    "    U_G = U_single**n\n",
    "    UG_list.append(U_G) # optional. for debugging\n",
    "    \n",
    "    Fid_n = (np.absolute((U_id.dag()*U_G).tr())**2 + 2) / 6\n",
    "    fid_list = np.append(fid_list,Fid_n)\n",
    "\n",
    "\n",
    "# plot\n",
    "fig,ax = plt.subplots(1,1, figsize=(10,6))\n",
    "ax.plot(N_list, 1 - fid_list)\n",
    "ax.set_xlabel('Number of SFQ pulses', fontsize=18)\n",
    "ax.set_ylabel('Gate error, 1-F', fontsize=18)\n",
    "ax.set_xscale('log')\n",
    "ax.set_yscale('log')\n",
    "ax.tick_params(axis='x', labelsize=18)\n",
    "ax.tick_params(axis='y', labelsize=18)\n",
    "ax.grid('on')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Two-level two-qubit gate with CR\n",
    "Easwar Eq. 3.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters\n",
    "w1 = 5.2 * 2 *pi # in GHz*2*pi \n",
    "w2 = 5 * 2 * pi  # target Q\n",
    "J = 0.003*2*pi\n",
    "Omega = 0.05*2*pi\n",
    "phi = 0 * 2 * pi\n",
    "\n",
    "b = destroy(2)\n",
    "\n",
    "# solve Lindblad master equation\n",
    "H0 = w1*tensor(b.dag()*b, qeye(2)) + w2*tensor(qeye(2), b.dag()*b) \\\n",
    "    + J*(tensor(b.dag(), b) + tensor(b, b.dag()))\n",
    "H1 = tensor((b.dag() + b), qeye(2))\n",
    "H1_coeff = lambda t, args: args['Omega'] * np.cos(args['w2']*t + args['phi'])\n",
    "H = [H0, [H1, H1_coeff]]\n",
    "\n",
    "psi0_g = tensor(basis(2,0), basis(2,0))\n",
    "psi0_e = tensor(basis(2,1), basis(2,0))\n",
    "# tlist = np.linspace(0, 1000, 5001)  # ns\n",
    "tlist = np.arange(0, 2500, 1)  # ns. if 0.5 is used, there is fast oscillation. why?\n",
    "args = {'Omega':Omega, 'w2': w2, 'phi':phi}\n",
    "\n",
    "# output_g = mesolve(H, psi0_g, tlist, [], [tensor(qeye(2), sigmaz()),tensor(qeye(2), sigmax()),tensor(qeye(2), sigmay())], args, progress_bar=True)  \n",
    "# output_e = mesolve(H, psi0_e, tlist, [], [tensor(qeye(2), sigmaz()),tensor(qeye(2), sigmax()),tensor(qeye(2), sigmay())], args, progress_bar=True)  \n",
    "U = propagator(H, tlist, [], args=args, progress_bar=True)\n",
    "x_g = expect(tensor(qeye(2), sigmax()), U*psi0_g)\n",
    "y_g = expect(tensor(qeye(2), sigmay()), U*psi0_g)\n",
    "z_g = expect(tensor(qeye(2), sigmaz()), U*psi0_g)\n",
    "\n",
    "x_e = expect(tensor(qeye(2), sigmax()), U*psi0_e)\n",
    "y_e = expect(tensor(qeye(2), sigmay()), U*psi0_e)\n",
    "z_e = expect(tensor(qeye(2), sigmaz()), U*psi0_e)\n",
    "\n",
    "# in rotating frame\n",
    "R = (-1j*(tensor(b.dag()*b, qeye(2)) + tensor(qeye(2), b.dag()*b))).expm()\n",
    "R*U*psi0_g\n",
    "\n",
    "# plot\n",
    "fig,ax = plt.subplots(1,1, figsize=(10,6))\n",
    "ax.plot(tlist, x_g, tlist,y_g, tlist, z_g)\n",
    "ax.plot(tlist, x_e, tlist,y_e, tlist, z_e)\n",
    "ax.set_xlabel('Time (ns)', fontsize=18)\n",
    "ax.set_ylabel(r'$<\\sigma_z>$', fontsize=18)\n",
    "ax.tick_params(axis='x', labelsize=18)\n",
    "ax.tick_params(axis='y', labelsize=18)\n",
    "ax.grid('on')\n",
    "\n",
    "# plot in Bloch sphere\n",
    "bloch = Bloch()\n",
    "bloch.add_points([x_g, y_g, z_g], meth='l')\n",
    "bloch.add_points([x_e, y_e, z_e], meth='l')\n",
    "bloch.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True\n",
       "Qobj data =\n",
       "[[1. 0. 0. 0.]\n",
       " [0. 1. 0. 0.]\n",
       " [0. 0. 1. 0.]\n",
       " [0. 0. 0. 1.]],\n",
       "       Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = False\n",
       "Qobj data =\n",
       "[[ 9.88799554e-01+4.56005726e-03j -1.31546012e-03+5.35383619e-04j\n",
       "  -8.76821917e-02-1.20683397e-01j  4.47211113e-05+6.15503756e-05j]\n",
       " [-1.31545943e-03+5.35383539e-04j  9.88644874e-01+4.62812187e-03j\n",
       "  -1.02750174e-02-1.41423240e-02j -8.76870078e-02-1.20680004e-01j]\n",
       " [-8.76821620e-02-1.20683319e-01j -1.02750163e-02-1.41423229e-02j\n",
       "   3.01105215e-01-9.41687529e-01j -9.15654659e-04+1.08567396e-03j]\n",
       " [ 4.47208202e-05+6.15496912e-05j -8.76870411e-02-1.20678407e-01j\n",
       "  -9.15643589e-04+1.08566806e-03j  3.01128459e-01-9.41842319e-01j]],\n",
       "       Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = False\n",
       "Qobj data =\n",
       "[[ 9.70828921e-01+3.01801081e-02j -3.40847213e-03+3.51586654e-03j\n",
       "  -2.26176477e-01-7.34876759e-02j  5.06651770e-04+1.64595549e-04j]\n",
       " [-3.40846634e-03+3.51586164e-03j  9.70430520e-01+3.06039949e-02j\n",
       "  -2.61131441e-02-8.48465161e-03j -2.26175625e-01-7.34679210e-02j]\n",
       " [-2.26176220e-01-7.34875575e-02j -2.61131527e-02-8.48465535e-03j\n",
       "  -8.03082142e-01-5.45641798e-01j  6.91259642e-04+4.84783624e-03j]\n",
       " [ 5.06639789e-04+1.64589005e-04j -2.26170535e-01-7.34644521e-02j\n",
       "   6.91270139e-04+4.84771775e-03j -8.03266510e-01-5.46063826e-01j]],\n",
       "       ...,\n",
       "       Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = False\n",
       "Qobj data =\n",
       "[[ 0.64164124-0.53321353j -0.34282982-0.40244146j -0.11087053-0.02526632j\n",
       "   0.10803095+0.02802051j]\n",
       " [-0.34289348-0.40251779j  0.60107302-0.58136651j  0.09532913+0.02514593j\n",
       "  -0.10810233-0.0264766j ]\n",
       " [-0.10852771-0.02469147j  0.09324146+0.02456276j -0.2550187 -0.7914282j\n",
       "   0.49797187-0.1903411j ]\n",
       " [ 0.10602391+0.02746522j -0.10612549-0.02593968j  0.49950508-0.19092124j\n",
       "  -0.29395408-0.7775232j ]],\n",
       "       Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = False\n",
       "Qobj data =\n",
       "[[ 0.64409197-0.50869939j -0.34284893-0.41182508j -0.18277089+0.07472229j\n",
       "   0.04025052-0.01528418j]\n",
       " [-0.34291192-0.41190353j  0.60316957-0.55812695j  0.01887791-0.00652412j\n",
       "  -0.18277412+0.07191035j]\n",
       " [-0.17883659+0.07315748j  0.01842543-0.00636599j -0.81627126+0.01655106j\n",
       "  -0.0335415 -0.53955355j]\n",
       " [ 0.03947307-0.01499016j -0.17935438+0.07062208j -0.03363916-0.54120989j\n",
       "  -0.81477993+0.05881682j]],\n",
       "       Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = False\n",
       "Qobj data =\n",
       "[[ 0.66438163-0.48453672j -0.34002214-0.40939252j -0.10748478+0.17010435j\n",
       "  -0.02214311+0.03618813j]\n",
       " [-0.34008463-0.40947079j  0.62354535-0.53338356j -0.03484704+0.05633135j\n",
       "  -0.10982081+0.16947811j]\n",
       " [-0.10514579+0.16644735j -0.03408873+0.05516326j -0.20482294+0.78942456j\n",
       "  -0.51999285-0.13348814j]\n",
       " [-0.02172662+0.03556566j -0.10773707+0.16632371j -0.52158457-0.13390281j\n",
       "  -0.16462829+0.80131301j]]], dtype=object)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     1
    ]
   },
   "outputs": [],
   "source": [
    "# Diagonalize and then apply CR drving\n",
    "# The result looks the same as without diagonalization\n",
    "def sorted_diag_op(eigen):\n",
    "    \"\"\"\n",
    "    Diagonalizing matrix by the order of eigenvector\n",
    "    \n",
    "    Args:\n",
    "        eigen : output of eigenstates() method in qutip\n",
    "       \n",
    "    Return: a new diagonalizing operator Qobj\n",
    "    \"\"\"\n",
    "       \n",
    "    index_list = [np.argmax(np.absolute(vec)) for vec in eigen[1]]\n",
    "    \n",
    "    X_array = np.column_stack([eigen[1][index_list.index(i)] for i in range(eigen[1].size)])\n",
    "    X = Qobj(X_array)\n",
    "    \n",
    "    return X\n",
    "\n",
    "# parameters\n",
    "w1 = 5.2 * 2 *pi # in GHz*2*pi \n",
    "w2 = 5 * 2 * pi  # target Q\n",
    "J = 0.003*2*pi\n",
    "Omega = 0.05*2*pi\n",
    "phi = 0 * 2 * pi\n",
    "\n",
    "b = destroy(2)\n",
    "\n",
    "# solve Lindblad master equation\n",
    "H0 = w1*tensor(b.dag()*b, qeye(2)) + w2*tensor(qeye(2), b.dag()*b) \\\n",
    "    + J*(tensor(b.dag(), b) + tensor(b, b.dag()))\n",
    "eigen = H0.eigenstates()\n",
    "\n",
    "# find U\n",
    "U = sorted_diag_op(eigen)\n",
    "U.dims = H0.dims\n",
    "\n",
    "# digonalize H0\n",
    "H0_d = U.dag() * H0 * U\n",
    "\n",
    "H1 = U.dag()*tensor((b.dag() + b), qeye(2))*U\n",
    "H1_coeff = lambda t, args: args['Omega'] * np.cos(args['w2']*t + args['phi'])\n",
    "H = [H0_d, [H1, H1_coeff]]\n",
    "\n",
    "psi0_g = tensor(basis(2,0), basis(2,0))\n",
    "psi0_e = tensor(basis(2,1), basis(2,0))\n",
    "# tlist = np.linspace(0, 1000, 5001)  # ns\n",
    "tlist = np.arange(0, 2500, 1)  # ns if 0.5 is used, there is fast oscillation. why?\n",
    "args = {'Omega':Omega, 'w2': w2, 'phi':phi}\n",
    "\n",
    "# output_g = mesolve(H, psi0_g, tlist, [], [tensor(qeye(2), sigmaz()),tensor(qeye(2), sigmax()),tensor(qeye(2), sigmay())], args, progress_bar=True)  \n",
    "# output_e = mesolve(H, psi0_e, tlist, [], [tensor(qeye(2), sigmaz()),tensor(qeye(2), sigmax()),tensor(qeye(2), sigmay())], args, progress_bar=True)  \n",
    "U = propagator(H, tlist, [], args=args, progress_bar=True)\n",
    "x_g = expect(tensor(qeye(2), sigmax()), U*psi0_g)\n",
    "y_g = expect(tensor(qeye(2), sigmay()), U*psi0_g)\n",
    "z_g = expect(tensor(qeye(2), sigmaz()), U*psi0_g)\n",
    "\n",
    "x_e = expect(tensor(qeye(2), sigmax()), U*psi0_e)\n",
    "y_e = expect(tensor(qeye(2), sigmay()), U*psi0_e)\n",
    "z_e = expect(tensor(qeye(2), sigmaz()), U*psi0_e)\n",
    "\n",
    "# plot\n",
    "fig,ax = plt.subplots(1,1, figsize=(10,6))\n",
    "ax.plot(tlist, x_g, tlist,y_g, tlist, z_g)\n",
    "ax.plot(tlist, x_e, tlist,y_e, tlist, z_e)\n",
    "ax.set_xlabel('Time (ns)', fontsize=18)\n",
    "ax.set_ylabel(r'$<\\sigma_z>$', fontsize=18)\n",
    "ax.tick_params(axis='x', labelsize=18)\n",
    "ax.tick_params(axis='y', labelsize=18)\n",
    "ax.grid('on')\n",
    "\n",
    "# plot in Bloch sphere\n",
    "bloch = Bloch()\n",
    "bloch.add_points([x_g, y_g, z_g])\n",
    "bloch.add_points([x_e, y_e, z_e])\n",
    "bloch.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Easwar paper, Eq.3.16, two-level qubit with CR drive. Ideal case.\n",
    "Delta = 0.2 *2*pi\n",
    "Omega = 0.02*2*pi\n",
    "J = 0.005*2*pi\n",
    "\n",
    "tlist = np.arange(0, 2000, 1)  # ns\n",
    "H0 = (Delta - np.sqrt(Delta**2 + Omega**2))/2*tensor(sigmaz(),qeye(2)) \\\n",
    "      - (J*Omega/np.sqrt(Delta**2 + Omega**2))/2*tensor(sigmaz(), sigmax())\n",
    "U = propagator(H0, tlist, [])\n",
    "\n",
    "psi0_g = tensor(basis(2,0), basis(2,0))\n",
    "psi0_e = tensor(basis(2,1), basis(2,0))\n",
    "\n",
    "x_g = expect(tensor(qeye(2), sigmax()), U*psi0_g)\n",
    "y_g = expect(tensor(qeye(2), sigmay()), U*psi0_g)\n",
    "z_g = expect(tensor(qeye(2), sigmaz()), U*psi0_g)\n",
    "\n",
    "x_e = expect(tensor(qeye(2), sigmax()), U*psi0_e)\n",
    "y_e = expect(tensor(qeye(2), sigmay()), U*psi0_e)\n",
    "z_e = expect(tensor(qeye(2), sigmaz()), U*psi0_e)\n",
    "\n",
    "bloch = Bloch()\n",
    "bloch.add_points([x_g, y_g, z_g])\n",
    "bloch.add_points([x_e, y_e, z_e])\n",
    "bloch.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(w1*tensor(b.dag()*b, qeye(2)) + w2*tensor(qeye(2), b.dag()*b) )/2/pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Two-level two-qubit gate with CR using sigmaz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters\n",
    "w1 = 5.2 * 2 *pi # in GHz*2*pi \n",
    "w2 = 5 * 2 * pi  # target Q\n",
    "J = 0.003*2*pi\n",
    "Omega = 0.05*2*pi\n",
    "phi = 0 * 2 * pi\n",
    "\n",
    "sp = sigmap()\n",
    "sm = sigmam()\n",
    "\n",
    "# solve Lindblad master equation\n",
    "H0 = w1/2*tensor(sigmaz(), qeye(2)) + w2/2*tensor(qeye(2), sigmaz()) \\\n",
    "    + J*(tensor(sp, sm) + tensor(sm, sp))\n",
    "H1 = tensor(sigmax(), qeye(2))\n",
    "H1_coeff = lambda t, args: args['Omega'] * np.cos(args['w2']*t + args['phi'])\n",
    "H = [H0, [H1, H1_coeff]]\n",
    "\n",
    "psi0_g = tensor(basis(2,0), basis(2,0))\n",
    "psi0_e = tensor(basis(2,1), basis(2,0))\n",
    "# tlist = np.linspace(0, 1000, 5001)  # ns\n",
    "tlist = np.arange(0, 2500, 1)  # ns if 0.5 is used, there is fast oscillation. why?\n",
    "args = {'Omega':Omega, 'w2': w2, 'phi':phi}\n",
    "\n",
    "# output_g = mesolve(H, psi0_g, tlist, [], [tensor(qeye(2), sigmaz()),tensor(qeye(2), sigmax()),tensor(qeye(2), sigmay())], args, progress_bar=True)  \n",
    "# output_e = mesolve(H, psi0_e, tlist, [], [tensor(qeye(2), sigmaz()),tensor(qeye(2), sigmax()),tensor(qeye(2), sigmay())], args, progress_bar=True)  \n",
    "U = propagator(H, tlist, [], args=args, progress_bar=True)\n",
    "x_g = expect(tensor(qeye(2), sigmax()), U*psi0_g)\n",
    "y_g = expect(tensor(qeye(2), sigmay()), U*psi0_g)\n",
    "z_g = expect(tensor(qeye(2), sigmaz()), U*psi0_g)\n",
    "\n",
    "x_e = expect(tensor(qeye(2), sigmax()), U*psi0_e)\n",
    "y_e = expect(tensor(qeye(2), sigmay()), U*psi0_e)\n",
    "z_e = expect(tensor(qeye(2), sigmaz()), U*psi0_e)\n",
    "\n",
    "# plot\n",
    "fig,ax = plt.subplots(1,1, figsize=(10,6))\n",
    "ax.plot(tlist, x_g, tlist,y_g, tlist, z_g)\n",
    "ax.plot(tlist, x_e, tlist,y_e, tlist, z_e)\n",
    "ax.set_xlabel('Time (ns)', fontsize=18)\n",
    "ax.set_ylabel(r'$<\\sigma_z>$', fontsize=18)\n",
    "ax.tick_params(axis='x', labelsize=18)\n",
    "ax.tick_params(axis='y', labelsize=18)\n",
    "ax.grid('on')\n",
    "\n",
    "# plot in Bloch sphere\n",
    "bloch = Bloch()\n",
    "bloch.add_points([x_g, y_g, z_g])\n",
    "bloch.add_points([x_e, y_e, z_e])\n",
    "bloch.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Two-level two-qubit with SFQ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def H1_coeff(t, args):\n",
    "    out = functools.reduce(operator.add, \n",
    "        [np.exp(-(t-args['tc']-i*2*pi/args['w2'])**2/2/args['tau']**2) \n",
    "         * (np.heaviside(t-i*2*pi/args['w2'],0) + np.heaviside(t+2*args['tc']-i*2*pi/args['w2'],0)) \n",
    "         for i in range(args['n'])])\n",
    "    \n",
    "    return out\n",
    "\n",
    "# parameters\n",
    "# w2 = 5 * 2 *pi # in GHz*2*pi \n",
    "tau = 2e-3 # in ns \n",
    "tc = 5*tau # in ns \n",
    "del_theta = 2*pi/360*2 # in rad\n",
    "n = 90\n",
    "\n",
    "# tlist_pulseTrain = np.linspace(0, n * 2*pi/w_01, 501)\n",
    "# args = {'tau':tau, 'tc': tc, 'w_01': w_01,'n':n}\n",
    "\n",
    "# parameters\n",
    "w1 = 5.2 * 2 *pi # in GHz*2*pi \n",
    "w2 = 5 * 2 * pi  # target Q\n",
    "# J = 0.003*2*pi\n",
    "# Omega = 0.05*2*pi\n",
    "# phi = 0 * 2 * pi\n",
    "\n",
    "b = destroy(2)\n",
    "\n",
    "# solve Lindblad master equation\n",
    "H0 = w1*tensor(b.dag()*b, qeye(2)) + w2*tensor(qeye(2), b.dag()*b) \\\n",
    "    + J*(tensor(b.dag(), b) + tensor(b, b.dag()))\n",
    "H1 =  del_theta/2/np.sqrt(2*pi)/tau*tensor((b.dag() + b), qeye(2))\n",
    "H = [H0, [H1, H1_coeff]]\n",
    "\n",
    "psi0_g = tensor(basis(2,0), basis(2,0))\n",
    "psi0_e = tensor(basis(2,1), basis(2,0))\n",
    "tlist = np.arange(0, np.round(n * 2*pi/w_01), 1)\n",
    "# tlist = np.arange(0, 2500, 1)  # ns if 0.5 is used, there is fast oscillation. why?\n",
    "args = {'tau':tau, 'tc': tc, 'w2': w2,'n':n}\n",
    "\n",
    "# output_g = mesolve(H, psi0_g, tlist, [], [tensor(qeye(2), sigmaz()),tensor(qeye(2), sigmax()),tensor(qeye(2), sigmay())], args, progress_bar=True)  \n",
    "# output_e = mesolve(H, psi0_e, tlist, [], [tensor(qeye(2), sigmaz()),tensor(qeye(2), sigmax()),tensor(qeye(2), sigmay())], args, progress_bar=True)  \n",
    "U = propagator(H, tlist, [], args=args, progress_bar=True)\n",
    "x_g = expect(tensor(qeye(2), sigmax()), U*psi0_g)\n",
    "y_g = expect(tensor(qeye(2), sigmay()), U*psi0_g)\n",
    "z_g = expect(tensor(qeye(2), sigmaz()), U*psi0_g)\n",
    "\n",
    "x_e = expect(tensor(qeye(2), sigmax()), U*psi0_e)\n",
    "y_e = expect(tensor(qeye(2), sigmay()), U*psi0_e)\n",
    "z_e = expect(tensor(qeye(2), sigmaz()), U*psi0_e)\n",
    "\n",
    "# plot\n",
    "fig,ax = plt.subplots(1,1, figsize=(10,6))\n",
    "ax.plot(tlist, x_g, tlist,y_g, tlist, z_g)\n",
    "ax.plot(tlist, x_e, tlist,y_e, tlist, z_e)\n",
    "ax.set_xlabel('Time (ns)', fontsize=18)\n",
    "ax.set_ylabel(r'$<\\sigma_z>$', fontsize=18)\n",
    "ax.tick_params(axis='x', labelsize=18)\n",
    "ax.tick_params(axis='y', labelsize=18)\n",
    "ax.grid('on')\n",
    "\n",
    "# plot in Bloch sphere\n",
    "bloch = Bloch()\n",
    "bloch.add_points([x_g, y_g, z_g])\n",
    "bloch.add_points([x_e, y_e, z_e])\n",
    "bloch.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.round(n * 2*pi/w_01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(H1_coeff(tlist,args))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "propagator?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Options?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.8.2"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "216.333px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}