Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
J Mol Graph Model. Author manuscript; available in PMC 2010 May 19.
Published in final edited form as:
PMCID: PMC2873194

Development and initial testing of an empirical forcefield for simulation of poly(alkylthiophenes)


Conductive polymers from the polythiophene (PT) family have attracted interest in numerous domains, including potential applications in biosensing. Despite this, atomistic simulations of PTs have tended to use general organic force fields without well-tuned PT parameters, and there exists no optimized and well-validated PT force field that is compatible and consistent with existing biomolecular simulation suites. We present here the development of a new PT forcefield following the AMBER approach, using the program ANTECHAMBER and ab initio calculations at the HF/6-31G* level of theory to assign partial charges and parameterize the critical backbone torsion potential. The optimized geometries and force field potentials match well with both empirical data and previous investigators' calculations. Initial testing of these parameters through a series of replica exchange simulations of two PT derivatives in aqueous and organic implicit solvents demonstrates that the parameters can match empirical expectations within the limits of an implicit solvent model. This new force field forms a framework for modeling of proposed PT-based devices and sensors, and is expected to accelerate device design and eventual deployment.

Keywords: Polythiophene, Conductive polymer, Molecular dynamics, Force field, Implicit solvent

1. Introduction

Polythiophenes (PTs) are well-known as highly conductive polymers that can be synthesized in a regioregular fashion and functionalized with a variety of side chains and end groups [13]. Because of these properties, many investigators have proposed the use of PTs as sensors to detect volatile organics [4] and biomolecules such as oligonucleotides [5]. In support of that line of research, we have recently demonstrated that a regioregular poly(alkylthiophene) (PAT) can form biocompatible self-assembled monolayers that have favorable properties as biomedical electrodes [6]. As the use of PTs in sensing and nanotechnological devices expands, there will be an increasing need to model devices by simulation, and particularly by atomistic simulation, prior to prototype construction.

Computational chemistry methods in general have been applied to PTs for nearly two decades [7,8]. Much of the recent work in this arena has focused on the use of quantum-level calculations to study conductivity and optical properties of PTs [5]; earlier work concerned itself more with structure and torsional potentials, and is referenced further below. However, molecular dynamics (MD) simulations of PTs are not unknown. Most investigators have used general organic force fields, and have achieved some success in multiple domains, including predicting solid-state structure [911], investigating conformations in solution [1214], and predicting the results of spectroscopic studies [1517]. One group has hand-built a PT-only force field that incorporates terms from a variety of the literature sources [18,19], and has used these parameters successfully to study phenomena related to PT crystal structure and the migration of dopant ions within films [2026]. These studies notwithstanding, it would be valuable to have a set of optimized PT parameters that are appropriately tuned for the charge models and underlying assumptions of common biomolecule-focused molecular dynamics engines such as AMBER [27], GROMACS [28], or CHARMM [29], and to have evidence that those parameters can predict empirical PT properties.

The work presented here addresses that need by developing, from first principles, a set of PT parameters that are consistent with the approach used to develop the current-generation AMBER force fields. Those parameters are then used for a series of replica exchange simulations of the conformation of PATs in implicit solvents. We have previously found that the addition of a branched alkyl side chain to a PT appears to confer solubility in a saturated hydrocarbon solvent, and that this solubility is lacking in PATs with unbranched side chains [30]. The implicit solvent simulations in this work produce results consistent with that finding. The same empirical data have been more thoroughly replicated in explicit-solvent simulations using these same parameters, and the results of those simulations are reported elsewhere [31,32]. In the following sections, we will detail the calculation of key parameters by ab initio methods, the translation/fitting of those parameters to AMBER functional forms, and the results of the initial replica exchange simulations. This new simulation framework, taken together with existing parameter sets for proteins, nucleic acids, and phospholipids, will allow the modeling of a wide range of hypothetical PT applications in the biological and chemical domains. One example of such modeling is our already-published work on the interaction of PATs with lipid bilayers, which has shown the feasibility of a new type of biosensor [31].

2. Methodology

2.1. Geometry and partial charges

Starting geometries, partial charges, and AMBER force field parameters were derived for undecorated PT and for two PAT derivatives: poly(3-hexylthiophene) (HPT), and poly(3-(2-ethylhexyl)thiophene) (EHPT). As noted above, EHPT's branched side chain gives it greater solubility in common organic solvents and greater disorder in the solid phase, whereas HPT is overall more ordered and less soluble. These differential solubilities can then be assessed in implicit solvents as a preliminary test of the force field's effectiveness.

The parameter development process was aided in part by the ANTECHAMBER component of the AMBER suite. Initial geometries were obtained by optimizing pentamers of plain PT, EHPT, and HPT at the Hartree-Fock level of theory with the 6-31G* basis set (HF/6-31G*). This and all other ab initio calculations described here were performed with GAUSSIAN03 [33]. Only the nearly planar all-anti conformation was optimized for each pentamer, based on others' findings that bond lengths and other parameters do not significantly change in other near-minimum conformations [34]. A pentamer is expected to be more than sufficient to capture the geometry and electron structure of the general polymer, as previous ab initio studies have found that geometries do not change substantially once the length exceeds three rings [35,36]. Partial charges were then derived from the HF/6-31G* electrostatic potential grid using the restrained electrostatic potential (RESP) method and program of the AMBER suite [37]. Finally, the optimized and charge-assigned pentamers were split into “N-terminal”, “C-terminal”, and “central” (derived from the center ring) monomer residues. Where necessary, charges were rounded off by hand to ensure that each monomeric unit had zero net charge.

Other investigators seeking to develop PT force fields have used higher levels of theory for their ab initio calculations, up to MP2/aug-cc-pVTZ [13]. However, we believe that HF/6-31G* is appropriate for this particular application. Perhaps the most important reason is that existing AMBER force fields were parameterized at this level of theory [38], and consistency with the existing parameters is a primary goal of our work. However, another deciding factor is that HF/6-31G* optimizations and rotations have been shown by other investigators to match well with experimental results [39], and the use of higher levels of theory (e.g., MP2) does not substantially change the relative energies of different conformations [40,41].

2.2. Force field atom types and parameters

