Skip to content

Commit 9ce766f

Browse files
nikosavolanikosavola
committed
Fixes to circulax notebook
Fix undefined transmon_frequency and bibtex URN; update notebook
1 parent 7e9f73e commit 9ce766f

5 files changed

Lines changed: 176 additions & 110 deletions

File tree

‎docs/bibliography.bib‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -478,7 +478,7 @@ @mastersthesis{Savola2023
478478
year = {2023},
479479
month = feb,
480480
pages = {68+6},
481-
urn = {http://urn.fi/URN:NBN:fi:aalto-202305213270},
481+
urn = {https://urn.fi/URN:NBN:fi:aalto-202305213270},
482482
langid = {english},
483483
school = {Aalto University},
484484
keywords = {Bayesian optimisation,electromagnetic simulation,energy participation ratio,finite element method,superconducting qubit}

‎notebooks/circulax_transmon_optimization.ipynb‎

Lines changed: 87 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -184,27 +184,37 @@
184184
"\n",
185185
"@component(ports=(\"p1\", \"p2\"), states=(\"phi\",))\n",
186186
"def JosephsonJunction( # noqa: N802\n",
187-
" signals: Signals, s: States, Ic: float = 50e-9, R_sub: float = 1e6\n",
187+
" signals: Signals,\n",
188+
" s: States,\n",
189+
" Ic: float = 50e-9,\n",
190+
" R_sub: float = 1e6,\n",
191+
" EJ2_EJ1_ratio: float = 0.0,\n",
188192
") -> tuple[dict, dict]:\n",
189-
" \"\"\"Josephson junction with sinusoidal current-phase relation.\n",
193+
" \"\"\"Josephson junction with multi-harmonic current-phase relation.\n",
190194
"\n",
191195
" The junction is modelled in the flux formulation:\n",
192196
" - The voltage across the junction equals d(flux)/dt, captured by the\n",
193197
" charge (storage) equation: q['phi'] = -(Φ₀/(2π)) * phi.\n",
194-
" - The supercurrent is I = Ic * sin(phi), entered as a flow equation.\n",
198+
" - The supercurrent is I = Ic * [sin(phi) + 2 * (EJ2/EJ1) * sin(2*phi)],\n",
199+
" accounting for higher-order harmonics in inhomogeneous tunnel barriers\n",
200+
" (Willsch et al. 2023).\n",
195201
" - A large parallel sub-gap resistance provides DC bias stability.\n",
196202
"\n",
197203
" Args:\n",
198204
" signals: Port voltages.\n",
199205
" s: Internal states; ``s.phi`` is the gauge-invariant phase.\n",
200-
" Ic: Critical current in amperes.\n",
206+
" Ic: Critical current (1st harmonic) in amperes.\n",
201207
" R_sub: Sub-gap (quasi-particle) resistance in ohms.\n",
208+
" EJ2_EJ1_ratio: Ratio of 2nd to 1st Josephson harmonic energy.\n",
209+
" Typically -0.02 to -0.11 for AlOx junctions.\n",
202210
"\n",
203211
" Returns:\n",
204212
" Tuple of (flow_dict, charge_dict) for the DAE formulation.\n",
205213
" \"\"\"\n",
206-
" # Supercurrent: I = Ic * sin(phi)\n",
207-
" i_jj = Ic * jnp.sin(s.phi)\n",
214+
" # Supercurrent with 2nd harmonic correction\n",
215+
" # I(phi) = (2e/hbar) * [EJ1 * sin(phi) + 2 * EJ2 * sin(2*phi)]\n",
216+
" # Ic = (2e/hbar) * EJ1\n",
217+
" i_jj = Ic * (jnp.sin(s.phi) + 2.0 * EJ2_EJ1_ratio * jnp.sin(2.0 * s.phi))\n",
208218
"\n",
209219
" # Small resistive shunt for numerical stability (sub-gap resistance)\n",
210220
" v_drop = signals.p1 - signals.p2\n",
@@ -246,6 +256,8 @@
246256
"metadata": {},
247257
"outputs": [],
248258
"source": [
259+
"from qpdk.models.capacitor import plate_capacitor_capacitance_analytical\n",
260+
"\n",
249261
"# Substrate properties\n",
250262
"ε_r_substrate = 11.45 # silicon relative permittivity\n",
251263
"\n",
@@ -256,29 +268,34 @@
256268
" pad_gap_um: float = 15.0,\n",
257269
" jj_area_um2: float = 0.04,\n",
258270
" Jc_A_per_cm2: float = 100.0,\n",
271+
" EJ2_EJ1_ratio: float = -0.05,\n",
259272
") -> dict:\n",
260-
" \"\"\"Convert layout dimensions to circuit parameters.\n",
273+
" \"\"\"Convert layout dimensions to circuit parameters using accurate analytical models.\n",
261274
"\n",
262275
" Args:\n",
263276
" pad_width_um: Capacitor pad width in μm.\n",
264277
" pad_height_um: Capacitor pad height in μm.\n",
265278
" pad_gap_um: Gap between pads in μm.\n",
266279
" jj_area_um2: Josephson junction area in μm².\n",
267280
" Jc_A_per_cm2: Critical current density in A/cm².\n",
281+
" EJ2_EJ1_ratio: Ratio of 2nd to 1st Josephson harmonic energy.\n",
268282
"\n",
269283
" Returns:\n",
270-
" Dictionary with Cs (shunt capacitance) and Ic (critical current).\n",
284+
" Dictionary with Cs (shunt capacitance), Ic (critical current), and EJ2_ratio.\n",
271285
" \"\"\"\n",
272-
" # Shunt capacitance (parallel-plate estimate for two pads)\n",
273-
" pad_area_m2 = (pad_width_um * 1e-6) * (pad_height_um * 1e-6)\n",
274-
" gap_m = pad_gap_um * 1e-6\n",
275-
" Cs = ε_0 * ε_r_substrate * pad_area_m2 / gap_m\n",
286+
" # Shunt capacitance using conformal mapping for coplanar pads\n",
287+
" Cs = plate_capacitor_capacitance_analytical(\n",
288+
" length=pad_height_um,\n",
289+
" width=pad_width_um,\n",
290+
" gap=pad_gap_um,\n",
291+
" ep_r=ε_r_substrate,\n",
292+
" )\n",
276293
"\n",
277294
" # Critical current from JJ area\n",
278295
" jj_area_cm2 = jj_area_um2 * 1e-8 # μm² to cm²\n",
279296
" Ic = Jc_A_per_cm2 * jj_area_cm2\n",
280297
"\n",
281-
" return {\"Cs\": Cs, \"Ic\": Ic}\n",
298+
" return {\"Cs\": Cs, \"Ic\": Ic, \"EJ2_ratio\": EJ2_EJ1_ratio}\n",
282299
"\n",
283300
"\n",
284301
"# Example: default transmon layout\n",
@@ -332,12 +349,13 @@
332349
"f_drive = f_target # drive at target frequency\n",
333350
"\n",
334351
"\n",
335-
"def build_transmon_netlist(Cs: float, Ic: float) -> dict:\n",
352+
"def build_transmon_netlist(Cs: float, Ic: float, EJ2_ratio: float = 0.0) -> dict:\n",
336353
" \"\"\"Build a driven transmon circuit netlist.\n",
337354
"\n",
338355
" Args:\n",
339356
" Cs: Shunt capacitance in farads.\n",
340357
" Ic: Junction critical current in amperes.\n",
358+
" EJ2_ratio: Ratio of 2nd to 1st Josephson harmonic energy.\n",
341359
"\n",
342360
" Returns:\n",
343361
" Circulax-compatible netlist dictionary.\n",
@@ -355,7 +373,7 @@
355373
" },\n",
356374
" \"JJ1\": {\n",
357375
" \"component\": \"josephson_junction\",\n",
358-
" \"settings\": {\"Ic\": Ic, \"R_sub\": 1e6},\n",
376+
" \"settings\": {\"Ic\": Ic, \"R_sub\": 1e6, \"EJ2_EJ1_ratio\": EJ2_ratio},\n",
359377
" },\n",
360378
" \"Cs1\": {\n",
361379
" \"component\": \"capacitor\",\n",
@@ -447,7 +465,9 @@
447465
"cell_type": "code",
448466
"execution_count": null,
449467
"id": "15",
450-
"metadata": {},
468+
"metadata": {
469+
"lines_to_next_cell": 2
470+
},
451471
"outputs": [],
452472
"source": [
453473
"# Set up harmonic balance at the drive frequency\n",
@@ -473,7 +493,9 @@
473493
{
474494
"cell_type": "markdown",
475495
"id": "16",
476-
"metadata": {},
496+
"metadata": {
497+
"lines_to_next_cell": 2
498+
},
477499
"source": [
478500
"### 1.7 Gradient-Based Optimization\n",
479501
"\n",
@@ -499,73 +521,84 @@
499521
"metadata": {},
500522
"outputs": [],
501523
"source": [
502-
"\n",
503-
"\n",
504-
"def transmon_frequency(Ic: float, Cs: float) -> float:\n",
505-
" \"\"\"Compute approximate transmon transition frequency.\n",
506-
"\n",
507-
" Uses the standard expression from Koch et al. (2007).\n",
524+
"def transmon_anharmonicity(Ic: float, Cs: float) -> float: # noqa: ARG001\n",
525+
" \"\"\"Compute transmon anharmonicity α ≈ -E_C/ℏ.\n",
508526
"\n",
509527
" Args:\n",
510-
" Ic: Critical current in amperes.\n",
528+
" Ic: Critical current in amperes (unused, kept for API symmetry).\n",
511529
" Cs: Total shunt capacitance in farads.\n",
512530
"\n",
513531
" Returns:\n",
514-
" Transition frequency in Hz.\n",
532+
" Anharmonicity in Hz (negative).\n",
515533
" \"\"\"\n",
516-
" E_J_val = Φ_0 * Ic / (2.0 * jnp.pi)\n",
517534
" E_C_val = e**2 / (2.0 * Cs)\n",
518-
" # Transmon frequency: f01 ≈ (√(8 E_J E_C) - E_C) / h\n",
519-
" return (jnp.sqrt(8.0 * E_J_val * E_C_val) - E_C_val) / h\n",
535+
" return -E_C_val / h\n",
520536
"\n",
521537
"\n",
522-
"def transmon_anharmonicity(Ic: float, Cs: float) -> float: # noqa: ARG001\n",
523-
" \"\"\"Compute transmon anharmonicity α ≈ -E_C/ℏ.\n",
538+
"def transmon_frequency(Ic: float, Cs: float) -> float:\n",
539+
" \"\"\"Compute transmon frequency f01 ≈ (sqrt(8*Ej*Ec) - Ec) / h.\n",
524540
"\n",
525541
" Args:\n",
526-
" Ic: Critical current in amperes (unused, kept for API symmetry).\n",
542+
" Ic: Critical current in amperes.\n",
527543
" Cs: Total shunt capacitance in farads.\n",
528544
"\n",
529545
" Returns:\n",
530-
" Anharmonicity in Hz (negative).\n",
546+
" Transition frequency in Hz.\n",
531547
" \"\"\"\n",
532548
" E_C_val = e**2 / (2.0 * Cs)\n",
533-
" return -E_C_val / h\n",
549+
" E_J_val = Φ_0 * Ic / (2.0 * jnp.pi)\n",
550+
" return (jnp.sqrt(8.0 * E_J_val * E_C_val) - E_C_val) / h\n",
534551
"\n",
535552
"\n",
536-
"def loss_fn(params_vec: jnp.ndarray) -> float:\n",
537-
" \"\"\"Loss function: squared error from target frequency and anharmonicity.\n",
553+
"def hb_loss_fn(params_vec: jnp.ndarray) -> float:\n",
554+
" \"\"\"Loss function: squared error from target frequency derived from HB simulation.\n",
538555
"\n",
539556
" Args:\n",
540557
" params_vec: Array [log(Ic), log(Cs)] (log-space for better conditioning).\n",
541558
"\n",
542559
" Returns:\n",
543-
" Scalar loss value.\n",
560+
" Scalar loss value based on the simulated 1st harmonic voltage.\n",
544561
" \"\"\"\n",
545562
" Ic = jnp.exp(params_vec[0])\n",
546563
" Cs = jnp.exp(params_vec[1])\n",
547564
"\n",
548-
" f01 = transmon_frequency(Ic, Cs)\n",
549-
" alpha = transmon_anharmonicity(Ic, Cs)\n",
565+
" # Rebuild netlist with current parameters\n",
566+
" net = build_transmon_netlist(Cs, Ic, EJ2_ratio=params[\"EJ2_ratio\"])\n",
567+
" grps, n_vars, pmap = compile_netlist(net, models)\n",
568+
"\n",
569+
" # We need a custom DC solver step here inside the JIT-able loss fn\n",
570+
" # For a simple LC circuit with no DC drive, the DC operating point is exactly zero\n",
571+
" y_dc_current = jnp.zeros(n_vars)\n",
572+
"\n",
573+
" # Set up harmonic balance at the target frequency\n",
574+
" # We drive at f_target, so the response at the 1st harmonic should be maximized\n",
575+
" # if the circuit is perfectly resonant at f_target.\n",
576+
" # To formulate this as a minimization, we penalize the inverse of the response.\n",
577+
" run_hb_current = setup_harmonic_balance(\n",
578+
" grps, n_vars, freq=f_target, num_harmonics=3\n",
579+
" )\n",
580+
"\n",
581+
" _, y_freq_current = run_hb_current(y_dc_current)\n",
550582
"\n",
551-
" # Normalized squared errors\n",
552-
" freq_error = ((f01 - f_target) / f_target) ** 2\n",
583+
" # Extract the 1st harmonic voltage magnitude at the junction node\n",
584+
" jj_node_idx = pmap.get(\"JJ1,p1\", 0)\n",
585+
" V_1st_harmonic = jnp.abs(y_freq_current[1, jj_node_idx])\n",
586+
"\n",
587+
" # Loss is inversely proportional to the resonance amplitude at the target frequency\n",
588+
" # We add a small epsilon to avoid division by zero, and a penalty for anharmonicity\n",
589+
" alpha = transmon_anharmonicity(Ic, Cs)\n",
553590
" alpha_error = ((alpha - alpha_target) / alpha_target) ** 2\n",
554591
"\n",
555-
" return freq_error + alpha_error\n",
592+
" resonance_loss = 1.0 / (V_1st_harmonic * 1e6 + 1e-12) # scale V to typical uV range\n",
593+
"\n",
594+
" return resonance_loss + 0.1 * alpha_error\n",
556595
"\n",
557596
"\n",
558597
"# Initial parameters in log-space\n",
559598
"params_init = jnp.array([jnp.log(Ic_init), jnp.log(Cs_init)])\n",
560599
"\n",
561600
"# Check initial loss\n",
562-
"loss_init = loss_fn(params_init)\n",
563-
"f01_init = transmon_frequency(Ic_init, Cs_init)\n",
564-
"alpha_init = transmon_anharmonicity(Ic_init, Cs_init)\n",
565-
"print(\n",
566-
" f\"Initial: f01 = {float(f01_init) / 1e9:.3f} GHz, α = {float(alpha_init) / 1e6:.1f} MHz\"\n",
567-
")\n",
568-
"print(f\"Target: f01 = {f_target / 1e9:.3f} GHz, α = {alpha_target / 1e6:.1f} MHz\")\n",
601+
"loss_init = hb_loss_fn(params_init)\n",
569602
"print(f\"Initial loss: {float(loss_init):.6f}\")"
570603
]
571604
},
@@ -589,19 +622,19 @@
589622
"outputs": [],
590623
"source": [
591624
"# Set up the optimizer\n",
592-
"optimizer = optax.adam(learning_rate=0.05)\n",
625+
"optimizer = optax.adam(learning_rate=0.02)\n",
593626
"opt_state = optimizer.init(params_init)\n",
594627
"\n",
595628
"params_opt = params_init\n",
596629
"losses = []\n",
597630
"freq_history = []\n",
598631
"\n",
599632
"# Compile the gradient function\n",
600-
"grad_fn = jax.jit(jax.grad(loss_fn))\n",
601-
"loss_jit = jax.jit(loss_fn)\n",
633+
"grad_fn = jax.jit(jax.grad(hb_loss_fn))\n",
634+
"loss_jit = jax.jit(hb_loss_fn)\n",
602635
"\n",
603636
"# Optimization loop\n",
604-
"n_steps = 200\n",
637+
"n_steps = 150\n",
605638
"for step in range(n_steps):\n",
606639
" grads = grad_fn(params_opt)\n",
607640
" updates, opt_state = optimizer.update(grads, opt_state)\n",
@@ -983,7 +1016,7 @@
9831016
" \"Cs1\": {\"component\": \"capacitor\", \"settings\": {\"C\": Cs_q1}},\n",
9841017
" \"JJ2\": {\"component\": \"josephson_junction\", \"settings\": {\"Ic\": Ic_q2}},\n",
9851018
" \"Cs2\": {\"component\": \"capacitor\", \"settings\": {\"C\": Cs_q2}},\n",
986-
" \"Cm\": {\"component\": \"coupling_cap\", \"settings\": {\"Cm\": float(Cm_val)}},\n",
1019+
" \"Cm\": {\"component\": \"coupling_cap\", \"settings\": {\"Cm\": Cm_val}},\n",
9871020
" },\n",
9881021
" \"connections\": {\n",
9891022
" \"GND,p1\": (\"Vpulse,p2\", \"JJ1,p2\", \"Cs1,p2\", \"JJ2,p2\", \"Cs2,p2\"),\n",
@@ -994,7 +1027,7 @@
9941027
" }\n",
9951028
"\n",
9961029
" grps, n_vars, pmap = compile_netlist(net, models_coupled)\n",
997-
" slvr = analyze_circuit(grps, n_vars)\n",
1030+
" slvr = analyze_circuit(grps, n_vars, backend=\"dense\")\n",
9981031
" y0 = slvr.solve_dc(grps, jnp.zeros(n_vars))\n",
9991032
"\n",
10001033
" sim_fn = setup_transient(groups=grps, linear_strategy=slvr)\n",

0 commit comments

Comments
 (0)