General Force-Field Parametrization Scheme for Molecular Dynamics Simulations of Conjugated Materials in Solution.

We describe a general scheme to obtain force-field parameters for classical molecular dynamics simulations of conjugated polymers. We identify a computationally inexpensive methodology for calculation of accurate intermonomer dihedral potentials and partial charges. Our findings indicate that the use of a two-step methodology of geometry optimization and single-point energy calculations using DFT methods produces potentials which compare favorably to high level theory calculation. We also report the effects of varying the conjugated backbone length and alkyl side-chain lengths on the dihedral profiles and partial charge distributions and determine the existence of converged lengths above which convergence is achieved in the force-field parameter sets. We thus determine which calculations are required for accurate parametrization and the scope of a given parameter set for variations to a given molecule. We perform simulations of long oligomers of dioctylfluorene and hexylthiophene in explicit solvent and find peristence lengths and end-length distributions consistent with experimental values.

Introduction

Semiconducting conjugated polymer materials offer a wealth of opportunity for the development of organic-based optoelectronics with several advantages including light weight, flexibility, low toxicity, and inexpensive fabrication for applications such as photovoltaic cells and light-emitting diodes.[1−3] There remains a challenge however to achieve sufficient power conversion efficiencies and durability for viable devices. One of the key reasons behind this is the important role of material morphology and conformation within such materials.[4]

While it is understood that many of the desirable properties, such as the control over solubility allowing solution processing and low bulk modulus providing for flexibility,[1] arise due to the polymeric nature of the constituent molecules and the inherently statistical nature of their mixing, it is often the case that this same principle can lead to the existence of phenomena such as deep-tail trap-states which can inhibit the conductivity.[5] This is primarily due to the delicate physics of both intramolecular conjugation - its sensitivity to local distortions along the backbone of a polymer[6] and the effect this has on the resulting optical absorption and emission dynamics[6−8] - and the intermolecular excitation transfer dynamics - the interplay of alignment and separation of conjugated segments or ‘chromophores’ and their spectral overlaps with the Förster-type[9] transfer of excitons and polarons.[10,11]

These questions of the interplay of the dynamics and statistical mechanics of conjugated polymer-based materials suggest the application of Molecular Dynamics (MD) to understand and predict macroscopic properties based on detailed knowledge on an atomic scale. In such simulations, in lieu of an explicit quantum mechanical treatment, the classical dynamics of a molecule are generated using force-fields, which appropriately describe the averaged effect of the molecular electrons on the covalent bonding and van der Waal type forces. This significantly reduces the computational expense of analyzing questions of conformational properties and, when coupled with modern computer power, allows for the simulation of reasonably long chains.[12−15]

It is imperative that the force-fields used to describe molecular motion are of sufficiently high accuracy to reproduce observed experimental behavior. Historically, MD methods are often designed with biochemical simulations in mind, e.g. dynamical simulations of large proteins, and, as such, there exists a number of available force-fields[16−24] - each parametrized in order to optimize their efficacy. While these force-fields contain many parameters which are transferable to the polymer system, certain aspects require a careful reparametrization due to the conjugated nature of the molecular physics involved. The challenge is to utilize as many transferable parameters as is possible while identifying which parameters require further attention and reparametrization.

Arguably, the most important term to be considered is the energetic profile governing the dihedral dynamics between monomers in a conjugated system. Accurate modeling of the dihedral profile and energy barriers between conformers is of utmost importance as the excited-state landscape of conjugated molecules is governed by the dihedral angles,[25] and, as such, the optical properties of the materials are crucially dependent on these torsions.

Another key element to be considered is how interatomic electrostatic interactions, which arise due to local deviations in electronic charge densities, are described. This is most commonly implemented in MD algorithms by assigning ‘partial charges’ to the atomic centers within the molecules. This requires fitting atomic charges to the calculated electrostatic potentials of the molecule. While there exist a number of such fitting schemes,[26−29] the RESP scheme is generally considered to be one of the more robust and accurate procedures for this task.[30]

The majority of methodological approaches[14,31,32] for generating parameters for MD simulations of conjugated polymers give strong weight to the importance of the acquisition of appropriate intermonomer dihedral angle profiles and partial charges in order to obtain accurate predictions from the resulting MD. However, many of these sources provide conflicting viewpoints on the appropriate levels of quantum chemical theory required. Furthermore, there is often a lack of systematic parameter development and benchmarking of given methodology, and, as such, there is a large degree of ambiguity in the accuracy of a given method as well as whether or not a prescribed computational route could be replaced by a significantly less computationally expensive one. Given the wide range of organic molecules which are of interest for applications the computational cost and complexity of obtaining parameters is a constraining barrier to progress through the exploitation of MD.

A further point that has, to date, been neglected is how sensitive a given parametrization scheme is to variation in length of molecules and lengths of their associated alkyl side chains. A typical solution-cast mixture of conjugated polymers requires the attachment of branched alkyl side chains to provide solubility. It has also been shown that these side chains and their interference with the intrinsic molecular motion of the conjugated backbone of the molecule have a key role in effects such as the inhibition of excitonic diffusion[5] as well as the emergence of exotic bulk behavior such as the well-studied β-phase of polyfluorene.[33]

In light of the above, the work presented here establishes a systematic approach to the force-field parametrization. In addition, our methodology also determines the domain of applicability of a given set of parameters for molecules of varying length and varying length of side chain. We also highlight cases in which the universality of a given parameter set to these variations is broken. In doing so, we put forward a parametrization protocol which conforms sufficiently to established benchmarks of accuracy; avoids unnecessarily computationally intensive calculation methods; and may be applied, in extension, to any type of conjugated system with alkyl side chains.

