%  LaTeX support: latex@mdpi.com 
%  In case you need support, please attach all files that are necessary for compiling as well as the log file, and specify the details of your LaTeX setup (which operating system and LaTeX version / tools you are using).

% You need to save the "mdpi.cls" and "mdpi.bst" files into the same folder as this template file.
%
%--------------------
% Class Options:
%--------------------
% journal
%----------
% Choose between the following MDPI journals:
% acoustics, actuators, addictions, admsci, aerospace, agriculture, agronomy, algorithms, animals, antibiotics, antibodies, antioxidants, applsci, arts, asi, atmosphere, atoms, axioms, batteries, bdcc, behavsci, beverages, bioengineering, biology, biomedicines, biomimetics, biomolecules, biosensors, brainsci, buildings, carbon, cancers, catalysts, cells, ceramics, challenges, chemengineering, chemosensors, children, cleantechnol, climate, clockssleep, cmd, coatings, colloids, computation, computers, condensedmatter, cosmetics, cryptography, crystals, cybersecurity, data, dentistry, designs, diagnostics, dairy, diseases, diversity, drones, econometrics, economies, education, electrochem, electrochemistry, electronics, energies, entropy, environments, epigenomes, est, fermentation, fibers, fire, fishes, fluids, foods, forecasting, forests, fractalfract, futureinternet, galaxies, games, gastrointestdisord, gels, genealogy, genes, geohazards, geosciences, geriatrics, hazardousmatters, healthcare, heritage, highthroughput, horticulturae, humanities, hydrology, informatics, information, infrastructures, inorganics, insects, instruments, ijerph, ijfs, ijms, ijgi, ijtpp, inventions, j, jcdd, jcm, jcs, jdb, jfb, jfmk, jimaging, jof, jintelligence, jlpea, jmmp, jmse, jpm, jrfm, jsan, land, languages, laws, life, literature, logistics, lubricants, machines, magnetochemistry, make, marinedrugs, materials, mathematics, mca, medsci, medicina, medicines, membranes, metabolites, metals, microarrays, micromachines, microorganisms, minerals, modelling, molbank, molecules, mps, mti, nanomaterials, ncrna, neonatalscreening, neuroglia, nitrogen, nutrients, ohbm, particles, pathogens, pharmaceuticals, pharmaceutics, pharmacy, philosophies, photonics, plants, plasma, polymers, polysaccharides, proceedings, processes, proteomes, publications, quaternary, qubs, reactions, recycling, religions, remotesensing, reports, resources, risks, robotics, safety, sci, scipharm, sensors, separations, sexes, sinusitis, smartcities, socsci, societies, soilsystems, sports, standards, stats, surfaces, surgeries, sustainability, symmetry, systems, technologies, toxics, toxins, tropicalmed, universe, urbansci, vaccines, vehicles, vetsci, vibration, viruses, vision, water, wem, wevj
%---------
% article
%---------
% The default type of manuscript is article, but can be replaced by: 
% abstract, addendum, article, benchmark, book, bookreview, briefreport, casereport, changes, comment, commentary, communication, conceptpaper, correction, conferenceproceedings, conferencereport, expressionofconcern, meetingreport, creative, datadescriptor, discussion, editorial, essay, erratum, hypothesis, interestingimages, letter, meetingreport, newbookreceived, opinion, obituary, projectreport, reply, reprint, retraction, review, perspective, protocol, shortnote, supfile, technicalnote, viewpoint
% supfile = supplementary materials
% protocol: If you are preparing a "Protocol" paper, please refer to http://www.mdpi.com/journal/mps/instructions for details on its expected structure and content.
%----------
% submit
%----------
% The class option "submit" will be changed to "accept" by the Editorial Office when the paper is accepted. This will only make changes to the frontpage (e.g. the logo of the journal will get visible), the headings, and the copyright information. Also, line numbering will be removed. Journal info and pagination for accepted papers will also be assigned by the Editorial Office.
%------------------
% moreauthors
%------------------
% If there is only one author the class option oneauthor should be used. Otherwise use the class option moreauthors.
%---------
% pdftex
%---------
% The option pdftex is for use with pdfLaTeX. If eps figures are used, remove the option pdftex and use LaTeX and dvi2pdf.

%=================================================================

