The influence of sorbitol doping on aggregation and electronic properties of PEDOT:PSS: a theoretical study

Many organic electronics applications such as organic solar cells or thermoelectric generators rely on PEDOT:PSS as a conductive polymer that is printable and transparent. It was found that doping PEDOT:PSS with sorbitol enhances the conductivity through morphological changes. However, the microscopic mechanism is not well understood. In this work, we combine computational tools with machine learning to investigate changes in morphological and electronic properties of PEDOT:PSS when doped with sorbitol. We find that sorbitol improves the alignment of PEDOT oligomers, leading to a reduction of energy disorder and an increase in electronic couplings between PEDOT chains. The high accuracy (r 2 > 0.9) and speed up of energy level predictions of neural networks compared to density functional theory enables us to analyze HOMO energies of PEDOT oligomers as a function of time. We find a surprisingly low degree of static energy disorder compared to other organic semiconductors. This finding might help to better understand the microscopic origin of the high charge carrier mobility of PEDOT:PSS in general and potentially help to design new conductive polymers.

Despite extensive experimental studies of PEDOT:PSS, the material is still not well understood from a microscopic perspective on the molecular level [15][16][17][18] While mesoscale domain formation is challenging to investigate on an atomistic level, our study aims to investigate the microscale arrangement and electronic properties of PEDOT oligomers and PSS chains and how the microscopic morphology changes when mixed with sorbitol, a widely used dopant for PEDOT:PSS that is known to enhance its conductivity.
In this work, we use a multiscale modelling approach [19] combining classical and quantum mechanical simulation tools with machine learning to microscopically investigate nanoscale morphology and electronic properties of PEDOT:PSS films. We start with the parameterization of a classical force field using density functional theory (DFT), followed by molecular dynamics (MD) simulations to generate simulation boxes of PEDOT:PSS with varying concentrations of D-sorbitol [20]. The boxes are then analyzed using DFT, yielding the energy disorder of the PEDOT oligomers caused by their individual conformations in the amorphous film. However, charge transport mostly depends on the static disorder, averaging over fluctuations of each single molecule. The analysis of a dynamical time series (MD trajectory) of each molecule in the film is computationally too costly, so we employed the data obtained by analyzing the last MD timestep to train a neural network model that can predict the HOMO energies of arbitrary conformations of PEDOT oligomers at a fraction of the cost of a DFT calculation.
We find that adding D-sorbitol to PEDOT:PSS films improves the intermolecular arrangement of PEDOT oligomers, increasing the electronic coupling and decreasing the energy disorder in the amorphous film. We furthermore reveal that-despite its flexibility-PEDOT oligomers show a surprisingly low amount of static disorder which correlates well with its high charge carrier mobility. Combined with the high level of doping and, potentially, doping induced disorder compensation effects [21], these findings can help explain the exceptionally high conductivity of PEDOT:PSS.
In section 2, we introduce the methodology developed in this work (see also figure 1), covering classical force field based MD simulations, quantum mechanical DFT calculations and neural networks architecture and training. Section 3 shows and discusses our findings, notably the influence of D-sorbitol on the microstructure and on electronic properties of PEDOT:PSS films. Finally, we draw our conclusions in section 4.

Molecular dynamics
The simulation reported herein were carried out with the LAMMPS software package [22]. The GAFF force-field has been chosen to describe the simulated systems [23]. Atomic charge parameters have been obtained from the results of the quantum chemical calculations (see subsection 'Quantum-mechanical calculations').
The MD simulations were carried out on PEDOT-PSS and PEDOT-PSS-sorbitol (10%, 20%, 40%) model systems in a cubic periodic box (see Table 1). PEDOT oligomers with eight repeat units and +2 charge were considered, while the PSS atactic chains consisted of 20 repeat units, with four deprotonated units randomly distributed in the sequence of each chain.
To ensure an adequate equilibration of the systems, the following approach was adopted. First, random starting configurations were generated with the PACKMOL program [24], and each one were allowed to relax with an energy minimization followed by a molecular dynamics simulation of 1 ns within the NVT ensemble at 300 K. After that, an NPT simulation of 1 ns duration at a temperature of 300 K and a pressure of 1 atm was performed to relax the dimensions of the box. Then, the temperature of the system was progressively increased to 700 K and decreased again to 300 K before the production run of 5 ns in the NPT ensemble at 300 K and 1 atm. Using that procedure, we generated 8 independent morphologies for each of the sorbitol concentrations using different random starting configurations. Boxes with various contents of sorbitol are shown in figures 1 and S13 (available online at stacks.iop.org/MLST/2/01LT01/mmedia).

