Morphology and ion diffusion in PEDOT:Tos. A coarse grained molecular dynamics simulation.

A Martini coarse-grained Molecular Dynamics (MD) model for the doped conducting polymer poly(3,4-ethylenedioxythiophene) (PEDOT) is developed. The morphology of PEDOT:Tos (i.e. PEDOT doped with molecular tosylate) and its crystallization in aqueous solution for different oxidation levels were calculated using the developed method and compared with corresponding all atomistic MD simulations. The diffusion coefficients of Na+ and Cl- ions in PEDOT:Tos are studied using the developed coarse-grained MD approach. It is shown that the diffusion coefficients decrease exponentially as the hydration level is reduced. It is also predicted that the diffusion coefficients decrease when the doping level of PEDOT is increased. The observed behavior is related to the evolution of water clusters and trapping of ions around the polymer matrix as the hydration level changes. The predicted behavior of the ionic diffusion coefficients can be tested experimentally, and we believe that molecular picture of ionic diffusion in PEDOT unraveled in the present study is instrumental for the design of polymeric materials and devices for better and enhanced performance.


Introduction
Conducting polymers provide an exciting possibility for cheap and easily processable printable and flexible electronic and bioelectronic devices.2][3][4][5][6][7][8][9] The PEDOT oxidation level can be tuned by a chemical reducing agent 10 or electrochemically. 112][3][4][5] The low work function (4.30 eV) and high conductivity make PEDOT an appropriate candidate for source and drain electrodes. 18The electrical conductivity, Seebeck coefficient and thermoelectric efficiency of PEDOT are dependent on its oxidation level. 10Because of the high ionic conductivity, PEDOT represents an excellent material for energy storage devices such as supercapacitors and fuel cells, 19 as well as for biolectronic devices such as sensors, electrochemical transistors 20 and ion pumps. 21Ionic conductivity enhancement with increased humidity level and the ionic Seebeck effect were recently reported for Na + ions in PEDOT:PSS. 5reviously the effect of hydration level was mainly studied experimentally 22,23 and theoretically 24,25 in proton transport through hydrated Nafion for application in fuel cells.
The ionic and electronic transport properties of PEDOT:Tos are strongly related to and determined by its complex morphology. 264]27 At the same time, theoretical work addressing the morphology of PEDOT:PSS is practically missing.][30] Recent All Atomistic (AA) Molecular Dynamics (MD) calculations for PEDOT doped with different molecular counterions (including Tos) reveal a complex morphology where small crystallites consisting of several p-p stacking PEDOT chains are embedded in a polymer matrix including PEDOT chains, Tos counterions and water molecules. 16,31Because of computational limitations, AA MD simulations do not always allow one to reach a lengthscale needed to study the morphology of a realistic material.For example, a formation of PEDOT-rich and PSS-rich phases in PEDOT:PSS that happens on the scale of B20-30 nm 32 is beyond the reach of AA MD.Theoretical investigation of the morphology of PEDOT on these scales requires coarse-grained (CG) approaches where the system is represented by reduced numbers of the degrees of freedom.In CG MD simulations the atoms are grouped in some specific beads to reduce complexity and computational cost for larger systems at longer time scales.
The Martini CG force field offers speed, accuracy, applicability, and versatility for MD simulations. 33So far, a Martini coarse-grained model for PEDOT is still missing, and therefore one of the aims of our study is to develop and validate the Martini coarse-grained model for doped PEDOT.A simplified CG model of PEDOT was used in ref. 34, where the entire EDOT monomer was treated as one bead.This model is however too crude and simplified and is not in a position to reproduce the most important morphological aspects of PEDOT on a subnanometer scale, such as its crystallization.On the contrary, our Martini CG model correctly reproduces the formation of p-p crystallites, as evidenced by the distributions of distance between chains.Our second aim is to apply the CG model to investigate and understand the ionic diffusion in PEDOT for different oxidation and water content.Such understanding is essential for improving device and material design because the diffusion of ions such as Na + and Cl À determines the performance of PEDOT-based devices with mixed electron-ionic conductivity, including supercapacitors, organic electrochemical transistors, and ion pumps.At the same time, theoretical understanding of ionic diffusion in PEDOT, in particular the dependence of the ionic diffusion coefficients on the water content, is still absent.
In the present study, we develop and validate a Martini CG MD model for PEDOT, and use it to simulate the polymer morphology and calculate the diffusion coefficients of Na + and Cl À ions in PEDOT:Tos.We find that the diffusion coefficient decreases exponentially as the hydration level is reduced.We also find that the diffusion coefficients decrease when the doping level of PEDOT is increased.The obtained behavior is analyzed based on the water cluster evolution and trapping of ions around the polymer matrix as the hydration level changes.The predicted behavior of the ionic diffusion coefficients can be tested experimentally, and we believe that the detailed molecular picture of ionic diffusion in PEDOT unraveled in the present study is instrumental for the design of polymeric materials and devices for better and enhanced performance.

