State Metrics and Solver Diagnostics

This tutorial shows how to compute common open-system quantities such as purity, entropy, fidelity, trace distance, populations, coherence, and Bloch-vector components. The same metric callbacks can be passed directly to solvers through state_observables; mcsolve aggregates the built-in linearizable metrics in the backend.

import matplotlib.pyplot as plt
import numpy as np

import openquantumsim as oqs

The metric functions accept kets or density matrices whenever the quantity is well defined.

up = np.array([1.0, 0.0], dtype=np.complex128)
down = np.array([0.0, 1.0], dtype=np.complex128)
plus = (up + down) / np.sqrt(2.0)
mixed = 0.5 * np.eye(2, dtype=np.complex128)

summary = {
    "purity(mixed)": oqs.purity(mixed),
    "entropy(mixed)": oqs.von_neumann_entropy(mixed),
    "fidelity(up, plus)": oqs.fidelity(up, plus),
    "trace_distance(up, down)": oqs.trace_distance(up, down),
    "bloch(plus)": oqs.bloch_vector(plus),
}
summary

state_metrics(...) builds a dictionary of scalar callbacks. Deterministic solvers and single trajectories evaluate callbacks at each requested time without returning full states unless Options(save_states=True) is also set. mcsolve can aggregate built-in metrics such as purity, entropy, pure-state fidelity, populations, and Bloch components without saving every trajectory state.

mesolve uses the Julia backend. Run python setup_julia.py once before executing this notebook locally.

gamma = 0.45
atom = oqs.SpinSpace(0.5, label="atom")

H = 0.0 * oqs.sigmaz(atom)
collapse = np.sqrt(gamma) * oqs.sigmam(atom)
excited = oqs.basis(atom, "up")
ground = oqs.basis(atom, "down")
rho0 = oqs.ket2dm(excited)
times = np.linspace(0.0, 8.0, 121)

metrics = oqs.state_metrics(
    purity=True,
    entropy=True,
    linear_entropy=True,
    participation_ratio=True,
    fidelity_to=excited,
    trace_distance_to=ground,
    population_indices=[0, 1],
    bloch_vector=True,
)

result = oqs.mesolve(
    H,
    rho0,
    times,
    c_ops=[collapse],
    state_observables=metrics,
    options=oqs.Options(rtol=1e-9, atol=1e-11, save_states=False),
)

sorted(result.state_observables)

obs = result.state_observables

fig, axes = plt.subplots(2, 1, figsize=(7, 6), sharex=True)
axes[0].plot(times, obs["population_0"].real, label="excited population")
axes[0].plot(times, obs["population_1"].real, label="ground population")
axes[0].plot(times, obs["fidelity"].real, "--", label="fidelity to initial")
axes[0].set_ylabel("probability")
axes[0].legend()

axes[1].plot(times, obs["purity"].real, label="purity")
axes[1].plot(times, obs["entropy"].real, label="entropy")
axes[1].plot(times, obs["linear_entropy"].real, label="linear entropy")
axes[1].set_xlabel("time")
axes[1].legend()
fig.tight_layout()

The same callback dictionary works for a single Monte Carlo wave-function trajectory. A trajectory stays pure, but the populations and fidelity change discontinuously when a jump occurs.

trajectory_metrics = oqs.state_metrics(
    purity=True,
    fidelity_to=excited,
    population_indices=[0, 1],
)

trajectory = oqs.single_trajectory(
    H,
    excited,
    times,
    c_ops=[collapse],
    state_observables=trajectory_metrics,
    options=oqs.Options(seed=2026, max_step=0.02, save_states=False),
)

traj_obs = trajectory.state_observables

fig, ax = plt.subplots(figsize=(7, 3.5))
ax.step(times, traj_obs["population_0"].real, where="post", label="excited")
ax.step(times, traj_obs["population_1"].real, where="post", label="ground")
ax.plot(times, traj_obs["purity"].real, "--", label="purity")
ax.set_xlabel("time")
ax.set_ylabel("trajectory metric")
ax.legend()

You can mix built-in callbacks with custom scalar diagnostics. This is useful for problem-specific return probabilities, subsystem entropies, or distances to reference states.

custom_metrics = oqs.state_metrics(purity=True, fidelity_to=excited)
custom_metrics["ground_infidelity"] = lambda rho: oqs.infidelity(rho, ground)

custom_result = oqs.mesolve(
    H,
    rho0,
    times,
    c_ops=[collapse],
    state_observables=custom_metrics,
    options=oqs.Options(rtol=1e-9, atol=1e-11),
)

custom_result.state_observables["ground_infidelity"].real[-1]