Quantum-mechanical calculations for force field parameterization
Quantum chemistry calculations were performed to parameterize the atomic forces for the molecular dynamic simulation. Calculations were performed at the HF/6-31 G(d,p)//B97-D [25]/6-31 G(d) level of theory as customary with the GAFF force field.
For the PEDOT-PSS complex, 200 starting geometries were generated using the PACKMOL software package. In a first step, they were optimized using a crude parameterized force field (FF) approach and were ranked according to their respective FF energy. The first twenty most stable structures were then relaxed at the B97-D/6-31 G(d,p) level of theory, using the Gaussian software package. For the sake of computational cost effectiveness, it is important to mention that the PSSpentamers were substituted by a single p-toluene-sulphonate unit. For D-sorbitol, three different starting structures were optimized. Note that relative stabilities of the different conformations were determined at the B97X-D/def2-TZVPD [26]//B97-D/6-31 G(d) level to reduce any artifacts arising from basis set superposition error. According to the Boltzmann distribution, the most stable conformation was found to have a population of 99%. As such, only the most stable conformation was used for analysis of the charges.

Quantum-mechanical calculations for morphology analysis
We performed single-point density functional theory (DFT) calculations (B3-LYP [27]/def2-SV(P) [28] level of theory as implemented in Turbomole [29]) of each PEDOT conformer extracted from the last snapshot of all MD trajectories. We chose a charge of 2+ for each PEDOT oligomer which allows us to interpret the energy of the highest occupied molecular orbital as the approximate energy of an additional third hole on the PEDOT chains, which is the typical charge of PEDOT oligomers when being ionized by PSS chains at the PEDOT:PSS ratio assumed in this work.
The electronic couplings were calculated using DFT (B3-LYP/def2-SV(P)) using a frontier (HOMO) orbital based diagonalization method [30,31]. Fock and overlap matrices were extracted from PEDOT dimer calculations and HOMO orbitals were extracted from calculations of isolated PEDOT molecules.

Machine learning methods
We used the neural network module as implemented in scikit-learn [32] and performed a Bayesian optimization based hyperparameter optimization [33]. The hyperparameters included the number of hidden layers, the size of the first hidden layer, the ratio between subsequent hidden layer sizes, the non-linear activation function, the L2 regularization parameter, the learning rate as well as the patience of the learning rate scheduler. The latter reduced the learning rate with a factor of 0.5. We used the adam optimization algorithm to train the model. The best performing model had one hidden layer with 125 neurons, an L2 regularization parameter of 10.53, a tanh activation function and a patience of the learning rate scheduler of 17 and an initial learning rate of 2.6 × 10 −4 . The conformation of the PEDOT oligomers consisted of all inverse bond lengths (short range distances), all inverse distances between sulfur and oxygen atoms (long range distances), the eigenvalues of the Coulomb matrix [34], the eigenvalues of modified Coulomb matrices using electronegativities (H: 2.2, O: 3.44, C: 2.55, S: 2.58) and ones instead of the atomic numbers Z.

Results and discussion
In this section, we discuss the impact of D-sorbitol on the packing motif as well as the electronic properties of PEDOT:PSS. We start with an analysis of PEDOT conformations as well as PEDOT:PEDOT and PEDOT:PSS arrangement, followed by a DFT and ML based analysis of the energy levels of PEDOT oligomers and the static and dynamic energy disorder in PEDOT:PSS.