Generating Molecular Dynamics Parameters

A minimal set of force-field parameters contains parameters of five types: three describing the energetics due to covalent bonding between atoms and two describing noncovalent interactions. The three covalent terms quantify bond-stretching, changes in the angle between three atoms (angle-bending), and changes in the dihedral angles between four atoms. These are modeled by functions ranging from quadratic (particularly for two-atom vibrations) to Fourier-based functions for the angular types. Noncovalent interactions are of the form of Lennard-Jones and Coulomb potentials which account for London dispersion, Pauli repulsion, and, in the case of the Coulomb potentials, electrostatic interactions between local variations in electronic density. We choose the OPLS[17−23] force-field parameter scheme (as implemented in Gromacs 4.6.5[34,35]) as our starting point due to its provision of parameters for many atoms in a multitude of different molecular frameworks as well as its use in previous works to parametrize conjugated polymers.[14,32,36]

To parametrize a given molecule, the first step is to build on the appropriate parameters found within OPLS for a monomer. For example, to build a molecule containing two connected units of thiophene (hereafter referred to as a 2mer and, for molecules with x connected units, as an xmer), the OPLS parameters for each atom of an individual thiophene unit are implemented for each unit. This leaves only the bonds, angles, and dihedrals associated with the linking bond between the two units and partial charges to be determined. This is illustrated in Figure .

Figure 1

Schematic highlighting the required new parameters in going from a parameter set for thiophene to dithiophene. The missing force-field parameters are the bond-stretching term between the two carbons circled in red; the angle-bending term between each combination of three blue-circled atoms (two from one unit and one from the other); and the dihedral angle term for each combination of four blue-circled atoms (two from one unit and two from the other).

We expect that bond-stretching and angle-bending within the monomer are reliably parametrized by the OPLS force-field terms thus the intermonomer dihedral profile is the primary unknown from the covalent terms. In the case of a nonbonding interaction, we can use existing Lennard-Jones terms without any modifications and focus on partial charge specification and validation.

The remainder of this Section outlines our general methodology for parametrization. Subsequent sections deal with the details of dihedral energy profile calculation (Section 3), partial charge generation (Section 4), implementation of modifications (Section 5), and preliminary MD calculations performed with our obtained parameter sets (Section 6). We also provide a further example of how our scheme might be applied for the copolymer F8BT in Section 7.

In order to generate sufficiently accurate dihedral potentials at low computational expense, we utilize the method of performing geometry optimizations over the span of a dihedral rotation (‘Scan’) using computationally inexpensive techniques and refine these results with a further single-point energy (‘SP’) calculation using more accurate methods. We find that the dihedral potentials for all lengths of conjugated backbone are accurately described by those of their corresponding 2mer.

We generate partial charges using the RESP[30] method with input electronic densities calculated using the above Scan - SP approach. With respect to molecules of varying backbone length, it is possible to generate a set of charges which are generalizable to any length of molecule. By considering the progression with length of the net charges on the internal monomer units, we can determinine the length at which the net charges go to zero. This indicates converged charge distributions which can be built into a general parameter set. With respect to variations of side-chains, there appears to be no means of generalizing such variations beyond recalculation.

Once the appropriate partial charges have been determined, these are directly implemented into the force-field. The dihedral profile, on the other hand, must be implemented by means of a ‘subtraction’ method so as to ensure that energetic terms already described by the force-field and partial charges are not double-counted. This is, in spirit, performing an analogous scan using the force-field parameters (without the required dihedral) in order to obtain the contribution already accounted for by the force-field and subtracting this from the calculated DFT profile. The resulting ‘subtracted’ profile is that which is implemented into the force-field. This is detailed in Section 5.

Determining Dihedral Energy Parameters

We choose fluorene and thiophene oligomers as our model systems with and without their alkyl side-chains (see Figure ). These molecules are experimentally well-characterized conjugated systems[33,37−39] against which to test our methodology.

Figure 2

Schematic of (a) fluorene and (b) thiophene 2mers with the associated dihedral angle, ϕ, highlighted in each case.

To obtain a reliable dihedral potential, we adopt a two-step ‘Scan - Single Point (SP)’ approach.[32,40] The ethos of the Scan - SP approach is to use lower levels of theory such as DFT for obtaining geometries and higher levels of theory, for example, local methods (MP2, CCSD, CCSD(T)) or DFT with a larger basis set for SP energy calculations.

The first step, the ‘Scan’, involves determining, by relaxed scan geometry optimization, a series of molecular geometries for different values of dihedral angle between the units of a 2mer. Geometry optimization can become a computationally cumbersome process, especially when large or exotic molecules are considered. However, in the case of many conjugated molecules, the calculation of reliable geometries for molecules of considerable length has been shown to be practically possible using moderate levels of theory.[32,36]

While it is possible to obtain accurate geometries in the Scan step, as we will go on to show, the energies obtained from this step are often quite inaccurate. This is why the second step, the ‘SP’ step, using higher levels of theory in order to obtain accurate energetic profiles is necessary.

For all DFT calculations, we employ the functional CAM-B3LYP[41] in Gaussian09.[42] We have determined that this initial choice performs comparably well to a number of DFT functionals as well as more involved methodology, such as MP2, CCSD(T), and density-fitted and local versions of each. Details of these tests and comparisons are provided in the Supporting Information (Section 1.1).

