Force Fields for Macromolecular Assemblies Containing Diketopyrrolopyrrole and Thiophene.

Utilising a force matching procedure, we parameterise new force fields systematically for large conjugated systems. We model both conjugated polymers and molecular crystals that contain diketopyrrolopyrrole, thiophene, and thieno[3,2-b]thiophene units. These systems have recently been found to have low band gaps, which exhibit high efficiency for photovoltaic devices. The equilibrium structures, forces and energies of the building block chromophores: diketopyrrolopyrrole thiophene, and thieno[3,2- b]thiophene computed using our parameters are comparable to those computed using the reference electronic structure method. We assess the suitability of this new force field for electronic property calculations by comparing the electronic excitation properties computed along classical and ab initio molecular dynamics trajectories. For both trajectories, we find similar distributions of TDDFT calculated excitation energies and oscillator strengths for the building block chromophore diketopyrrolopyrrole-thieno[3,2- b]thiophene. The structural, dynamical, and electronic properties of the macromolecular assemblies built upon these chromophores are characterised. For both polymers and molecular crystals, pronounced peaks around 0 degree or 180 degree are observed for the torsions between chromophores under ambient conditions. The high planarity in these systems can promote local ordering and π-π stacking, thereby potentially facilitating charge transport across these materials. For the model conducting polymers, we found that the fluctuations in the density of states per chain per monomer is negligibly small and does not vary significantly with chains comprising 20 to 40 monomers. Analysis of the electron-hole distributions and the transition density matrices indicates that the delocalised length is approximately four to six monomers, which is in good agreement with other theoretical and experimental studies of different conducting polymers. For the molecular crystals, our investigation of the characteristic timescale of the fluctuation in the excitonic couplings shows that low frequency vibration below 10 cm-1 is observed for the nearest neighbours. These observations are in line with previous studies on other molecular crystals, in which low frequency vibrations are believed to be responsible for the large modulation of the excitonic coupling. Thus, our approach and the new force fields provide a direct route for studying the structure-property relations and the molecular level origins of the high-efficiency of these classes of materials.


