API Reference

Mode{T}(omega, shift)

Stuct which purpose is to model vibrational LHO mode in Molecule.

$V(q) = \hbar \omega (q - q_0)^2, \quad \hbar = 1$

Arguments

• omega: The frequency of LHO ($\omega$).
• shift: The shift of the coordinate of LHO ($q_0$).

Molecule{T,C1,C2}(modes, Nvib, E)

Stuct which purpose is to model a molecule in Aggregate.

Arguments

• modes: Vector of modes (Mode).
• Nvib: Maximum number of vibrational states for all modes.
• E: Energy of ground and excited state of molecule (HOMO, LUMO).

electronic_indices(molCount)

Get the electric indices for all states on Aggregate.

Arguments

• molCount: Number of molecules in Aggregate.

franck_condon_factors(size, shift)

Get Franck-Condon factors for LHO mode calculated using ShiftOperator.

get_mol_frequencies(mol)

Get shifts of every mode on the Molecule state.

Arguments

• mol: Instance of Molecule.

get_mol_shifts(mol)

Get shifts of every mode on the Molecule state.

Arguments

• mol: Instance of Molecule.

get_mol_state_energy(mol, molElState, molVibState)

Get the energy of the Molecule state.

Arguments

• mol: Instance of Molecule.
• molElState: Electric state of the molecule in local basis (i.e. 1 or 2, where 1 is ground state).
• molVibState: Vibrational state of the molecule (e.g. [1, 5, 2, 2]).

get_mol_state_fc(mol, molElState1, molVibState1, molElState2, molVibState2)

Get the energy of the Molecule state.

Arguments

• mol: Instance of Molecule.
• molElState1, molElState2: Electric state of the molecule in local basis (i.e. 1 or 2, where 1 is ground state).
• molVibState1, molVibState2: Vibrational state of the molecule (e.g. [1, 5, 2, 2]).

vibrational_indices(maxInds)

Get the vibrational indices for all states on Molecule.

Arguments

• maxInds: Vector of maximum number of vibrational states.

AggregateCore{T,C1,C2}(molecules, coupling)
AggregateCore{T,C1,C2}(molecules)

Mutable stuct which purpose is to store important information about the whole aggregate.

Arguments

• molecules: Vector of molecules (Molecule).
• coupling: Matrix of couplings $J_{nm}$ between molecules.

get_agg_state_energy(agg, aggElState, aggVibState)

Get Hamiltonian of the Aggregate in a form of DenseOperator.

Arguments

• agg: Instance of Aggregate.
• aggElState: Aggregate electric state (e.g. [1, 1, 2]).
• aggVibState: Aggregate vibrational state (e.g. [[9, 1], [2, 5, 3], [10]]).

get_frequencies(aggCore)

Get frequencies for every mode on each molecule in the Aggregate.

get_nvib(aggCore)

Get maximum number of vibrational states of each molecule in the Aggregate.

get_shifts(aggCore)

Get shifts for every mode on each molecule in the Aggregate.

vibrational_indices(aggCore)

Get all electronic indices of the Aggregate.

Arguments

• aggCore: Instance of Aggregate.

getFCProd(agg, FCFact, aggIndices, vibindices)

Get product of Franck-Condon factors. This way the trace over bath will be faster

get_franck_condon_factors(aggCore, aggIndices)
get_franck_condon_factors(aggCore)

Get Frack-Condon factors of the Aggregate in a form of matrix.

Arguments

• aggCore: Instance of Aggregate.
• aggIndices: Aggregate indices generated by get_indices.

get_indices(aggCore)

Get all indices (electronic and vibrational) of the Aggregate.

Arguments

• aggCore: Instance of Aggregate.

getindicesmap(aggCore)

Get pointers (integers) to the indices of the Aggregate separated by electronic states (e.g. [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]]).

Arguments

• agg: Instance of Aggregate.

take_el_part(A, a, b, indicesMap)

Take electric part specified by electric indices a and b from the A (type of Array).

vibrational_indices(aggCore)

Get all vibrational indices of the Aggregate.

getagghaminteraction(agg, aggTools.indices, franckcondonfactors) getagghaminteraction(agg, aggTools.indices) getaggham_interaction(agg)

Get interation Hamiltonian of the Aggregate.

Arguments

• agg: Instance of Aggregate.
• aggTools.indices: Aggregate indices generated by get_indices.
• franck_condon_factors: Franck-Condon factors generated by get_franck_condon_factors.

getagghambathbig(agg)

Get Hamiltonian of the bath, $H_B$ only for $\mathcal{H}_S \otimes \mathcal{H}_B$.

getagghambathsmall(agg)

Get Hamiltonian of the bath, $H_B$ only for $\mathcal{H}_B$.

getagghamsystembath(agg, aggTools.indices, franckcondonfactors; groundEnergy = false) getagghamsystembath(agg, aggTools.indices; groundEnergy = false) getagghamsystembath(agg; groundEnergy = false)

Get system and bath Hamiltonian of the Aggregate in a form of DenseOperator.

Arguments

• agg: Instance of Aggregate.
• aggTools.indices: Aggregate indices generated by get_indices.
• franck_condon_factors: Franck-Condon factors generated by get_franck_condon_factors.

getagghamsystembig(agg)