To account for the possible contributions from dispersion interactions not natively incorporated in DFT methods, we have also repeated our DFT calculations utilizing the GD3BJ dispersion correction functional of Grimme et al.[43,44] In cases with small (up to ethyl) side-chains, we find that dispersion effects are negligible. As such, the following results do not include these correction terms, and a more thorough breakdown of our GD3BJ results is given in the Supporting Information (Section 1.3).

Following our Scan-SP approach, we have calculated 2mer dihedral profiles using various basis sets from the 6-31G[45−54] family with added polarization and diffuse functions for scans and the cc-pVTZ[55] basis set for SP calculations. Our results are also compared against the benchmark CBS-limit CCSD(T) calculations of Bloom et al.[40] for thiophene. This first step allows us to ascertain the appropriate level of theory required for further extending our investigation.

Throughout this work we utilize the ‘polymer convention’ for dihedral angle labeling. This convention casts the trans conformation at 0° and the cis conformation at 180°. Dihedral potentials taken from other works have been transformed so as to fit with this convention.

Figure offers a few notable examples of the basis sets chosen. A full breakdown of the various methodological choices made is given in the Supporting Information (Section 1.2). For both fluorene and thiophene (Figure (a) and (b) respectively), comparison of our choices is made against a scan performed using cc-pVTZ for both Scan and SP (herein referred to as the cc-pVTZ result), and further comparison is made to benchmark results, obtained by Bloom et al., performed with MP2/aug-cc-pVTZ for the Scan and CCSD(T) in the Complete Basis Set (CBS) limit for the SP (herein referred to as the CCSD(T)/CBS result). It can be seen from Figure (b) that the cc-pVTZ result of thiophene is in good agreement with the CCSD(T)/CBS result. As performing such a high-level calculation as that of Bloom et al. for as large a molecule as fluorene is computationally prohibitive, we use thecc-pVTZ result for fluorene as a benchmark for our results.

Figure 3

Dihedral profiles from geometry scans for CAM-B3LYP with various basis sets in (a) fluorene and (b) thiophene 2mers. Those with additional cc-pVTZ single-point calculations are labeled (SP). Full geometry optimizations using cc-pVTZ are also included and in (b) a comparison is made to the CCSD(T)/CBS thiophene result of Bloom et al.[40]

From the results for fluorene (Figure (a)), it is observed that the energetics of the geometry optimization alone are quite inaccurate when compared to the cc-pVTZ result: Those without diffuse functions considerably overestimate (∼2 kJ/mol) both the barriers at the planar conformations and the barrier at 90°; and those with diffuse functions overestimate the planar barrier but considerably underestimate the 90° barrier to a similar degree.

For thiophene (Figure (b)), those which include diffuse functions in the basis set provide considerably less accurate results than those without. Using the diffuse functions, the energetics of the conjugation barrier (∼90°) appear to be underestimated when compared to the cc-pVTZ and CCSD(T)/CBS results. This difference is ∼2–3 kJ/mol which at ∼20–30% of the barrier itself is a significant discrepancy. However, it is also observed that none of the 6-31G-type basis sets exhibit suitable accuracy as both overestimate the height of the 180° energy barrier in thiophene by ∼2 kJ/mol. This indicates that the SP calculation step is indeed required.

Note that while there exists significant deviation in the energetic barrier heights due to the inclusion of diffuse functions, upon performing SP calculations with cc-pVTZ this difference becomes negligible. This implies that all discussed choices of the basis set provide very similar geometries and that the final accuracy of the dihedral profile is more dependent on the choice of the basis set for the following SP calculation. It follows from this that the appropriate choice for geometry optimization is 6-31G* due to it being considerably less computationally expensive than the others.

Utilizing a further SP calculation is found to offer significant improvement to the resulting potentials. In thiophene, it is observed that the error in the SP curves is ∼50% of that of the initial 6-31G* calculation when compared to the CCSD(T) curve. For the largest deviation, at the 180° barrier, this translates to a reduction of the absolute error to <1 kJ/mol from 2 kJ/mol. When compared to the 6-31+G** curve, this error reduction is considerably greater.

Comparison is also made to the results of a full optimization using cc-pVTZ (i.e., performing the Scan step with cc-pVTZ without any further SP calculation). This comparison shows that the results of our Scan-SP approach are almost identical to that of using cc-pVTZ for geometry optimization. This lends further credence to the idea that reliable geometries are obtained from basis set choices as low as 6-31G*. With the above in mind, our choice of CAM-B3LYP/6-31G* Scans and CAM-B3LYP/cc-pVTZ SP’s is well justified.

Now that a computational methodology is established, we examine to what extent the dihedral profile is sensitive to the length of the conjugated backbone. Figure demonstrates the invariance of the calculated dihedral profile of the end-most intermonomer junction for both molecules over a backbone length range of 2–10 units. In fluorene, there is almost no dependence on chain length, while in thiophene there is a very slight deviation of <∼1 kJ/mol in the energetic barrier heights. We attribute this deviation to the dispersive contribution of the additional neighboring thiophene. As the thiophene units are considerably smaller than those of fluorene, there is a far smaller separation between the end dihedral and the third unit which, in turn, leads to long-range dispersive forces between the third unit and the first unit which are strongest at the 90° position.

Figure 4

Dihedral profiles for various backbone lengths of (a) fluorene and (b) thiophene. (Legend applies to both graphs.)

