Merge pull request #156 from BoxiLi/improve-circuitqed

Improve the SCQubits simulator

doc/pulse-paper/dj_algorithm.py

@@ -123,11 +123,6 @@
 ax3[9].set_xlabel(r"$t$")
 
 full_tlist = scqubits_processor.get_full_tlist()
-two_qubit_gate_region = np.array([[0, 6, -18], [0, 6, -4], [0, 9, -4], [0, 6, -18], [0, 6, -4], [0, 11, -4], [0, 6, -18], [0, 6, -4], [0, 6, -4]]) * 500
-for i, (point1, point2, point3) in enumerate(two_qubit_gate_region):
-    vmin, vmax = ax3[i].get_ylim()
-    ax3[i].fill_between([full_tlist[point2], full_tlist[point3]], [vmin ,vmin],  [vmax, vmax], color="lightgray", alpha=0.5)
-    ax3[i].vlines([full_tlist[point2], full_tlist[point3]], vmin, vmax, "gray", "--", linewidth=0.8, alpha=0.5)
 
 fig3.tight_layout()
 fig3.savefig("transmon_pulse.pdf")

doc/source/references.bib

@@ -62,4 +62,18 @@ @article{koch2007charge
   pages={042319},
   year={2007},
   publisher={APS}
-}
+}
+@article{motzoi2013,
+  title = {Improving frequency selection of driven pulses using derivative-based transition suppression},
+  author = {Motzoi, F. and Wilhelm, F. K.},
+  journal = {Phys. Rev. A},
+  volume = {88},
+  issue = {6},
+  pages = {062318},
+  numpages = {15},
+  year = {2013},
+  month = {Dec},
+  publisher = {American Physical Society},
+  doi = {10.1103/PhysRevA.88.062318},
+  url = {https://link.aps.org/doi/10.1103/PhysRevA.88.062318}
+}

src/qutip_qip/compiler/circuitqedcompiler.py

@@ -8,7 +8,7 @@
 
 
 class SCQubitsCompiler(GateCompiler):
-    """
+    r"""
     Compiler for :class:`.SCQubits`.
     Compiled pulse strength is in the unit of GHz.
 
@@ -27,6 +27,22 @@ class SCQubitsCompiler(GateCompiler):
         |``params``       | Hardware Parameters   |
         +-----------------+-----------------------+
 
+    For single-qubit gate, we apply the DRAG correction :cite:`motzoi2013`
+
+    .. math::
+
+        \Omega^{x} &= \Omega_0 - \frac{\Omega_0^3}{4 \alpha^2}
+
+        \Omega^{y} &= - \dot{\Omega}_0 / \alpha
+
+        \Omega^{z} &= - \Omega_0**2 / \alpha + \frac{2  \Omega_0^2)}{
+            4 \alpha
+        }
+
+    where :math:`\Omega_0` is the original shape of the pulse.
+    Notice that the :math:`\Omega_0` and its first derivative
+    should be 0 from the starts and the end.
+
     Parameters
     ----------
     num_qubits: int
@@ -79,30 +95,11 @@ def __init__(self, num_qubits, params):
         )
         self.args = {  # Default configuration
             "shape": "hann",
-            "num_samples": 1000,
+            "num_samples": 101,
             "params": self.params,
+            "DRAG": True,
         }
 