Get Hamiltonian of the system, $H_S$ for $\mathcal{H}_S \otimes \mathcal{H}_B$.

Arguments

• agg: Instance of Aggregate.

getagghamsystemsmall(agg)

Get Hamiltonian of the system, $H_S$ only for $\mathcal{H}_S$.

Arguments

• agg: Instance of Aggregate.

getagghamiltonian(agg, aggTools.indices, franckcondonfactors; , groundEnergy = true) getagghamiltonian(agg, aggTools.indices; groundEnergy = true) getagghamiltonian(agg; groundEnergy = true)

Get Hamiltonian of the Aggregate.

Arguments

• agg: Instance of Aggregate.
• aggTools.indices: Aggregate indices generated by get_indices.
• franck_condon_factors: Franck-Condon factors generated by get_franck_condon_factors.

setup_aggregate(agg; groundEnergy=true, verbose=false)

Generate all basic data from the Aggregate. Returns aggInds, vibindices, aggIndLen, basis, FCFact, FCProd, Ham, Ham_0, Ham_I.

setup_dimer(; E1=12500., E2=12700., J=100., mode_omega=200., mode_hr=0.02, nvib=3)

Create a fully set-up dimer Aggregate with one vibrational mode per molecule.

Returns an Aggregate ready for time-evolution calculations.

setup_linear_chain(; energies=[12500., 12700.], J_nearest=100., mode_omega=200., mode_hr=0.02, nvib=3)

Create a fully set-up linear chain Aggregate with nearest-neighbour coupling and one vibrational mode per molecule.

Returns an Aggregate ready for time-evolution calculations.

setup_trimer(; E=[12500., 12700., 12900.], J_matrix=zeros(3,3), mode_omega=200., mode_hr=0.02, nvib=3)

Create a fully set-up trimer Aggregate with one vibrational mode per molecule.

J_matrix is a 3x3 coupling matrix (only off-diagonal elements are used). Returns an Aggregate ready for time-evolution calculations.

evolutionApproximate!(op_array, op0, tspan, Hamiltonian)

Calculate approximate time evolution of the op0 inplace based on $U(t)$ that is calculated for $t_0$ and $t_\text{step}$. See evolution_operator. This method returns Vector of Arrays.

evolutionApproximate!(ket_array, ket0, tspan, Hamiltonian)

Calculate approximate time evolution of the ket0 inplace based on $U(t)$ that is calculated for $t_0$ and $t_\text{step}$. See evolution_operator. Argument ket_array is Vector of Arrays.

evolutionApproximate(ket0, tspan, Hamiltonian)

Calculate approximate time evolution of the ket0 based on $U(t)$ that is calculated for $t_0$ and $t_\text{step}$. See evolution_operator. This method returns Vector of Ket states.

evolutionApproximate(op0, tspan, Hamiltonian)

Calculate approximate time evolution of the op0 based on $U(t)$ that is calculated for $t_0$ and $t_\text{step}$. See evolution_operator. This method returns Vector of Operators.

evolutionExact!(op_array, op0, tspan, Hamiltonian; 
	diagonalize = true, approximate = false)

Calculate exact time evolution of the op0 state inplace see evolution_operator_iterator. The diagonalize argument decompose Hamiltonian into $\lambda_i, S, S^{-1}$ and calculate the exponential using eigenvalue decomposition. The approximate option asusmes that tspan is made of equidistant points, therefore $U(t)$ has to be calculated for $t_0$ and $t_\text{step}$.

evolutionExact!(ket_array, ket0, tspan, Hamiltonian; 
	diagonalize = true, approximate = false)

Calculate exact time evolution of the ket0 state inplace see evolution_operator_iterator. The diagonalize argument decompose Hamiltonian into $\lambda_i, S, S^{-1}$ and calculate the exponential using eigenvalue decomposition. The approximate option asusmes that tspan is made of equidistant points, therefore $U(t)$ has to be calculated for $t_0$ and $t_\text{step}$.

evolutionExact(ket0, tspan, Hamiltonian; diagonalize = true, approximate = false)

Calculate exact time evolution of the ket0 state see evolution_operator_iterator. The diagonalize argument decompose Hamiltonian into $\lambda_i, S, S^{-1}$ and calculate the exponential using eigenvalue decomposition. The approximate option asusmes that tspan is made of equidistant points, therefore $U(t)$ has to be calculated for $t_0$ and $t_\text{step}$.

evolutionExact(op0, tspan, Hamiltonian; diagonalize = true, approximate = false)

Calculate exact time evolution of the op0 state see evolution_operator_iterator. The diagonalize argument decompose Hamiltonian into $\lambda_i, S, S^{-1}$ and calculate the exponential using eigenvalue decomposition. The approximate option asusmes that tspan is made of equidistant points, therefore $U(t)$ has to be calculated for $t_0$ and $t_\text{step}$.

evolution_exact(rho0, tspan, Hamiltonian; diagonalize = false)

Calculate approximate time evolution of the rho0 based on $U(t)$ that is calculated for $t_0$ and $t_\text{step}$. The diagonalize argument decompose Hamiltonian into $\lambda_i, S, S^{-1}$ and calculate the exponential using eigenvalue decomposition. See evolution_operator. This method returns tspan and Vector of Operators.

evolution_exact(rho0, tspan, Hamiltonian; diagonalize = false)

