#!/usr/bin/env python
# coding: utf-8

# # Lecture 2A - Simulation of a two-qubit gate using a resonator as coupler
# 
# Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
# 
# This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
# 
# This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).

# In[1]:


import matplotlib.pyplot as plt
import numpy as np
from qutip import (about, basis, concurrence, destroy, expect, fidelity,
                   ket2dm, mesolve, phasegate, ptrace, qeye, sigmaz, sqrtiswap,
                   tensor)
from scipy.special import sici

get_ipython().run_line_magic('matplotlib', 'inline')


# ## Parameters

# In[2]:


N = 10

wc = 5.0 * 2 * np.pi
w1 = 3.0 * 2 * np.pi
w2 = 2.0 * 2 * np.pi

g1 = 0.01 * 2 * np.pi
g2 = 0.0125 * 2 * np.pi

tlist = np.linspace(0, 100, 500)

width = 0.5

# resonant SQRT iSWAP gate
T0_1 = 20
T_gate_1 = (1 * np.pi) / (4 * g1)

# resonant iSWAP gate
T0_2 = 60
T_gate_2 = (2 * np.pi) / (4 * g2)


# ### Operators, Hamiltonian and initial state 

# In[3]:


# cavity operators
a = tensor(destroy(N), qeye(2), qeye(2))
n = a.dag() * a

# operators for qubit 1
sm1 = tensor(qeye(N), destroy(2), qeye(2))
sz1 = tensor(qeye(N), sigmaz(), qeye(2))
n1 = sm1.dag() * sm1

# oeprators for qubit 2
sm2 = tensor(qeye(N), qeye(2), destroy(2))
sz2 = tensor(qeye(N), qeye(2), sigmaz())
n2 = sm2.dag() * sm2


# In[4]:


# Hamiltonian using QuTiP
Hc = a.dag() * a
H1 = -0.5 * sz1
H2 = -0.5 * sz2
Hc1 = g1 * (a.dag() * sm1 + a * sm1.dag())
Hc2 = g2 * (a.dag() * sm2 + a * sm2.dag())

H = wc * Hc + w1 * H1 + w2 * H2 + Hc1 + Hc2


# In[5]:


H


# In[6]:


# initial state: start with one of the qubits in its excited state
psi0 = tensor(basis(N, 0), basis(2, 1), basis(2, 0))


# # Ideal two-qubit iSWAP gate

# In[7]:


def step_t(w1, w2, t0, width, t):
    """
    Step function that goes from w1 to w2 at time t0
    as a function of t.
    """
    return w1 + (w2 - w1) * (t > t0)


fig, axes = plt.subplots(1, 1, figsize=(8, 2))
axes.plot(tlist, [step_t(0.5, 1.5, 50, 0.0, t) for t in tlist], "k")
axes.set_ylim(0, 2)
fig.tight_layout()


# In[8]:


def wc_t(t, args=None):
    return wc


def w1_t(t, args=None):
    return (
        w1
        + step_t(0.0, wc - w1, T0_1, width, t)
        - step_t(0.0, wc - w1, T0_1 + T_gate_1, width, t)
    )


def w2_t(t, args=None):
    return (
        w2
        + step_t(0.0, wc - w2, T0_2, width, t)
        - step_t(0.0, wc - w2, T0_2 + T_gate_2, width, t)
    )


H_t = [[Hc, wc_t], [H1, w1_t], [H2, w2_t], Hc1 + Hc2]


# ### Evolve the system

# In[9]:


res = mesolve(H_t, psi0, tlist, [], [])


# ### Plot the results

# In[10]:


fig, axes = plt.subplots(2, 1, sharex=True, figsize=(12, 8))

axes[0].plot(
    tlist,
    np.array(list(map(wc_t, tlist))) / (2 * np.pi),
    "r",
    linewidth=2,
    label="cavity",
)
axes[0].plot(
    tlist,
    np.array(list(map(w1_t, tlist))) / (2 * np.pi),
    "b",
    linewidth=2,
    label="qubit 1",
)
axes[0].plot(
    tlist,
    np.array(list(map(w2_t, tlist))) / (2 * np.pi),
    "g",
    linewidth=2,
    label="qubit 2",
)
axes[0].set_ylim(1, 6)
axes[0].set_ylabel("Energy (GHz)", fontsize=16)
axes[0].legend()