The majority of atom types in a PT and the non-bonded interactions between them are well-described by the existing generalized AMBER force field (GAFF). While the generalized aromatic parameters are not optimized for PTs, they are sufficient to accurately reproduce the molecular geometry and potentials, as will be shown below. The exception is the potential function describing the energy of torsion about the bond between two PT rings. This inter-ring torsion is known to be the most flexible degree of freedom among the polymer backbone atoms [35,34], and is the predominant determining factor for the polymer's conformational distribution in free-rotating gas and liquid phases [42]. In some cases, it has even been possible to predict oligomer structure solely from the potential of rotation about this bond [34]. Therefore, a new atom type and corresponding parameters must be added to GAFF to accurately model this critical potential.

All four carbons of the thiophene ring were initially assigned to the “ca” generic aromatic carbon type of the GAFF force field. To facilitate representation of the new torsional terms about the PT inter-ring bond, a new atom type “cp” was defined for the C2 and C5 carbons of each ring (see Fig. 1 for nomenclature). The majority of forcefield parameters for this new “cp” atom type were identical to the GAFF parameters for “ca”, but the equilibrium bond lengths and angles were altered to better match the geometry found from the HF/6-31G* optimizations.

Fig. 1
Diagram of ethylhexyl (EHPT, left) and hexyl (HPT, right) PT dimers to illustrate the ring and side chain atom nomenclature used in this paper. To enhance the clarity of the drawing, hydrogen atoms of the side chains are not shown. These side chain hydrogens ...

To properly parameterize the inter-ring torsion, the rotational potential of bithiophene (with no alkyl side chains) was calculated. The dihedral angle between the two rings (S1–C2–C5′–S1′ or C3–C2–C5′–C4′ in Fig. 1) was rotated over a full 360° in increments of 5°, with no geometry relaxation after the initial structure optimization. (For purposes of this work, the minimum-energy anti configuration is considered as 180°.) This “rigid rotor” method was used to facilitate the next step of matching the torsional potential between AMBER and the ab initio results. It should be noted that other investigators have found the use of a “relaxed rotor” in bithiophene to cause no substantial change in either the bond or angle parameters at each step or the observed energy profile [43,44]. The full 360° rotation (as opposed to the more commonly seen 180° rotation from planar anti to planar syn conformation) was undertaken based on previous reports that the rotor potential is not completely symmetric [45]. This appears to be the first instance where the bithiophene rotor potential has been calculated and reported in such detail; all previous reports used an increment of 10° or greater, and only show the potential for a 180° rotation.

Once the ab initio rotor potential was obtained, Fourier series decomposition was used to fit it to a series of cosine terms appropriate for use in an AMBER-compatible forcefield. Such decompositions have previously been shown to have good success in fitting thiophene rotor potentials with only a few terms [40,46,47].

The ab initio potential was decomposed into a set of terms in the AMBER-appropriate form:


for n = 1–4, using both cosines (ϕ = 180°) and sines (ϕ = 270°) for n = 3 and n =4. In this equation, Kn is the spring constant (magnitude of the Fourier coefficient), and θ represents the torsional angle. The obtained coefficients were then split between the S1–C2– C5′– S1′ and C3–C2– C5′– C4′ dihedral angles in the forcefield file, using four terms to describe each; the sine terms were assigned to the S1–C2– C5′– S1′ angle and the cosines to C3′C2′ C5′– C4′. For the exact coefficients and details of how to assign them to the dihedral terms, interested readers are referred to the parameter files available as supplemental online material.

In addition to the all-atom model just described, a united-atom version of this force field was also constructed, in anticipation of our use of these parameters with a united-atom lipid bilayer [31,32]. The hydrogens of the alkyl side chains (but not the ring backbone) were unified with their carbons, summing the partial charges together. Bond stretching, angle bending, torsion, and Van der Waals (VDW) parameters for the new united atoms were taken from the well-known Ryckaert-Bellemans potential for united alkanes [48]. Because these united-atom parameters require a different scaling of 1–4 VDW forces (a factor of 8 instead of the usual 2 for the AMBER forcefield) than do the standard AMBER parameters, the Fourier decomposition of the previous paragraph was repeated to properly parameterize the united-atom model.

2.3. Replica exchange MD in implicit solvents

Molecular dynamics simulations of the two PATs of interest were carried out using these newly derived parameters. EHPT, with its branched side chain, is known to show greater inherent disorder and behave differently from HPT in good and poor solvents. The chosen test case was therefore a comparison of the behavior of the two polymers in water (bad solvent) and hydrocarbon (presumed good solvent). These simulations were carried out using the revised generalized Born implicit solvent model described by Onufriev et al. [49], as implemented in AMBER 9. EHPT and HPT decamers were each simulated in implicit water (external dielectric constant 78.5) and a generic implicit saturated hydrocarbon (external dielectric constant 2.0). No surface area correction terms were used. It should be noted that simply changing the external dielectric constant does not fully capture the difference between the two solvents, as the underlying Born radii can also be solvent-dependent and the default AMBER parameters are designed for the aqueous case. Moreover, since neither model PT is water-soluble and EHPT is only somewhat hydrocarbon-soluble, we would expect only small differences between the runs. Nevertheless, this dielectric-only technique has recently been successfully used to compare the differential behavior of small peptides [50] and collagen [51] in aqueous and organic solvent. We therefore feel that it provides an adequate “sanity check” on the new parameters before they are deployed in more realistic simulations.

To accelerate configurational sampling, all four implicit solvent MD simulations were performed in a replica exchange ensemble. Replica exchange, also called parallel tempering, is a descendant of the simulated annealing method. Multiple copies (replicas) of the same system run in parallel at different temperatures. By occasionally exchanging system temperatures, each replica effectively undergoes several annealing cycles, and the simulation is more readily able to sample the full range of configuration space. For precise mathematical details of the setup and analysis of such simulations, readers are urged to consult one of the many excellent papers on the topic, including [52,53].The simulations described here used the built-in replica exchange capabilities of AMBER 9.

