\section{Macroscopic observables}
\label{sec:analysis}%
Spatial distributions of charge and current densities can provide a better insight in the microscopic mechanisms of charge transport.
%
If $O$ is an observable which has a value $O_\alpha$ in a state $\alpha$, its ensemble average at time $t$ is a sum over all states weighted by the probability $P_\alpha$ to be in a state $\alpha$ at time $t$
%
\begin{equation}
\label{equ:ensemble}
\left< O \right> = \sum_{\alpha} O_\alpha P_\alpha.
\end{equation}
%
If $O$ does not explicitly depend on time, the time evolution of $\left< O \right>$ can be calculated as
\begin{equation}
\begin{split}
\frac{d \left< O \right>}{dt} = \sum_{ \alpha, \beta}
\left[ P_\beta \Omega_{\beta \alpha} -
P_\alpha \Omega_{\alpha \beta} \right]
O_\alpha %\\
%
= \sum_{ \alpha, \beta}
P_\beta \Omega_{\beta \alpha}
\left[ O_\alpha - O_\beta \right] .
\end{split}
\end{equation}
%
If averages are obtained from KMC trajectories, $P_\alpha = s_\alpha / s$, where $s_\alpha$ is the number of Markov chains ending in the state $\alpha$ after time $t$, and $s$ is the total number of chains.
Alternatively, one can calculate time averages by analyzing a single Markov chain. If the total occupation time of the state $\alpha$ is $\tau_\alpha$ then
\begin{align}
\label{equ:time}
\overline{ O }
= \frac{1}{\tau} \sum_{\alpha} O_\alpha \tau_\alpha \,,
\end{align}
where $\tau = \sum_{\alpha} \tau_\alpha$ is the total time used for time averaging.
For ergodic systems and sufficient sampling times, ensemble and time averages should give identical results.
In many cases, the averaging procedure reflects a specific experimental technique. For example, an ensemble average over several KMC trajectories with different starting conditions corresponds to averaging over injected charge carriers in a time-of-flight experiment. In what follows, we focus on the single charge carrier (low concentration of charges) case.
\subsection{Charge density}
\label{sec:occupation}
For a specific type of particles, the microscopic charge density of a site $i$ is proportional to the occupation probability of the site, $p_i$
\begin{equation}
\rho_i = e p_i / V_i\, ,
\end{equation}
where, for an irregular lattice, the effective volume $V_i$ can be obtained from a Voronoi tessellation of space. For reasonably uniform lattices (uniform site densities) this volume is almost independent of the site and a constant volume per cite, $V_i = V/N$, can be assumed. In the macroscopic limit, the charge density can be calculated using a sxtpthing kernel function, i.e. a distance-weighted average over multiple sites. Site occupations $p_i$ can be obtained from \equ{ensemble} or \equ{time} by using the occupation of site $i$ in state $\alpha$ as an observable.
If the system is in thermodynamic equilibrium, that is without sources or sinks and without circular currents (and therefore no net flux) a condition, known as detailed balance, holds
%
\begin{equation}
\label{equ:detailed_balance}
p_j \omega_{ji} = p_i \omega_{ij},
\end{equation}
%
It can be used to test whether the system is ergodic or not by correlating $\log p_i$ and the site energy $E_i$. Indeed, if $\lambda_{ij} = \lambda_{ji}$ the ratios of forward and backward rates are determined solely by the energetic disorder, $\omega_{ji} / \omega_{ij} = \exp(-\Delta E_{ij} / k_\text{B} T)$ (see \equ{marcus}).
\subsection{Current}
\label{sec:vaverage}
\index{current}
If the position of the charge, $\vec{r}$, is an observable, the time evolution of its average $\left<\vec{r}\right>$ is the total current in the system
\begin{equation}
\vec{J} = e \left< \vec{v} \right> = e \frac{d \left< \vec{r}
\right>} {dt} = e \sum_{i, j} p_{j} \omega_{ji} ( \vec{r}_i -
\vec{r}_j ) .
\label{equ:current_def}
\end{equation}
Symmetrizing this expression we obtain
\begin{equation}
\vec{J} = \frac{1}{2} e \sum_{i, j} \left( p_{j}
\omega_{ji} - p_{i} \omega_{ij} \right) \vec{r}_{ij} ,
\label{equ:current}
\end{equation}
where $\vec{r}_{ij} = \vec{r}_{i} - \vec{r}_{j}$. Symmetrization ensures equal flux
splitting between neighboring sites and absence of local average fluxes in equilibrium. It allows to define a local current through site $i$ as\index{current!local}
\begin{equation}
\vec{J_i} = \frac{1}{2} e \sum_{ j} \left( p_{j} \omega_{ji} - p_{i} \omega_{ij} \right) \vec{r}_{ij} .
\label{equ:site_current}
\end{equation}
A large value of the local current indicates that the site contributes considerably to the total current. A collection of such sites thus represents most favorable charge pathways~\cite{van_der_holst_modeling_2009}.
\subsection{Mobility and diffusion constant}
\label{sec:mobility}
\index{mobility}
For a single particle, e.g. a charge or an exciton, a zero-field mobility can be determined by studying particle diffusion in the absence of external fields. Using the particle displacement squared, $\Delta {\bm r}_i^2$, as an observable we obtain
\begin{equation}
\begin{split}
2d D_{\gamma \delta} = \frac{d \left< \Delta{r}_{i, \gamma} \Delta{r}_{i, \delta} \right>}{dt}
= \sum_{\substack{i,j \\ i\ne j}} p_j\omega_{ji}
\left( \Delta r_{i,\gamma}\Delta r_{i,\delta} - \Delta r_{j,\gamma}\Delta r_{j,\delta} \right)
= \sum_{\substack{i,j\\ i\ne j}} p_j \omega_{ji} \left( r_{i,\gamma} r_{i,\delta} - r_{j,\gamma} r_{j,\delta} \right) \, .
\end{split}
\label{equ:diffusion}
\end{equation}
Here $\vec{r}_i$ is the coordinate of the site $i$, $D_{\gamma \delta}$ is the diffusion tensor, $\gamma, \delta = x,y,z$, and $d=3$ is the system dimension. Using the Einstein relation,
\begin{equation}
D_{\gamma \delta} = k_\text{B}T \mu_{\gamma \delta} \, ,
\end{equation}
one can, in principle, obtain the zero-field mobility tensor $\mu_{\gamma \delta}$. \Equ{diffusion}, however, does not take into account the use of periodic boundary conditions when simulating charge dynamics. In this case, the simulated occupation probabilities can be compared to the solution of the Smoluchowski equation with periodic boundary conditions (see the supporting information for details).
Alternatively, one can directly analyze time-evolution of the KMC trajectory and obtain the diffusion tensor from a linear fit to the mean square displacement, $\overline{ \Delta{r}_{i, \gamma} \Delta{r}_{i, \delta}} = 2d D_{\gamma \delta} t$.
The charge carrier mobility tensor, $\hat{\mu}$, for any value of the external field can be determined either from the average charge velocity defined in
\equ{current_def}
\begin{equation}
%\begin{split}
\langle \vec{v} \rangle = \sum_{i,j} p_j \omega_{ji} (\vec{r}_i - \vec{r}_j) = \hat{\mu} \vec{F} \, ,
%\end{split}
\end{equation}
or directly from the KMC trajectory. In the latter case the velocity is calculated from the unwrapped (if periodic boundary conditions are used) charge displacement vector divided by the total simulation time. Projecting this velocity on the direction of the field $\vec{F}$ yields the charge carrier mobility in this particular direction. In order to improve statistics, mobilities can be averaged over several KMC trajectories and MD snapshots.