axes[1].plot(tlist, np.real(expect(n, res.states)), "r",
             linewidth=2, label="cavity")
axes[1].plot(tlist, np.real(expect(n1, res.states)), "b",
             linewidth=2, label="qubit 1")
axes[1].plot(tlist, np.real(expect(n2, res.states)), "g",
             linewidth=2, label="qubit 2")
axes[1].set_ylim(0, 1)

axes[1].set_xlabel("Time (ns)", fontsize=16)
axes[1].set_ylabel("Occupation probability", fontsize=16)
axes[1].legend()

fig.tight_layout()


# ### Inspect the final state

# In[11]:


# extract the final state from the result of the simulation
rho_final = res.states[-1]


# In[12]:


# trace out the resonator mode and print the two-qubit density matrix
rho_qubits = ptrace(rho_final, [1, 2])
rho_qubits


# In[13]:


# compare to the ideal result of the sqrtiswap gate (plus phase correction)
# for the current initial state
rho_qubits_ideal = ket2dm(
    tensor(phasegate(0), phasegate(-np.pi / 2))
    * sqrtiswap()
    * tensor(basis(2, 1), basis(2, 0))
)
rho_qubits_ideal


# ### Fidelity and concurrence

# In[14]:


fidelity(rho_qubits, rho_qubits_ideal)


# In[15]:


concurrence(rho_qubits)


# # Dissipative two-qubit iSWAP gate

# ### Define collapse operators that describe dissipation
# 
# 

# In[16]:


kappa = 0.0001
gamma1 = 0.005
gamma2 = 0.005

c_ops = [np.sqrt(kappa) * a, np.sqrt(gamma1) * sm1, np.sqrt(gamma2) * sm2]


# ### Evolve the system

# In[17]:


res = mesolve(H_t, psi0, tlist, c_ops, [])


# ### Plot the results

# In[18]:


fig, axes = plt.subplots(2, 1, sharex=True, figsize=(12, 8))

axes[0].plot(
    tlist,
    np.array(list(map(wc_t, tlist))) / (2 * np.pi),
    "r",
    linewidth=2,
    label="cavity",
)
axes[0].plot(
    tlist,
    np.array(list(map(w1_t, tlist))) / (2 * np.pi),
    "b",
    linewidth=2,
    label="qubit 1",
)
axes[0].plot(
    tlist,
    np.array(list(map(w2_t, tlist))) / (2 * np.pi),
    "g",
    linewidth=2,
    label="qubit 2",
)
axes[0].set_ylim(1, 6)
axes[0].set_ylabel("Energy (GHz)", fontsize=16)
axes[0].legend()

axes[1].plot(tlist, np.real(expect(n, res.states)), "r", linewidth=2,
             label="cavity")
axes[1].plot(tlist, np.real(expect(n1, res.states)), "b", linewidth=2,
             label="qubit 1")
axes[1].plot(tlist, np.real(expect(n2, res.states)), "g", linewidth=2,
             label="qubit 2")
axes[1].set_ylim(0, 1)

axes[1].set_xlabel("Time (ns)", fontsize=16)
axes[1].set_ylabel("Occupation probability", fontsize=16)
axes[1].legend()

fig.tight_layout()


# ### Fidelity and concurrence

# In[19]:


rho_final = res.states[-1]
rho_qubits = ptrace(rho_final, [1, 2])


# In[20]:


fidelity(rho_qubits, rho_qubits_ideal)


# In[21]:


concurrence(rho_qubits)


# # Two-qubit iSWAP gate: Finite pulse rise time

# In[22]:


def step_t(w1, w2, t0, width, t):
    """
    Step function that goes from w1 to w2 at time t0
    as a function of t, with finite rise time defined
    by the parameter width.
    """
    return w1 + (w2 - w1) / (1 + np.exp(-(t - t0) / width))


fig, axes = plt.subplots(1, 1, figsize=(8, 2))
axes.plot(tlist, [step_t(0.5, 1.5, 50, width, t) for t in tlist], "k")
axes.set_ylim(0, 2)
fig.tight_layout()


# ### Evolve the system

# In[23]:


res = mesolve(H_t, psi0, tlist, [], [])


# ### Plot the results

# In[24]:


fig, axes = plt.subplots(2, 1, sharex=True, figsize=(12, 8))

axes[0].plot(
    tlist,
    np.array(list(map(wc_t, tlist))) / (2 * np.pi),
    "r",
    linewidth=2,
    label="cavity",
)
axes[0].plot(
    tlist,
    np.array(list(map(w1_t, tlist))) / (2 * np.pi),
    "b",
    linewidth=2,
    label="qubit 1",
)
axes[0].plot(
    tlist,
    np.array(list(map(w2_t, tlist))) / (2 * np.pi),
    "g",
    linewidth=2,
    label="qubit 2",
)
axes[0].set_ylim(1, 6)
axes[0].set_ylabel("Energy (GHz)", fontsize=16)
axes[0].legend()

axes[1].plot(tlist, np.real(expect(n, res.states)), "r", linewidth=2,
             label="cavity")
axes[1].plot(tlist, np.real(expect(n1, res.states)), "b", linewidth=2,
             label="qubit 1")
axes[1].plot(tlist, np.real(expect(n2, res.states)), "g", linewidth=2,
             label="qubit 2")
axes[1].set_ylim(0, 1)

axes[1].set_xlabel("Time (ns)", fontsize=16)
axes[1].set_ylabel("Occupation probability", fontsize=16)
axes[1].legend()

fig.tight_layout()


# ### Fidelity and concurrence

# In[25]:


rho_final = res.states[-1]
rho_qubits = ptrace(rho_final, [1, 2])


# In[26]:


fidelity(rho_qubits, rho_qubits_ideal)


# In[27]:


concurrence(rho_qubits)


# # Two-qubit iSWAP gate: Finite rise time with overshoot

# In[28]:


def step_t(w1, w2, t0, width, t):
    """
    Step function that goes from w1 to w2 at time t0
    as a function of t, with finite rise time and
    and overshoot defined by the parameter width.
    """

    return w1 + (w2 - w1) * (0.5 + sici((t - t0) / width)[0] / (np.pi))


fig, axes = plt.subplots(1, 1, figsize=(8, 2))
axes.plot(tlist, [step_t(0.5, 1.5, 50, width, t) for t in tlist], "k")
axes.set_ylim(0, 2)
fig.tight_layout()


# ### Evolve the system

# In[29]:


res = mesolve(H_t, psi0, tlist, [], [])


# ### Plot the results

# In[30]:


fig, axes = plt.subplots(2, 1, sharex=True, figsize=(12, 8))

axes[0].plot(
    tlist,
    np.array(list(map(wc_t, tlist))) / (2 * np.pi),
    "r",
    linewidth=2,
    label="cavity",
)
axes[0].plot(
    tlist,
    np.array(list(map(w1_t, tlist))) / (2 * np.pi),
    "b",
    linewidth=2,
    label="qubit 1",
)
axes[0].plot(
    tlist,
    np.array(list(map(w2_t, tlist))) / (2 * np.pi),
    "g",
    linewidth=2,
    label="qubit 2",
)
axes[0].set_ylim(1, 6)
axes[0].set_ylabel("Energy (GHz)", fontsize=16)
axes[0].legend()

axes[1].plot(tlist, np.real(expect(n, res.states)), "r", linewidth=2,
             label="cavity")
axes[1].plot(tlist, np.real(expect(n1, res.states)), "b", linewidth=2,
             label="qubit 1")
axes[1].plot(tlist, np.real(expect(n2, res.states)), "g", linewidth=2,
             label="qubit 2")
axes[1].set_ylim(0, 1)

axes[1].set_xlabel("Time (ns)", fontsize=16)
axes[1].set_ylabel("Occupation probability", fontsize=16)
axes[1].legend()

fig.tight_layout()


# ### Fidelity and concurrence

# In[31]:


rho_final = res.states[-1]
rho_qubits = ptrace(rho_final, [1, 2])


# In[32]:


fidelity(rho_qubits, rho_qubits_ideal)


# In[33]:


concurrence(rho_qubits)


# # Two-qubit iSWAP gate: Finite pulse rise time and dissipation

# In[34]:


# increase the pulse rise time a bit
width = 0.6

# high-Q resonator but dissipative qubits
kappa = 0.00001
gamma1 = 0.005
gamma2 = 0.005

