Simulations using COMSOL

Once the design has been completed using Qiskit-Metal, the design object can be used to run simulations in COMSOL. Currently, SQDMetal supports:

  • Capacitance matrix simulations

  • RF s-parameter simulations

The appropriate simulation parameters are generated by SQDMetal. The resulting simulation can either be run by opening COMSOL via the generated mph file or directly from the Python environment.

Preparing the COMSOL simulation

The API utilises the mph Python package. In order to run the commands below, COMSOL must be installed.

from SQDMetal.COMSOL.Model import COMSOL_Model
from SQDMetal.COMSOL.SimCapacitance import COMSOL_Simulation_CapMats
from SQDMetal.COMSOL.SimRFsParameter import COMSOL_Simulation_RFsParameters

#Initialise the COMSOL engine (needs to only be run once)
COMSOL_Model.init_engine()

#Instantiate a COMSOL model
cmsl = COMSOL_Model('leModel')

#Create simulations to setup - in this case capacitance matrix and RF s-parameter
sim_capMats = COMSOL_Simulation_CapMats(cmsl)
sim_sParams = COMSOL_Simulation_RFsParameters(cmsl, adaptive='Multiple', modal_min_freq_num_eigs=(5.2e9,25))

#(A) - Initialise model from Qiskit-Metal design object: design
cmsl.initialize_model(design, [sim_capMats, sim_sParams], bottom_grounded=True)

#(B) - Add metallic layers (use 0 if selecting ground plane)
cmsl.add_metallic(1, threshold=1e-10, fuse_threshold=1e-10)
cmsl.add_metallic(2, evap_mode='separate_delete_intersections', group_by_evaporations=True)
#Alternative syntax to combine multiple layers (via a list) in one go. In doing so, smooth_radius can be used to round off the corners
cmsl.add_metallic([3,4], threshold=2e-9, fuse_threshold=2e-9, smooth_radius=0.5e-6, multilayer_fuse=True)

#(C) - OPTIONAL: add ground plane
cmsl.add_ground_plane()

#(D) - Combine all contiguous metals into a single metallic electrodes
cmsl.fuse_all_metals()

#(E) - OPTIONAL: Fix the indices of the different electrodes via components on those electrodes
cmsl.reorder_conds_by_comps(['block0', 'launch0', 'launch1', 'launch2'])

#(F) - OPTIONAL: Set everything within a rectangle to be finely meshed
cmsl.fine_mesh_in_rectangle(-3e-3, -140e-6, -2.9e-3, 167e-6, 2e-7, 5e-6)
#(G) - OPTIONAL: Set conductors within a rectangle to be finely meshed
cmsl.fine_mesh_conductors_in_rectangle(-2.86e-3, -145e-6, -2.53e-3, 167e-6, 2e-7, 5e-6)

#(H) - Build model
cmsl.build_geom_mater_elec_mesh(skip_meshing=True, mesh_structure='Fine')

#(I) - Plot model
cmsl.plot()

#(J) - Save model to file
cmsl.save('Test')

The simulation treats the electrodes as 2D planar objects that live on top of a 3D block representing the dielectric subtrate chip. The chip is enclosed by a rectangular prism that represents the exterior boundaries of the simulation. The simulation automatically accounts for changes in the device geometry due to [shadow evaporations](PVD.md).