-    def _normalized_gauss_pulse(self):
-        """
-        Return a truncated and normalize Gaussian curve.
-        The returned pulse is truncated from a Gaussian distribution with
-        -3*sigma < t < 3*sigma.
-        The amplitude is shifted so that the pulse start from 0.
-        In addition, the pulse is normalized so that
-        the total integral area is 1.
-        """
-        #  td normalization so that the total integral area is 1
-        td = 2.4384880692912567
-        sigma = 1 / 6 * td  # 3 sigma
-        tlist = np.linspace(0, td, 1000)
-        max_pulse = 1 - np.exp(-((0 - td / 2) ** 2) / 2 / sigma**2)
-        coeff = (
-            np.exp(-((tlist - td / 2) ** 2) / 2 / sigma**2)
-            - np.exp(-((0 - td / 2) ** 2) / 2 / sigma**2)
-        ) / max_pulse
-        return tlist, coeff
-
     def _rotation_compiler(self, gate, op_label, param_label, args):
         """
         Single qubit rotation compiler.
@@ -133,9 +130,34 @@ def _rotation_compiler(self, gate, op_label, param_label, args):
             maximum=self.params[param_label][targets[0]],
             area=gate.arg_value / 2.0 / np.pi,
         )
-        pulse_info = [(op_label + str(targets[0]), coeff)]
+        if args["DRAG"]:
+            pulse_info = self._drag_pulse(op_label, coeff, tlist, targets[0])
+        else:
+            pulse_info = [(op_label + str(targets[0]), coeff)]
         return [Instruction(gate, tlist, pulse_info)]
 
+    def _drag_pulse(self, op_label, coeff, tlist, target):
+        dt_coeff = np.gradient(coeff, tlist[1] - tlist[0]) / 2 / np.pi
+        # Y-DRAG
+        alpha = self.params["alpha"][target]
+        y_drag = -dt_coeff / alpha
+        # Z-DRAG
+        z_drag = -(coeff**2) / alpha + (np.sqrt(2) ** 2 * coeff**2) / (
+            4 * alpha
+        )
+        # X-DRAG
+        coeff += -(coeff**3 / (4 * alpha**2))
+
+        pulse_info = [
+            (op_label + str(target), coeff),
+            ("sz" + str(target), z_drag),
+        ]
+        if op_label == "sx":
+            pulse_info.append(("sy" + str(target), y_drag))
+        elif op_label == "sy":
+            pulse_info.append(("sx" + str(target), -y_drag))
+        return pulse_info
+
     def ry_compiler(self, gate, args):
         """
         Compiler for the RZ gate
@@ -205,12 +227,7 @@ def cnot_compiler(self, gate, args):
         result += self.gate_compiler[gate1.name](gate1, args)
 
         zx_coeff = self.params["zx_coeff"][q1]
-        tlist, coeff = self._normalized_gauss_pulse()
-        amplitude = zx_coeff
         area = 1 / 2
-        sign = np.sign(amplitude) * np.sign(area)
-        tlist = tlist / amplitude * area * sign
-        coeff = coeff * amplitude * sign
         coeff, tlist = self.generate_pulse_shape(
             args["shape"], args["num_samples"], maximum=zx_coeff, area=area
         )

src/qutip_qip/compiler/gatecompiler.py

@@ -301,16 +301,16 @@ def _process_idling_tlist(
     ):
         idling_tlist = []
         if pulse_mode == "continuous":
-            # We add sufficient number of zeros at the begining
+            # We add sufficient number of zeros at the beginning
             # and the end of the idling to prevent wrong cubic spline.
             if start_time - last_pulse_time > 3 * step_size:
                 idling_tlist1 = np.linspace(
                     last_pulse_time + step_size / 5,
                     last_pulse_time + step_size,
-                    5,
+                    10,
                 )
                 idling_tlist2 = np.linspace(
-                    start_time - step_size, start_time, 5
+                    start_time - step_size, start_time, 10
                 )
                 idling_tlist.extend([idling_tlist1, idling_tlist2])
             else:
@@ -433,6 +433,14 @@ def generate_pulse_shape(cls, shape, num_samples, maximum=1.0, area=1.0):
     "cosine": np.pi / 2.0,
 }
 