Twelve replicas were used for each of the four simulations, with their temperatures selected according to the optimal spacing algorithm of Rathore et al. [54]. The temperature span for EHPT ran from 300 K to 709.27 K, while HPT ran from 300 K to 744.0 K, with all temperatures maintained by a Langevin thermostat with a coupling constant of 0.5 ps−1. Configuration information was written and exchanges attempted every 1 ps of simulated time. Representative samples of the raster and temperature history plots for the replica exchange simulations are shown in Fig. 2, and demonstrate that the simulations performed an adequate random walk in temperature space over the 5 ns of simulated time. The specific exchange rates varied from 20.8% to 38.6%. The latter value is slightly higher than optimal [55], but is not so high as to have compromised the quality of the simulation, especially given the shorter relaxation times of implicit solvent systems. The SHAKE algorithm was used to constrain bonds involving hydrogen, and an integration timestep of 1 fs was used due to the elevated temperatures applied to some replicas. The RESPA algorithm was used to only evaluate slowly varying forces at every other timestep. A non-bonded interaction cutoff of 20 Å was used, which allowed the side chain of each monomer unit to “see” as far as the side chains of its neighbors.

Fig. 2
Representative examples of raster plot (A, one dot per successful exchange, 500 ps sample) and temperature history plot (B, two randomly selected replicas, full 5 ns run) for the replica exchange implicit solvent simulations. All replicas show regular ...

Each replica-exchange simulation was run for 5 ns, giving a total simulation time of 240 ns. The output (trajectories, energies, and temperature histories) from these simulations were then used to calculate potentials of mean force (PMFs) using the weighted histogram analysis method (WHAM). All WHAM calculations and PMF error estimation used the methodology and code of Chodera et al. [56].

3. Results

3.1. Geometry and partial charges

Table 1 compares the ring geometry obtained from other groups' experimental and ab initio theoretical studies of oligothiophene to the geometry obtained from both the HF/6-31G* pentathiophene optimization and AMBER energy minimizations using our new force field. The agreement is generally good, with the root-mean-squared (RMS) error between empirical atom-atom distances and these calculated ab initio distances being 0.024 Å. RMS error between empirical and ab initio bending angles is 4.05° when three-atom bond angles (bending motions) about the interring C2– C5′ bond are included in the calculation, or 0.71° if they are not. This is not a surprising result; the empirical geometries for thiophene and terthiophene were measured in the solid phase, where intermolecular packing interactions would be expected to affect the bending and torsion angles about the inter-ring bond (the most flexible degrees of freedom). By contrast, the calculations reported here implicitly assume a gas phase molecule in vacuum, where those same parameters would be expected to be more relaxed.

Table 1
Optimized geometry of the thiophene ring in this work vs. selected studies

Because the torsion (and associated bending) about the interring bond is relatively free in polythiophene and is also sensitive to the molecule's environment, it is important to compare the values obtained in these calculations with empirical results that were also obtained in the gas phase. Two electron diffraction studies of gas-phase bithiophene have been carried out; the first reported an inter-ring dihedral angle of 146° [57], while the second reported 148 ± 3° [39]. The values reported here show good agreement with those experiments; the minimum of the HF/6-31G* rigid-rotor potential lies at 145.93°, and the same dimer minimized with the derived AMBER forcefield has a dihedral angle of 145.07° (Table 1). The center rings of pentamer structures from ab initio and AMBER minimizations show slightly more planar angles, at 150.15° and 151.56°, respectively. These are still within the experimental margin of error. Moreover, a pentamer should be more planarized than a dimer, since the increasing length should allow greater electron delocalization along the backbone, which will lower the energy of planar conformations [58].

It is also valuable to consider the effect of the side chains on the inter-ring dihedral angles. There are no known data on the gas-phase geometry of alkylated oligothiophenes. The closest match is a dynamic NMR study which found that alkylated dimers interconverted rapidly between syn-like and anti-like forms even at 140 K, and thus inferred that the barrier between configurations is very low [41]. Ab initio studies at levels of theory ranging from HF/3-21G* to HF/6-31G** also suggest a low barrier to ring rotation in oligo(alkylthiophenes), and that the addition of side chains should favor a geometry that is further deplanarized compared to plain oligothiophene [45,59,43]. These studies further indicate that longer alkyl chains should produce more lowering of the barrier and more deplanarization [45,59,41].

The results of adding all-atom or united-atom ethylhexyl or hexyl chains to pentathiophenes are shown in Table 2, and follow the same trend seen in others' calculations. The optimal inter-ring dihedral angles in HF/6-31G* optimizations of alkylated pentamers are slightly higher (more planar) than those observed in previous studies of dimers, possibly due once again to the increased electron delocalization of longer oligomers. AMBER minimizations (using the default steepest descent/conjugate gradient method) of those same molecules show substantially less deplanarization than the HF optimizations. There is no definitive explanation for this observation at this time. However, given what has just been noted about the broad, flat-bottomed potential well for alkylated bithiophene torsion (in contrast to the well-defined minima of the plain bithiophene rotor potential; see below), it is possible that the AMBER minimizations became trapped in local minima along the bottom of this potential well and were not able to reach the global minimum.

Table 2
S1–C2–C5′–S1′ dihedral angles for oligo(alkylthiophenes) in this work vs. others

Partial charges derived with the RESP procedure are shown in Table 3, and are generally unremarkable. As expected, the molecules are overall apolar, with the carbons and sulfurs of the thiophene rings showing higher partial charges (greater electron delocalization and polarizability) than the aliphatic side chains.

Table 3
Partial charges for ethylhexylthiophene and hexylthiophene monomers

3.2. AMBER fitting of ab initio torsion potential

Fig. 3 shows the HF/6-31G* rigid rotor potential for bithiophene, and compares it to the AMBER rotor potentials obtained using the all-atom and united-atom derived forcefields. The match is excellent, with the RMS error between AMBER and ab initio energies being 0.047 kcal/mol for the all-atom fit and 0.063 kcal/mol for the united-atom fit. The shape of the potential is also consistent with that observed by other groups using a self-consistent field approach and similar basis sets, and using both rigid and relaxed rotations [39,43,40,41]. The potential shows two sets of minima, corresponding to syn-like (0°) and anti-like (180°) configurations. As would be expected, the syn-like conformation has higher energy due to the close approach of the two sulfur atoms. The syn- and anti-like minima differ in energy by only 0.64 kcal/mol in the best case, and the barriers between them are relatively small, implying that PTs at equilibrium in solution should show a mixture of syn- and anti-backbone angles. This is consistent with prior NMR and gas-phase electron diffraction studies showing that bithiophene exists in a mixture between syn and anti conformations, with the anti conformation slightly favored (by 8 ± 4% in the gas phase) [60,41,39].

