wrapper.py

import pandas as pd
import numpy as np
import scipy as sci
import setup_helper_functions as helper_functions

# wrapper program
# calculates the numerical/analytical eigenmode solution for a general multi-qubit cavity system
# for which having an arbitrary number of qubits and cavity modes

# Output:
# Hamiltonian matrix
# calculates the Hamiltonian parameters and post-process

# Define the simulation parameters
# specify the paths that have relevant functions to be used in this .m file

# input:
# HFSS calculated cavity mode eigenfrequencies ( simulated for the same number of eigenmodes as in the EPR simulation ),
# HFSS field calculator output of V_terminal for each qubit [for each LJ and CJ]
# FEM calculated qubit eigenfrequencies, [calculated for the range of Fock_trunc]

# simP = simulation parameter. Will hold most of the data throughout calculations.
# Universal constants
simP = dict()
simP['epsilon0'] = 8.854e-12  # F/m
simP['mu_0'] = 4 * np.pi * 1e-7  # H/m
simP['electron_charge'] = -1.602e-19  # in Coulombs
simP['h_bar'] = 1.054571817e-34  # in J*s
simP['h_plank_const'] = simP['h_bar'] * 2 * np.pi
simP['c_freespace'] = 3e8  # m/s
simP['reduced_flux_quantum'] = simP['h_bar'] / (2 * simP['electron_charge'])

# # # input deck # # #
# components of our cQED
# # # cavity information # # #
# define the simulation parameters per subregion
simP['reg'] = np.array([])
# Region 1: cavity 1
simP['reg'] = np.append(simP['reg'], {'reg_Type': 0, 'num_mode': 5, 'omega': [1, 2, 3, 4, 5]})
# Referenced by simP['reg'][0]
# Region 2: port 1
simP['reg'] = np.append(simP['reg'], {'reg_Type': 1, 'num_mode': 1, 'omega': [1]})  # dummy frequency
# Region 3: port 2
simP['reg'] = np.append(simP['reg'], {'reg_Type': 1, 'num_mode': 1, 'omega': [1]})  # dummy frequency

# cavity dimension
simP['reg'][0]['x'] = 22.86e-3  # m
simP['reg'][0]['y'] = 10.16e-3  # m
simP['reg'][0]['z'] = 40e-3  # m

# info about the order of the cavity modes of our geometry
mode_string = np.array(['TE101', 'TE102', 'TE103', 'TE201', 'TE202'])

# probe 1 location
# coax x,z location in origin
simP['reg'][1]['x0'] = 0.5 * simP['reg'][0]['x']  # in m
simP['reg'][1]['y0'] = lambda probe_len: simP['reg'][0]['y'] - probe_len  # in m
simP['reg'][1]['z0'] = 0.25 * simP['reg'][0]['z']  # in m

# probe 2 location
# coax x,z location in origin
simP['reg'][2]['x0'] = 0.5 * simP['reg'][0]['x']  # in m
simP['reg'][2]['y0'] = lambda probe_len: simP['reg'][0]['y'] - probe_len  # in m
simP['reg'][2]['z0'] = simP['reg'][0]['z'] - simP['reg'][1]['z0']  # in m

# probe 1 dimension
simP['reg'][1]['inner_radius'] = 0.05e-3
simP['reg'][1]['outer_radius'] = 1e-3  # 0.9996e-3
simP['reg'][1]['data_radius'] = 2 * 1e-3

# probe 2 dimension
simP['reg'][2]['inner_radius'] = 0.05e-3
simP['reg'][2]['outer_radius'] = 1e-3  # 0.9996e-3
simP['reg'][2]['data_radius'] = 2 * 1e-3

# coaxial probes' length used for applying the microwave perturbation theory
simP['probe_len'] = 0.75e-3  # in m

# material properties in the subregions
# cavity
simP['eps_r'] = 1  # filled with vacuum material
simP['cavity_eps'] = simP['eps_r'] * simP['epsilon0']
simP['cavity_eta'] = np.sqrt(simP['mu_0'] / simP['cavity_eps'])

# constitutive parameters
simP['c0'] = 3e8  # in m/s
simP['mu_0'] = 4 * np.pi * 1e-7  # H/m
simP['eps_0'] = 8.854e-12  # F/m

simP['reg'][0]['eps'] = 1 * simP['eps_0']
simP['reg'][1]['eps'] = 12.92 * simP['eps_0']  # 2.1 * simP['eps_0']
simP['reg'][2]['eps'] = 12.92 * simP['eps_0']  # 2.1 * simP['eps_0']

simP['reg'][0]['eta'] = np.sqrt(simP['mu_0'] / simP['reg'][0]['eps'])
simP['reg'][1]['eta'] = np.sqrt(simP['mu_0'] / simP['reg'][1]['eps'])
simP['reg'][2]['eta'] = np.sqrt(simP['mu_0'] / simP['reg'][2]['eps'])

# # # transmon information # # #

# transmon location
simP['dipole'] = [{}, {}]
simP['dipole'][0]['x0'] = 0.5 * simP['reg'][0]['x']  # in m
simP['dipole'][0]['y0'] = 0.5 * simP['reg'][0]['y']  # in m
simP['dipole'][0]['z0'] = 0.25 * simP['reg'][0]['z']  # in m

