% typeset for review
\documentclass[journal = polymers, supfile, submit, moreauthors, pdftex, 10pt, a4paper]{mdpi}

%=================================================================
\firstpage{1} 
\makeatletter 
\setcounter{page}{\@firstpage} 
\makeatother
\pubvolume{xx}
\issuenum{1}
\articlenumber{1}
\pubyear{2018}
\copyrightyear{2018}
\externaleditor{Academic Editor: name}
\history{Received: date; Accepted: date; Published: date}
%\updates{yes} % If there is an update available, un-comment this line

%------------------------------------------------------------------
% The following line should be uncommented if the LaTeX file is uploaded to arXiv.org
%\pdfoutput=1

%=================================================================
% Add packages and commands here. The following packages are loaded in our class file: fontenc, calc, indentfirst, fancyhdr, graphicx, lastpage, ifthen, lineno, float, amsmath, setspace, enumitem, mathpazo, booktabs, titlesec, etoolbox, amsthm, hyphenat, natbib, hyperref, footmisc, geometry, caption, url, mdframed, tabto, soul, multirow, microtype, tikz

\graphicspath{{./images/}}
\usepackage{amssymb}
\usepackage{adjustbox}

\definecolor{darkblue}{rgb}{0,0.08,0.45}
\definecolor{darkgreen}{rgb}{0.2,0.6,0.2}