c_ops = [np.sqrt(kappa) * a, np.sqrt(gamma1) * sm1, np.sqrt(gamma2) * sm2]


# ### Evolve the system

# In[35]:


res = mesolve(H_t, psi0, tlist, c_ops, [])


# ### Plot results

# In[36]:


fig, axes = plt.subplots(2, 1, sharex=True, figsize=(12, 8))

axes[0].plot(
    tlist,
    np.array(list(map(wc_t, tlist))) / (2 * np.pi),
    "r",
    linewidth=2,
    label="cavity",
)
axes[0].plot(
    tlist,
    np.array(list(map(w1_t, tlist))) / (2 * np.pi),
    "b",
    linewidth=2,
    label="qubit 1",
)
axes[0].plot(
    tlist,
    np.array(list(map(w2_t, tlist))) / (2 * np.pi),
    "g",
    linewidth=2,
    label="qubit 2",
)
axes[0].set_ylim(1, 6)
axes[0].set_ylabel("Energy (GHz)", fontsize=16)
axes[0].legend()

axes[1].plot(tlist, np.real(expect(n, res.states)), "r", linewidth=2,
             label="cavity")
axes[1].plot(tlist, np.real(expect(n1, res.states)), "b", linewidth=2,
             label="qubit 1")
axes[1].plot(tlist, np.real(expect(n2, res.states)), "g", linewidth=2,
             label="qubit 2")
axes[1].set_ylim(0, 1)

axes[1].set_xlabel("Time (ns)", fontsize=16)
axes[1].set_ylabel("Occupation probability", fontsize=16)
axes[1].legend()

fig.tight_layout()


# ### Fidelity and concurrence

# In[37]:


rho_final = res.states[-1]
rho_qubits = ptrace(rho_final, [1, 2])


# In[38]:


fidelity(rho_qubits, rho_qubits_ideal)


# In[39]:


concurrence(rho_qubits)


# # Two-qubit iSWAP gate: Using tunable resonator and fixed-frequency qubits

# In[40]:


# reduce the rise time
width = 0.25


def wc_t(t, args=None):
    return (
        wc
        - step_t(0.0, wc - w1, T0_1, width, t)
        + step_t(0.0, wc - w1, T0_1 + T_gate_1, width, t)
        - step_t(0.0, wc - w2, T0_2, width, t)
        + step_t(0.0, wc - w2, T0_2 + T_gate_2, width, t)
    )


H_t = [[Hc, wc_t], H1 * w1 + H2 * w2 + Hc1 + Hc2]


# ### Evolve the system

# In[41]:


res = mesolve(H_t, psi0, tlist, c_ops, [])


# ### Plot the results

# In[42]:


fig, axes = plt.subplots(2, 1, sharex=True, figsize=(12, 8))

axes[0].plot(
    tlist,
    np.array(list(map(wc_t, tlist))) / (2 * np.pi),
    "r",
    linewidth=2,
    label="cavity",
)
axes[0].plot(
    tlist,
    np.array(list(map(w1_t, tlist))) / (2 * np.pi),
    "b",
    linewidth=2,
    label="qubit 1",
)
axes[0].plot(
    tlist,
    np.array(list(map(w2_t, tlist))) / (2 * np.pi),
    "g",
    linewidth=2,
    label="qubit 2",
)
axes[0].set_ylim(1, 6)
axes[0].set_ylabel("Energy (GHz)", fontsize=16)
axes[0].legend()

axes[1].plot(tlist, np.real(expect(n, res.states)), "r", linewidth=2,
             label="cavity")
axes[1].plot(tlist, np.real(expect(n1, res.states)), "b", linewidth=2,
             label="qubit 1")
axes[1].plot(tlist, np.real(expect(n2, res.states)), "g", linewidth=2,
             label="qubit 2")
axes[1].set_ylim(0, 1)

axes[1].set_xlabel("Time (ns)", fontsize=16)
axes[1].set_ylabel("Occupation probability", fontsize=16)
axes[1].legend()

fig.tight_layout()


# ### Fidelity and concurrence

# In[43]:


rho_final = res.states[-1]
rho_qubits = ptrace(rho_final, [1, 2])


# In[44]:


fidelity(rho_qubits, rho_qubits_ideal)


# In[45]:


concurrence(rho_qubits)


# ### Software versions

# In[46]:


about()