Introduction
Semiconductors based on small organic molecules and conjugated polymers offer the opportunity for large-area deposition on substrates and they can be chemically tuned or functionalized for many specific applications. 1 In the last two decades, they have emerged as excellent materials for displays, solid-state lighting, sensors and solar cells. [2][3][4][5] Their main advantages are the significant reduction in production cost compared to their organic counterparts and their processes can be easily scaled-up, since organic materials can be solution-processed.
Conjugated organic molecules exhibit a wealth of structures and properties, due to the range of possible molecular constituents. They can be completely disordered polymers formed by interlaced chains, partially ordered polymers, polycrystalline small molecule-based structures, or small molecule-based single crystals. 1,6 The detailed characterisation of the structure of conjugated organic assemblies remains a challenging task, often requiring simplified and qualitative morphological descriptions. In addition, the variety of these materials, their possible structures, and their complexity hinder the full elucidation of the relationship between the structures of conjugated organic assemblies and their electronic properties. 6,7 To investigate the fundamental processes that determine the efficiency of charge transport in organic semiconductors, many efforts have been devoted to understand their structureproperty relationships, [8][9][10][11] including the influence of molecular morphologies on optical and charge transfer events. However, the interplay between microstructure and electronic properties in these materials still remains unclear. 12,13 Amorphous materials or poly-crystalline materials especially are poorly understood. 6 In particular, a molecular level comprehension of the exciton transport and charge generation mechanisms is crucial, for example, to the rational design and enhancement of conjugated organic semiconductors for highly efficient photovoltaic devices. Theoretical and computational approaches can play a major role in the study of this important aspect, because they provide direct access to the microscopic and electronic properties of materials.
Molecular simulations offer a direct route to explore the relationship between molecular morphology and charge mobility. Amongst several morphology simulation techniques, 14 classical molecular dynamics (MD) is the most popular approach, as it gives atomic detail comparable with experimental observations. In the last decade MD simulations have penetrated into the study of optical and electronic properties of materials and biomolecules, where they are primarily used to sample the equilibrium structures explored by the system of interest. 7,[15][16][17] These equilibrium trajectories are utilised in large-scale quantum chemistry calculations to correlate the local structure with the observable electronic structure properties. 18 Consequently, it is vital that the force fields (FFs) used to describe molecular motion of the system are accurate enough to reproduce experiments. The main concern here is the mismatch in the FFs describing the interactions at the classical level and the electronic structure methods employed for subsequent quantum chemistry calculations.
Historically, FFs employed in classical MD simulations have been designed mainly for biochemical systems, e.g. protein simulations. Although these FFs contain parameters that are transferable to conjugated systems, many aspects require a careful parameterisation due to the conjugation. 19 For large conjugated systems such as polymers, FFs need to employ many atom types and parameters. Generally, the FFs used to describe the motions of molecular systems are seldom parameterised to reproduce electronic processes, such as electronic excitations, which are sensitive to the underlying description of the molecular structure. It still remains a challenge to find satisfactorily accurate FFs for applications where MD is the preliminary step toward the study of electronic structure properties. [19][20][21] Most existing procedures [22][23][24][25] are general tools for constructing FFs for small organic molecules, either in the gas phase or in a continuum solvent. When applied to specific problems in condensed phases, these FFs may not be directly transferable to the systems under investigation. For example, several recent studies of excitation energy transfer in light-harvesting (LH) complexes have shown that FFs used to describe the chromophores have a great influence on the resulting spectral density. [26][27][28] This problem is often referred to as the geometry mismatch, i.e. the mismatch between the classical and ab initio potential energy surfaces. As such, Prandi et al. proposed a new set of parameters for the different carotenoids found in the LH complexes of photosystem II, optimised to reproduce electronic excitation properties obtained by DFT. 29 To investigate this further, Andreussi et al. 20 systematically checked the feasibility of using classical FFs for quantum mechanics applications. They concluded that for photoinduced processes as an example, poor results are obtained from standard transferable FFs. Thus, in order to mitigate the mismatch in the geometry, the classical potential energy surface and/or its gradients of the system of interest should be mapped onto the ab initio ones by optimising the FFs. This would also implicitly incorporate many-body effects, which is crucial in the study of condensed phase systems. To address this, Ercolessi and Adams proposed the force-matching (FM) approach, 30 which was first employed to develop interatomic potentials that reproduce those computed from first principles. To overcome the computational cost incurred by using ab initio MD or QM/MM, Wang and coworkers introduced the adaptive force matching technique. 31,32 Instead of using ab initio MD, the adaptive force matching approach uses MM FFs to generate ensembles of equilibrium structures followed by decoupled QM/MM calculations on each configuration with the environment represented by point charges. The procedure starts with a set of guessed parameters and is repeated until convergence is reached. Since then, the method has been applied to systems including transition metal complexes, 33,34 anions, 35 ionic liquids, 36 and microporous materials. 37 Recently, the force-matching technique has also been employed to develop FFs for water, 38,39 water on graphene surface, 40 and explicit three-body interactions for molten carbon. 41 In this work, we extend the adaptive force matching technique to rapidly parameterize FFs for complex conjugated systems. We develop a strategy for fragmenting the "larger" molecule into smaller segments, obtain their parameters, and then assemble them together to derive the FFs for the "larger" molecule. This was a challenge, which was presented, but not achieved, in the previous work by Do and Troisi. 42 Do and Troisi 42 reported parameters for the monomer DPP and for the chromophore-triple DPP-TT-T, but for no other species, precluding the characterisation of the dynamics of the polymer and crystal systems that we study herein. Moreover, the suitability of the force field for electronic property calculations by comparing the electronic excitation properties computed along classical and ab initio molecular dynamics trajectories was not assessed. This is important and should be tested, at least, for the monomer. In the current study, we build on and go beyond the earlier work 42 to address three important aspects. (i) In the current work, we fragment the "larger" molecule into smaller segments, obtain their parameters, and then assemble them together to derive the FF for the"larger" molecule. This was a challenge, which was presented, but not achieved for the chromophore-quadruple (the monomer of the model polymer) in the previous work. 42 (ii) We validate the FF with the first MD simulation of the model polymer (for a single chain 1-mer) by comparing against counterpart ab initio MD simulations. (iii) We then predict the structural and the electronic excitation properties for the model polymer and also for a newly synthesised molecular crystal. The thrust of our approach is to obtain parameters consistent with any pre-defined electronic structure calculation method.
As an example, we parametrize FFs for conducting polymers and molecular crystals containing diketopyrrolopyrrole (DPP), thiophene (T), and thieno[3,2-b]thiophene (TT). The DPP units have strong intermolecular donor-acceptor interactions, which can promote selfassembly of the polymer and enhance the charge transport. 43 In addition, the incorporation of DPP into an oligothiophene or polythiophene backbone results in low band gap conjugated systems exhibiting high efficiency for photovoltaic devices. 44 Thus, FFs that generate accurate equilibrium structures via MD simulations of these systems are crucial for analyzing and understanding the structure-property relationships of conjugated semiconductors.