Conformational properties of the PEDOT and PSS chains
The conformational properties of the PEDOT chains were evaluated in terms of the end-to-end distance and the radius of gyration. Boxes without sorbitol as well as with 10% and 40% sorbitol concentration are shown in figures 1 and S1. Table 2 reports on the average values and standard deviation for the 64 PEDOT chains in each of the 8 boxes considered for each sorbitol concentration (table S1 shows the same information for PSS chains). Note that the PEDOT polymer chains have a very rigid structure, thus leading to end-to-end distances close to the value of the completely extended chain (i.e. approximately 30.07 Å). Nevertheless, it  can be seen that some of the PEDOT chains adopt a relatively bent geometry. Figure 2 displays the distribution of values of the end-to-end distance of PEDOT chains for the different sorbitol concentrations.
To further investigate the conformation of PEDOT oligomers, we quantified the extent of deviation from a planar conformation (using the residuals of a two component principal component analysis, see figure S2) and the RMSD between all PEDOT oligomers and the most elongated structure found (see figure S3). Both of these analysis confirmed an increase in order of PEDOT oligomers with increasing sorbitol concentration (more planar conformers and more distinct peaks for zero, one or more monomer flips).
The structural properties of the different systems, in particular the relative packing between the different components, were analysed in terms of atomic pair-pair radial distribution functions (RDF). Figure 3 shows the intermolecular radial distribution function for PEDOT oligomers in boxes with varying sorbitol concentrations. The plot shows an increase in the PEDOT-PEDOT contacts at short distances with increasing the sorbitol concentration. The four systems show a RDF maximum at approximately 4.2 Å (the experimental interlayer distance is about 3.5 Å), which is lower and broader for the system without sorbitol. As the sorbitol concentration increases, the peak becomes higher and narrower, suggesting that more PEDOT chains come closer to each other, but also at a more definite distance. It should be noted that the radial  distribution function displays a second peak at about 7.4 Å, which again is more pronounced at higher sorbitol concentration. We attribute this peak to the stacking of three PEDOT chains.
The increase in PEDOT:PEDOT contacts with the sorbitol concentration is linked to a decrease in the PEDOT:PSS and PSS:PSS contacts, as shown in figures 4 and S3, respectively. These results suggest that sorbitol molecules tend to aggregate close to the PSS chains, and, as a consequence, PEDOT oligomers are more prone to aggregate.

Energy disorder of PEDOT chains
To investigate the effect of sorbitol on the electronic structure and thus charge transport properties of PEDOT:PSS:sorbitol films, we computed the conformational energy disorder of the last snapshots of the MD trajectories presented above. The conformational disorder is the width of the distribution of energy levels arising due to conformational differences between molecules in amorphous thin films. It can be separated into static and dynamic disorder, which is discussed in the next section. As shown in figure 2, the fraction of elongated PEDOT chains is increasing with increasing sorbitol concentration while the fraction of strongly bent PEDOT oligomers with short end-to-end distances is vanishing. As we show in figure 5, the increasing homogeneity of PEDOT end-to-end distances leads to a decrease in conformational energy disorder with increasing sorbitol concentration from 77 meV at 0% sorbitol concentration to 64 meV at 40% sorbitol concentration. Compared to many other organic hole transport materials, the conformational disorder of PEDOT oligomers-despite their flexibility-is low, partially explaining the good hole transport properties of PEDOT:PSS(:sorbitol) films [35].