The coarse-grained model for PEDOT
In this section we present the CG model for PEDOT and outline our CG MD simulations using the Martini CG force field.8][39][40][41][42] The Martini model categorizes the beads into polar (P), nonpolar (N), apolar (C), and charged (Q) ones.For each kind of bead, there are subgroups which give more freedom for choosing the appropriate bead type.For example, the N-type beads are categorized based on the H-bonding capabilities: none (N0), donor (Nd), acceptor (Na) and donor-acceptor (Nda). 35There are also S-type beads with a smaller radius and weaker van der Waals interactions for modeling of ring structures. 35Recently, the Martini model was successfully applied to polymers like poly(3-hexyl-thiophene) (P3HT), [37][38][39]43 phenyl-C61-butyric acid methyl ester (PCBM), 37 polystyrene 40,41 and PSS. 42so, a polarizable CG model for water is developed based on one neutral and two charged beads. 44The validity of the Martini polarizable water model for application to systems with charges is confirmed for the ionic conductivity of Na + and Cl À ions in pure water. 45 Bad on the CG model for P3HT [37][38][39]43 we use the SC2 and SC1 bead types for the thiophene ring in PEDOT, where the SC2 bead is used for the sulfur atom and two less polarized SC1 beads are used for the four carbons, see Fig. 1(a). Th remaining dioxy group with hydrogen atoms is considered as a single SNda bead.The distance between nearest neighbor beads is constrained using the LINCS algorithm.46 A virtual site is placed in the center of mass of the sulfur and two neighboring carbon atoms to connect the thiophene rings and have more control of the PEDOT backbone.37,39,43 The connection between neighboring EDOT units is modeled by a harmonic force between virtual sites. Th virtual site has a zero mass and it does not interact with other beads via van der Waals interactions.We consider doped PEDOT with different oxidation levels C ox , where positive charges in the polymer backbone are compensated for by negative counterions.Ab initio calculations show that for different values of C ox , the distribution of positive and negative charges is practically the same for different monomers (with some deviations for monomers at the ends of the chain).16 Therefore, in our CG model, we assume that for each oxidation level C ox , positive charges are equally distributed on virtual sites.For example, for C ox = 33.3% each virtual site gets a charge of +0.333e.Our CG model of the tosylate molecule is based on the model for PSS recently developed by Vo ¨gele et al. 42 as shown in Fig. 1(b).The PSS model includes SCY and STY beads and considers the sulfonate group as a charged Q a bead.42 The water molecules are represented by the polarizable water model. 44The non-bonded parameters are described by Martini force field V2.2 alongside the required Martini beads and parameters for PSS and polarizable water.42,44,45,47 The bonded parameters, which include equilibrium distances, angles, dihedrals and force constants between PEDOT beads, were extracted from the comparison of different distributions in the CG and AA models.For the AA MD simulation we use the general AMBER force field (GAFF), 48 which was previously used to describe the morphology of PEDOT doped with various counterions including Tos 16,31 and a related self-doped polymer ETE-S. 49,50 A omparison between the bond lengths and angles calculated using the CG and AA approaches is given in section one of ESI, † SI1, see Fig. S1.The calculated parameters for the PEDOT CG model are listed in Table S1 (ESI †).Note that even though in the present study we focus on the morphology and ion diffusion for the case of molecular counterions such as Tos, our CG model for PEDOT can be used to study PEDOT crystallization in the presence of other polyanions such as PSS.
For the Tos molecule, the CG model from ref. 42 is validated using the AA model of Tos reported in ref. 16, see Fig. S2 in Section SI2 (ESI †).
In our simulation we use a computational cubic box filled with randomly distributed PEDOT chains, tosylate and polarized water beads under 3D periodic boundary conditions.(Typically the computational box contains 400 PEDOT chains of length N = 12 monomer units, and 52 000 water molecules; the number of Tos molecules depends on the oxidation level and varies between 1600 for C ox = 33.3% and zero for C ox = 0%).The MD simulations start with water equilibration, where three equilibration steps are performed as follows.First, a canonical NVT ensemble is run with a 5 fs time step for 10 ns using the Berendsen thermostat. 51Then, the equilibration continues with another NVT ensemble for 100 ns with a 20 fs time step.Lastly, the system is equilibrated in an isothermalisobaric NPT ensemble for 100 ns with a 20 fs time step by using the Berendsen barostat at 1 bar. 51During all water equilibration steps, we apply position restraints to both the PEDOT and tosylate molecules.Having equilibrated the water, we start the production run using the NPT ensemble and the Berendsen barostat for 200 ns and a 20 fs time step.The simulation temperature was T = 300 K during the equilibrations and the production run using the Berendsen thermostat.For all simulations, the long-range Coulomb electrostatic interactions are calculated using the Particle Mesh Ewald (PME) method. 52For the drying process, we remove the water beads in four steps and study the system with the hydration levels Hy = 80, 60, 40, 20 and 10 percent of water, where the hydration level (Hy) is defined as the water weight percent wt% with respect to the total solution weight.At each drying step, randomly chosen water molecules are removed from the box and the above MD equilibration routine is repeated.Our computation procedure for water evaporation is designed to mimic a corresponding experimental procedure performed by Palumbiny et al. for printed PEDOT:PSS. 536][57] For CG simulations, the van der Waals radius for CG beads and probe radius are set to the potential energy minimum distance.For the 6-12 Lennard-Jones potential in the Martini force field, the minimum is at r = 2 1/6 s, 58 where s = 0.47 nm and s = 0.43 nm for the ordinary and ring type beads, respectively.For the case of AA calculations, the probe radius is 0.14 nm, corresponding to the size of the water molecule, see ref. 16 for details.
We chose Na + and Cl À as diffusive ions because they are used as an electrolyte in a typical electrochemical setup. 32Also, Na + and Cl À represent biologically active ions providing the electrical activity needed to support muscle contractions and neuron activation and they are therefore utilized in various PEDOT-based bioelectronic devices. 59In CG calculations Na + and Cl À are respectively described by single Q d and Q a Martini CG beads. 35,44,45The diffusion constant is calculated in a standard way by using the mean square displacement and the Einstein relation. 45 where Dr 2 (t) is the mean square displacement of Na + and Cl À ions during time t.In MD calculations of ionic diffusion the concentration of Na + and Cl À ions was kept constant as water beads were evaporated, and the same two-step procedure of water equilibration following the production runs as described above was applied.The CG interactions are much smoother and the dynamics are faster compared to the atomistic model.1][62] The difference in the time between the Martini CG model and AA simulations was discussed in numerous publications and the speed up factor was estimated to be between 0.79-22 (depending on the particular parameterization and CG implementation). 63,646][37] In subsequent work, other physical and chemical quantities such as the mass density, 37,40,45,65,66 the dielectric constant, 42 the radius of gyration, 40,42,65,67,68 film morphology 37 and the diffusion coefficient 40 were used to validate the CG model.In the present work for the validation of our CG model of PEDOT, we focus at the formation of p-p stacking between PEDOT chains, the number of chains in crystallites, the mass density and the SASA, and compare the CG results with corresponding AA simulations as well as with experimental values (where available).