Methods
The idea behind the force-matching procedure is to find a set of parameters that minimizes the differences between the classical forces and the ab initio forces. Here, the latter can be computed using any reference electronic structure method. Generally, this is chosen to be consistent with the method that is employed in subsequent quantum chemistry calculations.
In addition, an appropriate FF that best represents the system of interest has to be determined in advance. The errors between the two types of forces are computed from a set of uncorrelated configurations generated from MD simulations during the production runs. In the approach we use, 42 a set of parameters {p j } is optimized by minimizing an objective function O({p j }), which is defined as the root mean square deviation (RMSD) of two types of forces: where N is the number of atoms in the molecule and M is the number of configurations

Computational Details
To start with, a set of M equilibrium configurations is generated using an initially "guessed" FF. The atom-type is carefully assigned to take into account the complex chemical structure of the system. As such, we tend to have more atom-types than standard FFs. A description of how the different atom types are chosen can be found in the SI. Apart from the point charges and equilibrium bonds and angles, all other "guessed" parameters are taken from OPLS. 24,[45][46][47] For missing parameters, we use analogy and quantum chemistry calculations to assign the parameters (e.g. equilibrium bonds and angles). Equilibrium structures are obtained from MD simulations under ambient conditions (i.e. 300 K and 1 bar). All MD simulations are carried out using the NAMD software. 48 After an energy minimization of 10000 steps, a 200 ps equilibration run is carried out, followed by a 1 ns production run, The algorithm proceeds by randomly selecting a parameter p j from the set and modifies it by adding a random number δp j uniformly distributed between −δp max j and δp max j . In other words, δp max j is the maximum "displacement" of the parameter δp j . This value is initialized to 0.3 × p j at the beginning of the optimization , and is adjusted in the course of the process to keep the acceptance rate of the MC moves at about 30%. The optimization is divided into MC "blocks". Each consists of 100 × N p attempts to randomly select and change a parameter p j . The idea here is to allow each parameter to be selected roughly 100 times in each MC block. This number can be adjusted and we found 100 times is sufficient. The acceptance rate is calculated at the end of each block and the maximum "displacement" δp max j is adjusted accordingly. The optimization is converged when δO < 10 −10 kcal mol where δO is the change in the objective function between two consecutive MC steps. Having The overlap function S r (r) between a hole and an electron can be defined as S r (r) = ρ hole (r)ρ elec (r), where ρ hole and ρ elec are the charge densities of the hole and electron, respectively. Thus, the S r -index is the integration of the S r (r) function over all space, which provides the extent of the overlap between hole and electron. The charge transfer (CT) length can be estimated by computing the distance between the centroid of hole and electron. To compute the X component of the centroid of the electron, for example, we integrate the x-coordinate of the position vector r weighted by the electron density i.e.
X elec = xρ elec (r)dr. The total magnitude of CT length is referred to as the D-index and can be computed as D-index = D 2