Note the following features: - (A) - To initialise the model, one must provide:

  • Qiskit-Metal design object

  • List of COMSOL simulation objects (they all implement COMSOL_Simulation_Base).

  • Optional pad_x, pad_y and pad_z (defaults to 0.5mm) that provide the padding from the chip to the exterior boundaries along the x, y and z axes respectively.

  • Optional bottom_grounded parameter sets the z-padding below the chip to be zero. That is, the chip will sit flush on the exterior boundary on the bottom; which, for all intents and purposes grounds the bottom of the chip.

      • The metallic polygons are added to the top of the chip as planar objects. They are added per layer in the Qiskit-Metal design object. Note that contiguous metallic objects are fused as single electrodes. One must provide:

    • Layer index corresponding to the Qiskit-Metal design object

    • Optional threshold (defaults to -1 where it is disabled) given as a positive length in metres. Idea is that:

      • After fusing all contiguous metallic objects, the resulting polygons may have adjacent points that are too close; this causes unnecessarily fine meshing in COMSOL.

      • If adjacent points in a given polygon electrode are within the given `threshold`, they are combined on the spot.

      • Use this to clean up a polygon by setting a size that is smaller than the smallest actual feature size.

    • Optional fuse_threshold (defaults to 1e-12) given as a positive length in metres. Idea is that:

      • The shapely package does not properly combine polygons in certain cases due to floating-point rounding errors.

      • If adjacent polygon electrodes are within the given `fuse_threshold` of each other, they are combined into a single polygon.

    • Optional evap_mode (defaults ‘separate_delete_below’) and evap_trim (defaults 20e-9) on how different electrodes across different are to be combined or separated. The options are the same as discussed in the usage of plot_layer in the [PVD_Shadows class](PVD.md).

    • Optional group_by_evaporations (defaults False) is relevant when using ‘separate_delete_intersections’ or ‘separate_delete_below’ for evap_mode. If True, all angled evaporation steps that are separated due to some deletion, are grouped together as a single metallic electrode - as discussed later in capacitance matrix simulations.

      • The ground plane is a metallic layer corresponding to the ground cut-outs/subtractions in Qiskit-Metal. Note that it currently does not account for any shadow evaporations for it is assumed that the fills are evaporated using single vertical evaporations.

      • Metallic electrods are added layer-by-layer. In cases of overlap between layers, it is important to treat them as single electrodes contiguous electrodes. Note that:

    • The function fuse_all_metals does not modify the underlying geometry - i.e. overlapping polygons across multiple layers are NOT combined into single polygons.

    • Instead fuse_all_metals combines the polygonal selections used to define the electrodes in the COMSOL simulations.

      • In cases where the simulation is swept, it is useful to fix the indices of different electrodes. This can be done deterministically via the reorder_conds_by_comps function. Note the following:

    • This function just requires a list of the names of Qiskit-Metal components that are attached to different electrodes. The electrodes are indexed according to the order given in this list. If the number of electrodes exceeds the number of supplied names, the remaining electrodes are automatically indexed.

    • The given component must be completely contiguous and must exist entirely within a given electrode. For example, a capacitor is not a good choice as it will split across two COMSOL electrodes (unless the two electrodes are electrically connected via other metallic polygons).

      • Everything on top of the chip within a given rectangular region is meshed via custom fine-meshing parameters. The arguments are:

    • First four parameters are x1, y1, x2 and y2 where (x1,y1) is the bottom-left corner while (x2,y2) is the top-right corner of the rectangular region.

    • Optional minElementSize (defaults 1e-7) sets the minimum mesh-element size of the triangular mesh in metres.

    • Optional maxElementSize (defaults 5e-6) sets the maximum mesh-element size of the triangular mesh in metres.

      • Instead of meshing everything finely, fine_mesh_conductors_in_rectangle is a counterpart to fine_mesh_in_rectangle (with the same arguments) in which only conductive regions falling within the rectangle are finely meshed with the given meshing parameters.

      • Build the geometry, setup the materials (i.e. the substrate dielectric and the metals), setup the individual simulation studies and build the resulting mesh in COMSOL. One may provide:

    • Optional skip_meshing (defaults False). If True, the mesh is setup but not explicitly built.

    • Optional mesh_structure (defults ‘Normal’) to define the general size of the meshing in COMSOL. Can be: ‘Extremely fine’, ‘Extra fine’, ‘Finer’, ‘Fine’, ‘Normal’, ‘Coarse’, ‘Coarser’, ‘Extra coarse’ or ‘Extremely coarse’.

    • Optional substrate_permittivity (defaults 11.7 for silicon) to set for the dielectric substrate chip.

      • The 2D model plot shows the different electrodes (and their selection indices) that have been setup for simulation. It should be used to ensure that the electrodes are properly mapped (e.g. when tuning the different threshold parameters in add_metallic). The plot function may take the following arguments:

    • Optional plot_chip_boundaries (defaults True) to include the chip boundaries if True. Useful to check if the different metallic elements fit within the chip in the simulation.

    • Optional plot_fine_mesh_bounds (defaults True) to include the fine mesh boundaries if True.

      • This command saves the current COMSOL model into an mph file; in this example, the file will be called ‘Test.mph’. Can include a valid folder path; do not include the mph extension as it will be automatically added.

Note that the order of operations given above must be generally satisfied. For example: - fuse_all_metals should be called after adding all metallic regions - reorder_conds_by_comps should be called after fuse_all_metals - Almost everything (including the fine-meshing) should be called before build_geom_mater_elec_mesh. - However, save may be called any time and the model in its current state will be saved.

The following sections will outline the different simulations, where one may directly run a given simulation from the Python environment.

Capacitance matrix simulations

Given a COMSOL model, the DC capacitances between different electrodes can be calculated via the COMSOL_Simulation_CapMats class. The indices of the different electrodes are assigned automatically or manually via the reorder_conds_by_comps command as outlined in the earlier section. The capacitance matrix can be extracted directly from the COMSOL_Simulation_CapMats as follows:

#Must call build_geom_mater_elec_mesh first and have sim_capMats registered when calling initialize_model
capMat = sim_capMats.run()

Note that the exterior boundaries are considered to be ground and held as a zero potential reference. Thus, no electrode should intersect with the exterior boundaries as setting the unit potential will fail if the conductor touches a region of zero potential.

RF s-parameter simulations

Given a COMSOL model, the s-parameters can be calculated between different ports. The COMSOL_Simulation_RFsParameters class is instantiated via the COMSOL model and some optional parameters:

  • adaptive (defaults ‘None’) specifies the kind of simulation to setup:

    • ‘None’ - the simulation is a frequency domain sweep where the solutions for every frequency sequentially. This gives the most accurate global phase, but it is the slowest simulation type.

    • ‘Single’ - the simulation finds the resonant frequency and width and interpolates the frequency scan given this resonant peak. This is fast and accurate if it is guaranteed that there is only a single resonant peak of interest in the s-parameters.

    • ‘Multiple’ - the simulation performs an eigenfrequency scan to find the resonant peaks and then interpolates the frequency scan given these peaks. This is fast and accurate if all resonant peaks within the given frequency range of interest are correctly found.

  • relative_tolerance (defaults 0.01) is only relevant when using adaptive=’Single’ and gives the relative simulation tolerance in the converged values regarding the resonant peak.

  • modal_min_freq_num_eigs (defaults (1e9,7)) is only relevant when using adaptive=’Multiple’:

    • Given as a tuple of the the minimum frequency (in Hertz) and number of expected eigenfrequency peaks.

    • Note that the minimum frequency should be less than or equal to the frequency range that is set for the simulation. All searched eigenfrequencies are assumed to be larger than this frequency.

    • The number of eigenfrequencies should be at least larger than the expected number of resonant peaks.

Before running `build_geom_mater_elec_mesh`, one must define the ports via any of the create_port_ commands:

#(A) - Create ports on LaunchpadWirebond objects
sim_sParams.create_port_CPW_on_Launcher('launch1')
sim_sParams.create_port_CPW_on_Launcher('launch2')

The supported create_port_ commands work as follows: - (A) - create_port_CPW_on_Launcher:

  • Creates a CPW-fed port on a Qiskit-Metal LaunchpadWirebond component at the start of the launcher. The CPW feed extends onto the ground plane.

  • The function requires the name of said LaunchpadWirebond component in the Qiskit-Metal design object.

  • Optionally provide len_launch (defaults 20e-6) to define the length of the CPW port inlet along the length of the launcher.

The order of the s-parameter ports are given in the order of calling the create_port_ commands.

Once the ports are defined and the simulation has been built via build_geom_mater_elec_mesh, one may set the frequency range upon which to calculate the s-parameters:

#Set frequency scan via an linspace-like construct (1000 points between 2GHz and 5GHz in this case):
sim_sParams.set_freq_range(2e9, 5e9, 1000)

#Set frequency scan via absolute frequency values (given as a list or numpy array):
sim_sParams.set_freq_values([2e9,2.5e9,3e9,5e9])

To run the simulation from the Python environment:

freq_s_params = sim_sParams.run()

The result is a numpy array where the rows are: frequencies, S11 values, S21 values, S31 values, …, SN1 values (for N defined ports). That is the excitation currently is set to the first port.

Finally, note that on running the simulation (as shown below), all field solutions are stored for all frequencies except for the individual frequency solutions (i.e. Multi-Modal) when using adaptive=’Multiple’, where only the fields on the ports (i.e. required for s-parameter calculations) are stored.

Magnetic field simulations

Given a COMSOL model, a magnetic field simulation can be done via the COMSOL_Simulation_MagneticField class. This simulation is specifically designed to send in a DC current, observe the resulting magnetic field and subsequently calculate the net flux through some unit area (e.g. a SQUID loop).

Before running `build_geom_mater_elec_mesh`, one must define the following:

sim_BField.set_current_feed_on_CPW_on_Route('flux_line')
sim_BField.set_Bfield_integration_area(coords)

The first command attaches a CPW-friendly current feeder (i.e. feeds the centre line while grounding both adjacent ground-planes simultaneously). The second command defines an area to which the $B_z$ component of the field is integrated. Note that a handy function in the PVD_Shadows class can be used to find the coordinates of the largest interior (i.e. hole) in a SQUID loop after shadow evaporation calculations:

pvd = PVD_Shadows(design)
# pvd.plot_layer(2)
coords, poly = pvd.get_shadow_largest_interior_for_component('jj_qubit')

To get the flux, simply run:

net_flux = sim_BField.run()

The units are in Webers (i.e. flux is the surface integral of $mathbf{B}$).

TODO: Update this upon documenting COMSOL classes