We also demonstrated that the dihedral energetic profile is invariant to its position along a chain longer than a 2mer and that a profile is invariant to the value of the neighboring dihedral angle. These are shown and discussed in detail in the Supporting Information (Section 1.4). Thus, we have shown that we can utilize the same dihedral energy profile for all intermonomer dihedrals within any length of molecule. This greatly reduces the amount of parametrization required for MD simulations.

We now turn to the influence of alkyl side chains on the dihedral energy profiles. Scans have been performed for side chains of lengths varying from 1 to 10 units. In the case of fluorene, initial geometries with side chains on opposing sides of the molecule have been chosen. There is no inclusion of any dispersion correction within these calculations. There are two reasons for this: First, the inclusion of dispersion interactions in DFT calculations of this type often lead to convergence problems as the relatively high flexibility of the side chains allows for many possible local minimal conformations. Due to this, the calculations both become increasingly computationally expensive and generate results which are strongly dependent on initial conditions, resulting in ambiguity in their translation into force-field parameters. Second, by effectively removing the side chain–side chain interactions, what is observed is the effect of the presence of the side chains on the electronic properties of the backbone. For the purposes of generating good force-fields for MD calculations, this is the crucial point as this determines whether or not, in principle, the bare dihedrals along the molecular backbone require modification due to the presence of the side chain.

With the above considered, Figure (a) gives the resulting dihedral profiles calculated for fluorene with various lengths of side-chain. We observe that increasing side-chain length has very little effect on the energetics of the dihedral. However, in the case of thiophene (Figure (b)), the problem becomes considerably more complex due to the close proximity of the side-chain to the dihedral itself. In this case, the dihedral energetics are strongly influenced by steric repulsion. This is illustrated in Figure . Fluorene has side-chains much further from the center of dihedral rotation than thiophene. As such, the notion of separating long-range interactions from those due to the electronic conjugation is considerably more complicated for thiophene, while it does not pose any issues in fluorene.

Figure 5

Dihedral profiles for 2mers of (a) fluorene, (b) thiophene, and (c) ‘mirrored’ thiophene for varying lengths of side-chain(s). The side-chain(s) are labeled by a chemical formula e.g. H1 and C5H11 refer to no side-chain(s) and pentyl side-chain(s), respectively. Each graph has a schematic of the associated molecule to the right. (Legend in (a) applies to all graphs.)

Figure 6

Schematic highlighting the difference in the location of the alkyl side-chains in 2mers of fluorene and thiophene. The blue area represents the location of the intermonomer dihedral; the green represents the alkyl side-chain; and, in thiophene, the red area highlights the steric conflict involved.

The issue of this unavoidable steric repulsion creates a problem in interpreting the dihedral profile. As our parametrization is built on the ability to transfer parameters for long-range dispersive forces, the steric contribution to the dihedral profile due to the side-chains is captured by the existing terms in the force-field. As such, we wish to determine the effect of the side-chain on the covalently bound component of the dihedral rotation i.e. the effect it has on orbital conjugation as it is this component which comprises the required parameters in the force-field.

In light of this, we invoke an approximation to the separable case by taking a ‘mirrored’ thiophene 2mer. (See schematics in Figure .) In the mirrored 2mer, the side-chain of one unit is moved from the 3 carbon position to the 2 carbon position so as to remove the steric conflict. The nature of the approximation is that the effect of the side-chain on the conjugation-dependent component of the dihedral energetics is the same in both cases. Given the nature of the delocalized orbitals across the molecule as well as evidence suggesting that side-chains do not have a strong effect on conjugated phenomena (which can be seen in the agreement of the fundamental transition energies and transition densities shown in Section 1.5 of the Supporting Information), we conclude that this approximation is valid. Figure (c) depicts the energetic profile of the mirrored case, and it is found there is only a slight deviation (∼1 kJ/mol) due to the presence of side-chains. We attribute this deviation, in a similar manner to those seen in Figure , to small remaining dispersive forces between the mirrored side-chain and the neighboring unit. These results strongly suggest that any variance in the dihedral energetics on account of the inclusion of side-chains is mediated entirely by known long-range, noncovalent interactions. This suggests that the dihedral profile of a force-field for these molecules is invariant to the addition of side-chains though this statement is analyzed quantitatively in Section 5.

Partial Charge Calculations

We obtain partial charge by utilizing the RESP scheme[30] using the Antechamber program within the AmberTools 14 suite.[56] Our choice of this scheme is based both on consistency with the recommended methodology for OPLS parametrization as well as the robustness of the method to slight perturbations in geometry as opposed to, for example, ESP charges. Not only does this result in more accurate partial charges for utilization within an MD force-field, the robustness of the method also allows us to accurately determine variations due to the changes in molecular environment that we are interested in. Input electronic densities for the RESP calculations are obtained from SP calculations with CAM-B3LYP/cc-pVTZ using geometries taken from CAM-B3LYP/6-31G* optimization.

When generating parameters for MD simulations of repeating structures of many units, such as polymers or oligomers, partial charges are often implemented using a ‘three-residue’ model in which three sets of charges are built: one for each end unit and one for the central unit. In order to build such a model that is fully scalable to a variety of lengths, a first prerequisite is the requirement that the net charge for each individual central residue and the sum of the end residues have a net zero charge. This ensures the overall charge neutrality of the molecule and that this neutrality will remain when generating longer molecules by inserting extra central units. As it is possible for the end residues to have a nonzero net charge, if one wishes for a model with which scaling to many lengths is possible, it is necessary to determine the length of molecule at which the above criteria are met. Only for molecule lengths greater than this length will the application of a three-residue model be valid.