PEDOT:Tos morphology
The aim of this section is to investigate the morphology of PEDOT:Tos within the CG model and compare it with the corresponding AA MD previously reported results.The oxidation level of PEDOT:Tos per monomer unit is considered in the range 0 o C ox o 33.3%; in most simulations the length of PEDOT chains is chosen to be N = 12 monomer units.(Note that the oxidation level of pristine, i.e. as synthesized, PEDOT typically corresponds to C ox = 33.3%, and the experimental value for the PEDOT chain length is N = 5-15 monomer units 8,69 ).
Fig. 2(a)-(d) shows snapshots of PEDOT:Tos for the hydration levels Hy = 80 and 40% and the oxidation levels C ox = 16.7 and 33.3%.All molecular configuration snapshots are prepared using the VMD package. 70A formation of small PEDOT crystallites containing several chains in p-p stacking and surrounded by Tos molecules is clearly visible.For higher oxidation levels, the Coulomb repulsion between PEDOT chains and the attraction between PEDOT-Tos separates PEDOT chains and some Tos molecules intercalate into crystallites as shown in the red boxes of Fig. 2(c).This is consistent with the results of recent AA MD simulations predicting the same phenomenon at high oxidation levels. 31Our results shows that PEDOT crystallization in the presence of Tos molecules occurs at the initial evaporation stage.During water evaporation, a denser polymer structure with more curvature along the PEDOT chains is obtained.The average end-to-end distance of the PEDOT chains decreases from 4.177 AE 0.009 for the solution with Hy = 80% to 3.846 AE 0.003 nm in the dry phase.(The length of a completely straight PEDOT chain for N = 12 is 4.27 nm).
To investigate PEDOT crystallization quantitatively, we calculate the distribution distance perpendicular to the PEDOT planes between virtual sites belonging to different chains with different oxidation levels and hydration levels and compare with present AA results. 16The distance distribution shows pronounced peaks at integer values of r/R p-p with R p-p = 4.6 Å being the p-p stacking distance between the adjacent chains, see Fig.

