QuTiP BRME DSL Example

QuTiP BRME DSL Example#

This example (example_qutip_brme_dsl.py) demonstrates a more specialized workflow (BRME-style notation) and shows how the DSL maps physics-style LaTeX into operator references consumed by the QuTiP backend.

What you’ll learn#

  • How to express BRME-style interaction terms in LaTeX.

  • How the DSL canonicalizes and maps those operators to QuTiP objects.

Source#

  1# flake8: noqa
  2r"""
  3Reimplementation of the QuTiP ``brmesolve`` tutorial using the LaTeX DSL.
  4
  5Original notebook steps (``qutip_example.ipynb``):
  61) Constant Hamiltonian H = a^\dagger a for a truncated oscillator.
  72) Add a time-dependent drive H(t) = H + sin(t)(a + a^\dagger).
  83) Add simple dissipation with a decaying envelope.
  94) Add a composite, explicitly time-dependent collapse operator.
 10
 11This example shows how each Hamiltonian/collapse operator can be expressed
 12in LaTeX, compiled with the parser, and then passed directly to QuTiP
 13``brmesolve`` while keeping the rest of the workflow unchanged.
 14"""
 15
 16from __future__ import annotations
 17
 18import sys
 19from pathlib import Path
 20
 21import numpy as np
 22import matplotlib.pyplot as plt
 23
 24from latex_parser.latex_api import compile_model, make_config
 25from qutip import about, basis, mesolve, plot_expectation_values
 26
 27
 28def _number_op(cutoff: int):
 29    """Use the DSL to build the bosonic number operator."""
 30    cfg = make_config(bosons=[cutoff])
 31    model = compile_model(
 32        H_latex=r"n_{1}", params={}, config=cfg, c_ops_latex=None, backend="qutip"
 33    )
 34    return model.H
 35
 36
 37def compile_constant_oscillator(cutoff: int):
 38    cfg = make_config(bosons=[cutoff])
 39    H_latex = r"a_{1}^{\dagger} a_{1}"
 40    model = compile_model(
 41        H_latex=H_latex, params={}, config=cfg, c_ops_latex=None, backend="qutip"
 42    )
 43    return model.H
 44
 45
 46def compile_driven_oscillator(cutoff: int):
 47    cfg = make_config(bosons=[cutoff])
 48    H_latex = r"a_{1}^{\dagger} a_{1} + \sin(t) \left(a_{1} + a_{1}^{\dagger}\right)"
 49    model = compile_model(
 50        H_latex=H_latex, params={}, config=cfg, c_ops_latex=None, backend="qutip"
 51    )
 52    return model.H
 53
 54
 55def compile_decay_collapse_ops(cutoff: int, kappa: float):
 56    cfg = make_config(bosons=[cutoff])
 57    model = compile_model(
 58        H_latex=r"a_{1}^{\dagger} a_{1}",
 59        params={"kappa": kappa},
 60        config=cfg,
 61        c_ops_latex=[
 62            r"\sqrt{\kappa} \exp(-t/2) \, a_{1}",
 63            r"\sqrt{\kappa} \exp(-t/2) \, a_{1}^{\dagger}",
 64        ],
 65        backend="qutip",
 66    )
 67    return model.H, model.c_ops, model.args
 68
 69
 70def compile_composite_collapse_ops(cutoff: int, kappa: float):
 71    cfg = make_config(bosons=[cutoff])
 72    model = compile_model(
 73        H_latex=r"a_{1}^{\dagger} a_{1}",
 74        params={"kappa": kappa},
 75        config=cfg,
 76        c_ops_latex=[
 77            r"\sqrt{\kappa} \, a_{1} \exp(I t)",
 78            r"\sqrt{\kappa} \, a_{1}^{\dagger} \exp(-I t)",
 79        ],
 80        backend="qutip",
 81    )
 82    return model.H, model.c_ops, model.args
 83
 84
 85def main() -> None:
 86    N = 2
 87    psi0 = basis(N, N - 1)
 88    times = np.linspace(0, 10, 100)
 89    n_op = _number_op(N)
 90
 91    # 1) Constant Hamiltonian
 92    H_const = compile_constant_oscillator(N)
 93    result_const = mesolve(H_const, psi0, times, e_ops=[n_op])
 94    plot_expectation_values(result_const, ylabels=["<n>"])
 95
 96    # 2) Driven Hamiltonian with sin(t)
 97    H_drive = compile_driven_oscillator(N)
 98    result_drive = mesolve(H_drive, psi0, times, e_ops=[n_op])
 99    plot_expectation_values(result_drive, ylabels=["<n>"])
100
101    # 3) Dissipation with decaying envelope (Bloch-Redfield a_ops form)
102    kappa = 0.2
103    H_decay, c_ops_decay, args_decay = compile_decay_collapse_ops(N, kappa)
104    result_decay = mesolve(
105        H_decay, psi0, times, c_ops=c_ops_decay, e_ops=[n_op], args=args_decay
106    )
107    plot_expectation_values(result_decay, ylabels=["<n>"])
108
109    # 4) Composite time-dependent collapse operator
110    H_comp, c_ops_comp, args_comp = compile_composite_collapse_ops(N, kappa)
111    result_comp = mesolve(
112        H_comp, psi0, times, c_ops=c_ops_comp, e_ops=[n_op], args=args_comp
113    )
114    plot_expectation_values(result_comp, ylabels=["<n>"])
115
116    # Parity checks mirroring the notebook assertions
117    assert np.allclose(result_const.expect[0], 1.0)
118    assert np.all(np.diff(result_comp.expect[0]) <= 0.0)
119    about()
120    plt.show()
121
122
123if __name__ == "__main__":
124    main()

Run#

python examples/example_qutip_brme_dsl.py

Notes#

  • This example is targeted at users integrating domain-specific LaTeX conventions (BRME) with the parser. Inspect IR output to verify mappings.

Contents