\newcommand*{\citen}[1]{%
\begingroup
\romannumeral-`\x % remove space at the beginning of \setcitestyle
\setcitestyle{numbers}%
\cite{#1}%
\endgroup   
}

\newenvironment{DIFnomarkup}{}{}
\newcommand{\todo}[1]{\textcolor{red}{\textit{#1} -- TODO}}
\newcommand{\evan}[1]{\textcolor{blue}{\textit{#1} -- EVAN}}
\newcommand{\ej}[1]{\textcolor{red}{~\textit{#1} -- EJ~}}
\newcommand{\mj}[1]{\textcolor{darkgreen}{\textit{#1} --MJ}}
\newcommand{\newtxt}[1]{\textcolor{red}{\textit{#1}-added}}
\newcommand{\deltxt}[1]{\textcolor{blue}{\textit{#1}-removed}}

\def\mobunits{cm$^{2}$/Vs}
%=================================================================

%% Please use the following mathematics environments: Theorem, Lemma, Corollary, Proposition, Characterization, Property, Problem, Example, ExamplesandDefinitions, Hypothesis, Remark, Definition
%% For proofs, please use the proof environment (the amsthm package is loaded by the MDPI class).

% Full title of the paper (Capitalized)
\Title{Tying together multiscale calculations for charge transport in P3HT: structural descriptors, morphology, and tie-chains}

% Author Orchid ID: enter ID or remove command
\newcommand{\orcidauthorA}{0000-0001-6831-6165} % Add \orcidA{} behind the author's name
\newcommand{\orcidauthorB}{0000-0003-2255-6686} % Add \orcidB{} behind the author's name
\newcommand{\orcidauthorC}{0000-0002-3267-1410} % Add \orcidC{} behind the author's name

% Authors, for the paper (add full first names)
\Author{Evan D. Miller $^{1}$\orcidA{}, Matthew L. Jones $^{1}$\orcidB{} and Eric Jankowski $^{1,}$*\orcidC{}}
%\dagger,\ddagger

% Authors, for metadata in PDF
\AuthorNames{Evan D. Miller, Matthew L. Jones and Eric Jankowski}

% Affiliations / Addresses (Add [1] after \address if there is only one affiliation.)
\address[1]{%
%$^{1}$ \quad Micron School of Materials Science and Engineering; ericjankowski@boisestate.edu \\
$^{1}$ \quad Micron School of Materials Science and Engineering; Boise State University, Boise 83705, United States of America; materials@boisestate.edu\\
}

% Contact information of the corresponding author
\corres{Correspondence: ericjankowski@boisestate.edu; Tel.: +1-208-426-5681}

% Current address and/or shared authorship
%\firstnote{Current address: Affiliation 3} 
%\secondnote{These authors contributed equally to this work.}

%\affiliation{$^{\dag}$Micron School of Materials Science and Engineering, Boise State University, 1910 University Dr., Boise, Idaho 83725, United States}



\abstract{
Evaluating new, promising organic molecules to make next-generation organic optoelectronic devices necessitates the evaluation of charge carrier transport performance through the semi-conducting medium.
In this work, we utilize quantum chemical calculations (QCC) and kinetic Monte Carlo (KMC) simulations to predict the zero-field hole mobilities of $\sim 100$ morphologies of the benchmark polymer poly(3-hexylthiophene), with varying simulation volume, structural order, and chain-length polydispersity.
Morphologies were generated previously using an optimized molecular dynamics force-field and represent a spectrum of nanostructured order.
We discover that a combined consideration of backbone clustering and system-wide disorder arising from side-chain conformations are correlated to hole mobility.
Furthermore, we show that strongly interconnected thiophene backbones are the most important requirement for efficient charge transport - not necessarily the highest degree of order in the system. %TODO: My only problem with this is the fact that amorphous p3ht really should have lower mobility than semi-crystalline p3ht, right?
By marrying QCC and KMC over multiple length- and time-scales, we demonstrate that it is now possible to routinely predict the relationship between molecular nanostructure and device performance.
}

\keyword{organic photovoltaics, charge transport, semi-empirical, kinetic Monte Carlo}


\begin{document}

\newpage

\section{ZINDO and DFT Comparison}

In this section, we compare the ZINDO calculation of electronic properties used in this investigation, to a more rigorous DFT method to determine the accuracy of semi-empirical frontier molecular orbital energy calculations for P3HT.
We use three representative P3HT chromophore pairs selected from an equilibrated, ordered test morphology, visualizations of which are depicted in \autoref{fig:QCCCompare}.
The DFT calculations were performed using the B3LYP functional \cite{Becke1993} and the 6311++g** basis set \cite{Krishnan1980}.
    
\begin{figure}[h!]
\centering
\includegraphics[width=\linewidth]{{chromophores}.pdf}
\caption{
    The three representative chromophore pairs used to investigate the accuracy of the ZINDO/S semiempirical method.
    a) \textbf{0469-3714}, b) \textbf{0841-1237}, c) \textbf{2032-2900}.
Terminating hydrogens were added for the QCCs based on position in the thiophene ring.
}
\label{fig:QCCCompare}
\end{figure}
    

\begin{table*}[h!]
\caption{A comparison of the HOMO splitting and calculated transfer integrals for three representative P3HT chromophore pairs}
\small
\centering
\begin{tabular}{ccc}
\toprule
\textbf{HOMO Splitting}&\textbf{DFT (eV)}&\textbf{ZINDO/S (eV)}\\
\midrule
0469-3714&0.196&0.095\\
0841-1237&0.199&0.058\\
2032-2900&0.086&0.008\\
\bottomrule
\end{tabular}
\label{table:DFT_ZINDO}
\end{table*}

The calculated electronic properties of the chromophore pairs are shown in \autoref{table:DFT_ZINDO}.
ZINDO appears to consistently underpredict the HOMO splitting, which would lead to lower transfer integrals and slower transport than expected from more rigorous DFT methods.
However, the ZINDO calculations provide good agreement with the DFT results to within $\sim 100$ meV, which is already the rough cutoff for DFT accuracy.
Furthermore, changes in transfer integral of factors of 2-3 are not expected to significantly affect the charge transport properties given that morphological changes can result in orders of magnitude differences.
ZINDO calculations can be performed within 5-10 seconds for a chromophore pair, compared to several minutes to half an hour in the case of more rigorous DFT calculations (depending on the DFT level desired).
As such, the computational throughput is significantly improved at the cost of the smaller reduction in accuracy - a critically important point given that each morphology can contain upwards of 40,000 chromophore pairs to be considered.

In summary, the vastly improved computational efficiency at the cost of a small reduction in accuracy of ZINDO/S justifies our use of the semi-empirical calculations for our charge transport properties instead of more rigorous DFT methodologies.

\section{Developing $\psi^{\prime}$ to Explicitly Consider Transfer Integrals}

One possible short-coming in our modified order parameter $\psi^{\prime}$ is that ir is an indirect predictor of the electronic transfer integral $J_{i,j}$ between two chromophores.
The hopping rate between two choromophores is strongly dependent on the electronic transfer integral $J_{i,j}$ (see Equation 2 of the main text) and the energy level difference $\Delta E_{i,j}$ of the two chromophores.
Additionally, $\Delta E_{i,j}$ is also partially encoded into $J_{i,j}$, in which chromophores with incompatible energy levels (for instance, a large $\Delta E_{i,j}$) will also reduce $J_{i,j}$ (see Equation 1 of the main text).
As such, the transfer integral seems like a good ``one-size-fits-all'' parameter to describe clustering.

\begin{figure}[h!]
\centering
%\begin{figure*}[t!]
    \includegraphics[width=\linewidth]{{cry_semi_dis_2a}.pdf}
%\begin{tabular}{ccc}
%    \includegraphics[width=0.33\linewidth]{{cryst_TI}.pdf}&
%    \includegraphics[width=0.33\linewidth]{{semicry_TI}.pdf}&
%    \includegraphics[width=0.33\linewidth]{{dis_TI}.pdf}\\
%\end{tabular}
\caption{
Distributions of chromophore Voronoi neighbor transfer integrals for the representative 1,000 molecule a) amorphous, b) semi-crystalline, and c) crystalline morphologies.
The red line shows the Gaussian filtered distribution shape that was used to determine the cluster cut-off criterion.
The black vertical line shows the value of the cut-off criterion, which was automatically determined to be at the minimum for each system - $J_{i,j} >$ 0.562, 0.549, and 0.457 eV for the crystalline, semi-crystalline, and amorphous morphologies respectively.
}
\label{fig:alg2a_TI}
\end{figure}

The transfer integral distributions for each representative system are shown in \autoref{fig:alg2a_TI}.
In all three cases, the distribution has a large spike at very low transfer integrals and a bump at high TI corresponding to pairs within the same P3HT chain.
Initially, we set the transfer integral cut-off to the location of the minimum for each morphology, such that only connections with transfer integrals greater than the cut-off are added to the same cluster.
It is convenient to set cut-offs to maxima and minima as these can be determined automatically, rather than being calibrated manually for each separate system.
For the crystalline, semi-crystalline, and amorphous morphologies, the cut-offs were set to $J_{i,j} > 0.562$, 0.549, and 0.457 eV respectively.

\begin{figure}[h!]
\centering
%\begin{figure*}[t!]
    \includegraphics[width=6in]{{clust_2a}.pdf}\\
\caption{
    Visualizations of the clusters in the a) amorphous, b) semi-crystalline, and c) crystalline systems with size > 6 monomer units.
    Clusters were determined based on an automatically-defined transfer integral cut-off for each system based on the distributions in \autoref{fig:alg2a_TI}.
}
\label{fig:alg2a_clust}
\end{figure}

\begin{table}[h]
\caption{
    Table of cluster statistics for the three systems, given the automatically-determined transfer integral cut-off criteria.
}
\small
\centering
\begin{adjustbox}{max width=\textwidth}
\begin{tabular}{cccc}
\toprule
\textbf{Property}       &\textbf{Amorphous}   &\textbf{Semi-Crystalline}&\textbf{Crystalline}       \\
\midrule
Mobility (\mobunits)    & $1.02 \times 10^{-1}$  & $1.63 \times 10^{-2}$  & $1.16 \times 10^{-1}$  \\
J$_{i,j}$ cut-off (eV)  & $0.457$                & $0.549$                & $0.562$                \\
Total clusters          & $1067$                 & $1065$                 & $972$                  \\
Large ($> 6$) clusters  & $964$                  & $941$                  & $873$                  \\
Largest cluster size    & $60$                   & $60$                   & $90$                   \\
\bottomrule
\end{tabular}
\end{adjustbox}
\label{table:alg2a}
\end{table}

The resultant cluster visualization in \autoref{fig:alg2a_clust} suggests that these cut-off values are too large - in all morphologies, hops with $J_{i,j} > \sim0.5$ eV are generally only intra-molecular hops (red region in \autoref{fig:alg2a_TI}).
This leads to nearly every chain in the system being considered an individual cluster, with few occurrences of clusters forming between multiple chains.
There is no significant difference in the cluster distribution throughout the morphology between the three systems, suggesting that a different transfer integral cut-off should be used.

\begin{figure}[h!]
\centering
%\begin{figure*}[t!]
    \includegraphics[width=\linewidth]{{cry_semi_dis_2b}.pdf}
% \begin{tabular}{ccc}
%     \includegraphics[width=0.33\linewidth]{{cryst_TI_2}.pdf}&
%     \includegraphics[width=0.33\linewidth]{{semicry_TI_2}.pdf}&
%     \includegraphics[width=0.33\linewidth]{{dis_TI_2}.pdf}\\
% \end{tabular}
\caption{
    Distributions of chromophore Voronoi neighbor transfer integrals for the representative 1,000 molecule a) amorphous, b) semi-crystalline, and c) crystalline morphologies.
    The red line shows the Gaussian filtered distribution shape that was used to determine the cluster cut-off criterion.
    The black vertical line shows the value of the cut-off criterion, $J_{i,j} > 0.2$ eV.
}
\label{fig:alg2b_TI}
\end{figure}

We can, for instance, reduce the cut-off to something smaller in order to include higher $J_{i,j}$ inter-molecular hops.
This however, has the short-coming in that such a selection will likely be artitrarily chosen, not rather than an automatically identified minimum.
Regardless, reducing the $J_{i,j}$ cut-off to 0.2 eV (\autoref{fig:alg2b_TI}) provides significantly improved results as now a non-negligible proportion of inter-molecular hops have $J_{i,j} >$ cut-off, thereby, allowing clusters to form between molecules.

\begin{figure}[h!]
\centering
    \includegraphics[width=6in]{{clust_2b}.pdf}\\
\caption{
    Visualizations of the clusters in the a) amorphous, b) semi-crystalline, and c) crystalline systems with size > 6 monomer units, given the following clustering criteria: transfer integral $> 0.2$ eV.
}
\label{fig:alg2b_clust}
\end{figure}


\begin{table}[h]
\caption{
    Table of cluster statistics for the three systems, given the following clustering criterion: transfer integral $J_{i,j} > 0.2$ eV.
}
\small
\centering
\begin{adjustbox}{max width=\textwidth}
\begin{tabular}{cccc}
\toprule
\textbf{Property}&\textbf{Amorphous}&\textbf{Semi-Crystalline}&\textbf{Crystalline}\\
\midrule
Mobility                & $1.02 \times 10^{-1}$  & $1.63 \times 10^{-2}$  & $1.16 \times 10^{-1}$   \\
J$_{i,j}$ cut-off (eV)  & $0.200$                & $0.200$                & $0.200$                 \\
Total clusters          & $289$                  & $163$                  & $64$                    \\
Large ($> 6$) clusters  & $273$                  & $151$                  & $51$                    \\
Largest cluster size    & $2564$                 & $9914$                 & $12837$                 \\
\bottomrule
\end{tabular}
\end{adjustbox}
\label{table:alg2b}
\end{table}

Now, we compare the clusters identified with the $J_{i,j}$ cut-off between the three systems.
The crystalline morphology shows one large cluster (shown in red) and a few smaller clusters with opposing grain orientations, indicating that the crystalline system will have a high connectivity.
Conversely, the amorphous morphology is predicted to have poor connectivity based on this clustering metric stemming from the larger number of small clusters. 
Howevever, the connectivity in the semi-crystalline morphology again shows a cluster arrangement intermediate between the other two.
This cluster arrangement would again predict a high mobility for the crystalline morphology, a low mobility for the amorphous morphology and an intermediate morphology in the semi-crystalline case, which is contrary to our mobility calculations.
Varying the cut-off to any consistent value between the three morphologies always results in this conclusion, suggesting that the transfer integral distribution is not an adequate way to identify clusters in the morphology.

\section{Clustering Based on Hops}

One short-coming of the previous clustering algorithms is that it considers charge transport between two chromophores in isolation.
However, in the KMC algorithm, hops to all neighboring chromophores are considered and the preferential hop (based on the hopping rate between $i$ and $j$ and the random number $x$) will be chosen.
As such, a ``good'' hop may not occur because there is a better hop.

\begin{figure}[h!]
\centering
%\begin{figure*}[t!]
     \includegraphics[width=\linewidth]{{cry_semi_dis_3}.pdf}
\caption{
    Distributions of the frequencies with which carriers hop between chromophore Voronoi neighbors for the representative 1,000 molecule a) amorphous, b) semi-crystalline, and c) crystalline morphologies.
    The red line shows the Gaussian filtered distribution shape that was used to determine the cluster cut-off criterion.
    The black vertical line shows the value of the cut-off criterion, which was automatically determined to be at the final minimum of the frequency distribution: a total of 3264, 1566, and 1635 hops for the crystalline, semi-crystalline, and amorphous systems respectively.
}
\label{fig:alg3_freq}
\end{figure}

As such, defining clusters based on regions in which charges will freely move is prudent, however, we must still identify a sensible cut-off in hopping frequency to separate these regions.
The distributions of total hole hops between chromophore pairs in the three representative systems are shown in \autoref{fig:alg3_freq}.
Note that the x-axis in these plots is logarithmic, leading to quantization of the bins on the left-hand side of the plot.
In all three systems, a second peak appears at high hop frequencies.
This leads to a local minimum at 3264 hops in the crystalline case, 1566 hops in the semi-crystalline case, and 1635 hops in the amorphous case.
We therefore use these values as the clustering criteria - only chromophores with connections that are used more than this number during the simulation will be added to the same cluster.
We note that the exact values of the cut-off criteria are strongly dependent on the duration of the KMC simulation; the value may change significantly if fewer carriers iterations are performed or if simulation times are reduced.
In this study, all three systems used the same simulation time-scales for KMC and the same number of carriers were averaged over in order to obtain the charge transport properties.

\begin{figure}[h!]
\centering
%\begin{figure*}[t!]
    \includegraphics[width=\linewidth]{{clust_3}.pdf}
%\begin{tabular}{ccc}
%    \includegraphics[width=0.33\linewidth]{{cryst_3_clust}.pdf}&
%    \includegraphics[width=0.33\linewidth]{{semicry_3_clust}.pdf}&
%    \includegraphics[width=0.33\linewidth]{{dis_3_clust}.pdf}\\
%\end{tabular}
\caption{
    Visualizations of the clusters in the a) amorphous, b) semi-crystalline, and c) crystalline systems with size > 6 monomer units.
    Clusters were determined based on an automatically-defined hopping frequency cut-off for each system based on the distributions in \autoref{fig:alg3_freq}.
}
\label{fig:alg3_clust}
\end{figure}
    
\begin{table}[h]
\caption{
    Table of cluster statistics for the three systems, given the automatically-defined total hop frequency cut-off criteria.
}
\small
\centering
\begin{adjustbox}{max width=\textwidth}
\begin{tabular}{cccc}
\toprule
\textbf{Property}   &\textbf{Amorphous} &\textbf{Semi-Crystalline}  &\textbf{Crystalline}           \\
\midrule
Mobility                & $1.02 \times 10^{-1}$  & $1.63 \times 10^{-2}$  & $1.16 \times 10^{-1}$  \\
Hop frequency cut-off   & $1635$                 & $1566$                 & $3264$                 \\
Total clusters          & $410$                  & $1376$                 & $418$                  \\
Large ($> 6$) clusters  & $134$                  & $209$                  & $73$                   \\
Largest cluster size    & $11867$                & $3164$                 & $10254$                \\
\bottomrule
\end{tabular}
\end{adjustbox}
\label{table:alg3}
\end{table}

The cluster visualizations using the hop frequency cut-off are shown in \autoref{fig:alg3_clust} are very promising.
The crystalline morphology shows different crystalline grains very clearly - although the majority of the simulation is a single cluster (red), a large cluster with a different grain orientation is clearly visible in the morphology (blue).
It is harder to distinguish the cluster distributions of the semi-crystalline and amorphous systems using the visualizations, although the crystals present in the semi-crystalline morphology are clearly resolvable from the amorphous matrix surrounding them.
However, \autoref{table:alg3} shows the first set of cluster properties where the semi-crystalline morphology is not intermediate between the crystalline and amorphous system, in terms of the number of total clusters and the largest cluster size.
These clusters describe regions of the morphology that carriers are frequently hopping within.
With this definition, hops within the regions are more common than those between clusters and so carriers are effectively trapped in this region - time is still progressing as they hop around, but their mean squared displacement is not significantly increasing.
Therefore, a small number of large clusters is advantageous, whereas a large number of small clusters will strongly restrict charge transport properties.


\section{Intra-cluster trapping}


\begin{figure*}[!h]
\centering
%\begin{figure*}[t!]
\includegraphics[width=\linewidth]{{NetworkGraphs}.pdf}
\caption{
The network diagrams for the (a) amorphous, (b) semi-crystalline, and (c) crystalline systems show carrier pathways between connected chromophores (insets: zoomed regions).
Connections are colored based on a perceptually uniform, logarithmic heatmap, where brighter zones correspond to more frequently used pathways.
}
\label{fig:network}
\end{figure*}

In our investigation, we record the location history of every carrier as it hops through the system.
Using the carrier hopping history, we can construct network connectivity diagrams (\autoref{fig:network}) to observe the most frequently travelled paths for charges through the morphology.
These network connectivity diagrams are constructed by identifying the frequency with which holes in the KMC simulation hop between pairs of chromophores.
The centre-of-mass locations of the chromophores then form the nodes of the network, and the shortest paths between each of the chromophore pairs become the edges.
The ``net hopping frequency'' is calculated by subtracting the frequency of forward hops from backward hops and taking the absolute value.
These values are normalized to the highest net hopping frequency in the system, and then assigned a color based on the logarithmic color map to highlight preferred carrier transport routes through the morphology.

The differences in structure between the three classes of morphology are clearly evident in \autoref{fig:network}.
The amorphous network graph (\autoref{fig:network}a) shows that no crystallites have formed in the system.
There are several high-traffic nodes spread homogeneously throughout the system, explaining the highly isotropic carrier trajectory presented in the main text.
The crystalline network graph (\autoref{fig:network}c) shows the lamellar structure of the system, with nearly all chains aligned in layers moving left-to-right across the morphology.
The most frequently used pathways are along chains, and there are many connections in the $\pi$-stacking direction between chains within the crystal.
The semi-crystalline network graph (\autoref{fig:network}b) exhibits behaviour intermediate between the other two - crystallites with varying grain orientations are clearly visible, within an amorphous matrix.

The insets in \autoref{fig:network} show a zoomed region in the corner of the morphology, to highlight an area of `cross-hatching' in the network, where carriers frequently loop around the same subset of chromophores (located at the vertices of the patterns), without increasing mean squared displacement (MSD) from their initial position.
In the amorphous morphology (\autoref{fig:network}a), no loops are observed, and therefore every hop (no matter how slow) is contributing to the MSD, increasing mobility.
In the semi-crystalline (\autoref{fig:network}b) and crystalline (\autoref{fig:network}c) systems, significant looping can be seen.
Carriers in these regions are becoming `trapped' by the loops - even though transport may be fast between the chromophores, it becomes more difficult for the carriers to leave the crystal along the slow transport routes due to the wealth of fast hops available within.
For the crystalline case, the morphology is dominated by one large crystallite that extends across the full simulation volume.
Therefore, carriers getting trapped inside this crystallite are still able to move long distances, and the penalty to the mobility from the trapping is lessened (reflected by a high mobility and a larger anisotropy in the main text).
However, in the semi-crystalline case, the morphology is composed of multiple crystallites with various orientations, with loops present across all three dimensions.
Trapping therefore has a more significant effect - carriers get stuck in the small loops and are unable to increase their mean squared displacement over time in a single direction, restricting the carrier mobility within the system.
These conclusions are supported by the cluster maps presented in \autoref{fig:network}, as well as the cluster properties presented in Table 1; the crystalline and amorphous systems are dominated by a single, well-connected cluster of chromophores permitting a high mobility, whereas the semi-crystalline system is composed of many clusters with differing grain orientations.
The visualizations of the network in \autoref{fig:network} serve to provide additional evidence as to why the clusters described in the main text form within these morphologies.


\section{Polydisperse Simulations}

\begin{figure*}[!h]
\centering
%\begin{figure*}[t!]
\includegraphics[width=6in]{{pdi-si}.pdf}
\caption{
(a) The scaled target distribution of chain lengths and the histogram of the actual chain lengths used in the polydisperse simulations.
(b) The equilibrated polydisperse systems is able to produce periodic features not seen in systems where all the chains are longer.
}
\label{fig:distribution}
\end{figure*}

Here we present how we generate polydisperse P3HT simulations.
This can be broken into two steps: first creating a dictionary of P3HT oligomers of varying lengths from 1 to 50 monomers long.
Second is using a distribution to determine the amount of each chain length to place into the simulation.
To produce chains of arbitrary length, we use the open-source program mBuild in which a polymer can be easily created using monomer building blocks. 
We limit the chain length used in this study at 50 monomers long as to avoid unphysical interactions of chains feeling themselves across periodic boundaries.
To generate the distribution of chain-lengths, we use the Schulz-Flory distribution which is a commonly used mathematical description for polymer lengths in the form \cite{Flory1942b}:

\begin{equation}
P_{L} = \alpha ^2 D_P ((1-\alpha)^{D_P-1}),
\end{equation}
in which $P_L$ is the probability of seeing a chain of a given length, $D_P$ is the degree of polymerization of a particular chain, and $\alpha$ is a tunable parameter which affects the shape of the distribution.
The value for $\alpha$ used in this study was 0.1 and was chosen as this value produces polydispersities of $\sim$ 1.8.

To create the actual distribution of chain lengths we utilize a simple Monte Carlo algorithm.
In this algorithm we select a random chain length between 1 and 50 and a random number ($x$) associated with this chain length between 0 and 1.
If $x$ is less than the probability of seeing a chain of that length $P(L)$ we accept the chain otherwise the chain is rejected.
In addition to this, to ensure that we have the same number of monomers as the other simulations (15,000 monomers) we keep track of the number of monomers which have been added to the simulation.
When the number of monomers added to the simulation via the Monte Carlo algorithm is less than 50 monomers from 15,000 monomers, we terminate the Monte Carlo algorithm and add the remaining monomers to the simulation via one more single chain so that the total number of monomers is 15,000 monomers.

After the Monte Carlo algorithm is finished and a distribution of the chain lengths suggested, we calculate the polydispersity of the simulation with:

\begin{equation}
PDI = \frac{M_w}{M_n}, 
\end{equation}
in which $M_w$ is the weight average molecule weight and $M_n$ is the number average molecule weight.
$M_w$ and $M_n$ can be calculated with:

\begin{align}
M_w = \frac{\sum N_i M_i^2}{\sum N_i M_i}
\end{align}
\begin{align}
M_n = \frac{\sum N_i M_i}{\sum N_i},
\end{align}
in which $N_i$ is the number of chains of that length and $M_i$ is the molecular weight of that chain length.
If the PDI of the stochastically generated distribution of chain lengths is below 1.8, we reject the distribution and regenerate the distribution until PDI $\geq$ 1.8.
A comparison between the target distribution and the histogram of chain lengths is shown in \autoref{fig:distribution}a.

The distribution of chains presented in \autoref{fig:distribution}a is able to produce ordered morphologies with periodic features (\autoref{fig:distribution}b) along $(100)$ and $(010)$ signifying $\pi$- and alkyl-stacking.
These features are seen in experimental and 15mer scattering patterns \cite{Miller2018}.
However, when simulating systems that contain only 50 membered chains, the system requires much longer to relax into these periodic structures.

\acknowledgments{
This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575 \cite{Towns2014}.
This material is based upon work supported by the National Science Foundation under Grant No. (1229709) and (1653954), and (1658076).
The authors thank Chris Muhich for his assistance in comparing ZINDO/S to DFT.
}

\authorcontributions{
Funding acquisition, Eric Jankowski; Investigation, Evan Miller and Matthew Jones; Software, Evan Miller and Matthew Jones; Supervision, Eric Jankowski; Writing – original draft, Evan Miller, Matthew Jones and Eric Jankowski; Writing – review \& editing, Evan Miller, Matthew Jones and Eric Jankowski.
}
 
\conflictsofinterest{The authors declare no conflict of interest.}

\abbreviations{

15mer - P3HT chain containing 15 monomers

50mer - P3HT chain containing 50 monomers

$\varepsilon_s$ - solvent quality

KMC - Kinetic Monte Carlo

MD - Molecular Dynamics

$\mu_{0}$ Zero-field Mobility

OPLS - Optimized Performance for Liquid Simulations

OPV - Organic Photovoltaic

P3HT - Poly(3-hexylthiophene)

$\psi$ - order parameter

$\psi\prime$ - modified order parameter

$\rho$ - density

$\sigma$ - standard deviation

$T$ - Temperature

VRH - Variable Range Hopping

QCC - Quantum Chemical Calculations
}

\bibliography{library}

\end{document}