E R exp
p-p = 3.5 Å. 2,3,16 The discrepancy between the CG and the atomistic p-p stacking distances is related to the parametrization of the Martini force field mapping four atoms to one CG bead, making the size of the bead larger than the distance between the polymer chains.This discrepancy is known to be inherent to the Martini model and cannot be easily amended by adjusting the Martini force-field parametrization. 37However, it has been argued by Alessandri et al. 37 that the discrepancy between CG and AA stacking distances does not impact the nanomorphology evolution.This conclusion is also supported by the agreement in the features of the calculated CG and AA morphologies as will be discussed below in the present section.
The number of peaks in g(r) apparently corresponds to the number of PEDOT chains in the crystallites.Experimentally the crystallite size is determined from the width of the p-p stacking peak in grazing incidence wide-angle X-ray scattering (GIWAXS) spectra using the Scherrer equation. 71The PEDOT experimental crystallite sizes vary in different studies, from 1.2 nm, 8 to 3 nm, 72 and 4.5 nm, 2 which corresponds to 3-12 PEDOT chains.In our CG calculations the number of chains in the crystallites depends on the hydration and oxidation levels and varies between 3-9.For a given oxidation level the size of crystallites decreases when water is evaporated; at the same time, for a given water content, the larger the oxidation level, the smaller the crystallite size.This behavior is fully consistent with the corresponding results of AA MD calculations in Fig. 2.
For further validation of the CG results, we calculate the PEDOT:Tos mass density and compare it with our AA calculations, as well as with the available experimental data 27 and density functional theory (DFT) calculations for an ideal crystal, 73 see Table 1.The mass density of the AA model is very close to the experimental one. 27The difference between the CG and AA mass densities is lower than 10%.The lower mass density in CG is related to a larger bead size and a greater p-p stacking distance.A similar difference between the calculated CG and AA densities was found for other polymers, such as P3HT and PCBM. 37Note that for the ideal crystal in the DFT calculations 73 and orthorhombic model with experimental lattice parameters, 27 the densities, as expected, are higher than those obtained from the MD simulations.
Using the CG and AA models we calculate the solventaccessible surface area (SASA) as described in the model and method section.The SASA corresponds to an area obtained by rolling a probe sphere with the radius of the CG water molecule around PEDOT chains, as illustrated in Fig. 3(a and b).Fig. 3(c) shows the evolution of the SASA during PEDOT crystallization for different oxidation levels at Hy = 80%.The SASA is apparently maximal at the starting time of the production run, t = 0, when all randomly distributed PEDOT chains are accessible by solvents.As PEDOT chains crystallize, the SASA decreases from the initial value by nearly 50, 45 and 30 percent for C ox = 8.3, 16.7, 33.3%, respectively, see Fig. 3(c).The decrease of the SASA is related to formation of p-p stacking between chains where the inner surfaces are not included in the SASA, see Fig. 3(b).The decrease of the SASA is almost sudden during the first 20-30 ns initial time interval, which means that the crystallization mostly takes place during this time interval.The variation of the SASA with the hydration level is shown in Fig. 3(d) for the CG and AA models for different oxidation levels.There is a fair agreement between the calculated values of the SASA in the AA and CG models.Because of the difference between beads and atoms in CG and AA MD respectively, we do not expect the same values for the SASA, but we obtain the same trends for the two models.The CG model successfully reproduces the AA SASA behavior when the hydration and oxidation levels are varied.According to expectations, the SASA decreases as water evaporates because the evaporation results in a denser structure and thus a smaller SASA.At the same time, for a given hydration level, the crystallites are smaller for higher C ox , see Fig. 2. Thus, the available surface of PEDOT crystallites is larger, which results in the larger SASA for higher C ox , see Fig. 3(d).