The remainder of this section provides an overview of the main results from our partial charge calculations and their implications. A more exhaustive breakdown of our partial charge analyses (such as the progression of total monomer charges and individual partial charges with varying backbone and side-chain length) is given in the Supporting Information (Section 2).

As an example, Figure presents the total charge of each monomer unit for 9mer of fluorene, 9,9-dioctylfluorene, thiophene, and 3-hexylthiophene. For the first three of these, the total charge for each monomer is close to zero across the entire molecule, whereas for hexylthiophene, there is a small but notable residual charge on the two end units. The key feature that distinguishes hexylthiophene from the other molecules is the breaking of reflection symmetry between each end of the molecule. At one end of the hexylthiophene molecule, the hexyl side-chain is on the side of the thiophene nearest the neighboring unit, while, at the other end, the side-chain is on the side further from the neighboring unit. This asymmetry of molecular structure forms an intrinsic asymmetry in the charge distribution of each end of the molecule i.e. a dipole exists. Given the small net charge difference (∼0.1 e) between each end, the dipole itself is not strong; however, the distribution of charges on each of the end units due to this asymmetry is notably distinct from other units. (See Supporting Information Figure S14.) Thus, if one wishes to build a three-residue model for hexylthiophene with accurate partial charges, the end residues must encompass two units of the molecule. In doing so, these two end residues will each have a nonzero net charge which cancel each other.

Figure 7

Total charge on each monomer unit of 9mers of fluorene, dioctylfluorene, thiophene, and hexylthiophene. The end monomers of hexylthiophene retain a considerably larger net charge than those of the other molecules.

To generate a model which is generalizable to all lengths of side-chain, two requirements must be met: First, there must exist a length of side-chain beyond which the charge distribution of the conjugated backbone is invariant; and, second, it must be possible to generate a standard methylene (CH2) group and a methyl (CH3) group which can be added in repeatedly to create longer structures. To analyze whether these requirements may be met, the charge distributions of a monomer with varying lengths of side chain have also been calculated. Figure displays the convergence of the total charge on the main conjugated component excluding the side chain. It is observed in both cases that such a convergence exists, and this convergence also implies that the charges on this part of the molecule are invariant. As such, the first requirement is satisfied. Analyzing the charge distributions of the side-chains in this case, it is observed that generating a standardized methylene group is not possible at the lengths of side chain we have considered (up to 10 carbons). This is because, while there are convergences in the end-most and innermost groups on a given side-chain, those in the middle of a long chain vary considerably as side-chain length is varied. An example of this is given for thiophene in Figure . It is seen here that the methylene group taken from the middle of the side-chain still fluctuates greatly in charge, primarily on the carbon atoms, even up to a 10 carbon long side-chain, while the two end groups have stabilized in charge. This is symptomatic of the strong asymmetry between each end of the side-chain - one is connected to a large molecule while the other is terminated only by a hydrogen. As such, if it were to be possible to generate a sufficiently accurate ‘standard’ methylene group charges for subsequent additions, a much longer side-chain would be required. As lengths of side chain greater than 12 are not often used in practical settings, it remains that if one wishes to include partial charges for use with side-chains, it is necessary to compute these charges for each length of side-chain required.

Figure 8

Relationship of the total charge, excluding side-chains, of monomers of dialkylfluorene and 3-alkylthiophene of various lengths of an alkyl chain. Each displays a convergence to near zero at around 6 carbons side-chain length.

Figure 9

Sample of beginning, middle, and end side-chain charges for various lengths of a side-chain in a monomer of alkylthiophene. It is observed that, in contrast to the end groups, the charge of the carbon on the inner methylene group (yellow triangles) does not converge with increased side-chain length.

Force-Field Implementation

Each molecular configuration is generated using pre-existing bond-stretching and angle-bending terms from the OPLS force-field. For both fluorene and thiophene, full sets exist for individual monomers, while only fluorene has parameters for the bonds and angles around the intermonomer junction. As they are lacking for oligo-thiophenes, the fluorene intermonomer parameters are used for the thiophene intermonomer junctions. Each molecule has improper dihedral terms at junctions with one heavy (i.e., non-hydrogen) atom joined to three other heavy atoms to re-enforce planarity across the monomer units as well as at the intermonomer junction. These may be more clearly observed from the Gromacs topology (.top) files provided along with the Supporting Information.

Partial charges are implemented in the force-field by means of generating a ‘three-residue’ model as described in the previous section and inputting these into the force-field parameter set. The dihedral potential requires a more careful implementation procedure. This is due to the pre-existing dispersive and electrostatic interactions. If one were to directly implement a dihedral potential, the presence of these interactions would lead to a double-counting of terms already included within the quantum chemistry calculation. In order to avoid this double-counting, it is necessary to isolate the effect of the extra interactions and subtract these from the dihedral profile.

The required force-field contribution to the dihedral potential is isolated by performing scans over intervals of 10° from 0° to 180° in a manner analogous to that of the scans performed using DFT. As we wish to isolate all interactions relevant in the dihedral rotation which are not the covalent interaction and also wish to restrain the dihedral at each value in the scan, the four covalent energetic functions at each intermonomer juncture are used to impose restraints.

In order to generate an effective restraint at a given angle, ϕ0, each of the four dihedral terms are placed under the influence of a periodic potential, V, given by