Fig. 3
Rigid rotor potential of plain bithiophene, as computed by ab initio HF/6-31G* (circles), compared to the same rigid rotor potential computed in AMBER using the all-atom (SCNB = 2, squares) and united-atom (SCNB = 8, diamonds) parameterizations. RMS error ...

The energy barriers in this rigid rotor potential are compared in Table 4 to those observed by other groups in experimental and computational studies. The height of the highest barrier has been quantified by NMR in a liquid crystal solvent as 5 ± 2 kcal/mol [61]; the value of approximately 2.2 kcal/mol observed in this work is slightly outside that range, but one would expect the rotational barrier to be lower in gas phase than in a highly ordered solvent. For the energy difference between the anti-like and syn-like conformers, there is little agreement in the literature–measurements with different techniques have ranged over almost an order of magnitude, from 0.18 to 1.16 kcal/mol [61,39,62]. Four values for this gap can be extracted from each curve in Fig. 3, depending on which minima one chooses to compare. All of the calculated values fall within the broad experimental range.

Table 4
Relative energies from the rigid rotor potentials of Fig. 3 vs. selected studies

The agreement between the bithiophene rotor potential found in this work and that observed in other investigators' calculations is qualitatively good, but the actual energy values differ slightly, as also seen in Table 4. Our calculated maximum barrier height of 2.2 kcal/mol agrees better with the HF/6-31G** calculations of Hernandez and Navarrete [43] than with the slightly lower values of prior HF/6-31G* calculations [39,40,46]. The gap between syn and anti minima reported here could be either higher or lower than previous ab initio calculations, depending on which of the four values one chooses. The majority of these differences are likely attributable to the use of different calculation techniques. Previous investigations have used relaxed-rotor potential calculations, rotating only over 180°, with a minimum spacing of 10° between sample points. By contrast, the potentials of Fig. 3 are rigid rotations over a full 360° with 5° spacing. This higher resolution reveals the aforementioned asymmetry of the rotor potential and may be able to sample deeper minima, while the lack of atomic position relaxation would be expected to produce higher energies at the peaks, where atoms are undesirably close together.

3.3. Replica exchange MD in implicit solvent

Through replica exchange simulations of EHPT and HPT, we sought to determine whether the new force field can produce results consistent with our prior experimental observations of different solubilities for these two molecules. The end-to-end distance of the simulated decamers is used here as a proxy for solubility. PATs are well-known to form aggregates, particularly in poor solvents [1,63], and it has been suggested that the first step in this process is the collapse of an individual chain upon itself to form “intramolecular aggregates” [64,14]. Such collapse brings the two ends of the molecule closer together, and thus the end-to-end distance is a rough estimate of the solvent's ability to solvate the polymer under study.

Fig. 4 shows the potentials of mean force (PMFs) for EHPT and HPT in implicit water and hydrocarbon. The PMF can be considered as the free energy profile of the system across an arbitrary reaction coordinate, in this case the end-to-end distance of the simulated decamer. Numeric values for the PMF minima are given in Table 5. Both PATs show a lower energy in water than in hydrocarbon for conformers with lower end-to-end distance. This is consistent with aggregation, as would be expected for hydrophobic molecules in an aqueous environment. The other notable feature of these PMFs is the difference in response when the solvent is changed from implicit water to implicit hydrocarbon. Both molecules show a shift towards more extended conformations, implying that hydrocarbon is a more favorable solvent than water, but EHPT shifts by approximately 1.5 Å, whereas HPT shifts by only 0.5 Å. The observed differences are small, and the curves of Fig. 4 largely overlap. Given this, only limited conclusions can be drawn from these results. Nevertheless, the observed trend does match our prior empirical observations that EHPT is more “at home” in alkane-like environments than HPT is. Moreover, the difference is not expected to be large –although EHPT does dissolve in saturated hydrocarbons, they are at best a weak solvent, and a full response cannot be expected.

Fig. 4
Potentials of mean force for the end-to-end distance of EHPT (A) and HPT (B) in implicit water (solid curves) and hydrocarbon (dot-dashed curves) solvents. The insets display detail of the potential well around 34–37 Å. Both polymers show ...
Table 5
Dihedral angle distribution peaks and other parameters from the PT simulations

The distribution of the backbone dihedral angles between thiophene rings should also be affected by solute/solvent compatibility. The bithiophene rotor potential of Section 3.2 describes the “natural” tendencies for the mutual orientation of thiophene rings in a polymer backbone. If different solvents change the molecule's propensity to adopt certain conformations, this can be seen in the dihedral angle distribution as a deviation from the “ideal” distribution predicted by the rotor potential. The base potential of Fig. 3 shows effectively two minima, with the lower-energy minimum corresponding to an anti-like configuration and the higher-energy minimum corresponding to a syn-like configuration; one would therefore expect a dihedral angle distribution from a PT simulation to show two peaks, with the peak at anti-like angles higher than that at syn-like angles. It should be kept in mind, however, that the simple presence of side chains can shift and flatten the torsional potential (as also discussed in Section 3.2), and that this phenomenon will interact with solvent effects to produce the final distribution.

The backbone dihedral angle distributions observed in these implicit solvent replica exchange simulations are shown in Fig. 5, and the angles corresponding to the peaks are in Table 5. All four simulations show the expected two peaks corresponding to syn and anti orientations, but there are also visible shifts in the relative peak magnitudes across polymer and solvent conditions. EHPT is consistently biased towards an anti configuration, whereas HPT is biased towards the syn configuration (although the two peaks are roughly equal height for HPT in hydrocarbon). Both polymers show an increase in the probability of anti-like angles when moving from water to alkane. Since the anti configuration is favored in the potential for rotation in vacuum, one can consider simulations that show more of an anti bias to represent polymers that are in a more “relaxed” or “free” environment. Moreover, the helical phase of PAT aggregation in poor solvent requires the chain to adopt mostly syn angles [65,64,14]. Following these lines of reasoning, the replica exchange simulation results are consistent with both EHPT and HPT being more soluble in alkanes than in water, and with EHPT being overall more soluble than HPT, all well-known properties of these PATs.

Fig. 5
Distribution of dihedral angles between the individual rings of the EHPT (A) and HPT (B) backbones during 5 ns of replica exchange molecular dynamics in implicit water (solid curves) and implicit hydrocarbon (dot-dashed curves). All 5000 frames of each ...