Ion diffusion through PEDOT:Tos
We first start with the configuration of different components of the system (i.e.PEDOT chains, tosylate molecules, ions and water molecules) relative to each other.Fig. 4(a) shows the radial distribution function between positively charged PEDOT chains and negative Tos molecules, exhibiting peaks at the distances 0.53 nm, 0.74 nm, and 0.93 nm.The positions of the calculated peaks are in good agreement with the corresponding AA calculations. 16The radial distribution functions (RDF) for positive Na + ions and negative Tos molecules, as well as for negative Cl À ions and positive PEDOT chains, are shown in Fig. 4(b).The RDF function for Na + -Tos À shows a strong peak at the distance 0.5 nm, indicating that the majority of positive Na + ions are situated close to negative Tos À counterions.The RDF function for Cl À -PEDOT + shows peaks at the same   positions as those for the RDF function PEDOT + -Tos À in Fig. 4(a).The presence of several peaks in the RDF for Cl À -PEDOT + (PEDOT + -Tos À ) is apparently related to the fact that PEDOT is an extended molecule, such that the first peak corresponds to the distance from the ion (Cl À or Tos À ) to the neighboring monomer, and the remaining peaks to more distant monomers.
Let us now study the diffusion of Na + and Cl À ions in PEDOT:Tos for different hydration and oxidation levels.We first calculate the diffusion coefficients of Na + and Cl À ions at a concentration of 100 mM using the developed Martini CG MD model in pure water and compare them the with corresponding AA MD simulations reported previously 74 and experimental values in the infinitely diluted solution, 75-77 see Table 2.
There is a fair agreement between the experimental, CG, and AA diffusion coefficients.Also, as in the case of the SASA (calculated in the previous section), we do expect exact quantitative agreement between the AA and CG results.The previous AA simulations showed a system-size dependency 80,81 and an effect of the truncation scheme of non-bonding interactions 82 on the diffusion coefficient.In our calculations the system-size dependence is negligible because of a large size of the computational box, see Table S2 in Section SI3 (ESI †).For the truncation of the electrostatic interactions we use a standard PME method which is shown to be the most accurate one for the MD simulations with the polarizable Martini water model 42,44,83 (for comparison of different truncation schemes see Table S3 in Section SI3, ESI †).
To calculate the diffusion coefficients D Na/Cl of Na + and Cl À ions we add them to the system at each hydration level, keeping their concentration of 100 mM constant and use the Einstein relation for the last 10 ns of the ion's trajectory.In the calculation of the diffusion coefficients, the simulation time is chosen to be long enough to ensure that the diffusive regime is established (i.e.hDr 2 (t)i/t in the Einstein relation, eqn (1), saturates and becomes independent of time, see Fig. S4 in Section SI4, ESI †).Here we should mention that AA study of the diffusion process for our large simulation box and long time scale is computationally very expensive.Fig. 5 presents evolution of the diffusion coefficients as the hydration level Hy changes.It exhibits exponential dependence of the diffusion coefficients on the hydration level; the value of D Na/Cl decreases by three orders of magnitude when the hydration level is reduced from 80% to 10%.The exponential dependence of the ionic diffusion coefficients on the water content represents one of the main findings of the present study.We also find that the diffusion coefficient depends on the oxidation level.For example, for Hy = 10%, the diffusion coefficient of both the Na + and Cl À ions decreases by a factor of B2.5 when the oxidation level increases from C ox = 0 to C ox = 33.3%.For a given oxidation level the calculated ion diffusion coefficients can be fitted to an exponential function, where D Na/Cl 0 is the diffusion coefficient of the Na + /Cl À ion in pure water, and l Na/Cl is the fitting parameter which depends on the oxidation level C ox .For oxidation levels C ox = 0-33.3%,l varies continuously between 5.1-5.4 and 5.8-8.4 for Na + and Cl À , respectively.The value of l increases for larger C ox which reflect the role of Coulomb attraction between the polymer and ions.The polymer matrix traps ions more effectively at a higher oxidation level due to the stronger Coulomb attraction, see a more detailed discussion below.For all oxidation levels, l Cl À 4 l Na + , which reflects a stronger interaction of Cl À ions with series of positive charged beads in the PEDOT chains (see Fig. 4) and consequently stronger Cl À localization with respect to Na + .The large PEDOT chains and spatially distributed positive charge on the monomers make PEDOT more effective than Tos molecules in ion attraction.
To understand the origin of the exponential behavior of the diffusion coefficient, let us focus on the evolution of the water clusters in the polymer matrix as the hydration level changes.We define the cluster as being composed of water molecules Fig. 5 The diffusion coefficients of (a) Na + and (b) Cl À ions in PEDOT:Tos with different hydration levels and oxidation levels.(The diffusion coefficients are normalized to the corresponding calculated CG Na + and Cl À ionic diffusion coefficients in pure water as given in Table 2).Dashed lines are the exponential fit based on eqn (2).
where for any molecule in the cluster one can find another molecule at the distance smaller than a critical distance; the critical distance is set to the diameter of the Martini CG water bead.(For pure water all water molecules apparently belong to a single cluster).As the hydration level decreases the clusters get fragmented into smaller clusters as illustrated in Fig. 6(a-c).
The size of water clusters decreases as the hydration level is lowered, and for charged polymers (i.e.C ox a 0) this decrease is exponential, see Fig. 6(d).At higher hydration level, there are wide diffusion paths of water molecules for ions to diffuse through the polymer matrix, Fig. 6(a).As water evaporates, the width and number of diffusion paths decrease and water clusters with weak connections are formed, and many clusters become disconnected, Fig. 6(b) and (c).Because ions primarily move in the water clusters, this leads to a decrease of the diffusion coefficient.The size of water clusters is also dependent on the PEDOT oxidation level.By adding more charges to PEDOT and Tos ions, the charged beads of polarizable water are attracted to the polymer matrix which decreases the water cluster size.
The radial distribution functions of water molecules around Na + and Cl À ions show peaks corresponding to formation of the hydration shells around the ions, see Fig. 6(e).The first and second shells include water molecules situated respectively at the distances 0.45 o r o 0.65 nm and 0.65 o r o 1.2 nm from the ions.The number of water molecules in each shell decreases as the hydration level is reduced, see Fig. 6(f).The decrease of the number of water molecules in the hydration shells is connected to the thinning of the percolative paths as discussed above.Also, this decrease leads to the reduction of electrostatic screening of ions by water molecules, which, in turn, leads to the enhancement of the electrostatic interaction between ions and the polymer matrix (i.e.PEDOT + chains and Tos À molecules).
Further insight into the behavior of the ionic diffusion coefficient can be obtained by investigating the relative position between the ions and polymer matrix (consisting of PEDOT + chains and Tos À molecules).Fig. 7  Tos À (or PEDOT + ).There is however a strong qualitative difference in the tails of these distributions between the cases of high and low hydration levels Hy.For higher hydration levels (e.g.Hy = 80%) the calculated distribution shows a broad tail, indicating that many ions move quite far away from the polymer matrix.However, for Hy = 20% the distribution tends to zero already at the distance E0.5 nm, indicating that practically all ions are confined in the vicinity of PEDOT + or Tos À , see Fig. 7(a) and (b).This is clearly seen in the animations presented in the ESI, † SI5 (see below for a more detailed discussion of the animations).
To understand the importance of this difference, let us consider the average mean square displacement hDr provides a molecular explanation of the pronounced decrease of the diffusion coefficient for lower water content presented in Fig. 5. Indeed, for small hydration levels, practically all ions are confined in the vicinity of the polymer matrix where the mean square displacement hDr d2 i is small.As a result, the corresponding diffusion coefficient, eqn (1), is also small.For high hydration levels, many ions move away from the polymer matrix where the mean square displacement hDr d 2 i is large.
This apparently results in the large diffusion coefficient.This molecular picture of the diffusion is illustrated in the trajectory snapshots for Hy = 80% (Fig. 7(d)) and Hy = 20% (Fig. 7(e)) and in two animations presented in the ESI.† Animation 1 correspond to the case of low water content (Hy = 20%).All ions in the system are confined close to the polymer matrix and execute a ''trembling'' motion with a very low displacement.In contrast, for high water content (Hy = 80%), an ion wanders around the whole system until it finally gets trapped in the vicinity of the polymer matrix, see Animation 2. (Note that after some time the ion will be ''released'' again and will continue its motion in the system).The ionic conductivity s is directly proportional to the diffusion coefficient of anions and cations in the solution through the Nernst-Einstein equation, s ¼ nq 2 D k B T , where q is the ionic charge and n the ion concentration. 45,84Our results predicting the exponential enhancement of ionic conductivity in the polymer as the water content is increased can be tested in the laboratory.The conductivity enhancement with increasing water content is reported experimentally for Na + ions in PEDOT:PSS,  and nanofibrillated cellulose-PSSH. 89In the low molecular weight polymer LiCF 3 SO 3 PEG 10 the anion and cation diffusion coefficients increase exponentially by two order of magnitude on adding water. 90he increase in ionic conductivity is even more intense for high molecular weight polymer Pb(CF 3 -SO 3 ) 2 PEO 16 . 91