Care must be taken in choosing the value of k so as to find a balance between forming an effective restraint without inducing any unwanted distortion to the molecule. For molecules with methyl or no side-chains, the choice of k = 5 × 104 kJ/mol is suitable. In the case of ethylthiophene, a large reduction is necessary k = 103 kJ/mol which we attribute to the prevalence of large forces in the side-chain - dihedral area. From this, the geometry is then optimized in vacuum using the conjugate-gradients minimization algorithm within Gromacs 4.6.5,[34,35] and the total energy of each point along the scan is calculated to form the corresponding profile.

With the force-field (FF) contribution isolated, the required dihedral profile is obtained by subtracting the FF contribution from the DFT scan. The resulting ‘subtracted’ profile is then fitted to a fifth-order Ryckaert-Bellmans function

The fit function described in eq yields two sets of parameters. This results from the difference of 180° between one pair of dihedral angles and another of the four used. For example, for a dihedral angle of 0° in the polymer convention, the dihedral angle of the two pairs of four atoms in the trans position is ϕ°, while the two pairs in the cis position have a (ϕ + 180)° dihedral angle. As such, the function cosine terms in eq must be modified to cos(ϕ + 180) = −cos(ϕ) in order to yield the appropriate energy.

Figure provides examples of the curves obtained in the subtraction process for 2mers of fluorene, thiophene, and ethylthiophene. In the cases with no side-chains (Figure (a) and (b)), it can be seen that this procedure and the fitting with the Ryckaert-Bellemans function result in a force-field which quantitatively mimics the DFT dihedral potential to a very high degree of accuracy.

Figure 10

Subtraction profiles for (a) fluorene, (b) thiophene, and (c) ethylthiophene. Each figure displays the calculated DFT profile; the profile obtained from the FF ‘scan’; the resulting subtracted profile; the fit of the subtracted profile to a fifth-order Ryckaert-Bellmans function; and the resulting ’effective’ profile given by the addition of the FF scan profile and the fitted profile. (Legend applies to all graphs.)

It is noted that this procedure results in a far less smooth fit for ethylthiophene (Figure (c)). As discussed with respect to the DFT scans, the introduction of side-chains is a potential source of inconsistency in the calculations due to the relative freedom of side-chains. To minimize this inconsistency, the starting geometries are taken from the DFT calculations so as to reproduce the side-chain conformations as well as possible. However, this does not, in all cases, lead to a perfect correspondence in the side-chain conformations between the DFT and FF scans.

As is shown in Figure , the subtraction procedure leads to an overall dihedral potential from the force-field which corresponds closely to that of the DFT. To further test this correspondence, we have calculated the difference in energy, Δ, between the two local-minima (one in the first quadrant of the dihedral and one in the second) of each molecule from DFT and our force-field and compare them in Table . In all instances, we find agreement to within 0.2 kJ/mol (∼0.08 RT at room temperature).

Table 1

Comparison in the Difference in Energy between the Cis and Trans Minima, Δ, of 2mers of Fluorene (F0), Methylfluorene (F1), Thiophene (T0), and Ethylthiophene (T1) Calculated from DFT and from the Force-Field (FF)

	ΔEm (kJ/mol)
	DFT	FF
F0	0.04	0.08
F1	–0.05	0.05
T0	2.23	2.09
T2	0.84	0.86

In order to achieve the agreement between the DFT and FF values of δE, we modified the equilibrium bond lengths and angles from the OPLS force-field while retaining the orginal force-constants. This measure gives a better agreement in the geometries and minimal energies obtained from the force-field. A full discussion of the implications of this modification is given in Section 3 of the Supporting Information.

As expected from the DFT results for alkylfluorenes (Figure (a)), addition of methyl side-chains to fluorene has no effect on the energies of the dihedral rotation. As such, the dihedral potential needs no further modification to accommodate for this. In Section 3, we argued that the steric interactions responsible for large changes in dihedral potential for alkylthiophenes should be incorporated by the force-field. However, as is shown in Table , utilizing the dihedral potential fitted from a thiophene with a shorter side-chain leads to drastically overestimated barriers at the planar positions. As such, reparametrization of the dihedral term must be performed to accommodate for this. Given the tendency for inconsistency observed in ethylthiophene due to the side-chain degrees of freedom and that the difference in DFT dihedral potential between ethyl- and propylthiophene is slight (∼0.5 kJ/mol at the planar barriers, Figure (b)), we use the ethylthiophene potential for thiophene molecules with longer side-chains.

Table 2

Height of the Energetic Barriers at 0° (E0) and 180° (E180) for 2mers Ethylthiophene, Methylthiophene, and Thiophene Using Fitted Profiles Obtained from Scans Using Different Side-Chain Lengthsa

	ΔE0 (kJ/mol)		ΔE180 (kJ/mol)
	DFT	FF	DFT	FF
T0 (T0)	0.59	0.33	2.02	1.75
T1 (T0)	1.20	5.35		7.21
T1 (T1)	1.20	1.33	4.65	5.00
T2 (T0)	4.02	11.07	6.71	12.00
T2 (T1)	4.02	6.67	6.71	9.18
T2 (T2)	4.02	3.90	6.71	7.66

Each barrier is calculated relative to the closest local minimum (i.e. the trans minimum for Δ0 and the cis minimum for Δ180). The labels Tx (Ty) denote the energies of a 2mer with an x-yl side-chain with energetic profile taken from a y-yl side-chain scan. The DFT values shown are those from the dihedral scan of the Tx molecule.

