Kenley M.
Pelzer
*ab,
Álvaro
Vázquez-Mayagoitia
c,
Laura E.
Ratcliff
c,
Sergei
Tretiak
d,
Raymond A.
Bair
efg,
Stephen K.
Gray
af,
Troy
Van Voorhis
h,
Ross E.
Larsen
i and
Seth B.
Darling
aj
aCenter for Nanoscale Materials, Argonne National Laboratory, 9700 Cass Ave., Lemont, IL 60439, USA. E-mail: kpelzer@anl.gov; Tel: +1-630-252-7020
bMaterials Science Division, Argonne National Laboratory, 9700 Cass Ave, Lemont, IL 60439, USA
cArgonne Leadership Computing Facility, Argonne National Laboratory, 9700 Cass Ave., Lemont, IL 60439, USA
dTheoretical Division, Center for Nonlinear Studies, Center for Integrated Nanotechnologies, Los Alamos National Laboratory, Los Alamos, NM 87545, USA
eMathematics and Computer Science Division, Argonne National Laboratory, 9700 Cass Ave., Argonne, IL 60439, USA
fComputation Institute, University of Chicago, 5735 S. Ellis Ave., Chicago, IL 60637, USA
gComputer, Environment, and Life Sciences, Argonne National Laboratory, 9700 Cass Ave., Lemont, IL 60439, USA
hDepartment of Chemistry, Massachusetts Institute of Technology, 77 Massachusetts Ave, Cambridge, MA 02139, USA
iComputational Science Center, National Renewable Energy Laboratory, 15301 Denver W. Parkway, Golden, CO 80401, USA
jInstitute for Molecular Engineering, University of Chicago, 5747 S. Ellis Ave., Chicago, IL 60637, USA
First published on 11th January 2017
Organic photovoltaics (OPVs) are a promising carbon-neutral energy conversion technology, with recent improvements pushing power conversion efficiencies over 10%. A major factor limiting OPV performance is inefficiency of charge transport in organic semiconducting materials (OSCs). Due to strong coupling with lattice degrees of freedom, the charges form polarons, localized quasi-particles comprised of charges dressed with phonons. These polarons can be conceptualized as pseudo-atoms with a greater effective mass than a bare charge. We propose that due to this increased mass, polarons can be modeled with Langevin molecular dynamics (LMD), a classical approach with a computational cost much lower than most quantum mechanical methods. Here we present LMD simulations of charge transfer between a pair of fullerene molecules, which commonly serve as electron acceptors in OSCs. We find transfer rates consistent with experimental measurements of charge mobility, suggesting that this method may provide quantitative predictions of efficiency when used to simulate materials on the device scale. Our approach also offers information that is not captured in the overall transfer rate or mobility: in the simulation data, we observe exactly when and why intermolecular transfer events occur. In addition, we demonstrate that these simulations can shed light on the properties of polarons in OSCs. Much remains to be learned about these quasi-particles, and there are no widely accepted methods for calculating properties such as effective mass and friction. Our model offers a promising approach to exploring mass and friction as well as providing insight into the details of polaron transport in OSCs.
Because of the relatively large size of OPV systems (with a typical active layer >100 nm in thickness), modeling of charge transport on the device scale presents a daunting problem in terms of computational cost. Although there have been major developments in linear-scaling density functional theory (DFT),17 the tens of millions of atoms in a representative volume are prohibitively costly for dynamical methods such as time-dependent DFT. Even the simpler ground state calculations would be prohibited by memory issues on most supercomputers. Here, we propose a simplified approach in which charge transport is modeled with Langevin molecular dynamics (LMD).18 LMD is fundamentally a classical method, while charges are traditionally assumed to require quantum mechanical treatment. However, we argue that in some organic semiconductors, the unique properties of charges in these materials permit an accurate classical treatment.
In many organic electronic materials, charges are coupled to surrounding nuclei. The Coulombic interaction between the negative charge of an excess electron and the positive charge of surrounding nuclei pulls nuclei toward the charge (or in the case of a positively charged hole, repels the nuclei). This interaction with the nuclei creates a “polaron,” a charge dressed with phonons, which then transfers between molecules via a thermally activated hopping process. Understanding the electronic structure of these charges and their coupling to molecular motions presents a major challenge in theoretical modeling.15,19–25 However, in this challenge there may also lie an opportunity for simplified modeling. We argue that because of the mass of the nuclei to which charges are coupled, polarons can be treated as localized quasi-particles with an effective mass greater than that of bare charges. With this greater mass, a classical treatment can be applied in which the polaron acts as a pseudo-atom obeying Newton's laws. Moreover, if the dynamics of the polaron pseudo-atoms can be modeled on a potential energy surface defined by the atoms in the OPV active layer (allowing the properties of the active layer atoms to be considered without explicit treatment of their dynamics), the scale of the problem is drastically reduced. The typical OPV charge density of 1021 to 1023 polarons per m3 indicates that a 100 nm thick region would only contain 1–100 polarons.16,26–28 Here, we present a multiscale approach grounded in constrained density functional theory (CDFT)29 calculations. We demonstrate that once a set of CDFT calculations are conducted to build a potential energy surface (PES) for LMD charge transfer, the transfer process itself can be modeled at an almost negligible computational cost.
In this work, we do not tackle the full challenge of modeling a large region of an OPV, but rather we present the building blocks of such a method by studying the charge transfer properties of a pair of fullerenes. We apply our model to phenyl-C61-butyric acid methyl ester (PCBM) molecules, with the crystal structure shown in the upper part of Fig. 1.30 Fullerenes are widely used in organic semiconductors, with PCBM being one of the most commonly used electron acceptors in OPVs. PCBM is also a convenient choice for developing our model due to recent work suggesting that excess charges are mostly localized to a single molecule in PCBM materials,31 allowing a straightforward application of our model in which charges are assumed to hop from one molecule to the next in localized states. Defining the extent of delocalization of a charge in a given material is a challenge in itself (and in fact there is conflicting literature regarding PCBM, with ref. 32 suggesting that the charge tends to delocalize over multiple molecules). If a polaron is delocalized over multiple molecules but still moves via thermally activated hopping, our model is still in principle applicable. However, the computational cost would increase, as a more delocalized state would require a larger number of atoms in each CDFT calculation. Here for convenience we work with a small polaron restricted to a single PCBM molecule.
![]() | ||
Fig. 1 Here we show the crystal structure of two PCBM molecules as obtained by ref. 29. Below we show the difference in α and β spin densities for the two constrained structures used to construct the reaction coordinate. We see that CDFT produces an anionic electronic state in which one excess charge with α spin is constrained to a chosen set of atoms. |
We emphasize that the model we present of a polaron traveling via hopping is not an accurate description of charges in every organic semiconductor. The nature of charge transfer in these materials is a very active area of research. Dynamics of polaron transport derived from a master equation approach demonstrate that polarons may exhibit hopping or band-like transport depending on parameters such as system-phonon coupling strength and temperature.33 Some materials have shown signs of band-like transport, while other literature suggests that charge transport in organic semiconductors is not fully described by a hopping or band-like transport model.34–36 Our description of the charge as a classical particle governed by the dynamics of LMD is not appropriate for these cases. Rather, our model provides an effective description for the case of polaronic charges that transfer between molecules via a hopping process and are effectively weighed down by their coupling to nuclei. We argue that in such cases, the polaron can be conceptualized a classical pseudo-atom with unique opportunities for theoretical modeling.
The use of CDFT to produce diabatic PESs, and, via diagonalization of the resulting Hamiltonian, produce adiabatic PESs, is well established. For example, this method has been used to calculate charge transfer rates via Marcus theory.29,37–42 In 2000, Kim et al. suggested that a charge could be treated via the Langevin equation.43 Beyond CDFT and LMD, many other models of charge transport in organic semiconductors have been proposed (see ref. 14, 19, 35 and 44–47 for thorough discussions of previous work). However, to our knowledge no one has conducted molecular dynamics simulations of polaron dynamics using CDFT PESs. Here we present such simulations and the resulting predictions of transfer rates. Going beyond transfer rates that describe charge transfer in a material with a single number, we also show how the simulation data provide insight into exactly when and how individual transfer events occur. We demonstrate that the LMD simulation approach is unique both in the information that can be incorporated into the model and the information that can be extracted from the results.
Our PES describes the energy of the pseudo-atom as it moves from one molecule to another (the process of polaron transfer). Because our pseudo-atom describes a charge coupled to nuclei, this transfer process is described by the gradual change in the positions of the PCBM atoms that occurs as the charge moves in space from the initial state to the final state. We can think of this gradual change in nuclear positions as the “reaction coordinate” for polaron transfer. The PES is created as a function of this reaction coordinate. The PES as a function of the reaction coordinate can then be transformed to a PES as a function of a spatial coordinate, which is then utilized in LMD simulations. This derivation of the PES is discussed in detail below.
The propagation of a particle in one dimension in the LMD formalism is described by the Langevin equation:18
mẍ = −∇U(x) − γmẋ + R(t) | (1) |
(1) Perform CDFT calculations, as well as unconstrained DFT calculations, to find the optimized geometries of the initial state, final state, and transition state. These three optimized structures are used to form a “reaction coordinate”, where the reaction coordinate describes the gradual change in nuclear positions as the system moves from the initial state optimized geometry to the final state optimized geometry.
(2) Perform single-point CDFT calculations for a series of geometries along the reaction coordinate. There are two sets of single-point CDFT calculations. In one, the charge is constrained to a single fullerene (the initial state). In the next, the charge is constrained to the other fullerene (the final state). By implementing these constraints while gradually changing the geometries, we obtain diabatic PESs for charge transfer.
(3) Use diabatic PESs to compute the adiabatic PESs. The lower adiabat is then used to create our PES for LMD.
(4) This adiabatic PES, like the diabatic PESs, is a function of the reaction coordinate. However, we wish to obtain a pseudo-atom PES as a function of the spatial coordinate: the movement of the pseudo-atom in space as it travels over a distance of several angstroms from one fullerene to the next. In this fourth step, we convert the PES from a reaction coordinate to a spatial coordinate.
(5) Effective mass and friction are assigned (for PCBM, we will explore a range of values).
(6) The spatial-coordinate PES, the effective mass, and the friction are then used as input for our LMD simulations. Because the LMD equation involves a random force term (representing the thermal effects of the environment), many simulations should be performed to obtain an appropriate sampling of trajectories. Here we perform 500 simulations for each point in our parameter space (where the effective mass and friction terms are varied).
(7) LMD results are then analyzed to compute charge transfer rates and examine the details of polaron movement through space.
We describe this procedure in detail in the following sections, with each section presenting a given step in formalism and then discussing its application to PCBM.
As discussed above, it is also possible to use initial and final states in which a more delocalized charge is constrained to a set of several molecules. In this article, for the sake of simplicity, our language will refer to bimolecular systems in which a charge is constrained to a single molecule.
The optimized geometries are represented with arrays R, where each array contains the xyz coordinates of each atom in the PCBM molecules. The geometries corresponding to the reaction coordinate Rq are then extracted from the optimized geometries of the initial state (R1), transition state (RT), and final state (R2) for the charge transfer event. This gradual change can be described by the following interpolation suggested by Wu et al.:40
![]() | (2) |
In implementing step #1 for the PCBM dimer, the crystal structure of two PCBM molecules shown in the upper part of Fig. 1 is used as the starting structure for all optimizations. Although the disordered nature of OPVs almost certainly leads to significant variability in the relative orientations of PCBM molecules, here for simplicity we work only with the crystal structure. Other work applying DFT to the PCBM crystal structure shows that use of this geometry produces efficiency results consistent with experimental observations of charge mobility in PCBM.48 The parameters of the DFT calculations are described in the ESI.† Due to the difficulty of defining an implicit solvent that represents a solid state OPV environment, the DFT calculations are conducted in the gas phase (with effects of the environment inherently included in the mass and friction of the LMD formalism).
A full-dimensional LMD calculation would obviously require a full-dimensional PES, which, although in principle possible, would be a major endeavor to derive from a set of CDFT calculations. Here we suggest a simpler approach: we perform LMD for a charge transfer event in which the charge is assumed to move in a straight line from one molecule to another, allowing a one-dimensional (1D) PES to be used. A 1D approach to modeling polaron transfer has been successfully applied in other work studying organic semiconductors, where a master equation approach to modeling polaron transfer demonstrated good agreement with experimental charge transfer rates when applied to a 1D model system.33
In applying steps #2 and #3 to PCBM, some experimentation (including the variation of density functional models, constraints, and omission of outlying points) was necessary to obtain physically reasonable PESs. We describe this process in detail in the ESI,† noting that such careful treatment of CDFT results may frequently be necessary to create appropriate PESs for this model. Ultimately we obtained the diabats and adiabats shown in Fig. 2, which give energies as a function of the reaction coordinate. These diabats are obtained from 17 points along the reaction coordinate; because a continuous adiabat is needed for LMD, all simulations employed a cubic spline interpolation to create a continuous PES. As discussed above, this reaction coordinate describes the change in positions of the atoms in the PCBM molecules as a polaronic charge moves from one molecule to another. As described in eqn (2), this reaction coordinate is a function of the interpolation parameter q, which is shown on the x-axis of Fig. 2. (While a range of −1 ≤ q ≤ 1 is used in the creation of the reaction coordinate, the omission of outlying points described in the ESI† leaves the range of −1 ≤ q ≤ 0.6 shown in Fig. 2.) This charge transfer event is illustrated by the CDFT spin density images shown below the plot.
![]() | ||
Fig. 2 Diabatic and adiabatic PESs obtained from CDFT calculations. These PESs are a function of reaction coordinate rather than space. The x-axis shows the value of the interpolation parameter q of eqn (2) that is used to derive the reaction coordinate. Diabat #1/2 corresponds to constraint of the charge to fullerene #1/2 as referenced in the text. Images of the initial and final spin densities are shown to represent the change in the system as it moves along the reaction coordinate. As indicated by the black arrows, transport in our model represents movement along the lower energy adiabat. |
We note that the PES used for the calculations is asymmetric, with an adiabatic final state energy 100.4 meV lower than the adiabatic initial state energy, as shown in Fig. 2. This deviation will lead to somewhat higher transfer rates relative to a symmetric PES. The asymmetry of the PES is partly an artifact of the removal of outlying points, as discussed in the ESI.† However, we also found that the CDFT calculations consistently produced different diabatic energies for the initial and final states, in spite of the dimer system having Ci symmetry. This deviation reflects the challenge of performing geometry optimizations for the large systems involved in OPVs, where geometries are prone to become trapped in local minima. One could force symmetry into the reaction coordinate by performing a single CDFT calculation and applying the result to both the initial and final states via symmetry operations; however, there is nothing especially desirable about a symmetric PES for methods development purposes. An asymmetric PES would actually be the norm in real organic semiconductors because they frequently are disordered, so we view the convergence-based asymmetry as an asset for testing the overall LMD approach. Of course, the disordered nature of OPVs must lead to great variation in the PESs between different pairs of molecules, and we wish to be clear that the goal of this work is not to explore the vast array of geometries and PESs that are present in a disordered OPV. We simply seek one reasonable PES for the purpose of developing and testing the LMD model.
The conversion to a spatial coordinate is straightforward, using the variables c1 and c2 that describe the contribution of each diabat to the adiabat at a given point on the reaction coordinate. The variables c1 and c2 are calculated when converting diabats to adiabats following ref. 40. With the positions of the fullerenes denoted as x1 and x2, the spatial coordinate of each adiabatic energy point is calculated as:
![]() | (3) |
We thus obtain an adiabatic PES describing the energy of our pseudo-atom as it moves through space.
In our PCBM calculations, x1 and x2 are set 10 Å apart in accordance with the center-to-center distance in PCBM crystals.30Fig. 3a shows the lower adiabat converted from a reaction coordinate to a spatial coordinate as described by eqn (3). C60 cages are added to show the approximate location and size of the fullerenes. To create a complete PES for our LMD simulations, the PES was mirrored around the starting point and then expanded in space with periodic boundary conditions, as shown in Fig. 3b. Our LMD simulations begin with the pseudo-atom at point B, and we consider transfer to have occurred if the pseudo-atom reaches point A or C.
The creation of the PES described by steps #1–4 is summarized by the flow chart in Fig. 4. Beginning at the top of the flow chart, we see that the geometries of the initial, final, and unconstrained states are calculated, where yellow ovals indicate the atoms over which the charge is allowed to delocalize. These geometries are then applied viaeqn (2) to create a reaction coordinate Rq, using the interpolation parameter q that varies from −1 to 1. The coloring of the arrays R1, RT, and R2 shown in cyan, magenta, and purple in eqn (2) correspond to the coloring of the optimized geometries in the top panel. We then compute diabatic and adiabatic PESs as a function of this reaction coordinate, shown in this figure by a simplified representation of Fig. 2. Eqn (3) is then applied to convert the lower adiabatic PES from a reaction coordinate to a spatial coordinate xsp, with A1 and A2 giving the contributions of diabats #1 and #2 to the adiabat at each point on the reaction coordinate. The coloring of the variables x1 and x2 shown in red and blue correspond to the coloring of the diabats in the representation of Fig. 2. The adiabat as a function of spatial coordinate is then shown by a simplified representation of Fig. 3a. This adiabat as a function of space is expanded in one dimension to provide a PES that can be used for LMD, as shown by a simplified representation of Fig. 3b.
![]() | ||
Fig. 4 From top to bottom, a flow chart summarizes the procedure (described by steps #1–4 in the text) for creating a PES suitable for LMD simulations. First, we compute CDFT geometries of the initial and final states for charge transfer, as well as an unconstrained DFT geometry. In this top panel, we show these geometries with yellow ovals indicating the atoms over which the charge is allowed to delocalize. Eqn (2) converts this information to a reaction coordinate Rq using the interpolation parameter q that varies from −1 to 1. We then create diabatic and adiabatic PESs as a function of this reaction coordinate. Eqn (3) converts the lower adiabat from a reaction coordinate to a spatial coordinate xsp, with A1 and A2 giving the contributions of diabats #1 and #2 to the adiabat at each point on the reaction coordinate. This adiabat as a function of space is then expanded in one dimension to create a full PES for LMD. |
In our PCBM calculations, for simplicity we hold meff constant throughout the simulation. However, there is no reason why meff could not be varied as the pseudo-atom moves along the spatial coordinate, and this may in fact be more accurate and offer an advantage to the LMD simulation approach. In regions corresponding to a steeper diabat (reflecting a sharp increase in energy with movement of the nuclei), an increased effective mass may be appropriate. Friction seems more likely to be constant for a given temperature because it is a manifestation of the disordered vibrations of many atoms, but within our formalism both mass and friction can easily be varied as a function of space or time.
For our simulations we choose a wide range of trial values of meff that vary over several orders of magnitude from 7.22 × 10−2 to 7.22 × 103 amu, where the mass of an electron is 5.49 × 10−4 amu. We emphasize that this effective mass is neither the mass of the electron nor the mass of the nuclei to which is it coupled, but rather describes the properties of the quasi-particle that is defined by the coupling between charges and nuclei. In the wide range that we use in our calculations, it is unlikely that the polaron effective mass is well described by either the highest or lowest values. However, exploring results over a wide range of masses is interesting for the purposes of model development and comparison to experimental benchmarks. The variation in γ is more complex, as γ is a function of meff, with a heavier particle influenced less by random fluctuations in the environment. For our calculations a friction γ = pγD was utilized, with γD as a function of meff derived as shown in the ESI.† This γD is proportional to the square root of meff, and friction was varied over three orders of magnitude with p values of 0.036, 0.36, and 3.6. We found that this method optimized transfer rates, as discussed below.
In our calculations of PCBM, all simulations are conducted at room temperature, with an initial random velocity assigned from a normal distribution with a mean of 0 and standard deviation of . For each point in the meff/γ parameter space, 500 calculations were performed, all with the pseudo-atom beginning at point B of Fig. 3b. For all values of meff and γ presented here, the simulation times tmax were long enough to allow all 500 simulations to achieve transfer to point A or C. We do not present data for meff/γ values for which all 500 runs did not achieve transfer within 5 ns—in such cases, it is not straightforward to calculate a transfer rate that is an appropriate comparison to the
rates described above. Regardless, a system in which no transfer occurs within 5 ns would, by the Einstein relation discussed immediately below, correspond to a charge mobility several orders of magnitude lower than that needed for a functional OPV device, and therefore is not of practical interest.
With the dashed line we show the Kramers' theory rate, a rate that is derived from the Langevin equation not via simulation but via an analytical derivation that is valid given the assumption of moderate to strong friction.74 We adapt the rate equation of Ensing:75
![]() | (4) |
The frequencies ωR/B are derived from a harmonic oscillator treatment of the PES in the region of the reactants and the barrier, with for the reactants and
for the barrier. γ is a function of meff as described in the ESI.†EB corresponds to the difference in U between the reactants and the energy peak of the adiabat and β = 1/kBT. In this work we seek a method that is robust over a broad range of friction values and provides trajectories in addition to just a single rate. However, the Kramers' theory rates are a useful check to ensure that our application of LMD gives proper functional dependence.
The simulation data and Kramers' theory results show almost exactly the same relationship between pseudo-atom effective mass and transfer rate. It is not obvious why the Kramers' theory rates are shifted a bit higher than the LMD simulation rates, and we suggest two explanations. One is the fact that Kramers' theory relies on assumptions about friction, while our simulations work directly with the LMD equation without simplifying premises. The other limitation of Kramers' theory is that it does not contain complete information on the PES; rather, the frequencies that enter the equation are derived from information on the PES only in the close proximity of the reactants and the barrier. In contrast, the complete shape of all regions of the PES is incorporated into the LMD simulations.
Another advantage of our simulation approach is that we can observe the distribution of rates in an ensemble of systems, rather than just observing a single rate that provides only averaged information. In Fig. 6, we present normalized histograms of transfer time tt in ps for all 500 calculations, where the y-axis gives the probability density of a particular transfer time. The mean and median of the transfer times are noted for each mass. For all masses the histograms show a long right tail and a mean transfer time that exceeds the median transfer time by about 20%. These results demonstrate that although the majority of the calculations achieve transfer in a relatively short period of time, a few outlying cases exhibit high transfer times. The ease of obtaining such information on ensemble distributions is a useful feature of our model.
In Fig. 7, we show the relationship between γ and k for three values of γ, varying γ by an order of magnitude with each data point. Separate lines are shown for varying values of meff. For values of meff from 7.22 × 10−1 to 7.22 × 103 amu, we see a non-monotonic relationship between friction and efficiency, with the highest efficiency found at the intermediate value of p = 0.36. This non-monotonic relationship (“Kramers turnover”) is expected, having been predicted as early as 1940 by Kramers for the movement over a barrier of a particle subject to friction.74 We again see results consistent with experimental estimates, with rates in the 109 to 1012 s−1 range for these boundaries. The case that varies slightly from the others, with higher rates and a monotonic decrease in this range of meff, is that with the lowest meff of 7.22 × 10−2 amu.
As mentioned above, CDFT diabats can be easily used to calculate rates of charge transfer at the level of Marcus theory.76 Comparing the LMD rates to Marcus theory sheds some light on the challenge of defining meff. In Marcus theory the rate is:
![]() | (5) |
While Marcus theory is more effective in modeling a system at low friction, the Kramers' rate of eqn (4), in contrast, is valid only in the limit of moderate to high friction. The fact that the Kramers' rates of Fig. 5 are very close to the results of the LMD simulations (and experimental predictions) suggests that the moderate-to-high friction assumption is appropriate for this problem of polaron transfer (at least in PCBM). A key advantage of our LMD simulation approach is the fact that we are not constrained to simulating systems in either the low- or high-friction limit; we are free to choose any values of friction and meff that are appropriate for the system of interest. Higher friction can be explicitly included, and a relatively high effective mass can incorporate the sluggish nature of reorganization, offering an advantage to Marcus theory. However, we are also free to simulate low friction systems, where simulations do not suffer from the invalidation of eqn (4) that occurs at lower frictions. This provides a major advantage in terms of the flexibility of our model in simulating a variety of systems.
Another key distinction between the Marcus and Kramers' theory models is the question of their appropriateness in modeling an adiabatic transfer process. As noted in Section 3, previous work suggests that systems in which V ≪ kBT should exhibit nonadiabatic charge transfer. V > kBT in our PCBM simulations, consistent with adiabatic charge transfer. Cao and Jung apply a two-state diffusion equation modeling charge transfer to demonstrate that the nonadiabatic limit of charge transfer can be used to derive the Marcus rate expression, while the adiabatic limit can be used to derive the strong-friction limit of the Kramers’ rate.49 Thus, the fact that our simulations provide V > kBT is consistent with the fact that our rates differ greatly from Marcus rates and are much closer to the Kramers' rate of eqn (4).
Here we have presented one-dimensional PESs. Derivation of a three-dimensional PES from a series of one-dimensional intermolecular PESs would be difficult. We suggest an alternative in which transport over a sizable, three-dimensional region is simulated by patching together a series of one-dimensional transport events such as those that we present in this paper. It would be possible and valuable to calculate different PESs for different intermolecular orientations, thus representing the disordered nature of an organic semiconductor. (There would be an added computational cost in performing CDFT calculations for distinct orientations, but not at a prohibitive level.) Because there would be discontinuities in the PES as different curves were patched together, it would not be possible to smoothly continue the LMD simulations in time. Instead, after the completion of a transfer event, the simulation would pause to recalculate velocity based on the change in direction of the pseudo-atom as it proceeds to the next transfer event. Then the process would begin anew with an LMD simulation of transport to the next molecule. Of course, there are multiple molecules adjacent to the molecule on which the charge resides, so there are multiple transfer events that are possible. However, this does not significantly complicate things: the code can simply perform a simulation of transfer to each of the possible acceptors utilizing the same set of random numbers R(t), and the accepting molecule that attains transfer most quickly will then be selected as the new position of the pseudo-atom. We note that consideration of re-crossing is inherent in this approach.
As mentioned above, the small polaron description of a charge is not appropriate for all organic semiconductors, and on a device scale one must carefully consider whether this model is suitable for a given morphology. One challenge faced by theoretical modeling of organic semiconductors is the variation of morphology within a given device: many devices contain some combination of disordered regions and more ordered, crystalline regions. Because our potential energy surface is constructed solely based on the orientation between two molecules, information on overall morphology of a region is not considered in the PES. Intuitively, effective mass and friction are influenced by a region's morphology, and exploration of this relationship is an important direction for future work.
A major advantage of this classical model is computational cost: by reducing a problem of millions of atoms to a calculation treating a small number of pseudo-atoms, the length and time scales of OPV charge transport are accessible. The greatest computational cost lies in the construction of an accurate PES, which requires atomistic CDFT calculations of large molecules. Many OPVs contain molecules larger than PCBM, and given the disordered nature of these systems, calculations would ideally use PESs built from CDFT calculations for many different inter-molecular orientations. Fortunately, recent improvements in linear-scaling DFT have dramatically decreased the cost of treating large molecules. With CDFT now available in linear-scaling DFT software packages,17,78–81 PES development for a few dozen possible orientations is feasible and the disordered nature of these systems can be captured.
The slow rate of charge transport in organic photovoltaics and other semiconductors presents a major challenge in the development of viable technologies. The charge transport model we present here is essentially a multiscale method in which atomistic ab initio calculations of a small bimolecular system serve as a foundation for classical dynamics over large distances. This method is not only computationally affordable but also unique in the information that can be incorporated into the model and the information that is provided by the results. We emphasize that although the low computational cost makes simulation of large regions tractable, we have demonstrated here that simulating even a single transfer event provides substantial insight into the polaron and its transfer process. The use of this model at any scale is expected to provide multiple insights into charge transfer dynamics that may facilitate the development of the next generation of organic electronics.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6sc04547b |
This journal is © The Royal Society of Chemistry 2017 |