simP['dipole'][1]['x0'] = 0.5 * simP['reg'][0]['x']  # in m
simP['dipole'][1]['y0'] = 0.5 * simP['reg'][0]['y']  # in m
simP['dipole'][1]['z0'] = simP['reg'][0]['z'] - simP['dipole'][0]['z0']  # in m

# dipole antenna geometry
simP['theta'] = np.pi / 2  # orientation
simP['dL'] = 1e-3  # the entire length in m
simP['dW'] = 0.25e-3  # the circumference of the cylinder in m
simP['a_radius'] = simP['dW'] / (2 * np.pi)  # cylinder radius in m

# load information in the dipole antenna
# load 1: dipole load capacitance
simP['Cap_load'] = 50e-15  # in F

# load 2: dipole load from Josephson Junction
# load 2 - 1: Josephson Junction linear capacitance
CJ = 0.34e-15  # in F (fF)
# load 2 - 2: Josephson Junction linear inductance
# it is the looping variable in this wrapper program

# total capacitive load in the dipole antenna
Cap_tot_handle = lambda Cap_Bal: CJ + simP['Cap_load'] + Cap_Bal  # total capacitive load

simP['num_t'] = 2  # number of transmons
simP['num_r'] = 3  # number of cavity modes
simP['num_mode'] = simP['num_t'] + simP['num_r']  # total number of modes considered in the system

# the cavity modes considered in our system
simP['cavity_mode_string'] = mode_string[0:simP['num_r']]

# the number of excitations/ fock number states considered per mode
Nfock_modes = np.array([[3, 3, 3, 3, 3]])  # list in the order of qubits then cavity modes
simP = helper_functions.set_nf_subsystems(simP, Nfock_modes)

# # # # # # # # #  wrapper program # # # # # # # # #

# # # # # # # # transmon information # # # # # # # #
# load 2 - 2: Josephson Junction linear inductance
# used as a looping variable in the wrapper program

# Use the precalculated master data pairs between frequencies and LJ
masterlj_mat = sci.io.loadmat(r"precalculated_values/MASTERLJ_SM_analy2_2.mat")
Eval_points_deck = np.array(masterlj_mat.get("Eval_points_deck"))
master_LJ1 = np.array(masterlj_mat.get("master_LJ1"))

# define the qubit frequencies to fix or sweep
qubit1_frequency = 11.5e9  # fix qubit 1
qubit2_f_start = 10.5e9  # sweep qubit 2
qubit2_f_end = 12.5e9

# find the index of the (closest or) exact frequency value as desired from the master frequencies deck
# and use it to find the corresponding LJ values for each qubit
qubit1_f_fix_idx = np.where(Eval_points_deck == qubit1_frequency)[1][0]  # can change code to find the closest frequency as desired from the master frequencies deck
qubit2_f_start_idx = np.where(Eval_points_deck == qubit2_f_start)[1][0]
qubit2_f_end_idx = np.where(Eval_points_deck == qubit2_f_end)[1][0]

Eval_index_LJ1 = qubit1_f_fix_idx  # qubit 1 (has LJ1) - index
Eval_index_LJ2 = np.array([i for i in range(qubit2_f_start_idx, qubit2_f_end_idx + 1, 5)])  # qubit 2 (had LJ2) - index % can change the sampling rate between the frequency range

# Now use index to access the closest frequency value to the desired qubit frequency
Eval_points1 = Eval_points_deck[0][Eval_index_LJ1]  # looping variable for qubit 1
Eval_points2 = Eval_points_deck[0][Eval_index_LJ2]  # looping variable for qubit 2

# Now use index to access the corresponding LJ value to the qubit frequency
master_LJ2 = master_LJ1  # second qubit    # accessing variable inside the loop for qubit 2
master_LJ1 = master_LJ1[0][Eval_index_LJ1]  # accessing variable inside the loop for qubit 1

# user specified peak threshold
threshold = np.array([0.4])  # default is 0.4
# comment out the default command (first command among the two commands) in the post_processing subfunction's peakvalues2 subfunction
# and choose either the second or the third commands to implement the new
# threshold by using the passed input parameter threshold of the
# post_processing subfunction

# analytical solution of the cavity mode eigenfrequencies and the qubit-cavity modes coupling strength (= normalized Vterminal accross the transmons )
# 1. cavity eigenfrequencies
# outputs: cavity resonant frequency in linear and angular
# simP['f_cavity_arr']
# simP['omega_r']
simP = helper_functions.cal_cavity_eigfreq(simP)

# capacitance of the antenna determination
# outputs:
# simP['Cap_Bal']
simP = helper_functions.Capacitance_of_dipole_ant(simP)  # F
Cap_tot = Cap_tot_handle(simP['Cap_Bal'])  # F

# compute Josephson Junction energies, Ec's for each qubits
simP['Cap_tot'] = np.tile(Cap_tot, (1, simP['num_t']))  # F
simP['Ec_arr'] = np.divide(simP['electron_charge']**2, (2 * (simP['Cap_tot'])))

# # # # # 2. the normalized transmon terminal voltages
# input
# dipole location,
# simP['Cap_Bal'], simP['Cap_load'], cavity = simP.['reg'][0], simP['cavity_mode_string']

# output
# simP['normVter_numereigen']: normalized Vterminal - each rows: qubits, each columns: each cavity mode
simP = helper_functions.calculate_normVter(simP)