diff --git a/pycqed/simulations/cz_unitary_simulation.py b/pycqed/simulations/cz_unitary_simulation.py index a448c5ba67..8f154935ff 100644 --- a/pycqed/simulations/cz_unitary_simulation.py +++ b/pycqed/simulations/cz_unitary_simulation.py @@ -5,18 +5,20 @@ import time import numpy as np import qutip as qtp +from pycqed.measurement import detector_functions as det from scipy.interpolate import interp1d +from pycqed.measurement.waveform_control_CC import waveform as wf alpha_q0 = 250e6 * 2*np.pi -w_q0 = 6e9 * 2 * np.pi # Lower frequency qubit -w_q1 = 7e9 * 2 * np.pi # Upper frequency (fluxing) qubit +w_q0 = 6e9 * 2 * np.pi # Lower frequency qubit +w_q1 = 7e9 * 2 * np.pi # Upper frequency (fluxing) qubit -J = 2.5e6 * 2 * np.pi # coupling strength +J = 2.5e6 * 2 * np.pi # coupling strength tlist = np.arange(0, 250e-9, .05e-9) # operators -b = qtp.tensor(qtp.destroy(2), qtp.qeye(3)) # LSB is static qubit +b = qtp.tensor(qtp.destroy(2), qtp.qeye(3)) # LSB is static qubit a = qtp.tensor(qtp.qeye(2), qtp.destroy(3)) n_q0 = a.dag() * a n_q1 = b.dag() * b @@ -38,27 +40,29 @@ def coupled_transmons_hamiltonian(w_q0, w_q1, alpha_q0, J): w_q1 > w_q1 """ - H_0 = w_q0 * n_q0 + w_q1 * n_q1 + \ - 1/2*alpha_q0*(a.dag()*a.dag()*a*a) +\ + H_0 = w_q0 * n_q0 + w_q1 * n_q1 + \ + 1/2*alpha_q0*(a.dag()*a.dag()*a*a) +\ J * (a.dag() + a) * (b + b.dag()) return H_0 -H_0 = coupled_transmons_hamiltonian(w_q0 =w_q0, w_q1=w_q1, alpha_q0=alpha_q0, + +H_0 = coupled_transmons_hamiltonian(w_q0=w_q0, w_q1=w_q1, alpha_q0=alpha_q0, J=J) # U_target = qtp.Qobj([[1, 0, 0, 0, 0, 0], - [0, 1, 0, 0, 0, 0], - [0, 0, -1, 0, 0, 0], - [0, 0, 0, 1, 0, 0], - [0, 0, 0, 0, -1, 0], - [0, 0, 0, 0, 0, 1]], - type='oper', - dims=[[2, 3], [3, 2]]) + [0, 1, 0, 0, 0, 0], + [0, 0, -1, 0, 0, 0], + [0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, -1, 0], + [0, 0, 0, 0, 0, 1]], + type='oper', + dims=[[2, 3], [3, 2]]) U_target._type = 'oper' + def rotating_frame_transformation(U, t: float, - w_q0: float=0, w_q1:float =0): + w_q0: float=0, w_q1: float =0): """ Transforms the frame of the unitary according to U' = U_{RF}*U*U_{RF}^dag @@ -72,27 +76,26 @@ def rotating_frame_transformation(U, t: float, w_q1 (float): freq of frame for q1 """ - U_RF = (1j*w_q0*n_q0*t).expm() * (1j*w_q1*n_q1*t).expm() + U_RF = (1j*w_q0*n_q0*t).expm() * (1j*w_q1*n_q1*t).expm() + + U_prime = U_RF * U # * U_RF(t=0).dag() + return U_prime - U_prime = U_RF * U #* U_RF(t=0).dag() - return U_prime def phases_from_unitary(U): """ Returns the phases from the unitary """ - phi_00 = np.rad2deg(np.angle(U[0, 0])) # expected to equal 0 + phi_00 = np.rad2deg(np.angle(U[0, 0])) # expected to equal 0 phi_01 = np.rad2deg(np.angle(U[1, 1])) phi_10 = np.rad2deg(np.angle(U[3, 3])) phi_11 = np.rad2deg(np.angle(U[4, 4])) - phi_cond = (phi_11 - phi_01 - phi_10 - phi_00)%360 + phi_cond = (phi_11 - phi_01 - phi_10 - phi_00) % 360 return phi_00, phi_01, phi_10, phi_11, phi_cond - - def pro_fid_unitary_compsubspace(U, U_target): # TODO: verify if this is correct (it seems correct) """ @@ -109,6 +112,7 @@ def pro_fid_unitary_compsubspace(U, U_target): return ptrace/dim + def leakage_from_unitary(U): """ Calculates leakage by summing over all in and output states in the @@ -119,15 +123,16 @@ def leakage_from_unitary(U): for i in range(4): for j in range(4): bra_i = qtp.tensor(qtp.ket([i//2], dim=[2]), - qtp.ket([i%2], dim=[3])).dag() + qtp.ket([i % 2], dim=[3])).dag() ket_j = qtp.tensor(qtp.ket([j//2], dim=[2]), - qtp.ket([j%2], dim=[3])) - p = np.abs((bra_i*U*ket_j).data[0,0])**2 + qtp.ket([j % 2], dim=[3])) + p = np.abs((bra_i*U*ket_j).data[0, 0])**2 sump += p - sump/=4 # divide by dimension of comp subspace + sump /= 4 # divide by dimension of comp subspace L1 = 1-sump return L1 + def seepage_from_unitary(U): """ Calculates leakage by summing over all in and output states in the @@ -141,12 +146,13 @@ def seepage_from_unitary(U): qtp.ket([2], dim=[3])).dag() ket_j = qtp.tensor(qtp.ket([j], dim=[2]), qtp.ket([2], dim=[3])) - p = np.abs((bra_i*U*ket_j).data[0,0])**2 + p = np.abs((bra_i*U*ket_j).data[0, 0])**2 sump += p - sump/=2 # divide by dimension of comp subspace + sump /= 2 # divide by dimension of comp subspace L1 = 1-sump return L1 + def pro_fid_from_unitary(U, U_target): """ i and j sum over the states that span the full computation subspace. @@ -168,10 +174,9 @@ def pro_fid_from_unitary(U, U_target): # return F - def simulate_quantities_of_interest(H_0, tlist, eps_vec, sim_step: float=0.1e-9, - verbose:bool=True): + verbose: bool=True): """ Calculates the quantities of interest from the propagator U @@ -193,7 +198,6 @@ def simulate_quantities_of_interest(H_0, tlist, eps_vec, for phase errors. """ - eps_interp = interp1d(tlist, eps_vec, fill_value='extrapolate') # function only exists to wrap @@ -206,14 +210,101 @@ def eps_t(t, args=None): tlist_sim = (np.arange(0, np.max(tlist), sim_step)) t0 = time.time() U_t = qtp.propagator(H_t, tlist_sim) - t1=time.time() + t1 = time.time() if verbose: print('simulation took {:.2f}s'.format(t1-t0)) U_final = U_t[-1] phases = phases_from_unitary(U_final) - phi_cond=phases[-1] + phi_cond = phases[-1] L1 = leakage_from_unitary(U_final) L2 = seepage_from_unitary(U_final) - return {'phi_cond': phi_cond, 'L1':L1, 'L2':L2} - + return {'phi_cond': phi_cond, 'L1': L1, 'L2': L2} + + +class CZ_trajectory(det.Soft_Detector): + def __init__(self, H_0, fluxlutman): + """ + Detector for simulating a CZ trajectory. + Args: + fluxlutman (instr): an instrument that contains the parameters + required to generate the waveform for the trajectory. + """ + super().__init__() + self.value_names = ['Cost func', 'Cond phase', 'L1', 'L2'] + self.value_units = ['a.u.', 'deg', '%', '%'] + self.fluxlutman = fluxlutman + self.H_0 = H_0 + + def acquire_data_point(self, **kw): + tlist = (np.arange(0, self.fluxlutman.cz_length(), + 1/self.fluxlutman.sampling_rate())) + if not self.fluxlutman.czd_double_sided(): + f_pulse = wf.martinis_flux_pulse( + length=self.fluxlutman.cz_length(), + lambda_2=self.fluxlutman.cz_lambda_2(), + lambda_3=self.fluxlutman.cz_lambda_3(), + theta_f=self.fluxlutman.cz_theta_f(), + f_01_max=self.fluxlutman.cz_freq_01_max(), + J2=self.fluxlutman.cz_J2(), + f_interaction=self.fluxlutman.cz_freq_interaction(), + sampling_rate=self.fluxlutman.sampling_rate(), + return_unit='f01') + else: + self.get_f_pulse_double_sided() + + # extract base frequency from the Hamiltonian + w_q1 = np.real(self.H_0[3,3]) + eps_vec = f_pulse - w_q1 + + qoi = simulate_quantities_of_interest( + H_0=self.H_0, + tlist=tlist, eps_vec=eps_vec, + sim_step=1e-9, verbose=False) + + cost_func_val = abs(qoi['phi_cond']-180) + qoi['L1']*100 * 5 + return cost_func_val, qoi['phi_cond'], qoi['L1']*100, qoi['L2']*100 + + def get_f_pulse_double_sided(self): + half_CZ_A = wf.martinis_flux_pulse( + length=self.fluxlutman.cz_length()*self.fluxlutman.czd_length_ratio(), + lambda_2=self.fluxlutman.cz_lambda_2(), + lambda_3=self.fluxlutman.cz_lambda_3(), + theta_f=self.fluxlutman.cz_theta_f(), + f_01_max=self.fluxlutman.cz_freq_01_max(), + J2=self.fluxlutman.cz_J2(), + # E_c=self.fluxlutman.cz_E_c(), + f_interaction=self.fluxlutman.cz_freq_interaction(), + sampling_rate=self.fluxlutman.sampling_rate(), + return_unit='f01') + + # Generate the second CZ pulse. If the params are np.nan, default + # to the main parameter + if not np.isnan(self.fluxlutman.czd_theta_f()): + d_theta_f = self.fluxlutman.czd_theta_f() + else: + d_theta_f = self.fluxlutman.cz_theta_f() + + if not np.isnan(self.fluxlutman.czd_lambda_2()): + d_lambda_2 = self.fluxlutman.czd_lambda_2() + else: + d_lambda_2 = self.fluxlutman.cz_lambda_2() + if not np.isnan(self.fluxlutman.czd_lambda_3()): + d_lambda_3 = self.fluxlutman.czd_lambda_3() + else: + d_lambda_3 = self.fluxlutman.cz_lambda_3() + + half_CZ_B = wf.martinis_flux_pulse( + length=self.fluxlutman.cz_length()*(1-self.fluxlutman.czd_length_ratio()), + lambda_2=d_lambda_2, + lambda_3=d_lambda_3, + theta_f=d_theta_f, + f_01_max=self.fluxlutman.cz_freq_01_max(), + J2=self.fluxlutman.cz_J2(), + f_interaction=self.fluxlutman.cz_freq_interaction(), + sampling_rate=self.fluxlutman.sampling_rate(), + return_unit='f01') + + # N.B. No amp scaling and offset present + f_pulse = np.concatenate([half_CZ_A, half_CZ_B]) + return f_pulse diff --git a/pycqed/tests/simulations/test_cz_simulations.py b/pycqed/tests/simulations/test_cz_simulations.py index 0f717b5390..03d65cd410 100644 --- a/pycqed/tests/simulations/test_cz_simulations.py +++ b/pycqed/tests/simulations/test_cz_simulations.py @@ -9,6 +9,7 @@ from pycqed.simulations import cz_unitary_simulation as czu from pycqed.measurement.waveform_control_CC.waveform import \ martinis_flux_pulse +from pycqed.instrument_drivers.meta_instrument.LutMans import flux_lutman as flm class Test_cz_unitary_simulation(unittest.TestCase): @@ -217,6 +218,41 @@ def test_simulate_quantities_of_interest(self): self.assertAlmostEqual(qoi['phi_cond'], 260.5987691809, places=0) self.assertAlmostEqual(qoi['L1'], 0.001272424, places=3) + def test_simulate_using_detector(self): + + fluxlutman = flm.AWG8_Flux_LutMan('fluxlutman') + + + # Hamiltonian pars + alpha_q0 = 285e6 * 2*np.pi + J = 2.9e6 * 2*np.pi + w_q0 = 4.11e9 * 2*np.pi + w_q1 = 5.42e9 * 2*np.pi + + H_0 = czu.coupled_transmons_hamiltonian(w_q0, w_q1, + alpha_q0=alpha_q0, J=J) + + fluxlutman.cz_J2(J*np.sqrt(2)) + fluxlutman.cz_freq_01_max(w_q1) + fluxlutman.cz_freq_interaction(w_q0+alpha_q0) + fluxlutman.cz_theta_f(80) + fluxlutman.cz_lambda_2(0) + fluxlutman.cz_lambda_3(0) + fluxlutman.sampling_rate(2.4e9) + fluxlutman.cz_length(220e-9) + + + + d = czu.CZ_trajectory(H_0=H_0, fluxlutman=fluxlutman) + vals = d.acquire_data_point() + self.assertAlmostEquals(vals[0], 13.588, places=1) + self.assertAlmostEquals(vals[1], 166.87, places=1) + self.assertAlmostEquals(vals[2], 0.0920, places=1) + self.assertAlmostEquals(vals[3], 0.1841, places=1) + + fluxlutman.close() + + class Test_CZ_single_trajectory_analysis(unittest.TestCase):