Solver Selection Guide

OpenQuantumSystems.jl provides several solver families for propagating the density matrix of molecular aggregates coupled to vibrational baths. This page helps you choose the right solver for your problem.

Decision Tree

Start here
  |
  |-- Do you need the full system+bath density matrix W(t)?
  |     |
  |     YES --> Use LvN (Liouville-von Neumann) on the full space
  |     |         |
  |     |         |-- Want the interaction picture?  --> LvN_SI
  |     |         |-- Want the Schroedinger picture? --> LvN_SS
  |     |
  |     NO  --> You need the reduced electronic density matrix rho(t)
  |               |
  |               |-- Is the system small enough for exact propagation
  |               |   (few molecules, few vibrational levels)?
  |               |     |
  |               |     YES --> Exact evolution with bath trace
  |               |     |         |
  |               |     |         |-- Interaction picture? --> Evolution_sI_exact
  |               |     |         |-- Schroedinger picture? --> Evolution_sS_exact
  |               |     |
  |               |     NO  --> You need a master equation
  |               |               |
  |               |               |-- Do you need high accuracy and can
  |               |               |   afford the cost of numerical integration?
  |               |               |     |
  |               |               |     YES --> QME exact (second-order, no approximations)
  |               |               |     |         |
  |               |               |     |         |-- Reduced density matrix?
  |               |               |     |         |     --> QME_sI_exact / QME_sS_exact
  |               |               |     |         |-- Full density matrix?
  |               |               |     |               --> QME_SI_exact / QME_SS_exact
  |               |               |     |
  |               |               |     NO  --> Choose an approximate method
  |               |               |               |
  |               |               |               |-- Is the system-bath coupling weak
  |               |               |               |   and Markovian-like?
  |               |               |               |     |
  |               |               |               |     YES --> QME_sI_Redfield
  |               |               |               |     |
  |               |               |               |     NO  --> Do you have a zeroth-order
  |               |               |               |             solution already?
  |               |               |               |               |
  |               |               |               |               YES --> QME_sI_iterative
  |               |               |               |               |
  |               |               |               |               NO  --> QME_sI_ansatz
  |               |               |               |                       (pick an ansatz
  |               |               |               |                        for the bath state)

Solver Families

1. Liouville-von Neumann (LvN)

Solves the standard Liouville-von Neumann equation for the full system+bath density matrix W(t):

\[\frac{d}{dt} W(t) = -i [H, W(t)]\]

These solvers propagate the entire Hilbert space without any bath tracing or perturbative expansion. They serve as the exact reference when the full space is computationally tractable.

2. Exact Evolution (non-ODE)

Direct matrix-exponential propagation using U(t) = exp(-iHt). These compute W(t) = U(t) W(0) U^dag(t) at each time point without solving a differential equation.

3. QME Exact (Second-Order Quantum Master Equation)

Solves the second-order quantum master equation without further approximations. The memory kernel is evaluated by numerical integration (QuadGK) at each time step, making these solvers accurate but expensive.

\[\frac{d}{dt} \rho^{(I)}(t) = -i [H_I^{(I)}(t), \rho(t_0)] - \int_0^t d\tau\; [H_I^{(I)}(t), [H_I^{(I)}(\tau), \rho^{(I)}(\tau)]]\]

4. QME Ansatz (Approximate Bath State)

Like QME exact, but the bath state at time tau inside the memory integral is replaced by an ansatz. This avoids propagating the full W(tau) and dramatically reduces the cost. Different ansatze offer different accuracy/cost trade-offs.

Use the unified interface:

QME_sI_ansatz(W0, tspan, agg; ansatz=:const_sch)

Available ansatze (passed as the ansatz keyword):

5. Redfield

Second-order master equation with the Redfield approximation: the density matrix inside the memory integral is evaluated at the current time t rather than at the integration variable tau, and the correlation function depends on t - tau.

QME_sI_Redfield(W0, tspan, agg)

This is appropriate when the bath correlation time is short compared to the system dynamics (weak-coupling, Markovian-like regime).

6. Iterative QME

A self-consistent iterative scheme. You first obtain a zeroth-order solution (e.g., from QME ansatz), then feed it back to refine the bath state and re-solve. This can systematically improve accuracy.

QME_sI_iterative(W0, rho_0_int_t, W_0_bath_t, tspan, agg)

Markov-level variants are available via the method keyword:

Comparison Table

Cost key: Low = single matrix exponential per step; Medium = numerical integration of memory kernel with approximate bath; High = numerical integration with exact or iterative bath; Very high = full-space delayed differential equation.

Common Usage Example

A typical workflow for a molecular dimer: set up the aggregate, choose a solver, and convert results to the Schroedinger picture if needed.

using OpenQuantumSystems

# --- Define the aggregate ---
HR = 0.01
shift = sqrt(2.0 * HR)
mode = Mode(300.0, shift)
mols = [
    Molecule([mode], 5, [12500.0, 12750.0]),
    Molecule([mode], 5, [12500.0, 12800.0]),
]
agg = Aggregate(mols)
agg.coupling[2, 3] = 100.0
agg.coupling[3, 2] = 100.0
setup_aggregate!(agg)

# --- Initial state ---
T = 300.0  # temperature in K
W0 = thermal_state_composite(T, [0.0, 0.5, 0.5], agg.operators.Ham, agg.tools)

# --- Time grid ---
t_max = 0.5
tspan = get_tspan(0.0, t_max, 500)

# --- Solve with QME ansatz (good default for production) ---
tspan_out, rho_int_t = QME_sI_ansatz(W0, tspan, agg; ansatz=:const_sch)

# --- Convert from interaction picture to Schroedinger picture ---
rho_sch_t = interaction_pic_to_schroedinger_pic(rho_int_t, agg.operators.Ham_0, tspan)

7. Förster Theory

Computes excitation energy transfer (EET) rates in the weak electronic coupling regime using the spectral overlap of donor emission and acceptor absorption lineshapes.

# Single pair
k = forster_rate(J, mol_donor, mol_acceptor; T=300.0)

# Full aggregate rate matrix
K = forster_rate_matrix(aggCore; T=300.0)

The Förster rate is: $k_{DA} = 2\pi |J_{DA}|^2 \int I_D(\omega) A_A(\omega)\, d\omega$

Use when the electronic coupling is weak relative to the reorganization energy and the bath correlation time.

8. Modified Redfield Theory

Goes beyond standard Redfield by treating diagonal (population) fluctuations non-perturbatively. Works in the exciton basis and computes time-independent population transfer rates.

# Compute rate matrix and exciton basis
rates, energies, coefficients = modified_redfield_rates(aggCore; T=300.0)

# Propagate populations
tspan = collect(0.0:0.01:10.0)
ts, pop = modified_redfield_dynamics(aggCore, [1.0, 0.0], tspan; T=300.0)

Use for intermediate coupling regimes where standard Redfield breaks down but full QME is too expensive.

Tips

• Start with Evolution_sI_exact or LvN_sI to establish a reference solution for your system size, then switch to an approximate solver for production.
• Interaction picture (suffix _sI or _SI) is generally more numerically stable for oscillatory dynamics because the fast free-evolution phase is factored out.
• Tolerance tuning: The int_reltol and int_abstol parameters control the accuracy of the memory-kernel quadrature. Tighten them for benchmarks, loosen for speed.
• Ansatz choice: Start with :const_sch (cheapest nontrivial ansatz). If accuracy is insufficient, move to :upart1_sch or :upart2_sch, or switch to the iterative solver.