Conclusion
The ionic and electronic transport properties of highly doped conducting polymers, in particular of PEDOT, are strongly dependent on the material morphology.Massive experimental efforts are currently underway in order to investigate the complex morphology of PEDOT.At the same time, the corresponding theoretical modeling and simulations that experiments can rely upon are practically missing.Recently, all atomistic MD simulations of the morphology of doped PEDOT:Tos and its crystallization in aqueous solution have been reported. 16For many systems the utilization of all atomistic MD becomes computationally prohibitive, and the theoretical investigation of the morphology of conducting polymers on a scale exceeding tens of nanometers requires coarsegrained approaches where the system is represented by reduced numbers of the degrees of freedom.In the present study, we develop a Martini coarse-grained MD model for the doped conducting polymer PEDOT.The coarse-grained model is validated using a comparison with corresponding all atomistic MD simulations.In particular, we calculate the size of PEDOT crystallites; the radial distribution functions between PEDOT and Tos counterions, PEDOT and Na + and Cl À ions, and Tos and Na + and Cl À ions; the solvent-accessible surface area; and the PEDOT:Tos mass density.
Using the developed coarse grained model we study diffusion of Na + and Cl À ions in PEDOT:Tos.We find that the ionic diffusion coefficients are decreased exponentially when the hydration level is reduced, exhibiting a drop of almost three orders of magnitude when the hydration level Hy is decreased from 80% to 10% wt.We relate this behavior to the evolution of the water clusters in the polymer matrix.At higher hydration level, there are wide percolation paths of water molecules for ions to diffuse through the polymer matrix.As the hydration level decreases, the water clusters become fragmented into smaller and disconnected clusters and their sizes decrease exponentially.Another theoretical prediction is that the diffusion coefficients decrease when the doping level of PEDOT is increased.For example, for the hydration level Hy = 10%, the ionic diffusion coefficients in fully reduced PEDOT are by a factor of B2.5 larger than the diffusion coefficients in fully oxidized PEDOT at C ox = 33.3%.
A complementary insight into the observed exponential behavior of the diffusion coefficient is obtained from the analysis of the relative position of ions and the polymer matrix and the average mean square displacement hDr d 2 i for an ion that starts its motion at the distance d from the polymer matrix.We demonstrate that at high hydration levels many ions move quite far away from the polymer matrix, and hDr d 2 i is close to the value corresponding to pure water.However, for low hydration levels, most ions ions are confined in the vicinity of the polymer matrix where the mean square displacement hDr d 2 i is negligible.As a result, the corresponding diffusion coefficient is decreased exponentially.
The predicted exponential dependence of the ionic diffusion coefficients on the hydration level and the predicted decrease of the diffusion coefficient with the increase of the oxidation level can be tested experimentally.We believe that our findings concerning the behavior of the ionic diffusion coefficients apply not only to PEDOT:Tos, but to a wider class of doped polymers with molecular counterions including polypyrrole, polyaniline, thiophene, etc., which exhibit similar character of disordered amorphous morphology with limited crystallinity.We also believe that the understanding of ionic diffusion on the molecular level provided in this study is important for designing devices with mixed electron-ion conductivity (such as supercapacitors, organic electrochemical transistors, and ion pumps) where the ionic diffusion determines the device performance.