Calculate exact time evolution of the rho0 inplace based on $U(t)$. The diagonalize argument decompose Hamiltonian into $\lambda_i, S, S^{-1}$ and calculate the exponential using eigenvalue decomposition. See evolution_operator. This method returns tspan and Vector of Operators.

evolution_operator(Hamiltonian, t)

Get evolution operator as Operator using the definition

$U(t) = e^{-i H t / \hbar}, \quad \hbar = 1$.

evolution_operator_a(H_lambda, S, Sinv, t)

Get evolution operator as Array using the definition

$U(t) = S e^{-i H_\lambda t / \hbar} S^{-1},$

where $\hbar = 1$, $H_{\lambda,ii} = \lambda_i$ are eigenvalues of $H$, so $H$ has to be non-singular, otherwise $H_{\lambda,ij} = 0, i \neq j$. $S$ is obtained from eigendecomposition of $H$, for example

H_lambda, S = eigen(Hamiltonian.data)
Sinv = inv(S)
H_lambda = diagm(H_lambda)

and arguments have to be Arrays.

evolution_operator_array(Hamiltonian, tspan)

Get evolution operators as Vector or Operators using the definition

$U(t) = e^{-i H t / \hbar}, \quad \hbar = 1$.

evolution_operator_iterator(Hamiltonian, tspan; diagonalize = true, approximate = false)

Resumable function that returns evolution operator as Operator type at the time t from tspan. See evolution_operator. The diagonalize argument decompose Hamiltonian into $\lambda_i, S, S^{-1}$ and calculate the exponential using eigenvalue decomposition. The approximate option asusmes that tspan is made of equidistant points, therefore $U(t)$ has to be calculated for $t_0$ and $t_\text{step}$.

evolution_super_operator(Hamiltonian, t)

Get evolution operator as SuperOperator using the definition

$\mathcal{U}(t) \:\cdot\: = U(t) \:\cdot\: U^\dagger(t) = e^{-i H t / \hbar} \:\cdot\: e^{i H t / \hbar}, \quad \hbar = 1$.

evolution_operator_array(Hamiltonian, tspan)

Get evolution superoperators as Vector or SuperOperators using the definition

$\mathcal{U}(t) \:\cdot\: = U(t) \:\cdot\: U^\dagger(t) = e^{-i H t / \hbar} \:\cdot\: e^{i H t / \hbar}, \quad \hbar = 1$.

evolution_super_operator_iterator(Hamiltonian, tspan; 
	diagonalize = true, approximate = false)

Resumable function that returns evolution operator as Operator type at the time t from tspan. See evolution_super_operator. The diagonalize argument decompose Hamiltonian into $\lambda_i, S, S^{-1}$ and calculate the exponential using eigenvalue decomposition. The approximate option asusmes that tspan is made of equidistant points, therefore $U(t)$ has to be calculated for $t_0$ and $t_\text{step}$.

schroedinger(psi0, tspan, H; 
	reltol=1.0e-12, abstol=1.0e-12, fout=nothing, alg=OrdinaryDiffEq.Tsit5())

Integrate Schroedinger equation to evolve states or compute propagators

$\frac{d}{d t}\vert\psi(t)\rangle = - \frac{i}{\hbar} \hat H\vert\psi(t)\rangle, \quad \hbar = 1$.

Arguments

• psi0: Initial state vector (can be a bra or a ket) or initial propagator.
• tspan: Vector specifying the points of time for which output should be displayed.
• H: Arbitrary operator specifying the Hamiltonian.
• reltol: Relative tolerance for OrdinaryDiffEq solver and its inner states.
• abstol: Absolute tolerance for OrdinaryDiffEq solver and its inner states.
• fout=nothing: If given, this function fout(t, psi) is called every time an output should be displayed. ATTENTION: The state psi is neither normalized nor permanent! It is still in use by the ode solver and therefore must not be changed.
• alg: Algorithm with which OrdinaryDiffEq will solve Schroedinger equation.

schroedinger_dynamic(psi0, tspan, f; 
	reltol=1.0e-12, abstol=1.0e-12, fout=nothing, alg=OrdinaryDiffEq.Tsit5())

Integrate time-dependent Schroedinger equation to evolve states or compute propagators

$\frac{d}{d t}\vert\psi(t)\rangle = - \frac{i}{\hbar} \hat H(t)\vert\rho(t)\rangle, \quad \hbar = 1$.

Arguments

• psi0: Initial state vector (can be a bra or a ket) or initial propagator.
• tspan: Vector specifying the points of time for which output should be displayed.
• f: Function f(t, psi) -> H returning the time and or state dependent Hamiltonian.
• reltol: Relative tolerance for OrdinaryDiffEq solver and its inner states.
• abstol: Absolute tolerance for OrdinaryDiffEq solver and its inner states.
• fout=nothing: If given, this function fout(t, psi) is called every time an output should be displayed. ATTENTION: The state psi is neither normalized nor permanent! It is still in use by the ode solver and therefore must not be changed.
• alg: Algorithm with which OrdinaryDiffEq will solve Schroedinger equation.

liouvilleVonNeumann(rho0, tspan, H; 
	reltol=1.0e-12, abstol=1.0e-12, fout=nothing, alg=OrdinaryDiffEq.Tsit5())

Integrate Liouville-von Neumann equation to evolve states or compute propagators