Machine learning based analysis of dynamic and static disorder
Energy disorder is known to strongly influence the charge carrier mobility of organic semiconductors as wider distribution of energy levels is trapping charge carriers in its tail states [35][36][37]. However, the conformational energy disorder discussed in the last section has two contributions: dynamic disorder and static disorder, and only the latter is considered to influence the charge carrier mobility [38][39][40]. The calculation of dynamic disorder is a computationally expensive task, as it requires the computation of energy levels over large time-scales to observe molecule and conformation specific site energy fluctuations. This typically necessitates a quantum mechanical analysis of MD trajectories over multiple time frames. In our case, the analysis of 32 morphologies with 64 PEDOT oligomers each over e.g. 500 time frames would sum up to 1 024 000 DFT single point calculations which is computationally barely feasible.
To reduce the computational cost of this method, we employed a machine learning based approach in which we trained a neural network on 80% of the 2048 PEDOT conformers (last snapshots of 32 boxes with 64 PEDOT molecules each) that we analysed using DFT to obtain the energy levels shown in figure 5. To represent the conformation of the PEDOT molecules, we used the inverse bond lengths, inverse distances between all sulfur and oxygen atoms, all bond angles as well as the eigenvalues of the Coulomb matrix and a modified Coulomb matrix based on the electronegativities of the elements (see discussion in SI). We compared a linear regression model to neural networks, a gradient boosting regression model, a random forest model and a Gaussian Process model. Among all methods, a neural network obtained using Bayesian optimization based hyperparameter optimization (426 sequential steps) reached the highest accuracy with a mean absolute error of MAE = 0.0167 eV (0.385 kcal mol −1 ) and a coefficient of determination of r 2 = 0.905. The performance of linear regression as a base-line model was r 2 = 0.818. The predictions of this model on the validation and test set are shown in figure 6(a). Further tests of the accuracy of the neural network model are shown in figures S11-S13. The observed speedup of a HOMO energy prediction using the neural network compared to the reference DFT computations was (5.48 ± 1.42)·10 7 . The energy level evaluation for a single conformer was 34.65 ± 7.38 µs for the neural network compared to 1897.55 ± 86.43 s for DFT, each on a single CPU. This makes the generation of the required input vectors for the neural network from 3D coordinates the computational bottleneck.
We then applied the trained model to predict the energy levels of all PEDOT molecules in the MD calculations throughout the entire trajectory (4 sorbitol concentrations, 8 boxes at each concentration, 64 PEDOT oligomers and 500 snapshots per box, in total 1 024 000 PEDOT molecules). The results for the pure PEDOT:PSS morphologies are shown in figure 6(b). We find that the total disorder agrees well with the disorder predicted for only one snapshot. Analysis of the trajectories of each of the 512 PEDOT molecules allows us to extract the dynamic disorder, which surprisingly accounts for a large fraction of the total disorder. A time-average of all trajectories eventually leads to the static disorder (see figures 6(c) and (d)). The dynamic disorder is almost independent of the sorbitol concentration (59 meV, 59 meV, 58 meV and 58 meV for 0%, 10%, 20% and 40% sorbitol concentration, respectively), while the total disorder decreases with increasing sorbitol concentration (72 meV, 67 meV, 65 meV and 64 meV). Consequently, the static disorder decreases from 42 meV at 0% sorbitol concentration to 32 meV, 27 meV and 28 meV at 10%, 20% and 40% sorbitol concentration. These disorder values are exceptionally low compared to other amorphous organic semiconductors [35,40], which is in good agreement with the high hole conductivity of PEDOT:PSS. A low static disorder increases the hopping rates which exponentially depend on the energy differences between neighboring molecules. Strong temporal fluctuations at the same time enhance the hopping attempt frequency.
In addition, we observe that the distribution of static HOMO energies strongly deviates from a Gaussian distribution. While the distribution has a long tail towards lower HOMO energies, it sharply cuts off at higher energies, which has important consequences for hole transport in PEDOT:PSS(:sorbitol) films as well as attempts to model charge transport in such systems. Holes prefer to be localized on PEDOT molecules Figure 7. Electronic couplings between PEDOT(2+) oligomers as a function of the center-of-geometry distance of the chains. The pairs were extracted from PEDOT:PSS:sorbitol morphologies with different sorbitol concentrations (0%, 10%, 20% and 40%, respectively). The distribution of electronic couplings is shifted to higher values with increasing sorbitol concentration. The horizontal lines denote the mean of the ten highest electronic couplings at each concentration.
with high HOMO energies, where, due to missing tail states (shallow traps), the effective energy disorder is even lower than that of the full distribution. As a consequence, the hole mobility of PEDOT will be higher than that of another material with the same energy disorder parameter but with a Gaussian distribution of energy levels.
We note that the total disorder values obtained using the ML model are slightly different than the values reported in figure 5, firstly because we are now analysing the entire trajectory and secondly due to a small systematic error of the ML model (see figure S13). Figure 3 shows that with increasing sorbitol concentration, the PEDOT-PEDOT pair distribution function shows a pronounced peak at low distances. To compute the effect of this short term order on charge transport, we computed the electronic couplings between pairs of PEDOT oligomers in morphologies with varying sorbitol concentration. The results shown in figure 7 indicate that improved alignment of PEDOT chains at high sorbitol concentration leads to an increase of high electronic couplings in the order of 0.1 eV and above. Electronic couplings of this value are in the same order of magnitude as typical reorganization energies of organic semiconductors which can lead to a delocalization of electronic states beyond the 'small/localized polaron' picture, in which Marcus theory is valid [41][42][43]. This delocalization can lead to a significant increase in charge carrier mobility as well as to changes of the temperature dependence, which is in particular observed in organic crystals such as rubrene [44].

Conclusion
We presented a multi-scale simulation study of PEDOT:PSS doped with D-sorbitol and analysed the effect of sorbitol on morphological and electronic properties of the PEDOT oligomers in PEDOT:PSS. Using molecular dynamics simulations we found that an increased sorbitol concentration leads to a larger end-to-end distance and thus less bending as well as an enhanced packing of PEDOT oligomers. A DFT based quantum mechanical analysis showed that this reduces the energy disorder and increases the electronic couplings between PEDOT oligomers. Energy disorder is typically split into static and dynamic contributions, where only the static disorder influences charge carrier mobility and thus the conductivity of a material. To further analyse charge transport properties of PEDOT:PSS, we used a machine learning based approach to compute the static and dynamic disorder of PEDOT oligomers. We trained a neural network to predict the HOMO energies of PEDOT chains throughout entire MD trajectories (more than one million conformers), which would be computationally intractable using DFT. We found that the static disorder of PEDOT oligomers is exceptionally small compared to other organic semiconductors, contradicting the widespread intuition that rigid molecules are desirable to decrease (static) energy disorder and thus increase the charge carrier mobility.