Fig. 1
Fig. 1 (a) The atomic structure and the CG model of (a) PEDOT and (b) TOS.The beads (SC1, SC2, and SNda for PEDOT and SCY, STY, and Q a for Tos) are indicated in grey; virtual sites are shown in red.
2(e)-(h).(Snapshots and the distance distribution chains of lengths N = 6, 12 and 18 are shown in ESI, † SI2).The value of R p-p = 4.6 Å obtained from CG simulations is somehow larger than the one obtained from both the AA MD simulations and experimental measurements, R AA MD p-p

Fig. 2
Fig. 2 (a-d) Snapshots of PEDOT:Tos for different oxidation and hydration levels.PEDOT is shown in blue and Tos in green; water molecules are not displayed.Red boxes in (c) outline Tos molecules intercalated between PEDOT chains.(e-h) The distance distribution between virtual sites of PEDOT for different oxidation and hydration levels.The solid lines are CG results and dashed lines correspond to the AA simulations reported in ref. 16.

Fig. 3 (
Fig. 3 (a and b) A schematic illustration of the SASA for an representative crystallite composed of three PEDOT chains seen from different directions; the SASA probe spheres are in gray, and the PEDOT chains are in blue.(c) Evolution of the SASA during PEDOT crystallization.(d) The SASA for PEDOT:Tos in the AA (dash lines) and CG (solid lines) models for different hydration levels and oxidation levels.In (c) and (d) the SASA is normalized to the number of PEDOT chains in the box for both the AA and CG models.
(a) and (b) show the distribution of the minimum distance d between Cl À (Na + ) ions and the closest bead of PEDOT + (Tos À ) molecules for different hydration levels (a definition of the minimum distance d is visualized in the inset of Fig. 7(b)).The minimum distance distribution has a sharp peak around 0.5 nm for all hydration levels, which shows that most positive (or negative) ions are localized close to the corresponding

Fig. 6
Fig.6 (a-c) Water clusters in the PEDOT:Tos polymer matrix for different hydration levels Hy = 40%, 20%, and 10% respectively.(d) Variation of the average water cluster size with the hydration level (normalized to the cluster size at Hy = 100%) for different oxidation levels.(e) The distribution function g(r) of water molecules around Na + ions for Hy = 80% and Hy = 20% (the corresponding distribution functions for Cl À ions are practically the same).The inset shows the water shells around a single Na + ion with four different (blue-green-orange-cyan) colors.(f) The number of water molecules in the first and second shells as a function of the hydration level Hy for Na + (solid line) and Cl À ions (dashed line).

d 2 i
for an ion that starts its motion at the distance d from the polymer matrix (a definition of hDr d 2 i is visualized in the inset of Fig. 7(b)).Fig. 7(c) shows a dependence of hDr d 2 i on d at Hy = 80, 60 and 40%.For Hy = 80% and at sufficiently high distances, d \ 3 nm, hDr d 2 i is close to the value corresponding to pure water, i.e. ions do not feel the presence of a polymer matrix.The closer ions are situated to the polymer matrix, the smaller hDr d 2 i is.This behavior of hDr d 2 i, combined with the distribution of the minimum distance d in Fig. 7(a) and (b),

Fig. 7 2 i
Fig. 7 The distribution of a minimum distance d between (a) Na + -Tos and (b) Cl À -PEDOT for C ox = 16.7% and different hydration levels.(c) The mean square displacement hDr d 2 i of Na + and Cl À ions as a function of the initial distance d from the polymer matrix for Hy = 80, 60 and 40%.The inset in (c) visualizes the definition of hDr d

Table 1
The comparison between experimental and theoretical AA and CG densities of PEDOT:Tos