%=================================================================
% 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{\mh}[1]{\textcolor{magenta}{\textit{#1} --MH}}
\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{Optimization and Validation of Efficient Models for Predicting Polythiophene Self-Assembly}

% 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
\newcommand{\orcidauthorD}{0000-0002-3870-9993} % Add \orcidD{} behind the author's name
\newcommand{\orcidauthorE}{0000-0003-2442-8055} % Add \orcidE{} behind the author's name
\newcommand{\orcidauthorF}{0000-0001-6083-6502} % Add \orcidF{} behind the author's name

% Authors, for the paper (add full first names)
\Author{Evan D. Miller $^{1}$\orcidA{}, Matthew L. Jones $^{1}$\orcidB{}, Michael M. Henry$^{1}$\orcidD, Paul Chery$^2$\orcidE{}, Kyle Miller$^3$\orcidF{} and Eric Jankowski $^{1,}$*\orcidC{}}
%\dagger,\ddagger

% Authors, for metadata in PDF
\AuthorNames{Evan D. Miller, Matthew L. Jones, Michael M. Henry, Paul Chery, Kyle Miller, 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\\
$^{2}$ \quad Physics; Macalester College, St Paul, 55105, United States of America; \\
$^{3}$ \quad Physics; University of Puget Sound, Tacoma, 98416, United States of America; \\
}

% Contact information of the corresponding author
\corres{Correspondence: ericjankowski@boisestate.edu; Tel.: +x-xxx-xxx-xxxx}

% 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}

%% Title, authors and addresses
\abstract{
%The self-assembled active-layer morphology of organic electronic devices determines their performance, so there is a need 
We develop an optimized OPLS-UA force-field for the benchmark donor polymer poly(3-hexylthiophene) (P3HT) and demonstrate its utility for predicting thermodynamic self-assembly.
In particular, we consider short oligomer chains, model electrostatics and solvent implicitly, and coarsely model solvent evaporation.
We quantify the performance of our model to determine what the optimal system sizes are for exploring self-assembly at combinations of state variables.
We perform molecular dynamic simulations to predict the self-assembly of P3HT at $\sim$ 350 combinations of temperature and solvent quality.
Our structural calculations predict that the highest degrees of order are obtained with good solvents just below the melting temperature.
We find our model produces the most accurate structural predictions to date, as measured by agreement with grazing incident X-ray scattering experiments.
}

\keyword{molecular dynamics, organic photovoltaics, coarse-graining}

\begin{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newpage

\section{Force-Field Constraints for Bonded Atoms}

Here we present the intra-molecular interactions, e.g. two-body bond, three-body angle, and four-body dihedral constraints governing our P3HT model.
These bonded interactions are adapted from the atomistic models developed by Huang and Bhatta \cite{Huang2010, Bhatta2013a}.
As we employ united-atom models in this work, we compare our force-field values with the OPLS-UA force-field to ensure that we are not introducing spurious interactions when adapting the atomistic model \cite{Jorgensen1984}.

Bond potentials are harmonic:

\begin{equation}
    E_{\text{bond}} = \frac{1}{2}k(l-l_{0})^{2},
\end{equation}
in which $k$ is the spring constant and $l_0$ is the equilibrium bond length.
The $k$ and $l_{0}$ values are presented in \autoref{tab:p3ht_bonds}.
It should be noted that the CA-CA bond parameters presented in \autoref{tab:p3ht_bonds} only apply to bonds between thiophene rings.
The thiophene rings themselves are represented as rigid bodies, and so bonds within the rings do not change during the simulation.
If a flexible model were to be used for the thiophene ring, then distinct bead types would be needed to distinguish bonds between and within rings.

\begin{table}[h]
    \caption{The bond constraints used in our force-field to simulate P3HT.
    $^{\dag}$: Bonds between two monoers, rather than within a single, rigid thiophene. 
    $^{\ddag}$: Bonds completely enclosed by a rigid body, which are therefore fixed at $l_{0}$.
    }
\centering
\small
\begin{tabular}{ccc}
\toprule
\textbf{Bond} & \textbf{$l_{0}$ (\AA)} & \textbf{k (kcal mol$^{-1}$ \AA$^{-2}$)}\\
\midrule
CA-CA$^{\dag}$ & 1.43 & 392.0 \\
CT-CT & 1.53 & 268.0 \\
CA-CT & 1.51 & 300.0 \\
CA-S$^{\ddag}$ & 1.71 & 291.0 \\
\bottomrule
\end{tabular}
\label{tab:p3ht_bonds}
\end{table}

Angles constraints are also harmonic:

\begin{equation}
    E_{\text{angle}} = \frac{1}{2}k(\theta-\theta_{0})^{2},
\end{equation}
in which $\theta_0$ is the equilibrium bond angle.
The angle constraints are presented in \autoref{tab:p3ht_angle}.
As before, the CA-CA-CA and CA-CA-S angles specified describe angle constraints between monomers, rather than those in the rigid thiophene ring.

\begin{table}[h]
    \caption{
        The angle constraints used in our force-field to simulate P3HT.
        $^{\dag}$: Angles between two monomers, rather than within a single, rigid thiophene.
        $^{\ddag}$: Angles completely enclosed by a rigid body, which are therefore fixed at $\theta_{0}$.
    }
\centering
\small
\begin{tabular}{ccc}
\toprule
\textbf{Angle} & \textbf{$\theta$ (rad)} & \textbf{k (kcal mol$^{-1}$ rad$^{-2}$)} \\
\midrule
CA-CA-CA$^{\dag}$ & 2.27 & 54.7 \\
CA-CA-S$^{\dag}$ & 2.09 & 41.7 \\
CA-S-CA$^{\ddag}$ & 1.62 & 86.0 \\
CA-CT-CT & 2.16 & 70.0 \\
CA-CA-CT & 2.15 & 70.0 \\
CT-CT-CT & 1.97 & 37.5 \\
\bottomrule
\end{tabular}
\label{tab:p3ht_angle}
\end{table}

The dihedral parameters are defined by a multi-harmonic function in the form:

\begin{equation}
    E_{\text{dihedral}} = \sum_{n=0}^{4} k_{n} cos^{n} \varphi,
\end{equation}
in which $\varphi$ is the dihedral angle.
The dihedral constraints are presented in \autoref{tab:p3ht_dihedral}.

\begin{table}[h]
    \caption{
        The angle constraints used in our force-field to simulate P3HT.
        $^{\dag}$: Dihedrals between two monomers, rather than within a single, rigid thiophene.
    }
\centering
\small
\begin{tabular}{cccccc}
\toprule
\textbf{Dihedral} & \textbf{k$_0$ (kcal mol$^{-1}$)} & \textbf{k$_1$} & \textbf{k$_2$} & \textbf{k$_3$} & \textbf{k$_4$} \\
\midrule
S-CA-CA-S & 2.9533 & 0.1571 & -4.2326 & 0.39979 & 1.8855 \\
CA-CA-CA-S$^{\dag}$ & 2.9533 & -0.1571 & -4.2326 & -0.39979 & 1.8855 \\
CA-CA-CT-CT & 0.3175 & 1.127 & 14.143 & -22.297 & 6.7188 \\
CA-CT-CT-CT & 2.4469 & -6.3946 & 10.747 & 30.695 & 11.139 \\
CT-CT-CT-CT & 1.8922 & -3.4904 & 1.4665 & 7.1418 & 0.2859 \\
\bottomrule
\end{tabular}
\label{tab:p3ht_dihedral} 
\end{table}

\begin{figure}[!h]
\centering
    \includegraphics[width=3.25in]{{dihedral_explained}.pdf}
    \caption{
The flexible dihedrals used in this investigation that are not completely enclosed by a rigid body.
Reducing the number of distinct atom types in this model (compared to the all-atom models from which they are derived) leads to some conflicts in the dihedral constraints.
To avoid this, the coefficients for some of these dihedrals (CA-CA-CA-CT and CA-CA-CA-CA) are set to zero. 
    }
    \label{fig:dihedralexplained}
\end{figure}

Due to the rigid bodies used in this investigation, some dihedrals become ambiguously defined if the force-field from the literature are used.
For instance, the CA-CA-CA-CT and the CA-CA-CA-CA dihedral constraints are already considered by the CA-CA-CT-CT and S-CA-CA-S dihedrals respectively.
Furthermore, in the case of CA-CA-CA-CT, two possible definitions are applicable depending on whether the CA is located in the same ring as the alkyl sidechain or not.
When the CA is in a neighboring ring to the monomer containing the CT, the dihedral is in the cis conformation, however, when CA is in the same ring as the specified side chain, the dihedral is in the trans conformation.
This ambiguity leads to instabilities in the system when the dihedral parameters are not set to 0. 
As such, we set the CA-CA-CA-CT and CA-CA-CA-CA dihedral coefficients to zero, effectively deactivating them in our force-field.
\autoref{fig:dihedralexplained} shows which dihedral constraints are considered (solid lines) and which are omitted (dashed lines), to highlight the redundancies that exist due to reduction to three atom types and the incorporation of the rigid body.
We therefore assume that the other carbon aromatic to chain (CA-CA-CT-CT and CA-CT-CT-CT) dihedrals are sufficient in describing the position of the chains relative to the rings and that the (S-CA-CA-S and CA-CA-CA-S) dihedrals are sufficient in describing the orientation of the thiophene rings along the backbone.

\newpage

\section{The Effect of Including Explicit Charges in the Model}

Modeling the long-range electrostatic interactions between P3HT chains is both computationally expensive and challenging to do accurately.
Including explicit partial charges results in a factor of $~3$ increase in simulation time for the systems studied here.
Additionally, small conformational changes in conjugated polymers such as P3HT give rise to significant changes in electron densities and therefore the partial charges associated with each atom fluctuate over time and space.
Furthermore, prior calculations for modeling P3HT electrostatics do not reach consensus on which partial charges are correct, or sufficient \cite{Bhatta2013a, Dag2010a, To2012}
In this section we compare our base case implicit charge model against one with explicit partial charges whose forces are calculated with the Fourier based particle-particle-particle-mesh Ewald summation method \cite{LeBard2012a}.
We quantify structural and performance differences between these approaches and show the implicit charge model sufficiently captures the relevant assembly physics.

The charges used in this study are determined through first-principle calculations with the NWChem software \cite{Valiev2010} with the Becke, three-parameter, Lee-Yang-Parr (B3YLP) hybrid functional \cite{Becke1993} with the 6-311++g** basis set \cite{Krishnan1980}.
The partial charges for the atoms within the thiophene ring from literature and those calculated in this study are shown in \autoref{table:atomidentities}.
\begin{table}[h]
\caption{The charges for atoms given in literature for all-atom simulations. 
*Indicates that this was a united-atom simulation.
\dag Indicates that these were calculated as part of this work by the NWChem program \cite{Valiev2010}.
}
\small
\centering
\begin{tabular}{cccccc}
\toprule
\textbf{Study} & \textbf{S} & \textbf{C1} & \textbf{C2} & \textbf{C3} & \textbf{C4}\\
\midrule
Bhatta \cite{Bhatta2013a} & -0.22 & 0.18 & -0.04 & -0.31 & 0.03 \\
Obata \cite{Obata2013} & -0.07 & 0.12 & -0.07 & -0.27 & 0.05 \\
Moreno \cite{Moreno2010} & -0.12 & -0.01 & -0.01 & -0.18 & -0.03 \\
Huang \cite{Huang2010} & -0.15 & -0.14 & 0.075 & -0.18 & -0.18 \\
Borzdun \cite{Borzdun2016} & 0.33 & -0.24 & -0.28 & -0.16 & 0.08 \\
D'Avino* \cite{DAvino2013} & -0.06 & -0.21 & 0.22 & 0.10 & -0.29 \\
Trimer\dag & -0.17 & 0.28 & -0.03 & -0.37 & -0.11 \\
Pentamer\dag & -0.16 & 0.24 & 0.01 & -0.37 & -0.14 \\
\bottomrule
\end{tabular}
\label{table:atomidentities}
\end{table}
From the values presented in \autoref{table:atomidentities} we utilize the pentamer values for our simulations (\autoref{fig:charged}a).
Further, we note that we also zero-out the charge in the simulation by subtracting the average per-particle charge from every simulation bead to account for rounding errors that may occur on individual atoms.

\begin{figure*}[h]
\centering
    \includegraphics[width=1.in]{{thiophene_atom_labels}.pdf}
    \caption{\label{fig:thiopheneatomlabels} Labels for identifying the atoms in the thiophene in \autoref{table:atomidentities}.
    }
\end{figure*}

Since a united-atom model is used in this work, we sum the partial charges of the hydrogen atoms into the atom to which they are directly bonded.
These values are shown in \autoref{table:chainidentities}.
The CT$N$ labels, in which $N$=1$\rightarrow$6, are used to identify the aliphatic carbons with CT1 being bonded to the thiophene ring and CT6 being the farthest away from the thiophene ring.
We note that the C3 presented in \autoref{table:chainidentities} is the same as that shown in \autoref{table:atomidentities}, but with the hydrogen considered.
Because the works of Moreno and Huang consider polythiophenes with no alkyl sidechains, we exclude these works from the list of the side-chain charges.

\begin{table}[h]
\caption{The charges from literature summed up so the hydrogens are considered with the carbons.
C3 is the same as above. CT1 is bonded to the thiophene ring and each CT$N=1\rightarrow6$ is extending from the ring.}
\small
\centering
\begin{tabular}{cccccccc}
\toprule
\textbf{Study} & \textbf{C3} & \textbf{CT1} & \textbf{CT2} & \textbf{CT3} & \textbf{CT4} & \textbf{CT5} & \textbf{CT6} \\
\midrule
Bhatta \cite{Bhatta2013a} & -0.08 & 0.14 & -0.04 & 0.03 & 0.04 & -0.03 & -0.01 \\
Obata \cite{Obata2013} & -0.11 & 0.07 & 0.03 & -0.04 & 0.01 & 0.09 & -0.08 \\
Borzdun \cite{Borzdun2016} & 0.06 & 0.02 & -0.01 & -0.01 & 0.0 & 0.0 & 0.0 \\
D'Avino \cite{DAvino2013} & 0.10 & 0.22 & 0.08 & -0.07 & 0.01 & 0.04 & -0.03 \\
Trimer & -0.12 & 0.19 & 0.06 & -0.16 & 0.11 & 0.07 & -0.10 \\
Pentamer & -0.13 & 0.19 & 0.07 & -0.17 & 0.10 & 0.09 & -0.10 \\
\bottomrule
\end{tabular}
\label{table:chainidentities}
\end{table}

We assign partial charges to simulation beads (\autoref{fig:charged}a).
We then employ scattering, RDF, and cluster analyis described above to to determine the differences in molecular packing.
In our optimized implicit-charge model, we reduce the thiophene ring bead diameters to $\sim$ 3.44 \AA~from the $\sim$ 3.7 \AA~used in OPLS-UA and Amber \cite{Jorgensen1984, Yang2006a}, which permits the thiophene rings to $\pi$-stack at the same length-scale as in experiments.
In the explicit charge model benchmarked here, we employ the partial charges, but do not modify the Lennard-Jones $\epsilon$ or $\sigma$ parameters from our optimized implicit-charge model.
These parameters \textit{should} be updated to create an optimized explicit-charge model, but we hypothesize that such efforts are not worth their cost and check here to see how the addition of explicit charges affects structure and performance.

\begin{figure*}[t!]
\centering
\includegraphics[width=5.5in]{{charged-morph}.pdf}
\caption{
(a) The partial charges applied to the thiophene ring in the explicit-charge model. 
These charges were generated from a pentamer chain using first-principle calculations.
Incorporating these into the equilibrated implicit-charge structure at the same state-point Figure 7a results in (b), which is visually indistinguishable (CA-dark blue, S-yellow, CT-cyan; state-point: explicit charges, $T \sim 600$ K, $\varepsilon_{s} = 0.8$, $\rho = 1.11$ g/cm$^{3}$).
Short- and long-range structural features observed in the (c) GIXS (averaged over 6 orientations) and (d) thiophene centroid RDF (compared to Figure 6a and Figure 7e respectively) show that the explicit-charge model with $\sigma_{LJ} = 3.7$ \AA~is consistent with the implicit-charge model with a reduced $\sigma_{LJ} = 3.44$.
(e) The explicit- (red circles) and implicit-charge (blue triangles) systems exhibit similar trends in order ($\psi$) over different temperatures (state-point: $\varepsilon_{s}$ = 0.2 and $\rho$ = 0.72 g/cm$^{3}$).
However, explicit-charge systems systematically obtain a lower degree of order.
}
\label{fig:charged}
\end{figure*}

To test whether explicit charges destabilize an already-equilibrated structure from an implicit-charge model (Figure 7a), we instantaneously ``turn on'' electrostatics with the above partial charges and re-equilibrate a 100 15-mer $T=600$ K, $\rho=1.11$ g/cm$^3$ Protocol (2) system over 0.3 $\mu$s. 
The resulting morphology (\autoref{fig:charged}b) is not significantly destabilized, and shows qualitative agreement with GIXS and RDF metrics of the implicit-charge case.
The primary quantitative difference arises from the electrostatic repulsion of identical beads when thiophene rings are aligned: With explicit charges, aligned rings are on average 4.0 \AA~apart, instead of 3.9 \AA~as in the implicit charge case.

To test differences in self-assembled structure we perform new simulations at $\varepsilon_{s} = 0.2$ with Protocol (2) and the base-case 100 15-mers, but with explicit electrostatic charges active throughout. 
\autoref{fig:charged}e shows that order parameter $\psi$ dependence on temperature is qualitatively the same for both implicit- and explicit-charge models.
Independent of charge consideration, we observe kinetic trapping at low temperatures ($<$ 250 K), increased order at moderate temperatures (250 K to 500 K), and decreasing order as temperatures increase ($>$ 750 K).
The explicit-charge model is systematically less ordered ($\psi_{explicit}\sim \psi_{implicit}+0.1$) than the implicit model across all temperatures studied here, but otherwise assembles the same structures: lamellae of $\pi$-stacked backbones below a transition temperature and disordered melts above.
On one hand, we interpret these results to justify the factor of 3 performance benefit of the implicit charge model.
On the other hand, we acknowledge that the self-assembled structures are not identical between these two cases, and this \textit{could} have significant ramifications, for example in predicting charge mobilities.
We therefore recommend using implicit charge models for screening large parameter spaces as performed here, and then subsequently performing more expensive simulations when they are warranted for additional insight.

\newpage

\section{The Effect of Considering Longer Chains}

In this section we compare performance and packing of ``long'' 50 monomer (50mer) chains relative to the 15mer base-case used throughout this work.
Commercially available P3HT chains (from e.g. Sigma-Aldrich, Reike Metals) are significantly longer, typically ranging from 125 to 625 monomers, and more polydisperse (PDI$\sim 2$) despite the feasibility of synthesizing narrow chain length distributions \cite{Bronstein2009,Koch2013a,Grubbs2017}.
However, 50mers are the longest chains that we can simulate at the same conditions as the 15mers, while ensuring a chain cannot interact with itself through the periodic simulation boundaries.
This is a conservative modeling choice, as we do not anticipate that such self-interactions would significantly impact self-assembled morphology.
Additionally, we expect the 50mer chains to provide useful information about the effect of simulating longer chains compared to the 15mer base case, despite being shorter than commercially available.
We simulate 300 50mer chains at $\rho =1.11$ g/cm$^3$ using Protocol (2) at $T \sim 600$ K and $\varepsilon_s=0.8$, and compare against 1000 15-mers at the same conditions as subsection 4.5 of the maintext.

The 50mer systems (\autoref{fig:polymer}a) demonstrate less ordering than their 15mer counterparts under the same conditions.
Using clustering analysis, we still see short-range order exists in the form of small crystallites that are randomly oriented with respect to each other (\autoref{fig:polymer}b).
The lack of long-range order for 50mers is apparent from GIXS (\autoref{fig:polymer}c), which lacks distinct reflections corresponding to lamellae of $pi$-stacked backbones as seen before.
The thiophene ring RDF (\autoref{fig:polymer}d) shows some short-range ordering, with anti-aligned thiophene rings (expected peak at $5.3$ \AA) less common in 50mer systems. 
The major structural difference between 50mers and 15mers is that the 50mers are observed to ``fold'' and form $\pi$-stacks with themselves (\autoref{fig:polymer}e)
Such self-stacking allows a single chain to form multiple layers in one crystallite, and is observed experimentally \cite{Ihn1993, Mena-Osteritz2000, Brinkmann2006a}.

\begin{figure*}[h]
\centering
\includegraphics[width=\linewidth]{{polymer}.pdf}
\caption{
(a) 300 $\times$ 50mer chains produce morphologies containing no visually discernible periodic order at $T \sim 600$ K and $\varepsilon_{s}$ = 0.8 (CA-dark blue, S-yellow, CT-cyan).
(b) However, short-range order is visible when thiophene rings are colored by cluster.
(c) Short-range order is confirmed through simulated GIXS (averaged over 60 orientations) and (d) the RDF of the thiophene centroids, which both show less defined peaks than in the 15mer case (Figure 6a and Figure 7e respectively).
(e) The 50mer chains are long enough to undergo $\pi$-stacking with themselves, which is observed experimentally \cite{Ihn1993, Mena-Osteritz2000, Brinkmann2006a} and not observed with the shorter chains.
}
\label{fig:polymer}
\end{figure*}


We estimate the time required for the system containing 300 $\times$ 50 monomer long chains (50mer) to relax to an equilibrated structure compared to the systems containing 1000, 15 monomers long chains (15mer) by comparing the evolution of the Lennard-Jones potential energies and also the structure factor over time.
As a first test, we assume that the systems have the same equilibrated structure and that these two structures will have the same non-bonded energies.
Therefore, we fit the 50mer and 15mer Lennard-Jones energy over time to equations and determine where the 50mer's equation will equal the final energy of the 15mer system.
In choosing what equation should be used to represent the energy evolution, we qualitatively observe that the P3HT systems undergoes three stages during relaxation: fast exponential decay, slow approximately linear decay and constant.
As such, we fit the per-particle potential decay over simulation time $t_{s}$ to:  

\begin{equation}\label{eq:fit}
f(t_{s}) = a \times exp(b \times t_{s}) + c \times t_{s} + d, 
\end{equation}
to get the coefficients $a$, $b$, $c$, and $d$ for the exponential and linear terms; the fitted curve is showed by the dashed red line in \autoref{fig:chain-relax}a.
With the coefficients known, we solve the equation to determine when it would equal the final energy of the 15mer system, which results in approximately 9,000 times longer for the 50mer system to run to equal the final energy as the 15mer system.

However, it is not necessarily valid to assume that the 50 and 15mer systems will have the same final energy.
For instance, it is possible that the 50 system (which is more prone to entanglements) will reach a metastable state and will require the lifetime of the universe to relax out of the metastable state.
Such behavior often occurs in real systems, in that P3HT devices have crystalline and amorphous regions rather than the more energetically favored perfect crystal.
As such, if we do not assume that the two systems will have the same equilibrium structure and hence same energy, we can instead utilize the rate of decay of the energy as our descriptor for system evolution.
Because the linear decay is the slowest process, we can compare the ratio of the linear coefficients between the 50mer and 15mer system to describe the difference in relaxation times.
As such, we find that the 50mer's coefficient is -1.94 $\times$ 10$^{-5}$ and the 15mer's is -3.6 $\times$ 10$^{-5}$.
This suggests a $\sim$ 2$\times$ faster decrease in energy in the 15mer system than the 50mer system.

\begin{figure*}[h]
\centering
\includegraphics[width=\linewidth]{{chain-relax}.pdf}
\caption{
(a) The potential energy in the (blue) 15mer system decays faster and achieves a final lower per-particle than the (red) 50mer system.
The structure factor - a radial average of the GIXS pattern - shows the growth of peaks corresponding to (red) $\pi$-stacking over time in the (b) 15mer and the (c) 50mer system.
(d) The evolution of the red peaks in b and c show that the 15mer system reaches a stable state faster than the 50mer system.
}
\label{fig:chain-relax}
\end{figure*}

Rather than using the energy, which is an indirect measure of the structure's evolution, we can instead utilize a direct measurement of the structure such as the structure factor.
We therefore conduct the scattering experiment to obtain the structure factor - a radial average of the GIXS pattern - every 25 ns of simulation run time (\autoref{fig:chain-relax}b and c).
We record the height of the peak in the structure factor located at 1.6 \AA$^{-1}$, corresponding to $\pi$-stacking, over time.
We again fit the evolution of the structure factor peak over time to an exponential function and compare the relative changes in structure evolution between the 50mer and 15mer systems (\autoref{fig:chain-relax}d).
From this, we again predict that the 50mer system will require twice as long to reach equivalent scattering intensities of that of the 15mer system.
Because the structure factor is the most direct calculation of structure and it agrees with the change of energy we conclude that the 50mer system requires twice as long as the 15mer system to order.
In summary, the longest chains that can be practically equilibrated should be used, 15mers are sufficient for predicting experimental GIXS patterns, and systems of 50mers can in principle be equilibrated, but the factor of two increase to simulation time precludes routine sweeps of large parameter spaces.

\newpage

\section{Linking System Evolution to Energy}

\begin{figure*}[h]
\centering
%\begin{figure*}[t!]
\includegraphics[width=\linewidth]{{pes-for-sizes}.pdf}
\caption{
The per-particle Lennard-Jones potential energies for the various systems shows three regions, highlighted with blue circles, cyan triangles, and green squares for the systems with various sizes.
}
\label{fig:pes-for-sizes}
\end{figure*}

In this section we present our analysis of the Lennard-Jones energy to infer structural evolution.
We detect three transients in the energy curve corresponding to different phases of morphology evolution: exponential decay representing formation of crystals from a disordered configuration, slow decay signifying the process of going from many small crystallites to few, highly ordered crystals and constant energy signifying equilibrium.
The regions corresponding to each transient are shown in \autoref{fig:pes-for-sizes} for simulations containing various amounts of 15mer chains.
We note that not all the final energies are exactly equal.
This is due to in part to systems reaching metastable states, which will require indefinite run times to relax out of.
Additionally, changes in the simulation volume between the simulations can lead to unfavorable box lengths.
For example, one system may be able to form a more perfect crystal when the simulation box length is commensurate with the crystal structure - resulting in a lower energy, whereas another system may never be able to form a perfect crystal because its periodicity is interrupted by the periodic boundary of the simulation volume.
To quantify the duration of the exponential decrease, we conduct a least-squares fit of the per-particle Lennard-Jones energy to \autoref{eq:fit}, which is also used to compare the 15mer and 50mer systems.
We consider a linear term in addition to an exponential term to better account for the slow transformation in Region 2, which is qualitatively linear.
We calculate the derivate of \autoref{eq:fit} to obtain the exponential and linear components of the slope, and identify the transition time $t_{s}$ as the point at which the linear component becomes a larger than the exponential component.

To determine when equilibrium is reached and Region 2 is ended, we must first assume that the final energy measurement is in an equilibrium state (which is only a valid assumption for the systems with fewer beads, $<$ 325 15mers).
We next assign the measured energies to bins containing 200 measurements each, which equates to 36 ns, and calculate the average $\bar{E}_{i}$ and standard deviation $\sigma_{i}$ of each bin. 
We then iterate through the bins in reverse order and compare $\bar{E}_{i}$ and $\sigma_{i}$ the final bin's energy $\bar{E}_{final}$ and deviation $\sigma_{final}$.
When $\bar{E}_{i}$ is calculated to be more than the sum of the standard deviations higher than $\bar{E}_{final}$, we consider the system as not equilibrated, i.e.:

\begin{equation}
\text{System is }
\begin{cases}
\text{Not Equilibrated}, & \text{ if } \bar{E}_{i} - \bar{E}_{final} > \sigma_{i} + \sigma_{final} \\
\text{Equilibrated}, & \text{otherwise}.
\end{cases}
\end{equation}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newpage

\section{Small and Large Comparison for Semiordered and Disordered Systems}

In the main text we consider size effect differences between the ordered large and small systems, here we expand this analysis to large and small semiordered and disordered systems (\autoref{fig:semi-dis}).
In the large semiordered system, some crystallites have aggregated to form alkyl-stacked lamellae (shown with colored, opaque beads), corresponding to the experimentally observed crystalline domains surrounded by an amorphous matrix (shown with gray, diffuse beads) \cite{Brinkmann2006}.
However, in the small system, one large crystallite interspersed with smaller crystallites dominates the morphology, and the alkyl-stacked lamellae are not observed in the smaller system.
The difference between these two structures is confirmed with GIXS patterns: periodic alkyl stackings are only observed in the large system. 
That said, the $\pi$-stacking feature along the (010) planes forms a diffuse band rather than a distinct peak - likely due to the wide range of crystallite orientations in the large morphology in the large system, whereas this feature is much stronger in the smaller system's pattern.
Therefore, there is some dependence on the system size when looking at semi-ordered systems: larger systems are able to have regions of high order and amorphous regions, which are missing in smaller systems that tend to form fewer, larger crystals.
This corresponds to crystallite sizes measured experimentally ($\sim$ 10-20 nm), in that the small system with a box length of $\sim$ 7 nm is smaller than the typical crystallite, whereas the larger systems with length 15 nm are the same as the average experimentally observed crystals \cite{Liu2014a}.
We therefore suggest using larger system sizes ($N=165,000$) to more accurately represent the long-range ordering characteristics of polymers.

In the case of disordered systems, both the 100 and 1000 chain systems are similar: they are comprised of many small, randomly-oriented regions of $\pi$-stacking, consisting of just a couple of chains in each case.
The GIXS patterns corroborate the visual observations: there is little periodic order within the system. 
As such there is minimal system size dependence when the structures are more disordered; smaller simulations can be used to accurately investigate high temperatures.
In summary, smaller system sizes are beneficial for fast investigations of the relationship between state-point and structure, whereas larger systems including more molecules are be needed to explore large-scale structural features.

\begin{figure}[h]
\centering
\includegraphics[width=6in]{{semi-dis-comb}.pdf}
\caption{
%TODO More to caption?
In cases where amorphous and structured regions exist, large volumes (a) provide more insight than their smaller counterparts (b) in that systems containing 100 15mers do not show periodic lamellae features seen in larger systems (c, d).
Conversely in disordered systems, large volumes do not provide significant structural insight as compared to small volumes (e-h).
}
\label{fig:semi-dis}
\end{figure}

\newpage

\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).
We would like to acknowledge high-performance computing support of the R2 compute cluster (DOI: 10.18122/B2S41H) provided by Boise State University's Research Computing Department.
}

\authorcontributions{Under the direction of Dr. Eric Jankowski, Evan Miller and Dr. Matthew Jones optimized the P3HT force-field and performed the molecular dynamic simulations. Michael Henry evaluated the MD performance on various GPUs. The code to organize and run the molecular dynamic simulations was written by Dr. Jankowski.}
 
\conflictsofinterest{The authors declare no conflict of interest.}

\abbreviations{

15mer - P3HT chain containing 15 monomers

50mer - P3HT chain containing 50 monomers

AA - All-Atom

CA - Aromatic Carbon

CT - Aliphatic Carbon

$\epsilon$ - Energy Parameter

$\varepsilon_s$ - solvent quality

GIXS - Grazing Incident X-ray Scattering

GPU - Graphical Processing Unit

HPC - High Performance Computing

$m$ - Mass

MD - Molecular Dynamics

OPLS - Optimized Performance for Liquid Simulations

OPV - Organic Photovoltaic

P3HT - Poly(3-hexylthiophene)

$\psi$ - order parameter

$\rho$ - density

RDF - Radial Distribution Function

S - Sulfur

$\sigma_{LJ}$ - Lennard-Jones van der Waals radius

$t$ - Simulation Unit of Time

$T$ - Temperature

TPS - Timesteps Per Second

UA - United-Atom

}

\bibliography{library}

\end{document}
