Conjugated polymers are an important class of electronic materials for organic photovoltaics (OPVs), light emitting diodes (OLEDs), field effect transistors (OFETs), and bioelectronics. [1][2][3][4] These materials are advantageous for their lower cost, light weight, and flexibility compared to the traditionally used inorganic materials. [1][2][3][4] The typical conjugated polymer structure includes a conjugated backbone (e.g. polymerized thiophene rings) with added alkyl side chains for increased solubility. This allows for charge transport either along the backbone (i.e. intra-chain charge transfer) or through the pi-orbital overlap of neighboring backbones (i.e. inter-chain charge transfer). It is known that the molecular morphology of conjugated polymers is connected to their macroscopic performance, and a thorough understanding and control of that morphology is crucial to the development of improved organic electronics. [5][6][7] Moreover, polymer dynamics have been explored for their effects on charge transport. 5,8 For example, Alberga and coworkers studied the effects of morphology fluctuation on charge mobility in amorphous and crystalline phases of poly(3-hexylthiophene) (P3HT) and poly(2,5-bis(3-alkylthiophen-2-yl)thieno [2,3-b]thiophene) (PBTTT). Their respective structures changed the chain rigidity and motions along the backbone at time and length scales relevant to inter-and intra-chain charge transport mechanisms. 5 Moreover, these polymer fluctuations occur at similar time scales (picosecond to nanosecond range) as the exciton thermalization and recombination in organic bulk heterojunction solar cells 9 and become increasingly relevant at anticipated device operating temperatures, which are usually well above the glass transition temperature for these materials. 8 The combination of quantum mechanics calculations with either molecular dynamics (MD) or Monte Carlo (MC) simulations has been one strategy used to study the effect of dynamic disorder on charge transport for these materials. 10 MD or MC simulations are used to first equilibrate a conjugated polymer system. Then charge transport properties are calculated using quantum mechanical (ab initio) methods on snapshots of the MD/MC systems. McMahon and coworkers studied the effects of dynamic disorder on crystalline P3HT, equilibrating the system via an MD simulation at 300K and 1 atm. 11 They noted that charge traps formed in regions of high backbone planarity with lifetimes on the order of ten picoseconds in crystalline P3HT. They also determined that the addition of defects in regioregularity to the crystalline system had no influence on the location or depth of the charge traps. Vukmirovic and coworkers also studied crystalline P3HT, using MC simulations to separately examine the effects of side chain and backbone disorder on the electronic density of states (DOS). 12,13 They concluded that both side chain and backbone disorder contributed to broadening of the DOS near the valence band edge and that backbone disorder played a larger role in increasing charge localization. For the case of amorphous P3HT, they found that the primary contributions to charge localization come from the electrostatic potentials of the backbone units, influencing their relative position and orientation, which in turn affects inter-chain pi-orbital overlap. 12 Poelking and coworkers also investigated the effect of regio-regularity on charge transport in P3HT using this approach. 10,14 They calculated time autocorrelation functions for two parameters that affect charge transport according to the Marcus theory of charge transport, which is related to charge transfer integrals and hopping site energies. They determined that the time scales for the decorrelation of these parameters was in the picosecond to nanosecond range for fully regio-regular P3HT (RR-P3HT).
Counterintuitively, their results indicated that decorrelation time scales for the 90% regio-regular system were longer than those of the fully regio-regular system. 14 More recently, Jones et al.
used combined coarse-grained molecular dynamics, quantum calculations and Monte Carlo techniques to evaluate the effects of molecular weight and annealing on the morphology and charge mobility of P3HT, making qualitative comparisons with experimental work in the literature. An increase in annealing temperature was determined to improve charge mobility due to a higher molecular ordering of the system. Higher molecular weights also increased charge mobility, but this was attributed to more chains connecting the crystalline regimes through amorphous domains. 15 Because all of these aforementioned methods use snapshots of trajectories obtained from MD and MC simulations, their ability to describe charge transport that is representative of real bulk materials depends on the accuracy of the simulation parameters.
However, there is limited agreement in the literature between the methods used for reparameterizing MD force fields (FFs) for conjugated polymers. This work aims to address these concerns by utilizing wide-angle X-ray scattering (WAXS) and quasi-elastic neutron scattering (QENS) techniques to critically evaluate the results of MD simulations for amorphous P3HT.
QENS experiments provide quantitative information about the dynamics of materials at molecular length scales. In comparison to other polymer systems, relatively few QENS experiments of conjugated polymers have been performed. Obrzut and Page related the dynamics of RR-P3HT to its conductivity by comparing QENS data with dielectric spectroscopy data. 16 They observed a "liberation" of the side chains at ~175 K, with an activation energy for these motions of ~9 kJ/mol and with a corresponding change in the conductivity of bulk P3HT.
Paterno and coworkers observed that the dynamics of P3HT-PCBM films were unaffected by the solvent used to make the film. 17 They additionally determined that increasing the PCBM concentration in the film increased the time-scales for the observed dynamics. They hypothesized that the longer time scales resulted from increased P3HT side chain confinement caused by steric interactions with PCBM. Etampawala and coworkers also performed QENS measurements of P3HT/PCBM solvent-cast films. 18 They verified the results of Paterno and coworkers, and further examined the Q dependence of the P3HT side chain dynamics. For length scales greater than 7 Å, they observed dynamics dominated by simple diffusion, and attributed these dynamics to collective motions of the side chain. The dynamics corresponding to length scales smaller than 6 Å were attributed to motions within the side chain (i.e. methyl rotations). 18 Guilbert and coworkers used time-of-flight QENS to study the dynamics of P3HT and PCBM blends above and below the glass transition temperature of P3HT. At higher temperatures, blending of the two materials caused a dampening of P3HT dynamics while enhancing those of PCBM. Moreover, the authors observed alkyl side chain motions of both polymers in the picosecond time range that are expected to impact the macroscopic performance of organic electronic devices. 19 QENS measurements and MD simulations are complementary techniques. QENS provides information about the dynamics of specific fragments of the experimental system (e.g. selfcorrelation of motions for hydrogen atoms). On the other hand, MD simulations provide detailed information about all aspects of the system, but this information is only meaningful if the simulations are experimentally validated and determined to be accurate. QENS also provides a powerful method for validating MD simulations since the measured scattering from an experiment can be directly compared to simulated scattering data calculated directly from the MD trajectories. Traditional methods of MD simulation validation usually involve the use of structural and thermodynamic properties such as heat capacity, surface tension, glass transition temperature, melting temperature or mass density to validate FF parameters. However, these benchmarks are calculated by time averaging and may not provide accurate information about localized molecular fluctuations. In contrast, simulations validated with QENS experiments provide ample quantitative information about the dynamics of the system at various length scales and time scales that correspond to specific molecular motions. Richter, Colmenero, and Arbe et al. proposed the use of QENS to directly validate MD simulations of polymers to obtain a clear understanding of polymer dynamics. 20,21 For example, they were able to validate MD simulations of poly(alkyl oxides) (PAOs) that used the COMPASS FF. 22,23 After validating the simulation, they could identify the contributions from different methyl groups to the overall dynamics and determine characteristic relaxation times for backbone and side chain motions separately. They have conducted similar analyses on poly(n-alkyl methacrylates), which showed reduced localized motions due to increased side chain confinement relative to PAOs. 24 Only a handful of combined QENS-MD studies have been conducted on conjugated polymers.
Sniechowski and coworkers adopted this method when investigating the dynamics of polyaniline (PANI) doped with the plasticizing dopant di-(2-butoxyethoxyethyl) ester of 4-sulfophthalic acid (DB3EPSA). [25][26][27] They investigated the change in the mean-squared displacement of protons in the dopant as a function of time and found that the time scales for relaxations of different dopant protons spanned several orders of magnitude. They observed that the protons closest to the PANI chains had the longest relaxation times, confirming that the polymer chains were effectively stationary with respect to the dopant. Guilbert and coworkers adopted the method used by Colmenero for RR-P3HT and RR-poly(3-octylthiophene). 28 They were able to demonstrate that QENS could be used to validate MD simulations of these materials and determined activation energies for different motions in these complex systems. Guilbert and coworkers also used MD to further understand their previously discussed work using QENS to capture dynamics in P3HT/PCBM blends. 19,29 They explained the dampening of P3HT motions and enhancement of PCBM motions by a wrapping of the P3HT chains around the PCBM molecules in the amorphous blend domains. Moreover, they explored a larger dynamics range in-silico (picoseconds to nanoseconds) and observed similar trends in the slower dynamics of the materials. 29 In this work, we further explore this new approach by performing multiple MD simulations of regio-random P3HT (RRa-P3HT) in order to quantitatively compare the impact of different simulation parameters on polymer dynamics in the solid-state. We do not focus on regio-regular P3HT, in spite of its superior electronic properties with respect to RRa-P3HT, because it contains coexisting crystalline and amorphous regions. This would require simulation of sufficiently large systems to include adequate sampling of both domains. 30 Moreover, if the dynamics of amorphous phases is not properly demonstrated, it would be confounding to analyze a system containing both amorphous and crystalline domains. Simulating a multi-phase system introduces other complexities into the simulation since there can be differing fractions and sizes of crystalline and amorphous domains between experiments and simulations. For these reasons, we first focus on RRa-P3HT, which is fully amorphous over a broad range of temperatures. This ensures that differences in the structure and dynamics of simulations arise primarily from the parameters that are used in each simulation (e.g. force-fields, equilibration).
Experimental QENS and WAXS data are used to assess the accuracy of individual MD simulations in order to develop guidelines for future computational analyses of conjugated polymer systems. Contrast variation of the RRa-P3HT is also used to highlight fluctuations in the side chains and in the backbone with complementary experiments. With these methods, five different MD FFs that have been parameterized specifically for poly(3-alkylthiophene)s using ab initio methods are compared for accuracy and the most impactful model parameters are identified. Finally, we investigate the effects of dynamic fluctuations on charge transport parameters (e.g. the charge transfer integral and conjugation length).
Regio-random P3HT (RRa-P3HT) (Mw = 48 kg/mol, Đ = 2.4) was purchased from Rieke Metals and is utilized in the first portion of this work. Another fully hydrogenated regio-random P3HT (P3HT-H14) (Mw = 90 kg/mol, Đ = 6.9) sample and a partially deuterated (side chain only) regio-random P3HT (P3HT-D13) (Mw = 83 kg/mol, Đ = 6.3) sample were synthesized by the procedure outlined below.
Nuclear magnetic resonance (NMR) spectroscopy (Figure S1 in the supplementary information), differential scanning calorimetry (DSC) (Figure S2 in the supplementary information), and size exclusion chromatography (SEC) (Figure S3 in the supplementary information for RRa-P3HT) are provided. NMR and SEC spectrum for RRa-P3HT were acquired from Rieke Metals. For synthesized P3HT-H14 and P3HT-D13, products were characterized by NMR on a Bruker 300 MHz AVance series instrument. Size exclusion chromatography (SEC) was performed in chlorobenzene for polymer samples (0.5 mg/mL) on a Malvern Viscotek TDA-305 instrument coupled with UV and RI detectors, at a flow rate of 0.5 mL/min referenced to polystyrene standards. All DSC measurements were performed with a TA Instruments Q200 using aluminum sample pans.
1-bromohexane-D13 was purchased from CDN Isotopes. Magnesium turnings, 1,3bis(diphenylphosphino)propanenickel(II) chloride (Ni(dppp)Cl2), 3-bromothiophene, iron(III) chloride, and anhydrous chloroform were purchased from Sigma-Aldrich and used without further purification. 3-hexylthiophene was purchased from TCI Chemicals. Anhydrous diethyl ether was purchased from Fisher Scientific and stored under activated 4Å molecular sieves. The syntheses of P3HT-H14 and P3HT-D13 were performed under nitrogen using standard Schlenk line techniques.
The synthesis of 3-hexyl-D13-thiophene was performed with modification to literature procedures. 31,32 Magnesium (0.8750 g, 36.00 mmol) was heat dried in an acid-washed 500 mL 3neck round bottom flask. Anhydrous diethyl ether (280 mL) and several small crystals of iodine were added to the flask. A portion of 1-bromohexane-D13 (5.012 g, 28.13 mmol total) was added via syringe. The reaction was initiated by gentle heating and crushing of magnesium with a glass stir rod under a nitrogen flow. Upon initiation, the remaining 1-bromohexane-D13 was gradually syringed and stirred for 2 hours at room temperature. In a separate 1 L flask, 3-bromothiophene (2.37 mL, 25.29 mmol) and Ni(dppp)Cl2, (0.0230 g, 0.0424 mmol) were stirred in anhydrous diethyl ether (280 mL) and cooled to 0° C in an ice bath. The aforementioned Grignard solution was carefully added to the cooled flask via syringe and stirred overnight. The contents were washed with 0.1 M HCl, and the organic layer was further washed with DI water and saturated aqueous NaCl solution. The organic layer was dried with anhydrous Na2SO4 followed by solvent removal with rotary evaporation. The product was purified by distillation under reduced pressure to collect a clear liquid in 57% yield. 1 H NMR (300 MHz, CDCl3): δ 7.23 (m, 1H), 6.92 (m, 2H).
Regiorandom poly(3-hexyl-D13-thiophene) (P3HT-D13) and poly(3-hexylthiophene) (P3HT-14) were synthesized according to literature procedures 33 in 65% and 62% yield, respectively. Figure 1 shows synthetic steps for regiorandom P3HT with deuterated and hydrogenated side chains. P3HT-D13: 1 H NMR (300 MHz, CDCl3): δ 6.97 (m, 1H). P3HT-H14: 1 H NMR (300 MHz, CDCl3): δ 6.98 (m, 1H), 2.79 (b, 2H, head-tail), 2.55 (b, 2H, head-head), 1.71 (m, 2H), 1.48-1.28 (m, 6H), 0.91 (m, 3H) RR = 49%. For P3HT-D13, the same RR was assumed because of the same synthetic procedure used. See more information in Figure S1 in the supplementary information.
WAXS measurements were performed with an Anton-Paar (Graz, Austria) SAXSess X-Ray scattering instrument with a Cu-Kα source (wavelength 1.54 Å) in line collimation mode. Samples were prepared by melting polymer powder into a ~1 mm thick film between two 10-μm thick Mylar sheets all held in place by an aluminum frame. The sample was held at the target temperature under vacuum, and a thermocouple put in contact with the polymer out of the direct beam path accurately recorded the temperature. Fujifilm (Greenwood, SC) image plates and a PerkinElmer Cyclone (Covina, CA) image plate reader were used for data acquisition. The SAXSQuant software was used to reduce the data from 2-D scattering images to 1-D plots of intensity (I) versus scattering vector (q).
All QENS samples were melt processed in aluminum foil pouches shaped cylindrically to fit along the inside wall of the standard aluminum sample canisters. The polymer was melted at approximately 538 K in an argon glovebox and spread evenly across the unfolded pouch. Enough material was used to form films between 0.1 mm thick (hydrogenated samples) and up to 0.2 mm thick (deuterated sample). For the P3HT-H14 and P3HT-D13 samples, a thin lead wire was used to seal the edges and prevent the material from leaking. The aluminum was then folded to enclose the material and a large heated aluminum block (approximately 2.3 kg) was placed on top of the pouch to reduce heterogeneity in film thickness. QENS measurements were carried out using three different instruments: the Disk Chopper Spectrometer (DCS) 34 at the NIST Center for Neutron Research (NCNR), the Backscattering Spectrometer (BASIS) 35 at the Spallation Neutron Source (SNS) at Oak Ridge National Laboratory (ORNL), and the High Flux Backscattering Spectrometer (HFBS) 36 at the NCNR. DCS is sensitive to energy transfers corresponding to time scales ranging from 0.1-10 ps, BASIS from 10-1000 ps, and HFBS from 100-2000 ps. Scattering vectors (Q) ranged from 0.3 to 1.9 Å -1 for DCS and BASIS, and ranged from 0.3 to 1.7 Å -1 for HFBS. Dynamics measurements were completed at temperatures between 273 K and 473 K (resolution taken at 50 K), and elastic scans (or fixed window scans, HFBS only) were completed between 50 K and 523 K. All measurements were performed in a top loading closed cycle refrigerator (CCR) sample environment, and samples were kept under vacuum to protect the material from degradation at high temperatures.
For standard dynamics measurements, data reduction and Fourier transformation of the data from the frequency domain to the time domain was performed using the Data Acquisition and Visualization Environment (DAVE) software 37 developed by the NCNR and the Mantid software 38 developed by the SNS. Data reduction from HFBS fixed window scan measurements to mean squared displacement was also performed utilizing the DAVE software. 37 A QENS experiment measures the neutrons scattered (𝝈) per energy change (ℏ𝝎) and per scattering angle (𝛀), which defines the double differential scattering cross section (
). 20,21,39 This is closely related to the scattering function 𝑺(𝑸, 𝝎), or the observed intensity as a function of scattering vector 𝑸 and frequency 𝝎. Additionally, a QENS fixed window scan collects scattering intensity at the elastic condition of ℏ𝝎 = 𝟎. 36 The double differential cross section can be further broken down into its coherent and incoherent components, ∑ 𝒃 𝜶 𝒃 𝜷 𝜶,𝜷 𝑺 𝒄𝒐𝒉 𝜶,𝜷 (𝑸, 𝝎) and ∑ 𝚫𝒃 𝜶 𝟐 𝜶 𝑺 𝒊𝒏𝒄 𝜶 (𝑸, 𝝎), respectively, weighted by scattering cross sections of isotopes in the samples. 20,21,39 𝒃 𝜶 and 𝒃 𝜷 are neutron scattering lengths for isotopes 𝜶 and 𝜷, and 𝚫𝒃 𝜶 𝟐 is the neutron scattering cross section for isotope 𝜶. Thus, the coherent and incoherent scattering functions provide information about pair-and self-motions of the atoms in the sample, respectively. Finally, a series of Fourier transforms recover the intermediate scattering functions,
(𝑸, 𝒕) and 𝑰 𝒊𝒏𝒄 𝜶 (𝑸, 𝒕), and the van Hove correlation functions, 𝑮 𝜶,𝜷 (𝒓, 𝒕) and 𝑮 𝜶 (𝒓, 𝒕). 20,21,39 Transforming the QENS data from the energy domain (collected) to the time domain with a Fourier transform not only allowed us to remove the effects of instrument resolution, but also allowed us to compare data across multiple instruments over a broad range of accessible time scales. This was achieved with the following equation 21 :
where 𝑺(𝑸, 𝝎) is the measured scattering signal and 𝑹(𝑸, 𝝎) is the instrument resolution function. Since the resolution differed across instruments, all experimental data was vertically shifted to align overlapping time regions. This equates to a change in magnitude of the flat background. This does not modify the captured dynamics as normalization or scaling the data would.
3 Computational Methods
Molecular dynamics simulations were carried out using the LAMMPS Molecular Dynamics Simulator. 40 A customized version of LAMMPS was built for systems utilizing the FF of Huang and coworkers 41 only to incorporate the additional dihedral parameters and is available upon a request to the authors. The system was built with 64 chains of regio-random P3HT with 60 monomers (RRa-P3HT) or 41 monomers (P3HT-H14/-D13) per chain to accurately represent the molecular weight of the corresponding polymer used in the experiments associated to each simulation (see Figures S3 andS4 in the supplementary information). In all of our systems, polydispersity of these materials were not accounted for and held at 1. In Figure S4 of the supplementary information, we explore the effects of chain lengths on the dynamics of the system. Above a chain length of 40 monomers, there is little change in the dynamics with increased chain length. Since the distribution of molecular weights in our sample is log-normal and biased towards longer chain lengths (Figure S3 in the supplementary information), we expect minimal effects due to polydispersity at the time and length scales of interest to this study.
Moreover, since MD simulations of polydisperse chains are very computationally expensive, we use the number average molecular weight to simulate monodisperse systems.
Chains were initialized in a large simulation box (corresponding to an initial density of 0.0671 g/cm 3 ) to allow for quick relaxation. Box dimensions then equilibrated to approximately 113 Å x 99 Å x 99 Å for RRa-P3HT and 91 Å x 92 Å x 91 Å for P3HT-H14/-D13. All chains were also generated such that the regio-regularity of the system reflected the experimental values determined from NMR (see Figure S1 in the supplementary information). Finally, a Metropolis algorithm was used to set the initial distribution of cis and trans monomer-monomer configurations to a population that would reflect the corresponding backbone torsion potential for each FF (see Figure S5 in the supplementary information).
All simulations in this work utilized the NPT ensemble. The velocity-Verlet algorithm was used with a timestep of 1 fs for integration of the equations of motion. The long-range Coulomb interactions were calculated using the particle-particle/particle-mesh (pppm) Ewald algorithm. 42 Pressure was set at 0.001 atm, which approximated the sample conditions during QENS and WAXS experiments. Temperatures were explored in the range of 13-523 K to again match the experimental conditions. Pressure and temperature were maintained using the Nose-Hoover barostat and thermostat, respectively, with characteristic time constants of 𝝉 𝒑𝒓𝒆𝒔𝒔𝒖𝒓𝒆 = 1 ps and
For standard dynamics simulations (RRa-P3HT 60mers), an initial equilibration at 600 K for 2 ns was performed followed by another equilibration at 473 K for 4 ns. Subsequently, the system was cooled with a 200 ps run every ten degrees until 13 K, providing starting places for simulations at a variety of temperatures. At each temperature point desired, another 4 ns equilibration was completed followed by another 5 ns production run from which all trajectory data was gathered. For those simulations used to calculate mean squared displacement curves (P3HT-H14/-D13 41mers), a similar procedure was used except that the second equilibration occurred at 523 K rather than 473 K (to increase the temperature range explored). Once again, the system was cooled with a 200 ps run every ten degrees until 13 K. At each temperature point desired, a 1.7 ns production run was completed to calculate the net mean squared displacement in the system.
To properly determine the effect that different FF parameters will have on the structure and dynamics of simulations of RRa-P3HT, it was first necessary to perform a sensitivity analysis on the FF independent parameters. We first investigated the impact that different equilibration methods have on the structural characteristics of the MD simulations. The timescales for chain diffusive motions (i.e. reptation) for conjugated polymers are much longer than is accessible with these simulations or with QENS experiments due to the comb-like structure of chains, which causes large molecular friction. Therefore, it is difficult for chains in MD simulations to reconfigure and achieve an equilibrium structure within a reasonable time-frame. As a result, it was essential to evaluate the effects of altering the equilibration process on the structure of the system. Common equilibration strategies for polymer systems (i.e. those without pi-conjugated orbitals), often involved placing chains within a simulation cell followed by heating to temperatures well above the glass transitions to facilitate relaxation and equilibration. 43 Because the energetic barriers for configurational changes of these polymers are relatively low, the polymers have enough energy to fully explore configurational space. The cell is then slowly cooled, with sufficient time between temperature steps for the chains to fully relax, until the desired temperature is reached. Previous simulations of P3HT have adopted this method. 44 A similar method involves placing individual chains in a very large simulation cell, corresponding to densities less than 0.1 g/cm 3 and then stepping down the volume of the cell while waiting for the chains to fully relax between steps. 45 Various combinations of these methods were tested to equilibrate the P3HT simulations and investigate their impact on the final polymer structure. Figure S6 of the supplementary information shows the static structure factor and radius of gyration for the RRa-P3HT system equilibrated with various ensembles, initial densities and temperatures between 473 and 600 K.
Evaluating radius of gyration is especially important for conjugated polymers as it affects the extent of conjugation of a chain and the intra-chain charge localization length. More discussion on the effects of these ensembles can be found in the supplementary information, but in the end, we chose an isotropic NPT equilibration at a temperature of 600 K with an initial density of 0.0671 g/cm 3 .
We also investigated the effect of equilibration at an extreme temperature (1500 K) on the structure and dynamics of our system. Following the work of Alexiadis et al. 46 , we tracked the chain end-to-end autocorrelation and system density during various equilibration procedures at temperatures of 600 and 1500 K for times as long as 30 ns. We then observed the effect these procedures had on the relaxations of the system (dynamic structure factor). This information can be found in Figure 2 in the main text and Figure S7 of the supplementary information. Although it was observed that a temperature of 1500 K was indeed required to fully relax the end-to-end autocorrelations, there was little sensitivity of the resulting chain relaxations measured at the relevant time and length scales to the longer equilibration periods. For example, no differences were observed apart from small changes at the longest times and highest temperatures (i.e. > 1 ns for MD simulations at 473 K). Because we are primarily concerned with the short-range structure and dynamics of these materials, as measured by QENS, we use equilibration procedures at 600 K for 2 ns for computational efficiency.
In this work, we utilized five different FFs that were published by Bhatta and coworkers 47 , Huang and coworkers 41 , and Moreno and coworkers 48 . They are all parameterized using different methods and validated against different criteria (e.g. crystal structure, density, surface tension, glass transition temperature and/or melting temperature). All FFs are based on the OPLS-AA FF 49- 51 with special modifications in backbone and side chain torsion potentials and atomic partial charges to account for the effect of conjugation in P3HT. The modified backbone torsion potential accounted for the increased rigidity of the polymer due to pi-conjugation between the neighboring monomers. The atomic partial charges are also modified to more accurately represent the distribution of charges that result from charge delocalization, as determined from first principles methods. All adjusted model parameters were derived from calculations using post-Hartree-Fock or hybrid-density-functional-theory methods. 52 For additional discussion regarding the relationship between backbone and side chain regimes during FF parameterization, see Figure S8 in the supplementary information.
After initial evaluation of the FFs, a large discrepancy was noticed between the side chain torsion potentials derived by Bhatta and coworkers 47 and the potentials used by Moreno 48 and
Huang 41 , especially for the first and second dihedrals (β1 and β2). Bhatta and coworkers had calculated the torsion barrier using a single point energy calculation from an optimized P3HT dodecamer, but it was noted that side chains of a rotating monomer could intersect with neighboring side chains at low dihedral angles (cis transformation of thiophene rings). This Additional differences were noted in the backbone torsion potentials applied to the bond connecting adjacent monomers (α) between all FFs. The Bhatta potential was derived from single point calculations on a P3HT decamer at the B3LYP/6-31+G(d,p) level. 47 The Huang potential was derived through a combined MD and DFT method from which the intrinsic torsion potential was extracted. 41 The Moreno potential (identical among FF1, FF2 and FF3) was derived at the B3LYP/6-311G(d,p) level for quarterthiophene with an added molecular mechanics step to ascertain the intrinsic torsion potential. 48 Figure 4 shows the results of these differing methods, including a range in the torsion barrier heights at 90˚ as well as conflicting preferences for cis and trans conformations along the thiophene backbone. Lastly, significant variation among the FFs was uncovered in atomic partial charges, especially for those atoms on the thiophene ring. Minimal agreement exists in either magnitude or sign of the partial charges, as shown in Table 1. Across all FFs, these parameters were developed with different thiophene oligomers and ab-initio basis sets. The partial charges of the modified Bhatta FF were developed at the MP2/6-31+G(d,p) level of theory with restrained electrostatic potential (RESP) across a P3HT decamer. 47 The authors allowed partial charges in all atoms of the P3HT monomer to vary, including those of the side chain. The Huang FF determined partial charges for the thiophene backbone from previous ab-initio calculations of tetrathiophene at the MP2/aug-cc-pVTZ level of theory with the electrostatic potential (ESP). 41,53- 55 Most of the partial chargers for the alkyl side chain were then taken from the OPLS-AA force field, except for the first alkyl carbon bonded to the thiophene ring which was adjusted to ensure neutrality of the monomer. 41 The three variations of the Moreno FF only differed in the atomic partial charges implemented. Moreno FF1 used general OPLS-AA parameters, Moreno FF2 was developed at the B3LYP/6-311G(d,p) level of theory with the electrostatic potential (ESP) from four different thiophene oligomers, and Moreno FF3 was an average of the FF1 and FF2 partial charges. 48 In the Moreno FF2 variant, most partial charges of the alkyl side chain were taken from the OPLS-AA force field except for the first alkyl carbon bonded to the thiophene ring. 48
The trajectories obtained from the production runs were used to calculate the static and dynamic structure factors using the following formulae 56,57 : (5) Equation ( 2) was used to calculate the static structure factor, weighted by atomic number of each atom (Z) to allow for comparison to experimental WAXS results. Equations ( 3) and ( 4) were used to calculate the incoherent and coherent dynamic structure factors, respectively, weighted by neutron scattering cross sections for each atom to allow for comparison to experimental QENS results. For hydrogenated samples, the QENS signal is dominated by incoherent contributions and we only utilize Equation (3). To efficiently run these calculations, the nMolDyn software was used 56 . For the static structure factor, 100 Q-vectors were generated from 0 to 2.5 Å -1 with a step size of 0.02 and a Q-shell width of 0.02 Å -1 . For the calculation of I(Q, t), 100 Q-vectors were generated from 0.3 to 1.9 Å -1 with a step size of 0.2 Å -1 and Q-shell width of 0.1 Å -1 .
Dihedral autocorrelation functions (DACFs) were calculated for the backbone (α) and the bond between the second and third carbon atoms in the side chain (β1) using the following formula 58 :
where ϕ(t) is the torsion angle at time t. In our system, a cis conformation of the thiophene rings equates to a 0˚ dihedral angle, while a trans conformation equates to a 180˚ dihedral angle.
For mean squared displacement (MSD) calculations, each simulation was cooled to 13 K utilizing the 1.7 ns cooling steps discussed above. At each temperature (every 10 degrees) 200 ps of run time was reserved for equilibration at the new temp. The remaining 1.5 ps was used to calculate a time-averaged (𝑡) MSD with the following formula:
where 𝑁 is total number of atoms in the system, and 𝜎 𝛼,𝑖𝑛𝑐 and 𝜎 𝛼,𝑐𝑜ℎ are the incoherent and coherent neutron scattering cross sections 59,60 , respectively. The characteristic time scale, 𝜏, was set at 1 ns to reach similar time scales to those of HFBS and balance with computational cost of the system.
To track changes in the distance between neighboring thiophene rings (e.g. pi-stacking distance), a mean squared relative displacement was calculated with the following formula:
where 𝑁 𝛼𝛽 is the total number of pairings between thiophene rings on separate chains or rings on opposite ends of the same chain. The distance to each ring (𝑟 𝛼 , 𝑟 𝛽 ) was calculated to the center of each ring. In this work, we refer to the net mean squared displacement of all atoms as MSD (Equation 7) and the mean squared relative displacement of thiophene-thiophene distances as MSRDthio (Equation 8).
To determine system ordering of the thiophene ring-to-ring distances, the radial distribution of the rings was calculated with the following formula for the radial distribution function:
where 𝑁 𝑟,𝑡ℎ𝑖𝑜 is the number of thiophene rings in the shell at distance 𝑟 with thickness 2Δ𝑟 and 𝑁 𝑡ℎ𝑖𝑜 is the total number of thiophene rings. Total system density of thiophene rings is represented by 𝜌 𝑡ℎ𝑖𝑜 .
The chain end-to-end autocorrelation function (EEACF) was calculated using Equation 9,
where 𝑢 is a chain's end-to-end vector measured at an instant in time 46 :
4 Results
MD simulations are increasingly used as the foundation for ab initio calculations to evaluate the electronic properties of conjugated polymers. Therefore, it is necessary to evaluate the accuracy of FFs used in MD simulations with a variety of experimental methods. Our results
show that there are still discrepancies between the experimental data and observables calculated from simulations when using different FF parameters for P3HT, a model semiconductive polymer. Moreover, certain FFs show better agreement with experiments depending on the specific benchmark that is being used to assess them.
The first assessment compares the static structure factors calculated from MD simulations to experimental WAXS data for the RRa-P3HT sample (Figure 5). All peaks were fit with Gaussian curves and their locations are shown below in Table 2. Interestingly, all FFs (Mod.
Bhatta, Huang, and Moreno FF1, FF2 and FF3) performed fairly similarly, capturing both peak locations reasonably well. Notably, the simulations show a lower intensity for the first peak in comparison to the experiments. This is most likely due to the limited simulation box (i.e. ~100 Å in all dimensions), which can effectively 'cap' the observed long-ranged correlations. We would expect the relative intensity of the first peak to increase if the simulation size was increased but this would significantly increase computational expense. Regardless, all MD simulations show reasonable agreement, as judged by peak positions and curve shape, with the WAXS profiles across all length scales. This suggests that the quantitative comparison of correlations (i.e.
structure factors) in MD simulations of amorphous P3HT to WAXS data is not a good metric, on its own, for assessing the quality of MD simulation FFs. In contrast, differences between the models become more apparent when comparing dynamic structure factors calculated from MD simulations with experimental QENS data. Figure 6 shows theoretical and experimental data at 473 K, while the remaining data at 273 K, 373 K and 423 K can be found in the supplementary information (Figure S10, S11 and S12, respectively).
Because hydrogen has an incoherent cross section that is almost two orders of magnitude higher than that of all other atoms, 39,59,60 the QENS signal for P3HT is most sensitive to side chain fluctuations. Moreover, the P3HT side chain contributes to 93% of the incoherent scattering, which dominates the signal at these q-values. While there is reasonable agreement to experimental QENS data across all FFs at the highest Q-values (Q > 1.1 Å -1 ), it becomes apparent the Mod. Bhatta FF is inadequate in describing the data. It fails to capture dynamics across all Qvalues and time scales for the lower temperatures (especially at 273 K in Figure S10 of the supplementary information) and drastically underestimates decays at lower Q-values and longer time scales at 473 K (see Figure 6). Moreno FF1 and FF3 do better at capturing the trends at low Q and long times, but then start to fail at higher Q-values and long times. While Moreno FF2 and Huang FF also underestimate the decay at low Q and long times, it is to a significantly lesser extent than the Mod. Bhatta FF. Moreover, they show similar reasonable agreement to all experimental QENS data across all temperatures. To further quantify these trends, Figure 7 shows average χ 2 values between the calculated and experimental QENS data from all temperatures, showing the validity of either the Huang or Moreno FFs in simulating amorphous P3HT.
Figure 8 shows a comparison of the mass density, as obtained from MD simulations for each FF, and the experimentally determined values. All simulations significantly underestimate the experimental density. This is observed both for MD simulations performed with sequential cooling and also during production runs at a fixed temperature. Here, the Huang FF is found to come closest to approximating density values that are obtained from neutral buoyancy experiments. Guilbert and coworkers have also reported an underestimation of the experimental density from simulations of amorphous RR-P3HT, and have ascribed the difference to the presence of crystalline domains in the experimental system. 28 We have shown here that this discrepancy exists even for purely amorphous systems across all investigated FFs, further increasing the need for better-designed FFs for these systems.
An advantage of MD simulations is that one can explore the molecular system beyond what may be directly extracted from experimental data, i.e. WAXS and QENS measurements.
Figure 9 shows the backbone dihedral population and dihedral autocorrelation functions for all five FFs. The peaks of the distribution of P3HT backbone configurations (Figure 9a) in the Huang FF simulation are centered at 0° and 180°, which correspond to completely cis and trans backbone configurations, respectively. This implies a higher planarity is preferred between thiophene rings when compared to other FFs, for which the peaks are shifted to 30° and 150°.
The discrepancies stem from the underlying backbone potential for each of the FFs (Figure 4). Figure 10 shows a similar analysis of the inter-chain thiophene ring relationships in the modeled polymer systems. The radial distribution function gives information about the spacing between thiophene rings on neighboring chains or on opposite ends of the same chain (Figure 10a). This is significantly important when considering the inter-chain charge transfer mechanisms due to the pi-orbital overlap as previously discussed. Besides the Moreno FF1 system that shows an overall closer ordering between thiophene rings, many of the modeled systems show similar results. However, the mean squared relative displacement (MSRDthio) of these distances (Figure 10b) shows clear differences in structural fluctuations. Since backbone partial charges change drastically for different MD FFs, it is logical to anticipate that this would have an effect on thiophene-thiophene separation distance between neighboring rings (i.e. pi-orbital overlap).
We have now shown that torsional parameters and partial charges have a significant impact on the underlying structure and dynamics, especially for the thiophene rings along the backbone. However, as previously discussed, it is not straightforward to decouple backbone and side chain fluctuations from QENS measurements using a single sample. Moreover, some important dynamic fluctuations (backbone) may be obscured by other dominant fluctuations (side chains) in the scattering data. Hydrogen/deuterium labelling can be used to alter the neutron contrast and to 'highlight' areas of interest in complex molecules. In this work specifically, we replaced all of the side chain hydrogens of P3HT with deuterium. This significantly reduced the side-chain contributions to the incoherent cross section from 93% to 25%, and biased the QENS signal towards measuring backbone motions, which are known to be key to charge transport in conjugated polymers. thiophene motions that occurs at a much higher temperatures. It is important to note that these two samples were prepared to be as identical as possible with the only major difference being the hydrogen / deuterium content in the side chains. When we compare experimental and MD results quantitatively, the simulations show larger system motions. This is not unexpected due to the underestimated density (i.e. more free volume) in the simulated systems, which could leave more free-volume for fluctuations to occur. In more qualitative comparisons, the simulations capture the MSD of the hydrogenated material very well with similar activation temperatures.
Side chain relaxations occur in all systems at approximately 100 K and the profiles compare well across the broad range of temperatures. More notably, simulations completely fail to capture MSD of the deuterated material as measured in QENS experiments. The delayed relaxation of the backbone motions above 400 K was not present and the material exhibits a similar relaxation to that of the hydrogenated material. Although we have shown that each FF has its own strengths and weaknesses, all models have proven to be inadequate at capturing the backbone motions of RRa P3HT, which is crucial to describing charge transport in this material.
Finally, we consider all of these results in terms of their relevance to charge transport mechanisms. The Marcus rate is a general expression commonly used to analytically describe charge transport in conjugated polymers using a hopping model that is temperature-dependent.
The equation is as follows: 14
where 𝑘 𝑖𝑗 is the charge hopping rate from site 𝑖 to site 𝑗, Δ𝐸 is the energy difference of the charge transfer, 𝜆 is the reorganization energy, ℏ is Planck's constant, 𝑘 𝐵 is Boltzmann's constant, and 𝐽 𝑖𝑗 is the transfer integral, or electronic coupling. 14 While most of these parameters can only be evaluated for P3HT by means of computationally expensive ab initio calculations, the dependence of the charge transfer integral on the relative positions of the conjugated rings is well established. 61 Figure 12 shows values of the charge transfer integral calculated using quantum mechanical methods for both inter-chain and intra-chain charge hopping in P3HT along with fits. 61 It is evident that the magnitude and sign of the transfer integral for intra-chain charge transport (Figure 12b) are highly dependent on the backbone torsion angle (α). Similarly, the strength of the inter-chain transfer integral exponentially decays as the distance between adjacent thiophene rings increases. It is therefore possible to qualitatively assess how charge transport properties calculated from snapshots of bulk MD simulations of P3HT could be affected by the choice of FF in a particular MD simulation. When examining Figure 9 and Figure 10 through the Marcus theory lens, it is immediately clear that different FF parameters (particularly those parameterized from ab initio methods) will affect the accuracy of an ab initio calculation made from a snapshot of an MD simulation. For example, Figure 9 suggests that the Mod. Bhatta FF would represent a system with facilitated intra-chain charge transport along the backbone when compared to Moreno FF1 because more trans backbone configurations would be present. The converse would be true when considering inter-chain charge transport because the probability of 'seeing' a ring from another chain within the range required for inter-chain charge hopping is greater for Moreno FF1 (Figure 10a). Finally, we reconsider the results presented in Figure 11.
The sensitivity of the charge transfer integral to conformational changes along the backbone reinforces the importance of accurately capturing the backbone motions. This points to the need for more rigorous scrutiny when calculating FF parameters using ab initio methods and in the selection of FFs for the estimation of bulk electronic properties from simulations.
). This is not a physically meaningful method to determine partial charges. We further consider the parts of the FF that were not re-parameterized for conjugated systems, such as the non-bonded Lennard Jones (LJ) parameters. For all FF investigated in this work, LJ parameters were taken directly from the OPLS-AA FF. 51 For example, carbon and hydrogen parameters were borrowed from those of a benzene molecule. However, these classical non-bonded values were determined based on liquid-phase experimental data of these simple molecules. 62 The effects of long-range conjugation, such as stronger dispersion forces 63 , between polymer chains would not be properly accounted for in these FFs. For example, the work of Schmit et al. investigated the binding energy between two approaching poly(acetylene) chains and found a strong attraction due to electron tunneling in the system. 64 This non-bonded binding energy between different chains is almost four times stronger than what the LJ interaction would have predicted between carbon atoms. Therefore, there is reason to believe that there is also a significant underestimation of non-bonded attraction in P3HT when adopting LJ parameters from classical FFs (e.g OPLS). Importantly, stronger non-bonded interactions would directly impact the system's density and the backbone fluctuations, which are the two aspects of greatest discrepancy between experiments and simulations in this work. Thus, we emphasize the need for more thorough evaluation and re-parameterization of non-bonded parameters of new FFs for conjugated materials.
In this work, we have also gained valuable insight into the limitations of experimental methods, especially their lack of sensitivity to key parameters that could only be probed in silico (Figure 9 and Figure 10). For example, none of the x-ray or neutron scattering data was determined to be sensitive to differences in the torsion populations or the radial distribution functions. Furthermore, comparisons between experimental QENS data and simulated dynamics for fully hydrogenated P3HT was not sufficient, on its own, to select the optimum FF from Moreno FF2, Moreno FF3, and Huang FF. Rather, it was through the use of isotope-labelled samples and QENS that a significant weakness of all FFs was identified. This is likely due to the fact that side-chains, which dominate QENS signals in hydrogenated P3HT, have 'classical'
interactions that are adequately captured by most FFs. In contrast, the backbone dynamics, as measured with samples with deuterated side-chains, are greatly affected by electronic conjugation and require the development of new FFs for adequate representation.
Finally, we make note of the challenges ahead in the development of accurate FFs for conjugated polymers with future in-silico material design in mind. Although P3HT is a model material, our work shows that we still fail to accurately capture molecular structure and dynamics relevant to charge transport mechanisms. We must first focus on correcting these faults if we are to develop accurate FFs for more complex and higher performing CPs. A key component to this process will be the use of meaningful experimental data for a critical assessment of FF parameters. To facilitate this, we provide all experimental QENS data and relevant meta-data that is presented in this work in the supplemental information. As we move forward towards open science, we encourage others to use this data for the development of improved modelling methods for these materials. We believe wide access to this type of experimental data within the community is essential for advancement in the field. Moreover, this takes us one step closer towards development of accurate, high-throughput methods for FF parameterization with validation to enable in-silico design of novel materials. A bottleneck in this work is the exploration of the simulation parameter space with many MD simulations with high uncertainty. Streamlining this process is critical if we hope to achieve efficient in silico material design that informs chemical synthesis for better performing materials.
In this work, we have extensively explored parameters involved in running MD simulations of conjugated polymers, specifically regio-random P3HT. We have quantitatively and qualitatively assessed the performance of five different FFs from Bhatta and coworkers, 47 Huang and coworkers, 41 and Moreno and coworkers. 48 WAXS alone was insufficient in quantitatively evaluating different models. However, when comparing simulations to QENS measurements, Moreno FF2, Moreno FF3 and Huang FF showed better performance in capturing system motions. Of these three FFs, the Huang FF was closest to matching the experimental density but there were still significant discrepancies. We also show that both the magnitude and the sign of the atomic partial charges play a significant role in determining the dynamics and the structure of the system. By comparing the three Moreno FFs, we saw that the sign of the partial charges