Force Field Development
Recently, a series of interesting thienothiophene-flanked diketopyrrolopyrrole (DPPTT)based copolymers with various branched alkyl side-chains have been synthesized. 5 In this section we shall describe the procedure to construct the FFs for the above mentioned conducting polymers and molecular crystals that contain diketopyrrolopyrrole (DPP), thiophene (T) and thieno[3,2-b]thiophene (TT) chromophores. Instead of optimising the parameters for the whole system in one go as in the work of Do and Troisi, 42 our strategy is to fragment the "larger" molecule into smaller segments, obtain their parameters, and then assemble them together to derive the FFs for the "larger" molecule. The idea is to parameterize the FFs for each single chromophore (e.g. DPP, TT and T) and utilize them to construct the FFs for larger segments, e.g. segments that contain two or more chromophores. This can be achieved by running the optimization procedure for the "larger" molecules and only allow the missing parameters to be optimized while keeping the rest fixed at those obtained from the single chromophores. In practice, we perform a further optimization for the whole "larger" molecules (including sidechains if there is any) to obtain better sets of parameters. These calculations are fairly quick, as the original "guessed" parameters are near optimum. Thus, our fragmentation approach is particularly suitable for large systems with many parameters, for which the previous technique 42 is less well suited.  TC5   TC6   TC6   TS1   TS1   TC5   TC4 TC2   TC3   TO1   TO1   TC1   TC3   TC1   TC2   TC5   TC4  TC6   TN1   TS1   TC6   TS1   TC5   TC4 TC4   TC8   TS1   TC8   TC7   TH0   TH1   TH1   THN   THN   TH1   TH1   TH1   TH0  The procedure for obtaining the parameters for the chromophore-quadruple is similar to that of the chromophore-pairs. In the first stage, all parameters apart from those of the "soft" dihedrals are obtained using the planar structures. In these geometries, the forces of the "soft" dihedrals are essentially zero. The initial "guessed" parameters are taken from the single chromophores. Missing parameters due to the combination of two chromophores (e.g. the bond parameters for C6-C7 or the angle parameters for C6-C7-S1) are taken from the chromophore-pairs. In the optimization process, it is possible to allow only the missing parameters to be optimized while the rest remain fixed at those obtained from the parameterization of the single chromophore. This approach is particular suitable for very large molecules with many parameters that require a long optimization time. For the case of the chromophore-quadruple (and also the chromophore-pair), it is possible to carry out the optimization processes for the whole molecule, and since the initial "guessed" parameters are reasonably good, it only required two iterations to converge. For chromophore-quadruple (and larger molecules e.g. oligomers), the "soft" dihedrals are not parameterized individually, and their parameters are taken directly from those in the chromophore-pairs. Figure 2 shows that the inclusion of the T ring slightly alters the torsion potential of the DPP-TT torsion.

Chromophore-quadruple
However, this effect is small enough to be safely neglected.

Results and Discussion
Validation of the Parameters for the Chromophore-quadruple A comparison of the force and energy RMSD is shown in Table 1. Similar to the cases for single chromophores and chromophore-dimers (See tables S1 and S2 in the SI), AM1 gives the largest deviation for forces followed by HF/6-31G * . B3LYP/3-21G* gives the smallest deviation from the reference method B3LYP/6-31G * for both forces and energies. The force RMSD is comparable to that of B3LYP/3-21G * , while the energy RMSD is comparable to HF/6-31G * .
To demonstrate that our parameters are reliable for generating structures for quantum chem- Table 1: Force and energy RMSD between calculations at B3LYP/6-31G * levels and the optimized FF or other levels of theory, for 100 arbitrary geometries of the chromophorequadruple TT-DPP-TT-T not used in the FF optimization process.

Method
Force RMSD (kcal mol −1 Å −1 ) Energy RMSD (kcal mol −1 ) Finally, Figure 4 compares the structural properties of the chromophore-quadruple computed for 100 snapshots along the classical and ab initio trajectories. As an example, we only show distributions for a bond, an angle and a soft dihedral of the chromophore in this figure.
Other bonds, angles and dihedrals show the same trend. The good agreement observed for the structural properties in Figure 4 for the chromophore-quadruple underpin the agreement of the electronic excitation properties observed in Figure 3. Hence, this provides further confidence that our FFs can reproduce the ab initio potential energy surfaces, to which their parameters are fitted.  Figure 4: Comparison of the structural properties of chromophore-quadruple computed using 100 uncorrelated frames taken from the classical (black) and the ab initio (red) trajectories: (a) bond between C6 and C7, (b) angle between C3, C6 and C7, and (c) dihedral between C3, C6, C7 and S1.