$\frac{d}{d t} \rho(t) = - \frac{i}{\hbar} [ \hat H, \rho(t) ], \quad \hbar = 1$.

Arguments

• rho0: Initial state vector (can be a bra or a ket) or initial propagator.
• tspan: Vector specifying the points of time for which output should be displayed.
• H: Arbitrary operator specifying the Hamiltonian.
• reltol: Relative tolerance for OrdinaryDiffEq solver and its inner states.
• abstol: Absolute tolerance for OrdinaryDiffEq solver and its inner states.
• fout=nothing: If given, this function fout(t, rho) is called every time an output should be displayed. ATTENTION: The state rho is neither normalized nor permanent! It is still in use by the ode solver and therefore must not be changed.
• alg: Algorithm with which OrdinaryDiffEq will solve LvN equation.

master(W0, tspan, Ham; 
	reltol=1.0e-12, abstol=1.0e-12, int_reltol=1.0e-8, int_abstol=0.0, 
	fout=nothing, alg=DelayDiffEq.MethodOfSteps(DelayDiffEq.Vern6()))

Integrate Quantum Master equation

$\frac{d}{d t} \rho(t) = - \frac{i}{\hbar} [ \hat{H}, \rho(t_0) ] -\frac{1}{\hbar^2} \int_{t_0}^{t_1} \text{d} \tau \: [ \hat{H}, [ \hat{H}, \rho(\tau) ]] ,\quad \hbar = 1.$

Arguments

• W0: Initial state vector (can be a bra or a ket) or initial propagator.
• tspan: Vector specifying the points of time for which output should be displayed.s
• Ham: Arbitrary operator specifying the Hamiltonian.
• reltol: Relative tolerance for DiffEqCallbacks solver and its inner states.
• abstol: Absolute tolerance for DiffEqCallbacks solver and its inner states.
• int_reltol: Relative tolerance for QuadGK solver and its inner states.
• int_abstol: Absolute tolerance for QuadGK solver and its inner states.
• fout=nothing: If given, this function fout(t, rho) is called every time an output should be displayed. ATTENTION: The state rho is neither normalized nor permanent! It is still in use by the ode solver and therefore must not be changed.
• alg: Algorithm with which DiffEqCallbacks will solve QME equation.

master_int(W0, tspan, Ham_0, Ham_I; 
	reltol=1.0e-12, abstol=1.0e-12, int_reltol=1.0e-8, int_abstol=0.0, 
	fout=nothing, alg=DelayDiffEq.MethodOfSteps(DelayDiffEq.Vern6()))

Integrate Quantum Master equation

$\frac{d}{d t} \rho^{(I)}(t) = - \frac{i}{\hbar} [ \hat{H}_I^{(I)}(t), \rho^{(I)}(t_0) ] -\frac{1}{\hbar^2} \int_{t_0}^{t_1} \text{d} \tau \: [ \hat{H}_I^{(I)}(t), [ \hat{H}_I^{(I)}(\tau), \rho^{(I)}(\tau) ]]$

$H = H_S + H_B + H_I = H_0 + H_I, \quad \hbar = 1.$

Arguments

• W0: Initial state vector (can be a bra or a ket) or initial propagator.
• tspan: Vector specifying the points of time for which output should be displayed.s
• Ham_0: System and bath Hamiltonian as Operator.
• Ham_I: Interaction Hamiltonian as Operator.
• reltol: Relative tolerance for DiffEqCallbacks solver and its inner states.
• abstol: Absolute tolerance for DiffEqCallbacks solver and its inner states.
• int_reltol: Relative tolerance for QuadGK solver and its inner states.
• int_abstol: Absolute tolerance for QuadGK solver and its inner states.
• fout=nothing: If given, this function fout(t, rho) is called every time an output should be displayed. ATTENTION: The state rho is neither normalized nor permanent! It is still in use by the ode solver and therefore must not be changed.
• alg: Algorithm with which DiffEqCallbacks will solve QME equation.

LineshapeResult

Container for a normalized molecular lineshape (absorption or emission spectrum).

Fields

• freqs: frequency grid (same units as molecule energies, typically cm⁻¹).
• intensities: normalized spectral intensities (∫ intensities dν ≈ 1).

absorption_spectrum(mol; T=0.0, sigma=50.0, n_points=2000)

Compute the normalized absorption spectrum of a Molecule.

Transitions from thermally-weighted ground-electronic vibronic states to all excited-electronic vibronic states, broadened by a Gaussian of width sigma.

Arguments

• mol: Molecule instance.
• T: temperature in Kelvin (energies assumed in cm⁻¹; uses BOLTZMANN_CM).
• sigma: Gaussian broadening standard deviation (same units as energies; default 50 cm⁻¹).
• n_points: frequency grid resolution (default 2000).

emission_spectrum(mol; T=0.0, sigma=50.0, n_points=2000)

Compute the normalized emission spectrum of a Molecule.