The force-field parameters obtained from the above procedure are available, in a form suitable for use with the Gromacs MD package, along with the Supporting Information.

Molecular Dynamics Results

All MD simulations were carried out using Gromacs 4.6.5[34,35] with an integrator step-size of 2 fs and system coordinates sampled every 10 ps. Each simulation was performed at ambient temperature and pressure (273.15 K, 1.01325 bar). In all cases, the following measures have been taken prior to the MD run: steepest-descent minimization; followed by 0.5 ns of both NVT and NPT ensemble equilibration with position restraints on the heavy atoms of the molecule followed by 5 ns of unrestrained NPT equilibration.

We perform our calculations in a fully solvated manner i.e. in a periodic cubic box large enough to avoid any possible interactions of periodic images. In the largest molecules we have looked at (32 units in length), this condition leads to the majority of computational effort being spent on calculating the solvent dynamics. Each simulation was performed with chloroform as a solvent, and solvent parameters were acquired from the molecule database at virtualchemistry.org.[57,58]

Long-range electrostatics in all simulations have been treated with Reaction-Field with a dielectric constant of 4.81 for chloroform. We have found from preliminary simulations that this choice gives results comparable to those obtained using the Particle-Mesh Ewald method with a significant reduction in the simulation run-time. Full details of our tests are provided in Section 4 of the Supporting Information.

As a test of our parameter sets, we performed simulations of 32mers of fluorene with octyl side-chains (PF8) and thiophene with hexyl side-chains (P3HT) over the course of 50 and 100 ns, respectively, in chloroform and calculated persistence lengths, n (in number of monomer units) and l (in nm). We perform this calculation by generating vectors, v, across the first and last carbon of unit i and generating a correlation function, A(θ), of the angles between each pair of vectors, θ (eq )

The persistence length[59,60] is defined by the e–1-point of A(θ) i.e.

The correlation curves for fluorene and thiophene with their respective side-chains are shown in Figure . By fitting these to the exponential decay function given in eq , the persistence lengths are n = 12.9 and n = 8.5 units, respectively. In order to recast these into nm, we calculate the average unit length, l, from the MD simulation to be l = 0.832 nm for PF8 and l = 0.397 nm for P3HT which gives l = 10.8 nm and l = 3.4 nm, respectively. As a means of comparison, persistence lengths for the polymer in each case have been reported as l = 8 ± 1 nm for PF8[37,38] and l = 2.4 ± 0.3 nm for P3HT;[39] both measured by a combination of gel permeation chromatography and light scattering in THF solution.

Figure 11

Angle correlation functions, A(θ), of 32mers of dioctylfluorene (F8) and hexylthiophene (T6) simulated in chloroform. Each calculated correlation function is given along with their respective fitted curves (dashed lines) (given by eq ). From the fitted curves, n, indicated by arrows corresponding to the crossing of the fitted curve and the e–1 line, is found to be 12.9 and 8.5, respectively.

The remarkable agreement found between the MD and measured persistence lengths should not be overplayed and is, to an extent, deceptive. Our simulations are performed in a different solvent (chloroform) from the experiments (THF). We make this choice as simulations with THF are ∼5 times more expensive computationally than those with chloroform. As we are performing our calculations in a fully solvated manner i.e. in a periodic cubic box large enough to avoid any possible interactions of periodic images, this means we have ∼2 × 105 and 2 × 104 solvent molecules per simulation for fluorene and thiophene, respectively. While we are currently working on means of simulating fully solvated molecules in a more efficient manner, as it currently stands, this method restricts us computationally to smaller solvent molecules such as chloroform. Furthermore, we have found it is necessary to use molecules considerably larger than the persistence length in order to obtain a converged result. Particularly, performing a similar simulation on a 16mer of fluorene resulted in the value l = 19.5 units which corresponds to 16.24 nm. This values is ∼2 times that of the 32mer and of the polymer experiment. Similar behavior was observed for thiophenes of near persistence-length lengths although this was to a considerably lesser extent (e.g., l = 9.5 units for a 5mer).

Distributions of end-to-end length have also been computed from our simulations. In Figure , these distributions are given for 16mers of dioctylfluorene and 16mers and 32mers of hexylthiophene in chloroform. The distributions for 32mers of dioctylfluorene have not been shown as, while the tangent correlation curves were converged, after 50 ns the end-to-end length distributions were far from converged. In each example, the end-to-end length is given as a fraction, r, of the straight length of the molecule, Nl, where N is the number of units and l is the mean unit length given above.

Figure 12

End-to-end length distributions (solid lines) for a 16mer of dioctylfluorene (F8) and a 16mer and 32mer hexylthiophene (T6) in chloroform. The end-to-end length is scaled to give each length as a fraction of the fully extended length of each molecule. Comparison is given for each distribution to the distribution r2Pr obtained using eq (dashed lines)[61,62] with a curve for an F8 32mer given based on the calculated persistence length.

In comparing each length of thiophene, it can be seen that the 32mer has a much wider distribution which peaks at a lower length fraction (0.7) than that of the 16mer which is narrower and peaked at ∼0.85. Given the shorter persistence length calculated, the distribution for the 16mer of fluorene is less spread and peaked at a higher length fraction ∼0.9 than the 16mer and 32mer of thiophene. The distribution obtained for the 16mer of fluorene is consistent with measurements performed by Muls et al.[63] Using end-marked hexylfluorenes of a similar length scale (∼42 monomer units with a polydispersity of 1.8) in an inert Zeonex matrix, they measured end-to-end length distributions which are centered at a length fraction (based on the fully extended 42mer) of ∼0.89.