Both peaks in all four backbone dihedral distributions are shifted by 50–60° from the fully planar syn and anti configurations. This is consistent with the effect of alkyl side chains described above, and the same phenomenon has also been observed in other investigators' PT simulations [13,16]. Deplanarized configurations have even been observed for hydrophilic PTs in water, demonstrating that this is not a solvent-driven effect per se [66]. Ab initio studies of PT polarizability suggest that even with this substantial deplanarization, there should still be modest π-orbital overlap and electron conductivity [58]. That prediction can be evaluated by examining the persistence length of the polymers. The persistence length is defined as half the Kuhn length, which in turn can be considered as the average length over which the polymer does not bend. Since bends interrupt electronic conjugation and decrease conductivity, longer persistence lengths imply greater conductivity, and comparisons of the persistence length between two simulations can provide a rough estimate of the relative conductivity of the polymers in those simulations. For the short molecules simulated here, which contain only a few Kuhn lengths, the persistence length λ is best estimated by a modification of the method given by Heffner et al. [67]:


where R02 is the mean squared end-to-end distance, n is the degree of polymerization, and L0 is the length of a monomeric unit (taken to be 3.9 Å in this case). For light-scattering measurements, the formula is usually couched in terms of molecular weight and monomer mass, since the precise degree of polymerization is not directly known. In the case of simulations, n is known, and the weight terms can be canceled out.