Transitions from thermally-weighted excited-electronic vibronic states (Kasha's rule) to all ground-electronic vibronic states, broadened by a Gaussian of width sigma.

Arguments

• mol: Molecule instance.
• T: temperature in Kelvin.
• sigma: Gaussian broadening standard deviation (default 50 cm⁻¹).
• n_points: frequency grid resolution (default 2000).

forster_rate(J, mol_donor, mol_acceptor; T=0.0, sigma=50.0, n_points=2000)

Compute the Förster excitation energy transfer rate:

$k_{DA} = 2\pi |J_{DA}|^2 \int I_D(\omega) A_A(\omega)\, d\omega$

Arguments

• J: electronic coupling (cm⁻¹).
• mol_donor: donor Molecule.
• mol_acceptor: acceptor Molecule.
• T: temperature in Kelvin (default 0).
• sigma: Gaussian broadening width in cm⁻¹ (default 50).
• n_points: frequency grid resolution (default 2000).

Returns

Rate in cm⁻¹. To convert to 1/fs: multiply by 2π × c where c = 3×10⁻⁵ cm/fs.

forster_rate_matrix(aggCore; T=0.0, sigma=50.0, n_points=2000)

Compute the matrix of pairwise Förster excitation energy transfer rates for all molecules in an AggregateCore.

Returns an N×N matrix where entry (i,j) is the Förster rate from molecule i to molecule j. Diagonal entries are zero.

The electronic coupling $J_{ij}$ is read from aggCore.coupling[i+1, j+1] (the +1 offset accounts for the ground state row/column in the padded coupling matrix).

Arguments

• aggCore: AggregateCore instance.
• T: temperature in Kelvin (default 0).
• sigma: Gaussian broadening width in cm⁻¹ (default 50).
• n_points: frequency grid resolution (default 2000).

spectral_overlap(donor_emission, acceptor_absorption)

Compute the spectral overlap integral:

$J = \int I_D(\omega) A_A(\omega)\, d\omega$

Both spectra must be normalized (∫ dν = 1). Returns the overlap in units of 1/[frequency], where frequency is in the same units as the molecule energies.

Uses linear interpolation onto the donor emission frequency grid.

SpectralDensity(frequencies, huang_rhys)

Discrete spectral density defined by a set of vibrational mode frequencies and their Huang-Rhys factors. Used to compute lineshape functions for Förster and modified Redfield theories.

Fields

• frequencies: mode frequencies (cm⁻¹).
• huang_rhys: Huang-Rhys factors $S_\alpha$ for each mode.

lineshape_derivative(sd, t, T)

Compute the first derivative $\dot{g}(t)$ of the lineshape function.

$\dot{g}(t) = \sum_\alpha S_\alpha \omega_\alpha \left[(2\bar{n}_\alpha + 1)\sin\omega_\alpha t + i(\cos\omega_\alpha t - 1)\right]$

lineshape_function(sd, t, T)

Compute the lineshape function $g(t)$ for a discrete spectral density at temperature T (Kelvin):

$g(t) = \sum_\alpha S_\alpha \left[(2\bar{n}_\alpha + 1)(1 - \cos\omega_\alpha t) + i(\sin\omega_\alpha t - \omega_\alpha t)\right]$

Returns a complex number.

lineshape_second_derivative(sd, t, T)

Compute the second derivative $\ddot{g}(t)$ of the lineshape function.

$\ddot{g}(t) = \sum_\alpha S_\alpha \omega_\alpha^2 \left[(2\bar{n}_\alpha + 1)\cos\omega_\alpha t - i\sin\omega_\alpha t\right]$

reorganization_energy(sd::SpectralDensity)

Compute the reorganization energy $\lambda = \sum_\alpha S_\alpha \omega_\alpha$.

spectral_density(mol::Molecule)

Extract a SpectralDensity from a Molecule by reading the frequency and Huang-Rhys factor (shift/2) from each mode.

exciton_basis(aggCore, aggTools)

Diagonalize the system Hamiltonian (excited-state block) to obtain exciton energies and transformation coefficients.

Returns (energies, coefficients) where:

• energies: N-element vector of exciton energies (cm⁻¹).
• coefficients: N×N matrix where coefficients[m, k] is the contribution of site k to exciton state m.

exciton_spectral_density(site_sds, coefficients, m, n)

Compute the spectral density in the exciton basis for exciton pair (m, n):

$\tilde{S}_{mn,\alpha,k} = c_{mk}^2 c_{nk}^2 S_{k\alpha}$

Returns a combined SpectralDensity with all site contributions.

modified_redfield_rates(aggCore; T=300.0, t_max=1000.0, gamma=0.0, rtol=1e-6, atol=1e-8)

Compute the modified Redfield population transfer rate matrix in the exciton basis.

Returns (rates, energies, coefficients) where:

• rates: N×N rate matrix. rates[m, n] is the rate from exciton n to exciton m. Diagonal elements are set so that columns sum to zero (population conservation).
• energies: exciton energies.
• coefficients: exciton transformation coefficients.

The rate from exciton state n to exciton state m is:

$R_{mn} = 2 \operatorname{Re} \int_0^\infty dt\, \left[\ddot{g}_{mn}(t) + \left(\dot{g}_{mm}(t) - \dot{g}_{nn}(t)\right) \dot{g}_{mn}(t)\right] e^{+i\varepsilon_{mn}t - f_{mn}(t) - \gamma t}$

where $f_{mn}(t) = g_{mm}(t) + g_{nn}(t) - 2g_{mn}(t)$ and $\varepsilon_{mn} = \varepsilon_m - \varepsilon_n$.

Arguments

• aggCore: AggregateCore instance.
• T: temperature in Kelvin (default 300).
• t_max: upper integration limit (default 1000, increase for slow dynamics).
• gamma: phenomenological damping rate for convergence (default 0, set > 0 for weakly damped systems).
• rtol, atol: tolerances for numerical integration.

modified_redfield_dynamics(aggCore, rho0_exciton, tspan; T=300.0, t_max=1000.0, gamma=0.0)

Propagate exciton populations using modified Redfield rates.

Arguments

• aggCore: AggregateCore instance.
• rho0_exciton: initial population vector in the exciton basis (N-element vector).
• tspan: time points for output.
• T: temperature in Kelvin.
• t_max: integration limit for rate calculation.
• gamma: phenomenological damping for rate convergence.

Returns (tspan, populations) where populations is a Matrix{Float64} of size (length(tspan), N).

AggregateCore(molecules, dipoles)

Construct an AggregateCore where the coupling matrix is computed automatically from point-dipole interactions.

Arguments

• molecules: Vector of Molecule instances.
• dipoles: Vector of TransitionDipole instances (one per molecule).

TransitionDipole(position, dipole)

Stores the 3D position and transition dipole vector of a molecule, used to compute point-dipole couplings between molecules in an aggregate.

Arguments

• position: 3-element vector giving the molecular position in Å.
• dipole: 3-element transition dipole vector in Debye (magnitude encoded in norm).

coupling_from_dipoles(dipoles)

Build an N×N coupling matrix (in cm⁻¹) from a vector of TransitionDipoles. Diagonal entries are zero. Off-diagonal entry (i,j) is the point-dipole coupling between molecule i and molecule j.

Arguments

• dipoles: Vector of TransitionDipole instances, one per molecule.

dipole_dipole_coupling(td1, td2)

Compute the point-dipole coupling (in cm⁻¹) between two transition dipoles.

Uses the formula:

J = DIPOLE_COUPLING_FACTOR × (μ₁·μ₂ - 3(μ₁·R̂)(μ₂·R̂)) / R³

where R is the inter-molecular distance in Å, μ in Debye, and the prefactor converts D²/ų to cm⁻¹.

Arguments

• td1, td2: TransitionDipole instances.

ad(rho, rho_bath, agg, FCProd, aggIndices, vibindices)

This is the inverse operation to the trace over bath trace_bath and get_rho_bath defined as follows

$\rho = \operatorname{ad}\{\rho_\text{tr}, \rho_\text{bath} \}$

correlation_function(t, rho0_bath, Ham_0, Ham_I, agg, FCProd, aggInds, vibindices)

Get time dependent correlation function for a specified time t using following definition

$C(t) = \operatorname{tr}_B \{ \hat{H}_I^{(I)}(t) \hat{H}_I \rho_\text{bath}(0) \}$

get_rho_bath(rho, agg, FCProd, aggIndices, vibindices; justCopy=false)

This method will return the bath part of rho knowing the result of trace_bath defined as follows

$\rho_\text{bath} = \operatorname{tr}_S \{\rho\}$

$\rho_{\text{bath}, ab} = \rho_{ab} / \langle a \vert \operatorname{tr}_B \{ \rho \}\vert b \rangle$

trace_bath(rho, a, b, agg, FCProd, aggIndices, vibindices)

Trace out bath degrees of freedom from rho without the product of Franck-Condon factors. The trace will be done only on the Hilber space for electric bra part a and ket part b. Input density matrix rho is for the whole Hilber space. This method returns number.

trace_bath(rho, agg, FCProd, aggIndices, vibindices)

Trace out bath degrees of freedom from rho

$\rho_\text{tr} = \operatorname{tr}_B \{\rho\} = \sum_{k} \langle k \vert \left( \sum_{ab} \rho_{am, bn} \vert am \rangle \langle bn \vert \right)\vert k \rangle$

trace_bath_part(rho, a, b, agg, FCProd, aggIndices, vibindices)

Trace out bath degrees of freedom from rho without the product of Franck-Condon factors. The trace will be done only on the Hilber space for electric bra part a and ket part b. Input density matrix rho is only for the subspace. This method returns number.

trace_bath_slow(rho, agg, FCFact, aggIndices, vibindices)

Trace out bath degrees of freedom from rho without the product of Franck-Condon factors.

thermal_state(T, mu_array, Ham, aggIndices; 
	boltzmann_const = 0.69503476, diagonalize = false, diagonal = false)

Get initial state as thermal state excited with ultra-fast laser pulse. In this version we suppose that after the thermal state is excited with laser pulse, the whole population of ground state is distributed over electric states in mu_array. We assume Condon approximation.

$\rho_\text{thermal} = \exp( -\frac{i}{\hbar} H ), \quad \hbar = 1$`.

Arguments

• T: Temperature of the initial thermal state.
• mu_array: Vector of electric states in local basis, see electronic_indices. The first index is for ground state of the aggregate the rest are first excited states in local basis.
• Ham: Arbitrary operator specifying the Hamiltonian.
• aggIndices: Aggregate indices, see get_indices.
• boltzmann_const: Boltzmann const in $\mathrm{cm^{-1}}$.
• diagonalize: Decompose Hamiltonian into $\lambda_i, S, S^{-1}$ and calculate the exponential using eigenvalue decomposition.
• diagonal: Return only the diagonal part of the excited density matrix.

thermal_state_composite(T, mu_weighted, Ham, aggIndices; 
	boltzmann_const::Float64 = 0.69503476, diagonalize::Bool=false, diagonal=false)

Functionality of this method is similar to thermal_state, but the final state is constructed from partial thermal_states with weight specified in mu_weighted. For example

julia thermal_state_composite(T, [0.0, 0.8, 0.2], ...) = 0.8 * thermal_state(T, [1, 2, 1], ...) + 0.2 * thermal_state(T, [1, 1, 2], ...)`

$\rho_\text{thermal} = \exp( -\frac{i}{\hbar} H ), \quad \hbar = 1$.

Arguments

• T: Temperature of the initial thermal state.
• mu_weighted: Vector of weights of electric states in local basis, see electronic_indices.
• Ham: Arbitrary operator specifying the Hamiltonian.
• aggIndices: Aggregate indices, see get_indices.
• boltzmann_const: Boltzmann const in $\mathrm{cm^{-1}}$.
• diagonalize: Decompose Hamiltonian into $\lambda_i, S, S^{-1}$ and calculate the exponential using eigenvalue decomposition.
• diagonal: Return only the diagonal part of the excited density matrix.

MemoryKernel_1_traced(H_II_t, H_II_tau, W_bath, agg, FCProd, aggIndices, indicesMap)

Calculate the first part of Memory Kernel with the definition

$\mathcal{M}_1(t, \tau) = \operatorname{tr}_B \{ \hat{H}_I^{(I)}(t) \hat{H}_I^{(I)}(\tau) W_\text{bath} \}$.

Arguments

• H_II_t: Interaction Hamiltonian in interaction picutre at the time t, $\hat{H}_I^{(I)}(t)$.
• H_II_tau: Interaction Hamiltonian in interaction picutre at the time tau, $\hat{H}_I^{(I)}(\tau)$.
• W_bath: Density matrix representing bath part of the density matrix, see get_rho_bath.
• agg: Aggregate of molecules, see Aggregate.
• aggIndices: Aggregate indices, see get_indices.
• indicesMap: Aggregate vibrational indices, see get_indices_map.

MemoryKernel_2_traced(H_II_t, H_II_tau, W_bath, agg, FCProd, aggIndices, indicesMap)

Calculate the second part of Memory Kernel with the definition

$\mathcal{M}_2(t, \tau) = \operatorname{tr}_B \{ \hat{H}_I^{(I)}(t) W_\text{bath} \hat{H}_I^{(I)}(\tau) \}$.

Arguments

• H_II_t: Interaction Hamiltonian in interaction picutre at the time t, $\hat{H}_I^{(I)}(t)$.
• H_II_tau: Interaction Hamiltonian in interaction picutre at the time tau, $\hat{H}_I^{(I)}(\tau)$.
• W_bath: Density matrix representing bath part of the density matrix, see get_rho_bath.
• agg: Aggregate of molecules, see Aggregate.
• aggIndices: Aggregate indices, see get_indices.
• indicesMap: Aggregate vibrational indices, see get_indices_map.

MemoryKernel_3_traced(H_II_t, H_II_tau, W_bath, agg, FCProd, aggIndices, indicesMap)

Calculate the third part of Memory Kernel with the definition

$\mathcal{M}_3(t, \tau) = \operatorname{tr}_B \{ \hat{H}_I^{(I)}(\tau) W_\text{bath} \hat{H}_I^{(I)}(t) \}$.

Arguments

• H_II_t: Interaction Hamiltonian in interaction picutre at the time t, $\hat{H}_I^{(I)}(t)$.
• H_II_tau: Interaction Hamiltonian in interaction picutre at the time tau, $\hat{H}_I^{(I)}(\tau)$.
• W_bath: Density matrix representing bath part of the density matrix, see get_rho_bath.
• agg: Aggregate of molecules, see Aggregate.
• aggIndices: Aggregate indices, see get_indices.
• indicesMap: Aggregate vibrational indices, see get_indices_map.

MemoryKernel_4_traced(H_II_t, H_II_tau, W_bath, agg, FCProd, aggIndices, indicesMap)

Calculate the fourth part of Memory Kernel with the definition

$\mathcal{M}_4(t, \tau) = \operatorname{tr}_B \{ W_\text{bath} \hat{H}_I^{(I)}(\tau) \hat{H}_I^{(I)}(t) \}$.

Arguments

• H_II_t: Interaction Hamiltonian in interaction picutre at the time t, $\hat{H}_I^{(I)}(t)$.
• H_II_tau: Interaction Hamiltonian in interaction picutre at the time tau, $\hat{H}_I^{(I)}(\tau)$.
• W_bath: Density matrix representing bath part of the density matrix, see get_rho_bath.
• agg: Aggregate of molecules, see Aggregate.
• aggIndices: Aggregate indices, see get_indices.
• indicesMap: Aggregate vibrational indices, see get_indices_map.

MemoryKernel_traced(H_II_t, H_II_tau, W_bath, agg, FCProd, aggIndices, indicesMap)

Calculate Memory Kernel with the definition

$\mathcal{M}(t, \tau) = \operatorname{tr}_B \{ [ \hat{H}_I^{(I)}(t), [ \hat{H}_I^{(I)}(\tau), W_\text{bath} ]]\}$.

Arguments

• H_II_t: Interaction Hamiltonian in interaction picutre at the time t, $\hat{H}_I^{(I)}(t)$.
• H_II_tau: Interaction Hamiltonian in interaction picutre at the time tau, $\hat{H}_I^{(I)}(\tau)$.
• W_bath: Density matrix representing bath part of the density matrix, see get_rho_bath.
• agg: Aggregate of molecules, see Aggregate.
• aggIndices: Aggregate indices, see get_indices.
• indicesMap: Aggregate vibrational indices, see get_indices_map.

validate_state(rho; tol=1e-6) -> Bool

Check a density matrix (Operator) for physical validity. Warns (does not error) if trace deviates from 1 or if Inf/NaN values are found. Returns true if valid, false otherwise.

validate_trajectory(rho_t; tol=1e-6) -> Bool

Check a vector of density matrices for physical validity. Applies validate_state to each element. Returns true if all valid, false otherwise.

SimulationResult

Wrapper around solver output that provides a uniform interface regardless of which solver was used.

Fields

• tspan::Vector{Float64} — time points at which states were recorded.
• states::Vector{<:Operator} — density matrices at each time point.
• method::Symbol — solver method that produced this result.
• extra::Dict{Symbol,Any} — optional extra data (e.g. :W_1_bath_t).

populations(r::SimulationResult) -> Matrix{Float64}

Return a length(r) x dim real matrix whose (t, k) entry is the k-th diagonal element of the density matrix at time step t.

solve(W0, tspan, agg; method=:qme_exact_si, kwargs...)

Unified entry point for all solvers in OpenQuantumSystems.

Dispatches to the appropriate solver function based on method:

All keyword arguments except method and iterative_args are forwarded to the underlying solver.

For :qme_iterative, pass the required positional arguments (rho_0_int_t, W_0_bath_t) as a tuple via the iterative_args keyword:

solve(W0, tspan, agg; method=:qme_iterative,
      iterative_args=(rho_0_int_t, W_0_bath_t))

AnnihilationOperator{BL,BR}(basis_l, basis_r)

Dense annihilation operator as a mutable struct.

CreationOperator{BL,BR}(basis_l, basis_r)

Dense creation operator as a mutable struct.

MomentumOperator{BL,BR}(basisl, basisr)

Dense momentum operator as a mutable struct without mass and frequency.

PositionOperator{BL,BR}(basis_l, basis_r)

Dense position operator as a mutable struct without mass and frequency.

ShiftOperator{BL,BR}(basis_l, basis_r, shift)

Dense shift operator as a mutable struct using the definition

$D(\alpha) = \exp(\alpha a^\dagger - \alpha^* a)$.

Arguments

• basis_l: Bra basis.
• basis_r: Ket basis.
• shift: Shift or $\alpha$ parameter, can be complex number.

OneDenseOperator(basis_l, basis_r)

DenseOperator with ones on the diagonal.

Commutator(A)

Create commutator in a form of superoperator from a given operator.

$\text{Commutator}(A) \:\cdot\: = [A, \:\cdot\:]$

tracedistance(rho, sigma)

Trace distance between rho and sigma. It is defined as

\[T(ρ,σ) = \frac{1}{2} \operatorname{tr}\{\sqrt{(ρ - σ)^† (ρ - σ)}\}.\]

It calls tracenorm which in turn either uses tracenorm_h or tracenorm_nh depending if $ρ-σ$ is hermitian or not.

tracedistance_h(rho, sigma)

Trace distance between rho and sigma. It uses the identity

\[T(ρ,σ) = \frac{1}{2} \operatorname{tr}\{\sqrt{(ρ - σ)^† (ρ - σ)}\} = \frac{1}{2} \sum_i |λ_i|\]

where $λ_i$ are the eigenvalues of rho - sigma.

tracedistance_nh(rho, sigma)

Trace distance between rho and sigma. Note that in this case rho and sigma don't have to be represented by square matrices (i.e. they can have different left-hand and right-hand bases). It uses the identity

\[ T(ρ,σ) = \frac{1}{2} \operatorname{tr}\{\sqrt{(ρ - σ)^† (ρ - σ)}\} = \frac{1}{2} \sum_i σ_i\]

where $σ_i$ are the singular values of rho - sigma.

tracenorm(rho)

Trace norm of rho. It is defined as

\[T(ρ) = \operatorname{tr}\{\sqrt{ρ^† ρ}\}.\]

Depending if rho is hermitian either tracenorm_h or tracenorm_nh is called.

tracenorm_h(rho)

Trace norm of rho. It uses the identity

\[T(ρ) = \operatorname{tr}\{\sqrt{ρ^† ρ}\} = \sum_i |λ_i|\]

where $λ_i$ are the eigenvalues of rho.

tracenorm_nh(rho)

Trace norm of rho. Note that in this case rho doesn't have to be represented by a square matrix (i.e. it can have different left-hand and right-hand bases). It uses the identity

\[ T(ρ) = \operatorname{tr}\{\sqrt{ρ^† ρ}\} = \sum_i σ_i\]

where $σ_i$ are the singular values of rho.

integrate(tspan, df::Function, x0::Vector{ComplexF64},
        state::T, dstate::T, fout::Function; kwargs...)

Integrate using OrdinaryDiffEq