Molecular Dynamics Simulations
MD simulations are carried out with NAMD 48 for 20 40-mers of the DPPTBT bulk conducting polymer. The monomer of DPPTBT is the chromophore-quadruple described in the previous section. Thus, the topology and parameters for DPPTBT oligomers are built upon the chromophore-quadruple. The polymer is initially placed into a cubic box of 109×109×109 Å 3 at a density of 0.5 g/cm 3 without overlaps. To achieve this, polymer chains are added into the box by growing segment by segment using the configurational bias Monte Carlo technique 63 To avoid steric clashes, the growing process takes into account interaction with all atoms already positioned, whilst continuously monitoring the single chain conformations.
The result of this procedure is that low energy sites are preferred over high energy sites.
Thus, bulk disordered systems containing chain molecules in realistic equilibrium conformations are created.
Three-dimensional periodic boundary conditions are implemented. The system is minimized for 5000 steps and gradually heated to 300 K using 50 ps in the NVT ensemble. Then, the NPT ensemble at 1 atm and 300 K is employed for 200 ns to equilibrate the system. Finally, a further 100 ns is simulated in the NVT ensemble for calculating ensemble average properties. The Langevin thermostat 64 with a damping coefficient equal to 1 ps −1 is employed to keep the temperature constant. The period and decay parameters of the Langevin piston are set to 100 and 50, respectively, to maintain the pressure. A time step of 1 fs for integration of the equations of motion is used throughout the simulation. A cutoff of 12 Å is used for nonbonded interactions. The particle mesh Ewald algorithm 65 is used to calculate long-range electrostatic interactions. Figure 5 gives a visualization of the system and Figure   6 shows the density variation with respect to the simulation time. Starting from 0.5 g/cm 3 , the system reaches about 1 g/cm 3 after 100 ns, and eventually converges to about 1.035 g/cm 3 . A density of around 1 g/cm 3 is a typical value for conducting polymers. The same final density is also observed for different starting densities and different chain lengths.
First of all, we inspect the radius of gyration of the polymer under ambient conditions. The range of gyration radii ( Figure S28) is between 125 to 160 Å, which indicates these  In general, the flexible dihedral angles play an important role in the electronic properties of conducting polymers. A planar configuration leads to high hole mobility, due to the formation of the π orbitals along the backbone. 66 To examine the planarity of this polymer, we compute the distribution of the soft dihedrals along the chain backbone. Figure 7 shows the distribution between DPP and TT of the first monomer in chain one. A strong peak near 0 • is observed. This distribution is also found for the torsion DPP-TT along chain number one and in other chains. Similarly, for the distribution of the torsions between TT and T and between DPP and T, we also observe peaks around 180 • or 0 • . This suggests that dihedral angles near 0 • and 180 • are favoured thermodynamically, which indicates planar configurations for DPPTBT at the temperature and density studied.  Figure 7: Distribution of dihedrals of polymer during 100 ns NVT simulation. Figure 8 shows the RDFs of all N1, S1 and O1 atoms in the system, respectively. These RDFs give quantitative information about the separation (or packing) of surrounding atoms from any given atom. There are three distinctive peaks observed for all three atom types at roughly 4 , 8 and 20 Å. These are the radii for which other atoms most likely to be found.
The first two peaks in the RDF for sulfur atoms (Figure 8b) are similar to that observed for the one chain simulations of the PCDTBT conducting polymer, where sulfur peaks at about 5 and 9 Å were observed. 66 Thus, the first two peaks at short distances in our RDFs can be associated with the distributions of atoms along the chain (i.e. intra-chain packing) and the last peak (at about 20 Å) gives the inter-chain distributions.

Electronic Structure Properties
To gain insights into the electronic structure of the polymer, we first study the bulk density of states (DOS). Figure 9 shows the bulk DOS for the system of 20 chains of 40-mers averaged over ten snapshots from 10 to 100 ns each separated by 10 ns. Here, the background charges are not included to save the computational time, due to the effect of electrostatic disorder, which is relatively weak for semiconducting polymer. 7 The long-range corrected functional ωB97X-D is employed along with the 3-21G * basis set for this task. We found that the fluctuations in the DOS per chain per monomer for this system are almost negligible. Similarly, the DOS does not vary significantly with the chain lengths between 20 to 40 monomers, which agrees with a similar study on amorphous MEH-PPV polymers. 7 Calculations (only for one snapshot) using a larger basis set, i.e. 6-31G * , indicate that the discrepancy in the computed DOS between the two basis sets is negligible, especially around the HOMO-LUMO gap. Thus, it is sufficient to employ the smaller basis set (3-21G * ) for characterising the electronic excitation in this type of conducting polymer using TDDFT.  Figure 9: Bulk density of states computed for the 40-mer chains using the long-range corrected functional ωB97X-D. Solid red line is the result using 3-21G * basis set averaged over 10 snapshots from 10 to 100 ns separated by 10 ns. Dashed black line is the result using 6-31G * basis set for a single chain at 100 ns.
To investigate the nature of the electronic excitation of the polymer, we analyse the e-h distributions and the transition densities for the lowest three singlet excitations of a single chain. Among several other descriptors, the Λ-index proposed by Peach et al. 54 has gained much attention, and has been widely employed as a diagnostic tool for monitoring TDDFT results. 67 This method is based on the overlap of molecular orbital moduli. However, in some difficult cases, this index may not be able to locate problematic excitations. 68 To provide a robust tool for exploring the excited state metric, another index called the ∆r-index has been introduced. 55 This index is based on the charge centroids of the orbitals involved in the transition. It is easy to interpret, because it measures the average e-h distance upon excitation. The ∆r-index works particularly well for the CT excitations for which reliable charge transfer length can be obtained.
To characterise the electronic excitation of the conducting polymer, we computed all the descriptors discussed previously and we present the results for the three lowest singlet excitations for the 40-mer chain in Table 2. All TDDFT calculations are performed at the ωB97X-D/3-21G * level of theory. Table 2 shows that the S 0 →S 1 is clearly a local excitation with a relatively small D-index and ∆r-index, which indicate that the hole and electron are quite close. This is also confirmed by, for example, the high value in the exciton binding energy and the overlap S r -index. Further investigation of this transition using the contributions of the molecular orbitals of atoms and chromophores in the polymer chain to the e-h ( Figure 10) shows that both the hole and electron originate from monomer number 30 and the e-h also resides on this monomer (mainly on the DPP chromophore). This is confirmed by a plot of the fragment transition density matrix (Figure 11), which shows the distribution of charges along the polymer chain. According to the heat map shown in Figure 11, we can also observe that the electron and hole are mainly distributed on monomer 30. Since there is no clear off-diagonal element, this transition does not cause significant electron transfer between monomers. Although there are some small fractions of charge on the nearest neighbours, as can be observed from Figures 10 and 11, the main character of this excitation is its localization on monomer 30.

