SCQubits Example

SCQubits Example#

This example demonstrates interoperability with the scqubits package (if available) and shows how to compose models that include non-standard qubit implementations.

What you’ll learn#

  • How to adapt latex_parser outputs to other quantum libraries.

  • How to inspect and adapt operator shapes when integrating with third-party packages.

Source#

  1#!/usr/bin/env python
  2# flake8: noqa
  3"""
  4Examples: integrating scqubits transmon/fluxonium with the LaTeX DSL.
  5
  6These examples show how to:
  71) Build scqubits devices.
  82) Register their operators as CustomSpec entries.
  93) Use operator-valued functions like cos(phi) directly in the DSL.
 104) Verify compiled Hamiltonians against scqubits-provided matrices.
 11
 12No long simulations are run by default; the examples stop after compiling
 13to QuTiP objects and, optionally, validating against scqubits to keep
 14runtime light. Run this file directly to execute both examples if
 15scqubits is installed.
 16"""
 17
 18from __future__ import annotations
 19
 20import sys
 21from pathlib import Path
 22import logging
 23from typing import Tuple
 24
 25import numpy as np
 26from qutip import Qobj, destroy, qeye, tensor  # type: ignore
 27from scipy.linalg import cosm, sinm  # type: ignore
 28
 29ROOT = Path(__file__).resolve().parents[1]
 30if str(ROOT) not in sys.path:
 31    sys.path.insert(0, str(ROOT))
 32
 33from latex_parser.dsl import CustomSpec, HilbertConfig
 34from latex_parser.latex_api import compile_model as compile_latex_model, make_config
 35
 36logger = logging.getLogger(__name__)
 37
 38
 39def _require_scqubits():
 40    try:
 41        import scqubits as scq  # type: ignore
 42    except ImportError as exc:  # pragma: no cover - optional dependency
 43        raise RuntimeError(
 44            "scqubits is required for these examples. Install it with "
 45            "'pip install scqubits'."
 46        ) from exc
 47    return scq
 48
 49
 50def fluxonium_model(
 51    EJ: float = 8.0,
 52    EC: float = 1.0,
 53    EL: float = 0.5,
 54    cutoff: int = 35,
 55    n_g: float = 0.0,
 56    flux: float = 0.0,
 57) -> Tuple[object, HilbertConfig]:
 58    """
 59    Build a fluxonium Hamiltonian
 60
 61        H = 4 EC (n - n_g)^2 + 0.5 EL phi^2 - EJ cos(phi)
 62
 63    using scqubits operators mapped into a CustomSpec. Returns the compiled
 64    model and the HilbertConfig.
 65    """
 66    scq = _require_scqubits()
 67
 68    flux_dev = scq.Fluxonium(EJ=EJ, EC=EC, EL=EL, cutoff=cutoff, flux=flux)
 69    phi_op = flux_dev.phi_operator()
 70    n_op = flux_dev.n_operator()
 71
 72    # scqubits >=3 returns NumPy arrays; wrap into Qobj for the DSL.
 73    if not hasattr(phi_op, "dims"):
 74        phi_op = Qobj(phi_op, dims=[[cutoff], [cutoff]])
 75    if not hasattr(n_op, "dims"):
 76        n_op = Qobj(n_op, dims=[[cutoff], [cutoff]])
 77
 78    # Register custom subsystem with phi and n operators.
 79    spec = CustomSpec(
 80        label="c",
 81        index=1,
 82        dim=phi_op.dims[0][0],
 83        operators={"phi": phi_op, "nphi": n_op},
 84    )
 85    cfg = HilbertConfig(qubits=[], bosons=[], customs=[spec])
 86
 87    H_latex = r"4 EC (nphi_{1} - n_g)^2 + \frac{1}{2} EL \phi_{1}^2 - EJ \cos(\phi_{1})"
 88    params = {"EC": EC, "EL": EL, "EJ": EJ, "n_g": n_g}
 89
 90    model = compile_latex_model(
 91        H_latex=H_latex,
 92        params=params,
 93        config=cfg,
 94        t_name="t",
 95    )
 96    return model, cfg
 97
 98
 99def fluxonium_model_with_flux_offset(
100    EJ: float = 8.0,
101    EC: float = 1.0,
102    EL: float = 0.5,
103    cutoff: int = 35,
104    n_g: float = 0.0,
105    phi_ext: float = 0.2,
106) -> Tuple[object, HilbertConfig]:
107    """
108    Fluxonium with external flux offset via the shifted cosine
109
110        -EJ cos(phi - phi_ext) = -EJ [cos(phi) cos(phi_ext) + sin(phi) sin(phi_ext)].
111
112    Uses operator-valued cos/sin on the phi operator and scalar parameters
113    for phi_ext.
114    """
115    scq = _require_scqubits()
116
117    flux_dev = scq.Fluxonium(EJ=EJ, EC=EC, EL=EL, cutoff=cutoff, flux=phi_ext)
118    phi_op = flux_dev.phi_operator()
119    n_op = flux_dev.n_operator()
120
121    if not hasattr(phi_op, "dims"):
122        phi_op = Qobj(phi_op, dims=[[cutoff], [cutoff]])
123    if not hasattr(n_op, "dims"):
124        n_op = Qobj(n_op, dims=[[cutoff], [cutoff]])
125
126    spec = CustomSpec(
127        label="c",
128        index=1,
129        dim=phi_op.dims[0][0],
130        operators={"phi": phi_op, "nphi": n_op},
131    )
132    cfg = HilbertConfig(qubits=[], bosons=[], customs=[spec])
133
134    H_latex = r"""
135        4 EC (nphi_{1} - n_g)^2
136        + \frac{1}{2} EL \phi_{1}^2
137        - EJ \left(
138            \cos(\phi_{1}) \cos(\phi_{ex})
139            + \sin(\phi_{1}) \sin(\phi_{ex})
140        \right)
141    """
142    params = {
143        "EC": EC,
144        "EL": EL,
145        "EJ": EJ,
146        "n_g": n_g,
147        "phi_ex": phi_ext,
148    }
149
150    model = compile_latex_model(
151        H_latex=H_latex,
152        params=params,
153        config=cfg,
154        t_name="t",
155    )
156    return model, cfg
157
158
159def verify_fluxonium_against_scqubits(
160    EJ: float = 8.0,
161    EC: float = 1.0,
162    EL: float = 0.5,
163    cutoff: int = 15,
164    n_g: float = 0.0,
165    flux: float = 0.0,
166) -> float:
167    """
168    Compile the fluxonium Hamiltonian from LaTeX and compare to an explicit
169    matrix built from the same scqubits ``phi`` and ``n`` operators. Returns
170    the operator norm of the difference (expected to be ~0).
171    """
172    scq = _require_scqubits()
173    model, cfg = fluxonium_model(EJ=EJ, EC=EC, EL=EL, cutoff=cutoff, n_g=n_g, flux=flux)
174    H_dsl = model.H if isinstance(model.H, Qobj) else model.H[0]
175    if not isinstance(H_dsl, Qobj) and hasattr(H_dsl, "shape"):
176        dim = H_dsl.shape[0]
177        H_dsl = Qobj(H_dsl, dims=[[dim], [dim]])
178
179    flux_dev = scq.Fluxonium(EJ=EJ, EC=EC, EL=EL, cutoff=cutoff, flux=flux)
180    phi = flux_dev.phi_operator()
181    n_op = flux_dev.n_operator()
182    if not hasattr(phi, "dims"):
183        phi = Qobj(phi, dims=[[cutoff], [cutoff]])
184    if not hasattr(n_op, "dims"):
185        n_op = Qobj(n_op, dims=[[cutoff], [cutoff]])
186
187    cos_phi = Qobj(cosm(phi.full()), dims=phi.dims)
188    H_ref = 4 * EC * (n_op - n_g) ** 2 + 0.5 * EL * (phi**2) - EJ * cos_phi
189
190    H_sq = flux_dev.hamiltonian()
191    if not isinstance(H_sq, Qobj) and hasattr(H_sq, "shape"):
192        dim = H_sq.shape[0]
193        H_sq = Qobj(H_sq, dims=[[dim], [dim]])
194
195    diff = (H_dsl - H_ref).norm()
196    diff_sq = (H_ref - H_sq).norm()
197    if diff_sq > 0:  # pragma: no cover - informational
198        logger.info(
199            "scqubits Hamiltonian differs from reference by %.3e (truncation effects?)",
200            diff_sq,
201        )
202    return diff
203
204
205def verify_fluxonium_with_flux_offset_against_scqubits(
206    EJ: float = 8.0,
207    EC: float = 1.0,
208    EL: float = 0.5,
209    cutoff: int = 15,
210    n_g: float = 0.0,
211    phi_ext: float = 0.2,
212) -> float:
213    """
214    Verify flux-offset fluxonium by comparing the DSL-compiled Hamiltonian
215    to an explicit reference built from scqubits ``phi``/``n`` operators and
216    matrix cosine/sine. Returns the operator norm of the difference
217    (expected ~0).
218    """
219    scq = _require_scqubits()
220    model, cfg = fluxonium_model_with_flux_offset(
221        EJ=EJ, EC=EC, EL=EL, cutoff=cutoff, n_g=n_g, phi_ext=phi_ext
222    )
223    H_dsl = model.H if isinstance(model.H, Qobj) else model.H[0]
224    if not isinstance(H_dsl, Qobj) and hasattr(H_dsl, "shape"):
225        dim = H_dsl.shape[0]
226        H_dsl = Qobj(H_dsl, dims=[[dim], [dim]])
227
228    flux_dev = scq.Fluxonium(EJ=EJ, EC=EC, EL=EL, cutoff=cutoff, flux=phi_ext)
229    phi = flux_dev.phi_operator()
230    n_op = flux_dev.n_operator()
231    if not hasattr(phi, "dims"):
232        phi = Qobj(phi, dims=[[cutoff], [cutoff]])
233    if not hasattr(n_op, "dims"):
234        n_op = Qobj(n_op, dims=[[cutoff], [cutoff]])
235
236    cos_phi = Qobj(cosm(phi.full()), dims=phi.dims)
237    sin_phi = Qobj(sinm(phi.full()), dims=phi.dims)
238    H_ref = (
239        4 * EC * (n_op - n_g) ** 2
240        + 0.5 * EL * (phi**2)
241        - EJ * (cos_phi * np.cos(phi_ext) + sin_phi * np.sin(phi_ext))
242    )
243
244    H_sq = flux_dev.hamiltonian()
245    if not isinstance(H_sq, Qobj) and hasattr(H_sq, "shape"):
246        dim = H_sq.shape[0]
247        H_sq = Qobj(H_sq, dims=[[dim], [dim]])
248
249    diff = (H_dsl - H_ref).norm()
250    diff_sq = (H_ref - H_sq).norm()
251    if diff_sq > 0:  # pragma: no cover - informational
252        logger.info(
253            "scqubits Hamiltonian with flux offset differs from reference by %.3e",
254            diff_sq,
255        )
256    return diff
257
258
259def transmon_cavity_model(
260    EJ: float = 20.0,
261    EC: float = 0.25,
262    cutoff: int = 10,
263    cav_cutoff: int = 6,
264    g: float = 0.05,
265) -> Tuple[object, HilbertConfig]:
266    r"""
267    Build a simple transmon + cavity charge-coupled Hamiltonian:
268
269        H = H_transmon + omega_c n_cav + g n (a + a^\dagger)
270
271    where H_transmon = 4 EC n^2 - EJ cos(phi). The transmon operators are
272    built internally using harmonic-oscillator zero-point fluctuations
273    (no scqubits operators are used in the compiled model); the cavity
274    is a boson in the DSL.
275    """
276    H_tmon, phi_op, n_op = _build_transmon_ops_internal(EJ=EJ, EC=EC, cutoff=cutoff)
277
278    spec = CustomSpec(
279        label="c",
280        index=1,
281        dim=H_tmon.dims[0][0],
282        operators={"Htm": H_tmon, "phi": phi_op, "ntm": n_op},
283    )
284    # One cavity boson a_1 with cutoff cav_cutoff
285    cfg = make_config(qubits=0, bosons=[(cav_cutoff, "a")], customs=[spec])
286
287    H_latex = r"Htm_{1} + \omega_c n_{1} + g \, ntm_{1} (a_{1} + a_{1}^{\dagger})"
288    params = {"omega_c": 5.0, "g": g}
289
290    model = compile_latex_model(
291        H_latex=H_latex,
292        params=params,
293        config=cfg,
294        t_name="t",
295    )
296    return model, cfg
297
298
299def verify_transmon_cavity_against_scqubits(
300    EJ: float = 20.0,
301    EC: float = 0.25,
302    cutoff: int = 8,
303    cav_cutoff: int = 4,
304    g: float = 0.05,
305    omega_c: float = 5.0,
306) -> float:
307    """
308    Compile the transmon-cavity Hamiltonian from LaTeX and compare to an
309    explicit construction built from the same internal transmon operators
310    used in :func:`transmon_cavity_model`. Returns the operator norm of the
311    difference (expected ~0). Any deviation between this reference and
312    scqubits' transmon Hamiltonian is logged for information.
313    """
314    scq = _require_scqubits()
315    model, cfg = transmon_cavity_model(
316        EJ=EJ, EC=EC, cutoff=cutoff, cav_cutoff=cav_cutoff, g=g
317    )
318    H_dsl: Qobj = model.H if isinstance(model.H, Qobj) else model.H[0]
319    H_tm, _, n_tm = _build_transmon_ops_internal(EJ=EJ, EC=EC, cutoff=cutoff)
320
321    # Cavity operators (boson is first subsystem in cfg ordering, custom next)
322    a = destroy(cav_cutoff)
323    n_cav = a.dag() * a
324
325    # Tensor order: [boson, custom]
326    I_tm = qeye(H_tm.dims[0])
327    I_cav = qeye([cav_cutoff])
328
329    H_cav = omega_c * tensor(n_cav, I_tm)
330    H_tm_full = tensor(I_cav, H_tm)
331    n_tm_full = tensor(I_cav, n_tm)
332    a_full = tensor(a, I_tm)
333    adag_full = tensor(a.dag(), I_tm)
334
335    H_expected = H_tm_full + H_cav + g * n_tm_full * (a_full + adag_full)
336    diff = (H_dsl - H_expected).norm()
337
338    # Informational check against scqubits' transmon implementation
339    try:  # pragma: no cover - optional dependency precision differences
340        tmon = scq.Transmon(EJ=EJ, EC=EC, ng=0.0, ncut=cutoff, truncated_dim=cutoff)
341        H_tm_scq = tmon.hamiltonian()
342        n_tm_scq = tmon.n_operator()
343        if not hasattr(H_tm_scq, "dims"):
344            dim_tm = H_tm_scq.shape[0]
345            H_tm_scq = Qobj(H_tm_scq, dims=[[dim_tm], [dim_tm]])
346        if not hasattr(n_tm_scq, "dims"):
347            dim_tm = n_tm_scq.shape[0]
348            n_tm_scq = Qobj(n_tm_scq, dims=[[dim_tm], [dim_tm]])
349
350        H_tm_full_scq = tensor(I_cav, H_tm_scq)
351        n_tm_full_scq = tensor(I_cav, n_tm_scq)
352        H_expected_scq = (
353            H_tm_full_scq + H_cav + g * n_tm_full_scq * (a_full + adag_full)
354        )
355        diff_scq = (H_expected - H_expected_scq).norm()
356        if diff_scq > 0:
357            logger.info(
358                "Transmon reference differs from scqubits by %.3e "
359                "(basis/truncation effects).",
360                diff_scq,
361            )
362    except Exception:
363        pass
364
365    return float(diff)
366
367
368def _build_transmon_ops_internal(
369    EJ: float, EC: float, cutoff: int
370) -> Tuple[Qobj, Qobj, Qobj]:
371    """
372    Build transmon phi, n, and Hamiltonian using harmonic-oscillator
373    zero-point fluctuations (Koch et al.).
374    """
375    dim = 2 * cutoff + 1  # match scqubits Transmon dimension (2*ncut + 1)
376    a = destroy(dim)
377    adag = a.dag()
378    # Zero-point fluctuations
379    phi_zpf = (2 * EC / EJ) ** 0.25
380    n_zpf = (EJ / (32 * EC)) ** 0.25
381
382    phi_op = phi_zpf * (a + adag)
383    n_op = 1j * n_zpf * (adag - a)
384
385    # cos(phi) via matrix cosine
386    cos_phi = Qobj(cosm(phi_op.full()), dims=phi_op.dims)
387    H_tm = 4 * EC * (n_op * n_op) - EJ * cos_phi
388    return H_tm, phi_op, n_op
389
390
391if __name__ == "__main__":  # pragma: no cover - example execution
392    try:
393        flux_model, flux_cfg = fluxonium_model()
394        print("Fluxonium model compiled. time_dependent:", flux_model.time_dependent)
395        diff_flux = verify_fluxonium_against_scqubits()
396        print(f"Fluxonium H difference norm (vs reference): {diff_flux:.3e}")
397
398        flux_model_off, _ = fluxonium_model_with_flux_offset()
399        diff_flux_off = verify_fluxonium_with_flux_offset_against_scqubits()
400        print(
401            "Fluxonium (flux offset) H difference norm (vs reference): "
402            f"{diff_flux_off:.3e}"
403        )
404    except RuntimeError as exc:
405        print(exc)
406
407    try:
408        tmon_model, tmon_cfg = transmon_cavity_model()
409        print(
410            "Transmon-cavity model compiled. time_dependent:",
411            tmon_model.time_dependent,
412        )
413        diff_tmon = verify_transmon_cavity_against_scqubits()
414        print(
415            f"Transmon-cavity H difference norm: {diff_tmon:.3e}",
416        )
417    except RuntimeError as exc:
418        print(exc)

Run#

python examples/example_scqubits.py

Notes#

  • This example may require scqubits to be installed; it demonstrates integration patterns rather than a turnkey workflow.

Contents