Both the progression of the distributions and the difference in overall spread between fluorene and largely spread thiophene can be understood qualitatively by considering the increase in conformational entropy with increasing length and is consistent with the persistence lengths calculated previously. For a quantitative insight, we have also calculated the end-to-end length distributions, P(r), using the expression derived by a path-integral approach for semiflexible polymers by Wilhelm and Frey[61,62]where r is the end-to-end length fraction, ξ = n/N is the persistence length fraction, and is the normalization constant such that ∫01drr2P(r) = 1. The values of ξ used are ξ = 1.219, 0.403, 0.531, and 0.266 for the 16mer and 32mer of fluorene and the 16mer and 32mer of thiophene, respectively. For each fluorene, the corresponding persistence length was used due to the large difference for both lengths (as described above). The value calculated for the 32mer was used for both thiophenes. For all the distributions calculated, each MD distribution agrees well with the corresponding calculated result, and we expect that this behavior would be replicated by a fully converged fluorene 32mer distribution.

A Parametrization Example

In order to illustrate how our scheme would be applied in practice, this section provides a working example of the calculations required to parametrize the copolymer 9,9-dioctylfluorene-alt-benzothiadiazole, commonly known as F8BT (Figure ). The example is taken as an instance where one wishes to perform MD simulations on an F8BT molecule of greater than 8 units in length. Before carrying out the main body of the scheme suitable parameters must be obtained (either directly from the OPLS parameter set or from other sources) for the basic force terms (e.g., bond-stretching, L-J terms). One must then generate the lists of atoms and bound pairs within the molecule based on the OPLS atom naming specifications (which are intrinsically tied to L-J parameters).

Figure 13

Base monomer of the copolymer 9,9-dioctylfluorene-alt-benzothiadiazole (F8BT). The carbons circled in red are an example of where a substitution of a covalent force term (in this case, bond-stretching between two benzene carbons) may need to be made.

The next stage is to generate lists of the covalently bound force terms in the molecule. These are often generated by the simulation suite (such as is the case in Gromacs 4.6.5). The bonds, angles, and dihedrals which do not have a direct specification must be input from choices of similar bonds. For example, the carbon–carbon bond between the two central atoms of the BT unit (marked in Figure ) could be taken from the parameters for bond-stretching between two carbons of a benzene ring. All internal dihedrals which are not specified are to be replaced with the dihedral potential for 4 aromatic carbons. This choice works in general as this potential serves primarily to ensure the rigidity of the monomeric structure.

The first stage of the quantum chemistry calculations is obtaining a dihedral potential by performing a scan of the rotation between the fluorene (F8) and benzothiadiazole (BT) subunits of the monomer. As shown in Section 3, the inclusion of side-chains is unnecessary and may lead to complications. With this in mind, the octyl side-chains of the F8 subunit should be replaced with hydrogens and the scan performed on the monomer shown in Figure (hereafter referred to as the F0BT monomer). The scan is performed using CAM-B3LYP/6-31G* to obtain geometries and CAM-B3LYP/cc-pVTZ for corrected energies by SP calculation.

Figure 14

Base monomer of the copolymer fluorene-alt-benzothiadiazole (F0BT) which would be used for calculations of dihedral profile and in the subtraction procedure.

To obtain charges, one would take an 8mer and 9mer of the F8BT molecule optimized with CAM-B3LYP/6-31G* and electron density calculated by SP calculation using CAM-B3LYP/cc-pVTZ. Partial charges would then be calculated using the RESP scheme. If there is no notable difference in the total charges of the monomers between 8mer and 9mer, then the charges are converged and will be suitable to build charge sets. If not, longer molecules would need to be considered. The above also applies if one wishes to match the end units (e.g., put an extra F8 or BT on either end) to symmetrize it. As discussed, the convergence of the charges would occur at smaller lengths than in standard F8BT. However, it is likely that using an 8mer and 9mer will be sufficient.

The final step is performing the subtraction using the MD algorithm. This would be performed with the F0BT monomer in the manner described in Section 5. This concludes the parametrization of the F8BT molecule by our scheme.

Conclusions

We have determined that it is possible to replicate dihedral potentials of high levels of theory using an economical, two-step method of CAM-B3LYP/6-31G* and CAM-B3LYP/cc-pVTZ for geometry optimization scans and single-point energy calculations, respectively. We find that dihedral potentials are invariant to the length of the molecular conjugated backbone and that any effect due to the inclusion of side-chains seems to stem from long-range interactions alone.

Furthermore, we find that partial charge distributions converge quickly with varying length of molecule, but one must take into consideration deviations in end monomer distributions for molecules which do not possess end-to-end reflection symmetry. In terms of varying side-chains, it seems that there is no simple way of generalizing the inclusion of partial charges to this variation.

With these points in mind, we have developed a parametrization scheme which is applicable to a wide range of conjugated molecules. We have tested the transferability of parameters expected in molecules with high steric contributions to the dihedral energetics from close side-chains (such as in 3-alkylthiophenes) and found that delicate reparametrization is necessary, while, in molecules with side-chains far from the intermonomer junction, this careful approach is not required. Our preliminary tests indicate that our MD simulations are consistent with experimental measurements of persistence lengths and end-to-end lengths with long thiophene oligomers displaying far greater flexibility than their fluorene counterparts.