The persistence lengths for the four implicit solvent simulations were calculated using 5 ns of trajectory snapshots, taking at each timestep the configuration of the replica that was visiting the 300 K commanded temperature. These results are shown in Table 5. As would be expected from the end-to-end distance PMFs, both polymers show longer persistence lengths in alkanes than they do in water, and EHPT has a higher persistence length in hydrocarbon than does HPT. All of the reported persistence lengths are significantly different from each other (p < 0:0001 by Student's t-test). However, it is worth noting that these low p-values are due to the ability of MD simulations to generate an arbitrarily large number of datapoints, and that there is a great deal of overlap between the persistence length distributions of the four simulations. As we stated above, the difference between the end-to-end distances from which we derive the persistence lengths is not large and likely cannot be large in this particular polymer/solvent model; as such, the most that can be said is that a trend is observed and may be significant.

The more important point is that all four observed lengths hover around roughly 3.5 repeat units, giving a Kuhn length of approximately seven thiophene rings. This suggests a polymer that is behaving largely as a rigid “twisted rod” and that, even in the worst case of HPT in water (mean persistence length of 3.34 repeat units) is not collapsed into a ring or partial-helix conformation. It is also noteworthy that the persistence length varies substantially over the course of each simulation, with the tails of the distributions for all four polymers (as seen in Table 5) encompassing an extremely coiled form (persistence length roughly 1.5 units) and a fully extended form (persistence length roughly 5 units). At the same time, empirical measurements of the persistence length of regioregular PTs in good solvent (chloroform) using UV-vis absorption spectroscopy tend to observe lengths of 10 repeat units or greater, substantially above what is seen here [68,69].

One property not clearly reflected in these data is the well-known “solvatochromic” effect, where some regioregular PATs change colors (due to changes in backbone angle and resulting changes in electron delocalization) when solvent quality is altered from good to poor. This property is known to hold for HPT, and has been theorized to occur through the formation of the previously mentioned helical intramolecular aggregates [64,14]. Solvatochromism has not been reported in the literature for EHPT, and we have never observed it in chloroform/ethanol dilution tests of EHPT samples.3 If the helical aggregation mechanism does indeed drive HPT solvatochromism, one could argue that it is seen in the results of Fig. 5; a PT helix should show predominantly syn angles [14], and HPT does show more syn angles than EHPT, particularly when in water (which should be a non-solvent). However, solvatochromism should also manifest as changes in the electronic conjugation length (as reflected in the persistence length). Since HPT solvatochromism takes the form of a shift to purple colors in poor solvents, one would expect greater conjugation (longer wavelength absorption) in those solvents. This effect is not seen for EHPT or HPT in these replica exchange simulations—as noted above, both polymers show longer persistence lengths in hydrocarbon than in water, opposite to the expected trend. The evidence for solvatochromism in these simulations is therefore ambivalent at best.

4. Discussion

We have presented a new set of force field parameters specifically optimized for simulation of PATs using the most popular molecular dynamics engines. Replica exchange simulations of two model PATs for which we have empirical data produced theoretical results that are mostly consistent with experiment, despite the many approximations inherent in the implicit solvent model and the single-decamer system. Although this paper presents studies only of EHPT and HPT, the same parameters could be used to study any PT as long as the partial charges were properly assigned.

The AMBER parameters derived from ab initio computations show good geometrical agreement with both empirical and previous ab initio results, indicating that they can be relied on to predict the configuration of other polymers. The HF/6-31G* rotor potential and its AMBER sine/cosine fit also match the limited empirical data available; since this term is the most important part of the PT force field, its match to the empirical data bodes well for the ability of the force field to accurately predict conformations. The agreement with the calculated potentials of other groups is not perfect, but the differences in results can largely be explained by differences in methodology. It should also be noted that the discrepancies in energies are extremely small, about 0.2 kcal/mol in the worst case. This small barrier is unlikely to have a detectable effect on the behavior of a simulated system, as it is of the same order of magnitude as the thermally induced potential energy fluctuations of a system at room temperature.

Finally, we note that the ab initio rotor potential is not symmetric, despite the molecule itself being symmetric. One potential explanation is the way the curve was derived; we chose a rigid rotation (without energy-minimizing relaxation) to increase the speed of calculation. In such a rotation, it is possible for the effective +θ and −θ configurations to have slightly different atomic coordinates, thus allowing steric effects to break the expected symmetry. It is also possible that the effect captures some not-well-understood property of the actual bithiophene molecule; as we noted above, other groups have reported an asymmetry in the rotor potential [45] and that the use of a relaxed rotation does not appear to change the energy profile [43,44]. Regardless of the source, this asymmetry should not substantially impact the quality or outcomes of simulations based on the derived force field. Whether from simple steric effects or some other source, the asymmetry captures some actual behavior of the electron orbitals of the bithiophene system. That behavior is appropriately modeled in AMBER, since the force field tuning curve matches the ab initio potential, and thus the overall system should behave in a manner approximating an actual PT chain (modulo solvent effects, which are further considered below).

In considering the results of the replica exchange simulations, it is important to keep in mind that they represent a simplified first test of the new force field. Within the scope of that limitation, it is encouraging that the simulations show multiple trends that match previous experimental findings. Both EHPT and HPT are seen to adopt more extended conformations in hydrocarbon and more collapsed (fully or partly aggregated) conformations in water. EHPT is able to achieve a more extended state than HPT, although their potentials of mean force overlap substantially. Both polymers appear to favor a “twisted rod” conformation, with EHPT showing a predominance of anti-like angles and HPT a predominance of syn-like angles. This may imply that HPT is attempting to form a helical aggregate, particularly in water. PT helices are predicted to have at least 12 units per turn, so the full helix is not observable in these decamer simulations [14]. The relatively short length of the simulated molecules (compared to the high degrees of polymerization obtainable in real syntheses) may also explain why these data show only a slight collapse on moving from hydrocarbon to water solvation. These short rod-like molecules do not have much capacity for collapse, whereas longer molecules would have more flexibility to fold back upon themselves and stabilize bent conformations.

Not all of the results from the replica exchange simulations agree with expectations. The persistence lengths observed in all four simulations are shorter than those empirically measured for PTs in good solvents (which should be the shortest lengths possible). The persistence lengths also show an inverse solvatochromic effect, with both polymers showing more conjugation in alkanes than in water. These deviations from expected trends are most likely due to the many physical variables not captured in the system reported here. While variations in the solvent dielectric constant do have a small effect on PT conformations (as evidenced by the present results), a more accurate implicit solvent model would also change the Born radius and screening parameters when moving from aqueous to organic solvation. A real solvent would also alter the molecular configuration through electrostatic and VDW interactions between solvent and solute. Hydrocarbon solvents in particular will exhibit a high viscosity, since it is energetically unfavorable for the polymer to displace many long alkyl chains. The presence of any or all of these neglected phenomena would likely widen the gap between the observed properties of HPT and EHPT. It would also help to have a truly “good” solvent in which to test both model PTs; we compared implicit alkane to implicit water, but in point of fact, even EHPT is not nearly as soluble in hydrocarbons as it is in chloroform or toluene. Unfortunately, there is no good way to capture these solvents' properties within the generalized Born model used here. In work reported elsewhere, we have used these same parameters to simulate EHPT and HPT interacting with explicit water and lipid systems, and have seen substantial differences in their behavior [31]. Finally, there is the fact that the empirical data being used as a reference involve long-chain PTs interacting with each other in solution, whereas the data presented here represent short-chain PTs that are unable to experience intermolecular interactions with neighboring polymers. Multi-chain aggregation may be particularly important for the solvatochromic effect that was not clearly seen in these results [64]. Therefore, without the possibility of aggregation, it is unsurprising that these results are not able to fully reproduce the range of known PT properties.

In summary, we have demonstrated a new force field and framework for theoretical simulation of PTs that produces results consistent with known experimental data. This framework includes a straightforward and computationally reasonable method for parameterizing other related polymers and deriving partial charges. Disagreement between our simulations and empirical expectations is largely attributable to the use of overly simplified models, and was corrected when more realistic (but more computationally intensive) ensembles were employed. This new force field can be used to predict the behaviors and functionality of PTs in proximity to biomolecular systems, an area of increasing scientific and engineering relevance. We hope that these parameters and their accompanying results will be of value to the biosensing community, and we eagerly await the new developments in sensor technology that should result.


ASW was funded during this research by a Ruth L. Kirchstein National Research Service Award (grant number F30 NS051866-01) from the National Institutes of Health. The simulations reported here were carried out using the Cray XT3 machine of the Pittsburgh Supercomputing Center, and we thank PSC for its generous grants of computer time to MK (NSF PACI CHE030007P) and YM. MK's research is also supported by NIH (GM067962-01).


3ASW, unpublished observations.


[1] McCullough RD, Ewbank PC. Regioregular head-to-tail coupled poly(3-alkylthiophene) and its derivatives. In: Skotheim TA, Elsenbaumer RL, Reynolds JR, editors. Handbook of Conducting Polymers. Marcel Dekker; New York, NY: 1998. pp. 225–258. Chapter 9.
[2] Jeffries-El M, Sauvé G, McCullough RD. In-situ end group functionalization of regioregular poly-3-alkylthiophenes using the Grignard metathesis polymerization reaction. Adv. Mater. 2004;16:1017–1019.
[3] Jeffries-El M, Sauvé G, McCullough RD. Facile synthesis of end-functionalized regioregular poly(3-alkylthiophene)s via modified Grignard metathesis reaction. Macromolecules. 2005;38:10346–10352.
[4] Mabeck J, Malliaras G. Chemical and biological sensors based on organic thin-film transistors. Anal. Bioanal. Chem. 2006;384:343–353. [PubMed]
[5] Barbarella G, Melucci M, Sotgiu G. The versatile thiophene: an overview of recent research on thiophene-based materials. Adv. Mater. 2005;17:1581–1593.
[6] Widge AS, Jeffries-El M, Cui X, Lagenaur CF, Matsuoka Y. Self-assembled monolayers of polythiophene conductive polymers improve biocompatibility and electrical impedance of neural electrodes. Biosens. Bioelectron. 2007;22:1723–1732. [PubMed]
[7] Simandiras E, Handy N, Amos R. Correlated ab initio harmonic frequencies and infrared intensities for furan, pyrrole, and thiophene. J. Phys. Chem. 1988;92:1739–1742.
[8] van Bolhuis F, Wynberg H, Havinga E, Meijer E, Staring EG. The X-ray structure and MNDO calculations of α-terthienyl: a model for polythiophenes. Synth. Met. 1989;30:381–389.
[9] Shibaev P, Schaumburg K. Conformational transitions in chiral polythiophenes. Synth. Met. 2001;124:291–294.
[10] Vouyovitch L, Brown D, Neyertz S, Gallot B. Prediction of the crystalline structure of a novel polythiophene using molecular dynamics simulations. Soft Mater. 2003;1:93–114.
[11] Shibaev P, Arzhakov M, Schaumburg K, Vinokur R, Bjornholm T, Greve D, Komolov A, Norgaard K. The conformation of poly (3-dodecyl thiophene) within methacrylic polymer matrix. J. Polym. Sci., Part B: Polym. Phys. 1999;37:2909–2917.
[12] Shibaev P, Schaumburg K, Bjornholm T, Norgaard K. Conformation of polythiophene derivatives in solution. Synth. Met. 1998;97:97–104.
[13] Marcon V, Raos G, Allegra G. Tetrathiophene on graphite: molecular dynamics simulations. Macromol. Theory Simul. 2004;13:497–505.
[14] Kiriy N, Jähne E, Adler HJ, Schneider M, Kiriy A, Minko S, Jehnichen D, Simon P, Fokin AA, Stamm M. One-dimensional aggregation of regioregular polyalkylthiophenes. Nano Lett. 2003;3:707–712.
[15] Cirák J, Tomčik P, Barančok D, Bolognesi A, Ragazzi M. Dipole moment of a modified poly(3-alkylthiophene) at the air/water interface. Thin Solid Films. 2002;402:190–194.
[16] Diaz-Quijada GA, Weinberg N, Holdcroft S, Pinto BM. Conformational analysis of oligothiophenes and oligo(thienyl)furans by use of a combined molecular dynamics/NMR spectroscopic protocol. J. Phys. Chem. A. 2002;106:1277–1285.
[17] Gesquière A, Abdel-Mottaleb M, De Feyter S, De Schryver F, Schoonbeek F, van Esch J, Kellogg R, Feringa B, Calderone A, Lazzaroni R, Breédas J. Molecular organization of bis-urea substituted thiophene derivatives at the liquid/solid interface studied by scanning tunneling microscopy. Langmuir. 2000;16:10385–10391.
[18] Corish J, Morton-Blake D, Veluri K, Bénière F. Atomistic simulations of the structures of the pristine and doped lattices of polypyrrole and polythiophene. J. Mol. Struct. 1993;283:121–134.
[19] Corish J, Morton-Blake D. An atomistic investigation of helical polythiophene. Mol. Simul. 1995;14:381.
[20] Corish J, Morton-Blake D, Bénière F, Lantoine M. Interaction of side-chains in poly (3-alkylthiophene) lattices. J. Chem. Soc., Faraday Trans. 1996;92:671–677.
[21] O'Dwyer S, Xie H, Corish J, Morton-Blake D. An atomistic simulation of the effect of pressure on conductive polymers. J. Phys.: Condens. Matter. 2001;13:2395–2410.
[22] Corish J, Feeley D, Morton-Blake D, Beniere F, Marchetti M. Atomistic investigation of thermochromism in a poly (3-alkylthiophene) J. Phys. Chem. B. 1997;101:10075–10085.
[23] Veluri K, Corish J, Morton-Blake D, Beniere F. A lattice simulation of the migration of the BF4 ion in polythiophene and polypyrrole lattice. J. Mol. Struct. 1996;365:13–19.
[24] Corish J, Morton-Blake D. The migration of anions through poly(3-octylthiophene) matrices. Phys. Status Solidi C. 2005;2:681–684.
[25] O'Farrell R, O'Dwyer S, Morton-Blake D. The transport of an ion through a channel formed by a helical electroactive polymer. Mol. Simul. 2004;30:649–659.
[26] Morton-Blake D, Yurtsever E. The entry of molecular species into the lattice of an electroactive polymer during its dissolution. J. Mol. Liq. 2006;124:106–113.
[27] Case DA, et al. AMBER 9. University of California; San Francisco: 2006. 2006.
[28] Lindahl E, Hess B, van der Spoel D. GROMACS 3.0: a package for molecular simulation and trajectory analysis. J. Mol. Modell. 2001;7:306–317.
[29] MacKerell AD, Brooks B, Brooks CL, Nilsson L, Roux B, Won Y, Karplus M., III . The Encyclopedia of Computational Chemistry. vol. 1. John Wiley & Sons; Chichester: 1998. CHARMM: the energy function and its parameterization with an overview of the program; pp. 271–277.
[30] Widge A, Matsuoka Y. Conductive polymer “molecular wires” increase electrical conductance across artificial cell membranes. Proceedings of the 26th Annual International Conference of the IEEE Engineering in Medicine and Biology Society.2004. pp. 4330–4333.
[31] Widge AS, Matsuoka Y, Kurnikova M. In Silico insertion of poly(alkylthiophene) conductive polymers into phospholipid bilayers. Langmuir. 2007;23:10672–10681. [PMC free article] [PubMed]
[32] Widge AS. Ph.D. Thesis. Carnegie Mellon University; Pittsburgh, PA: 2007. Self-assembled monolayers of polythiophene “molecular wires”: a new electrode technology for neuro-robotic interfaces.
[33] Frisch MJ, et al. Gaussian 03, Revision B. 05. Gaussian, Inc.; Wallingford, CT: 2004.
[34] Ciofalo M, La Manna G. Ab initio conformational study of 2, 2′ : 5′, 2′′-terthiophene. Chem. Phys. Lett. 1996;263:73–78.
[35] Kofranek M, Ková[r with breve] T, Lischka H, Karpfen A. Ab initio studies on heterocyclic conductive polymers: structure and vibrational spectra of thiophene, oligothiophenes, and polythiophene. J. Mol. Struct. 1992;259:181–198.
[36] Ramirez F, Hernandez V, Navarrete J. Transferable semiempirical quadratic force fields: the case of polythiophene and shorter oligomers. J. Comput. Chem. 1994;15:405–423.
[37] Wang J, Wang W, Kollman P. Antechamber, an accessory software package for molecular mechanical calculations. J. Comput. Chem. 2005;25:1157–1174. [PubMed]
[38] Yang L, Tan C, Hsieh M, Wang J, Duan Y, Cieplak P, Caldwell J, Kollman PA, Luo R. New-generation Amber united-atom force field. J. Phys. Chem. B. 2006;110:13166–13176. [PubMed]
[39] Samdal S, Samuelsen E, Volden H. Molecular conformation of 2, 2′-bithiophene determined by gas phase electron diffraction and ab initio calculations. Synth. Met. 1993;59:259–265.
[40] Ortí E, Viruela PM, Sánches-Marín J, Tomás F. Ab initio determination of the geometric structure and internal rotation potential of 2, 2′-bithiophene. J. Phys. Chem. 1995;99:4955–4963.
[41] Diaz-Quijada GA, Weinberg N, Holdcroft S, Pinto BM. Investigation of barriers to conformational interchange in oligothiophenes and oligo(thienyl)furans. J. Phys. Chem. A. 2002;106:1266–1276.
[42] Westenhoff S, Beenken WJ, Yartsev A, Greenham NC. Conformational disorder of conjugated polymers. J. Chem. Phys. 2006;125:154903. [PubMed]
[43] Hernandez V, Navarrete J. Ab initio study of torsional potentials in 2, 2′-bithiophene and 3, 4′- and 3, 3′-dimethyl-2, 2′-bithiophene as models of the backbone flexibility in polythiophene and poly (3-methylthiophene) J. Chem. Phys. 1994;101:1369.
[44] Quattrocchi C, Lazzaroni R, Brédas J. Theoretical investigation of the conformational behavior of 2, 2′-bithiophene. Chem. Phys. Lett. 1993;208:120–124.
[45] Di Césare N, Belletête M, Leclerc M, Durocher G. A conformational study of ethyl-substituted bithiophenes. Semi-empirical versus ab initio methods. Synth. Met. 1998;94:291–298.
[46] Karpfen A, Choi CH, Kertesz M. Single-bond torsional potentials in conjugated systems: a comparison of ab initio and density functional results. J. Phys. Chem. A. 1997;101:7426–7433.
[47] Raos G, Famulari A, Marcon V. Computational reinvestigation of the bithiophene torsion potential. Chem. Phys. Lett. 2003;379:364–372.
[48] Ryckaert J, Bellemans A. Molecular dynamics of liquid alkanes. Faraday Discuss. 1978;66:95–106.
[49] Onufriev A, Bashford D, Case D. Exploring protein native states and large-scale conformational changes with a modified generalized Born model. Proteins. 2004;55:383–394. [PubMed]
[50] Xiang T, Anderson BD. Conformational structure, dynamics, and solvation energies of small alanine peptides in water and carbon tetrachloride. J. Pharm. Sci. 2006;95:1269–1287. [PubMed]
[51] Radmer RJ, Klein TE. Triple helical structure and stabilization of collagen-like molecules with 4(R)-hydroxyproline in the Xaa position. Biophys. J. 2006;90:578–588. [PubMed]
[52] Earl DJ, Deem MW. Parallel tempering: theory, applications, and new perspectives. Phys. Chem. Chem. Phys. 2005;7:3910–3916. [PubMed]
[53] Sugita Y, Okamoto Y. Replica-exchange molecular dynamics method for protein folding. Chem. Phys. Lett. 1999;314:141–151.
[54] Rathore N, Chopra M, DePablo JJ. Optimal allocation of replicas in parallel tempering simulations. J. Chem. Phys. 2005;122:02411. [PubMed]
[55] Kone A, Kofke DA. Selection of temperature intervals for parallel-tempering simulations. J. Chem. Phys. 2005;122:206101. [PubMed]
[56] Chodera JJ, Swope WC, Pitera JW, Seok C, Dill KW. Use of the weighted histogram analysis method for the analysis of simulated and parallel tempering simulations. J. Chem. Theor. Comput. 2006;3:26–41. [PubMed]
[57] Almenningen A, Bastiansen O, Svendsas P. Electron diffraction studies of 2,2′-diethienyl vapour. Acta. Chem. Scand. 1958;12:1671–1674.
[58] Mosley DH, Fripiat JG, Champagne B, André JM. Ab initio investigation of the static polarizability of planar and twisted infinite polythiophene chains. Int. J. Quantum Chem. 1994;28:451–467.
[59] Di Césare N, Belletête M, Durocher G, Leclerc M. Towards a theoretical design of thermochromic polythiophenes. Chem. Phys. Lett. 1997;275:533–539.
[60] Khetrapal C, Kunwar A. The conformation of 2, 2′-bithiophene as determined from PMR studies in a nematic solvent. Mol. Phys. 1974;28:4416.
[61] Bucci P, Longeri M, Veracini C, Lunazzi L. Nematic phase nuclear magnetic resonance investigations of rotational isomerism. III. Conformational preferences and interconversion barrier of 2, 2′-bithienyl. J. Am. Chem. Soc. 1974;96:13059.
[62] Chadwick J, Kohler B. Optical spectra of isolated s-cis- and s-trans-bithiophene: torsional potential in the ground and excited states. J. Phys. Chem. 1994;98:3631–3637.
[63] Yue S, Berry G, McCullough R. Intermolecular association and supramolecular organization in dilute solution. 1. Regioregular poly(3-dodecylthiophene) Macromolecules. 1996;29:933–939.
[64] Kiriy N, Jähne E, Kiriy A, Adler HJ. Conformational transitions and aggregations of regioregular polyalkylthiophenes. Macromol. Symp. 2004;210:359–367.
[65] Yang R, Evans D, Christensen L, Hendrickson W. Scanning tunneling microscopy evidence of semicrystalline and helical conducting polymer structures. J. Phys. Chem. 1990;94:6117–6122.
[66] Leith D, Morton-Blake D. A molecular dynamics simulation of the structure and properties of a self assembled monolayer formed from an amphiphilic polymer on a water surface. Mol. Simul. 2005;31:987–997.
[67] Heffner GW, Pearson DS, Gettinger C. Characterization of poly(3-octylthiophene). I. Molecular characterization in dilute solution. Polym. Eng. Sci. 1995;35:860–867.
[68] Sandstedt C, Rieke R, Eckhardt C. Solid-state and solvatochromic spectra of a highly regular polythiophene. Chem. Mater. 1995;7:1057–1059.
[69] Bidan G, De Nicola A, Enée V, Guillerez S. Synthesis and UV–visible properties of soluble regioregular oligo(3-octylthiophenes), monomer to hexamer. Chem. Mater. 1998;10:1052–1058.
[70] Bak B, Christensen D, Hansen-Nygaard L, Rastrup-Andersen J. The structure of thiophene. J. Mol. Spectrosc. 1961;7:58–63.
[71] Barbarella G, Zambianchi M, Bongini A, Antolini L. Conformational chirality of oligothiophenes in the solid state. x-ray structure of 3; 4′, 4′′-trimethyl-2, 2′ : 5′, 2′′-terthiophene. Adv. Mater. 1994;6:561–564.