Molecular Dynamics Simulations
As another example, we study the dynamics of the molecular crystal UBEQOK, which has been synthesised recently ( Figure S34). 59 The equilibrium structures generated from MD simulations will be valuable for understanding the exciton dynamics in this new material.
First of all, a 3 × 8 × 2 supercell, which comprises 192 molecules is created. Then, the system is minimized using 5000 steepest descent steps followed by 0.5 ns heating in the NPT ensemble (at 300 K and 1 atm). It is equilibriated for another 100 ns and continued for another 10 ns for production. All other simulation parameters are kept the same as those in the simulations of the polymeric system.
Since the simulations are performed under ambient conditions, it is expected that the molecules do not move very far from the original crystal structure and this can be observed from the root mean squared displacement of all molecules in the system (see Figure S27  Here, only the couplings of nearest neighbors are considered, as they have been shown to be adequate to describe quantitatively the excited state properties of molecular crystals. 71 In addition, the coupling between molecular pairs at longer distances is dominated by the electrostatic term, which is not sensitive to small thermal fluctuations. 70 Two types of dimers with closest intermolecular contact for the exciton transport in UBQOK molecular crystal are identified. Figure 13 shows the arrangement of these two   Figure 14 shows that the magnitude of the couplings for UBEQOK falls into the same range as those computed for other molecular crystals e.g. DCVSN5 72 and H 2 − OBP c . 73 In addition, the largest excitonic coupling is also observed for the parallel structure for this system.

Conclusion and Outlook
A robust approach for generating accurate FFs for conjugated systems has been applied and extended for applications where subsequent quantum chemical calculations are employed.
Built upon the force-matching technique, the method provides a consistency between the potential energy surface generated by empirical FFs and that generated by any ab initio method of choice. The approach is particularly suitable for applications that involve excited state calculations where accurate FFs to describe the excited state potential energy surface are needed. As an example, we parametrize the FFs and study the dynamics as well as the electronic excited properties of conjugated conducting polymers and molecular crystals containing diketopyrrolopyrrole, thiophene, and thieno [3,2-b]thiophene. The good agreement observed for both the structural and electronic excitation properties for all building block chromophores provides confidence that our FFs can reproduce the ab initio potential energy surfaces, for which their parameters are fitted to.
In this work, we employ OPLS FFs and the DFT at the B3LYP/6-31G * level of theory to illustrate how the desired optimised parameters can be obtained. However, the procedure is general and can be applied to any other combination of FF functional form and electronic structure method. The force RMSD between ab initio and empirical FFs can be improved by including more fitting parameters, such as introducing more atom types, including coupling between degrees of freedom and anharmonicity, or modifying the FF functional form. Our optimised FFs, by construction, are not transferable by atom types. However, the conjugated core such as DPP, TT, T, DPP-TT, DPP-T and TT-T could be treated as transferable. FFs for larger molecules can be constructed from these basic units, and thus allow access to the dynamics of conjugated macro-assemblies.

Supporting Information
A brief description of the OPLS force field; discussion on choosing different atom types; force field development for single chromophores and chromophore pairs; additional properties of the conducting polymer and molecular crystal; force field parameters. This information is available free of charge via the Internet at http://pubs.acs.org