+# Analytically implementing the pulse shape because the Scipy version
+# subjects more to the finite sampling error under interpolation.
+# More analytical shape can be added here.
+_analytical_window = {
+    "hann": lambda t: 1 / 2 - 1 / 2 * np.cos(2 * np.pi * t),
+    "hamming": lambda t: 0.54 - 0.46 * np.cos(2 * np.pi * t),
+}
+
 
 def _normalized_window(shape, num_samples):
     """
@@ -445,6 +453,9 @@ def _normalized_window(shape, num_samples):
     t_max = _default_window_t_max.get(shape, None)
     if t_max is None:
         raise ValueError(f"Window function {shape} is not supported.")
-    coeff = signal.windows.get_window(shape, num_samples)
     tlist = np.linspace(0, t_max, num_samples)
+    if shape in _analytical_window:
+        coeff = _analytical_window["hann"](tlist)
+    else:
+        coeff = signal.windows.get_window(shape, num_samples)
     return coeff, tlist

src/qutip_qip/device/circuitqed.py

@@ -27,6 +27,9 @@ class SCQubits(ModelProcessor):
     for simplicity, we only use a ZX Hamiltonian for
     the two-qubit interaction.
 
+    See the mathematical details in
+    :obj:`.SCQubitsCompiler` and :obj:`.SCQubitsModel`.
+
     Parameters
     ----------
     num_qubits: int
@@ -195,7 +198,7 @@ def _old_index_label_map(self):
         num_qubits = self.num_qubits
         return (
             ["sx" + str(i) for i in range(num_qubits)]
-            + ["sz" + str(i) for i in range(num_qubits)]
+            + ["sy" + str(i) for i in range(num_qubits)]
             + ["zx" + str(i) + str(i + 1) for i in range(num_qubits)]
             + ["zx" + str(i + 1) + str(i) for i in range(num_qubits)]
         )
@@ -220,6 +223,11 @@ def _set_up_controls(self):
             op = destroy_op * (-1.0j) + destroy_op.dag() * 1.0j
             controls["sy" + str(m)] = (2 * np.pi / 2 * op, [m])
 
+        for m in range(num_qubits):
+            destroy_op = destroy(dims[m])
+            op = destroy_op.dag() * destroy_op
+            controls["sz" + str(m)] = (2 * np.pi * op, [m])
+
         for m in range(num_qubits - 1):
             # For simplicity, we neglect leakage in two-qubit gates.
             d1 = dims[m]
@@ -334,6 +342,7 @@ def get_control_latex(self):
         labels = [
             {f"sx{n}": r"$\sigma_x" + f"^{n}$" for n in range(num_qubits)},
             {f"sy{n}": r"$\sigma_y" + f"^{n}$" for n in range(num_qubits)},
+            {f"sz{n}": r"$\sigma_z" + f"^{n}$" for n in range(num_qubits)},
         ]
         label_zx = {}
         for m in range(num_qubits - 1):

tests/test_device.py

@@ -161,7 +161,7 @@ def test_numerical_evolution(
     pytest.param(circuit2, LinearSpinChain, {}, id = "LinearSpinChain"),
     pytest.param(circuit2, CircularSpinChain, {}, id = "CircularSpinChain"),
     # The length of circuit is limited for SCQubits due to leakage
-    pytest.param(circuit, SCQubits, {"omega_single":[0.003]*3}, id = "SCQubits"),
+    pytest.param(circuit, SCQubits, {"omega_single":[0.02]*3}, id = "SCQubits"),
 ])
 @pytest.mark.parametrize(("schedule_mode"), ["ASAP", "ALAP", None])
 def test_numerical_circuit(circuit, device_class, kwargs, schedule_mode):
@@ -210,3 +210,58 @@ def test_pulse_plotting(processor_class):
     processor.load_circuit(qc)
     fig, ax = processor.plot_pulses()
     plt.close(fig)
+
+
+def _compute_propagator(processor, circuit):
+    qevo, _ = processor.get_qobjevo(noisy=True)
+    if parse_version(qutip.__version__) < parse_version("5.dev"):
+        qevo = qevo.to_list()
+        result = qutip.propagator(qevo, t=processor.get_full_tlist())[-1]
+    else:
+        result = qutip.propagator(
+            qevo,
+            t=processor.get_full_tlist(),
+            parallel=False
+            )[-1]
+    return result
+
+
+def test_scqubits_single_qubit_gate():
+    # Check the accuracy of the single-qubit gate for SCQubits.
+    circuit = QubitCircuit(1)
+    circuit.add_gate("X", targets=[0])
+    processor = SCQubits(1, omega_single=0.04)
+    processor.load_circuit(circuit)
+    U = _compute_propagator(processor, circuit)
+    fid = qutip.average_gate_fidelity(
+        qutip.Qobj(U.full()[:2, :2]), qutip.sigmax()
+    )
+    assert pytest.approx(fid, rel=1.0e-6) == 1
+
+
+def test_idling_accuracy():
+    """
+    Check if the switch-on and off of the pulse is implemented correctly.
+    More sampling points may be needed to suppress the interpolated pulse
+    during the idling period.
+    """
+    processor = SCQubits(2, omega_single=0.04)
+    circuit = QubitCircuit(1)
+    circuit.add_gate("X", targets=[0])
+    processor.load_circuit(circuit)
+    U = _compute_propagator(processor, circuit)
+    error_1_gate = 1 - qutip.average_gate_fidelity(
+        qutip.Qobj(U.full()[:2, :2]), qutip.sigmax()
+    )
+
+    circuit = QubitCircuit(2)
+    circuit.add_gate("X", targets=[0])
+    circuit.add_gate("X", targets=[1])
+    # Turning off scheduling to keep the idling.
+    processor.load_circuit(circuit, schedule_mode=False)
+    U = _compute_propagator(processor, circuit)
+    error_2_gate = 1 - qutip.average_gate_fidelity(
+        qutip.Qobj(U.full()[:2, :2]), qutip.sigmax()
+    )
+
+    assert error_2_gate < 2 * error_1_gate