PCBM: A

Subscriber access provided by Northern Illinois University

Article

Understanding how Isotopes Affect Charge Transfer in P3HT/PCBM: A Quantum Trajectory-Electronic Structure Study with Nonlinear Quantum Corrections Lei Wang, Jacek Jakowski, Sophya Garashchuk, and Bobby G. Sumpter J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00126 • Publication Date (Web): 09 Aug 2016 Downloaded from http://pubs.acs.org on August 10, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Understanding how Isotopes Affect Charge Transfer in P3HT/PCBM: A Quantum Trajectory-Electronic Structure Study with Nonlinear Quantum Corrections Lei Wang,† Jacek Jakowski,∗,‡ Sophya Garashchuk,∗,† and Bobby G. Sumpter‡ Department of Chemistry and Biochemistry, University of South Carolina, Columbia, SC 29208, and Center for Nanophase Materials Sciences and Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831 E-mail: [email protected]; [email protected]

Abstract The experimentally observed effect of selective deuterium substitution on the open circuit voltage for a blend of poly(3-hexylthiophene)(P3HT) and [6,6]-phenyl-C61 butyric acid methyl ester (PCBM) (Nat. Commun. 5:3180, 2014) is explored using a 221-atom model of a polymer-wrapped PCBM molecule. The protonic and deuteronic wavefunctions for the H/D isotopologues of the hexyl side chains are described within a Quantum Trajectory/Electronic Structure approach where the dynamics is performed with newly developed nonlinear corrections to the quantum forces, necessary to describe the nuclear wavefunctions; the classical forces are generated with a Density Functional Tight Binding method. The resulting protonic and deuteronic time-dependent ∗

To whom correspondence should be addressed Department of Chemistry and Biochemistry, University of South Carolina, Columbia, SC 29208 ‡ Center for Nanophase Materials Sciences and Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831 †

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

wavefunctions are used to assess the effects of isotopic substitution (deuteration) on the energy gaps relevant to the charge transfer for the donor and acceptor electronic states. While the isotope effect on the electronic energy levels is found negligible, the quantum-induced fluctuations of the energy gap between the charge transfer and charge separated states due to nuclear wavefunctions may account for experimental trends by promoting charge transfer in P3HT/PCBM and increasing charge recombination on the donor in the deuterium substituted P3HT/PCBM.

1

Introduction

Over the last two decades conducting polymers have continued to receive considerable attention from researchers. Their potential application include, among others, optoelectronic and spintronic device technology 1,2 as well as solar energy conversion. In general, the intra- and inter-polymer chain interactions are difficult to control and manipulate, in part due to coupling to a large density of molecular vibrational modes. Interestingly, it was recently demonstrated that deuterium substitution can influence electronic properties, the phase separation kinetics and morphology in some polymer blend systems. 3–5 This opens a new opportunity for tuning opto- and spin-electronic properties of these materials by isotopic substitution. A growing number of experimental and theoretical studies points to the significance of the Nuclear Quantum Effects (NQE) in photovoltaic cells 5–12 and to the need for understanding NQEs at the fundamental level. The most common organic photovoltaic cells are based on organic heterojunctions in which difference in the electrochemical potential of a chromophore and of some other material is used to convert light into electric current. The chromophore side of the heterojunction works as an electron donor. It absorbs visible light, which leads to the electronic excitation and formation of an exciton (an electron-hole pair); after electron-hole separation across the heterojunction the positive charge (holes) is carried away from the heterojunction through the donor. Some of the earlier polymeric solar cells used poly(3-hexylthiophene) (P3HT) as 2


Page 2 of 44

Page 3 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the donor material. 13 A typical electron acceptor is a molecule or nanomaterial, consisting of a functionalized fullerene, such as [6,6]-phenyl-C61 -butyric acid methyl ester (PCBM), which can transport electrons from the heterojunction after charge separation. The efficiency of solar cells composed of conducting polymers is fundamentally determined by the difference in electrochemical potentials, charge transfer integrals and electron/hole mobilities. These quantities are defined by the electronic structure such as the character of the frontier molecular orbitals. It is commonly assumed that the difference of the energy levels between the highest occupied molecular orbital (HOMO) of an electron donor and the lowest unoccupied molecular orbital (LUMO) of an electron acceptor in the solar cells determines its output voltage and that it correlates with open circuit voltage of the solar cells. 14,15 But the HOMO and LUMO are not dependent on the isotopes used. Thus the isotope dependence of the open circuit voltage is unexpected. 5 However, the difference in the protonic and deuteronic dynamics (including the nuclear wavefunction and electronic properties), can indirectly affect optoelectronic properties and can be used to explain experimental findings of Ref. 5 Light absorption, charge transfer and transport of the charge carriers (electrons and holes) in organic semiconductors are intrinsically time-dependent quantum-mechanical processes, and may be coupled to nuclear motion. In fact, it has been shown that: (a) isotopic substitution may have great impact on thermal and optical properties of materials, 16,17 and (b) isotopic substitution causes nuclear tunneling in conducting polymers that strongly modifies the optoelectronic properties. 6–8 Also, according to recent theories and modeling of Bittner and Silva 18 and Tamura and Burghardt 19 fluctuations or vibrational excitations may contribute (i) to resonant tunneling which ’couples photo excitations directly to photocurrent producing charge-transfer states’ 18 and (ii) to a ’vibronically hot nature of the charge-transfer state promotes charge dissociation’ 19 beyond the electron-hole Coulomb barrier on time scales of less than 100 fs. In the present work we explore the effects of the side chain deuterium substitution of

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

P3HT within a single molecule model of P3HT/PCBM shown in Fig. 1. Our focus on the side chain alkyl arms of the P3HT polymer is due to prior experimental results showing how they influence the open-circuit voltage. We use a computational approach based on quantum dynamics capable of estimating the NQE in a relatively large molecular system. This approach, quantum trajectory/electronic structure (QTES) dynamics, 20 uses the forces computed on-the-fly based on a Density Functional Tight Binding (DFTB) method. 21–24 The quantum corrections of QTES dynamics are computed from a nuclear wavefunction represented by an ensemble of trajectories. QTES-DFTB dynamics is practical for several hundred atoms 25 and has been used to model adsorption of a single H/D atom on a carbon flake 20 and to simulate substrate motion during a proton transfer in soybean lipoxygenase1. 26 The goals of this paper are the following: (i) to develop a method of quantum dynamics of multiple quantum nuclei, enabling studies of isotope effects in conducting polymers, photovoltaic cells and other organic materials with spintronic and optoelectronic applications in mind, (ii) to expand the QTES methodology to include higher-order quantum corrections for dynamics, 27 necessary to describe the CH bonds without the Zero Point Energy (ZPE) ’leakage’, 28,29 and (iii) to investigate the role of quantum nuclear effects (QNE) and of the dynamics of the hydrogenic/deuteronic wavefunction on the open-circuit voltage in a model for the P3HT/PCBM heterojunction. The paper is organized as follows: The methodology and a test application to CH4 are described in Section 2. Section 3 discusses the simulation of hydrogen/deuterium dynamics in a model P3HT-PCBM heterojunction. We compare properties of pristine, unsubstituted P3HT with those of its isotopologue in which the side chain alkyl arms are fully deuterated (Fig. 1). In this paper we use the term pristine isotopologue to mean, P3HT, while the side-chain fully deuterated form is referred to as SD-P3HT. The computational model consists of P3HT and PCBM molecules for a total of 221 atoms with 39 substituted H/D atoms treated as quantum nuclei. In Section 4 the opto-electronic properties, excitation and

4


Page 4 of 44

Page 5 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


charge separation in P3HT/PCBM heterojunction are analyzed/discussed based on density functional theories for the set of molecular structure snapshots collected from the QTES dynamics. These snapshots of molecular structures from QTES collectively represent a protonic/deuteronic nuclear wavefunction. Section 5 concludes.

Figure 1: The P3HT/PCBM model of 221 atoms: carbon (gray), hydrogen (white), sulfur (yellow) and oxygen (red) spheres, respectively. 39 hydrogens, highlighted in blue, are treated as quantum particles.

2

Dynamics with quantum corrections for selected atoms

The goal of this section is to describe the simulation methodology, used to study the dynamics of multiple quantum atoms (of hydrogen and/or deuterium) in molecules, including conducting polymers and organic heterojunctions. We use the Mixed Quantum/Classical Bohmian (MQCB) dynamics based on quantum trajectories (QT) to selectively account for the nuclear quantum effects (NQE) associated with H/D in P3HT and SD-P3HT. This dynamics approach, introduced in Ref. 25 and interfaced with the Density Functional Tight Binding (DFTB) method of electronic structure (ES) calculations, 30–33 has been used previously to explore NQE in adsorption of H/D on planar graphene ’flakes’, C37 H15 and C87 H23 ,. 20,34 In 5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(a) HOMO

(b) LUMO+3

(c) LUMO

(d) LUMO+2

(e) LUMO+1

Figure 2: Frontier orbitals of P3HT-PCBM from B3LYP/6-31G(d): HOMO (panel a) and LUMO+3 (panel b) are localized on the donor (P3HT). First three lowest unoccupied orbitals (LUMO, LUMO+1, LUMO+2) are localized on the acceptor (PCBM). Panels (c) to (e) show respectively LUMO, LUMO+1 and LUMO+2. Exciton dissociation involves charge transfer between primarily LUMO+3 (donor) and LUMO+2 (acceptor)”

6


Page 6 of 44

Page 7 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the QT approach, quantum corrections to dynamics of selected nuclei are introduced through (a) an ensemble of trajectories representing the quantum wavefunction evolving in time and (b) a quantum potential that, in addition to classical potential, governs motion of quantum nuclei. The quantum potential is a purely quantum mechanical object with no classical analogue and, in principle, it is responsible for tunneling, interference and spreading of the wavefunction of a freely moving particle. Mathematically the quantum potential is obtained from a corresponding wavefunction’s amplitude. We start this section by briefly outlining the MQCB framework. More details can be found elsewhere. 25 Next we describe evaluation of the quantum potential and the approximations used. In Ref. 20 we used the linear quantum force approximation. Here we go beyond the linear quantum force approximation and introduce quadratic fitting terms into the quantum potential expression in a practical manner. Benchmarking, validation and discussion of this approach for a test case of a methane molecule is presented at the end of this section.

2.1

Mixed Quantum/Classical dynamics

The MQCB approach is defined in the Madelung-de Broglie-Bohm 35,36 or QT formulation of the time-dependent Schrödinger equation (TDSE). Within the QT formulation it is implemented using the Approximate Quantum Potential approach. 37 However, for a few quantum particles the MQCB dynamics can be implemented with any conventional quantum or semiclassical wavefunction propagation technique. In the QT formulation the quantum effects arise from the so-called quantum potential in the Hamilton-Jacobi equation. This framework allows including NQE into the degrees of freedom (DOFs) describing light particles and to neglect the quantum potential in the DOFs describing heavy quasi-classical nuclei comprising a molecular system. The superscripts q and c will be used to denote the quantum and quasi-classical/classical subsystems henceforth. Consider a system of N p particles described in three-dimensional Cartesian space by 3N p coordinates. Out of N p particles, N q are expected to exhibit NQE and N c = N p − N q are expected to behave nearly classically. 7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 44

A column-vector ~x of dimensionality 3N q will represent the quantum DOFs (QDOFs), and a column-vector ~y of dimensionality 3N c , will represent the classical DOFs (CDOFs). In MQCB a wavefunction is factorized into slow and fast components describing CDOFs and QDOFs, respectively, ψ(~x, ~y , t) = ψ q (~x, t; ~y ) ∗ ψ c (~y , t),

(1)

and each component is expressed in the polar form: ı q ψ (~x, t; ~y ) = A (~x, t; ~y ) exp S (~x, t; ~y ) h ¯ ı c S (~y , t) . ψ c (~y , t) = Ac (~y , t) exp h ¯ q

q

(2)

The fast component ψ q is assumed to depend on ~y parametrically. To complete the matrix/vector notations, the inverse particle masses are arranged into the diagonal matrices Mq and Mc for the quantum and classical subspaces: q q M11 = (m1 )−1 , . . . , M3N q

3N q

c M11 = (mN q +1 )−1 , . . . , M3N c

= (mN q )−1 3N c

= (mN p )−1 .

(3) (4)

Substitution of Eqs (1) and (2) into the usual TDSE

ı¯h

∂ ˆ x, ~y , t), ψ(~x, ~y , t) = Hψ(~ ∂t

(5)

and identification of the trajectory momentum with the gradient of the phase,

p~x = ∇x S q ,

p~y = ∇y S c

(6)

leads to the following quantum/classical Hamilton-Jacobi equation: 1 1 ∂S q ∂S c + = − (~px )T · Mq · p~x − (~py )T · Mc · p~y − V (~x, ~y ) − U (~x, ~y , t). ∂t ∂t 2 2 8


(7)

Page 9 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The non-local quantity U is the quantum potential,

U (~x, ~y , t) = −

h ¯2 (∇T · Mq · ∇x + ∇Ty · Mc · ∇y )|ψ(~x, ~y , t)|. 2|ψ(~x, ~y , t)| x

(8)

Neglecting the gradients of Aq and S q and Ac with respect to ~y and averaging the classical force acting on CDOFs, ~y , over QDOFs, ~x, one obtains the time-evolution equations in the Lagrangian frame-of-reference for the quantum and classical subspaces from Eq. (7): d~yt = Mc · p~yt , dt d~xt = Mq · p~xt , dt

d~pyt = − h∇y V iq |~y=~yt dt d~pxt = − ∇x (V (~x; ~yt ) + U q (~x, t; ~yt ))|~x=~xt . dt

(9) (10)

Subscript t is used to indicate quantities associated with the trajectories, such as p~xt as opposed to functions of coordinates, such as p~(x, y). In Eq. (9) specifying the dynamics of a trajectory (~yt , p~yt ) in the CDOFs, the averaging of the classical force is performed over the QDOFs for a ’sub-wavefunction’ ψ q (~x, t; ~yt ), guided by that trajectory. Eq. (10) describes evolution of ψ q (~x, t; ~yt ) in the QT framework with the quantum effects arising from the quantum potential U q approximating the full potential U of Eq. (7),

U ≈ U q (~x, t; ~y ) ≡ −

h ¯2 ∇T · Mq · ∇x Aq (~x, t; ~y ). 2Aq (~x, t; ~y ) x

(11)

In the QT formulation the wavefunction phase along the QT ~xt in the quantum subspace, representing a sub-wavefunction ψq (~x, t; ~yt ) which is guided by the mean-field trajectory ~yt in the CDOFs, is defined as 1 1 d q (S + S c ) = (~pxt )T · Mq · p~xt + (~pyt )T · Mc · p~yt − V (~xt , ~yt ) − U q (~xt , ~yt ). dt 2 2

(12)

The wavefunction amplitude evolves according to the continuity equation for the probability density. However, when computing expectation values, we simply use the following conserva-

9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 44

tion property: the probability density within the volume element associated each trajectory, referred to as trajectory weights 37 is constant in time. The trajectory weight wkl – of the k th trajectory in QDOFs guided by the lth trajectory in CDOFs – are defined by the initial wavefunction, (k)

(l)

wkl = |ψ(~x , ~y , 0)|

ˆ q hΩi

2.2

(l) y=yt

=

X

2

q 3N Y

(k) dxi

(l)

dyj ,

(13)

j=1

i=1

(k)

c 3N Y

(l)

Ω(~xt , ~yt )wkl .

(14)

k

The non-linear approximation to the quantum force

The quantum potential, in general, is a complicated singular function. For a practical highdimensional implementation the quantum potential has to be computed approximately. The simplest, physically meaningful approximation is the Linearized Quantum Force (LQF). 37 LQF is a globally-defined mean-field approximation, which is exact for correlated timedependent Gaussian wavefunctions and is capable of describing the dominant quantum effects, such as wavefunction bifurcation, moderate tunneling, ZPE and the energy of wavefunction localization on a time-scale dependent on the anharmonicity of the potential. The LQF parameters are defined variationally which ensures conservation of the total energy in a system. The LQF idea is extended here to allow a non-linear approximation of the quantum forces, which we find to be necessary (discussed below) for the study of P3HT. We construct the approximate quantum potential by fitting each component of the socalled nonclassical momentum defined for the QDOFs,

~r =

∇ x Aq . Aq

(15)

The quantum potential of Eq. (11) in terms of ~r becomes:

Uq = −

1 T ~r · Mq · ~r + ∇Tx · Mq · ~r . 2

10


(16)

Page 11 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The fitting of ~r is done in a basis f~ of Nb unspecified polynomials of the QDOF coordinates, f~ = [x1 , x2 , . . . , 1], and the expansion coefficients of the nth component is a vector ~cn of the size Nb , n = [1, 2, . . . , 3N q ]. The least squares fit of ~r is equivalent to the minimization of hU q i with respect to ~cn . Arranging the expansions coefficients as a matrix C, the optimal coefficients solve the linear matrix equation,

SC = B, where C = [~c1 , ~c2 , . . . , ~c3N q ].

(17)

In Eq. (17) the overlap matrix S and the right-hand-side matrix B are 1 S = hf~T ⊗ f~iq , B = − h∇T ⊗ f~iq . 2

(18)

The expectation values in Eq. (18) are computed using the trajectory weights according to Eq. (14). In the LQF method the basis is linear, f~ = [~x, 1] with a basis size of Nb = 3N q + 1, which makes it efficient for large systems. If all the quadratic terms in xi are included in a basis, then the basis size will scale with the system size as (N q )2 , which is expensive, Nbquad = 2 × (3N q ) + 3N q × (3N q − 1)/2 + 1. However, not all the DOFs are correlated and we find that for P3HT, a basis with quadratic terms for each hydrogen – without the nonlinear terms involving different atoms – reproduces the zero point energy (ZPE) over multiple oscillation periods typical for the CH bond. Thus, we use a basis of 10 functions for each quantum nucleus to generate their respective quantum potentials according to Eq. (16), Unq , n = [1, N q ]. For the first nucleus, n = 1, the basis is f~1T = [x1 , x2 , x3 , x21 , x22 , x23 , x1 x2 , x1 x3 , x2 x3 , 1].

(19)

Three vectors of the expansion coefficients {~c1 , ~c2 , ~c3 } (the vector length is Nb = 10) are determined by Eqs (17) and (18) describing the three components of the nonclassical momentum

11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 44

associated with the first nucleus. Using the fitting functions, ~ri = ~ci · f~i for i = {1, 2, 3} in Eq. (16), one obtains an approximate analytical U1q . The remaining quantum particles are treated similarly, one at a time. The total quantum potential is a simple sum over the quantum particles, q

q

U (x1 , x2 , . . . , x3Nq ) =

N X

Unq (xna+1 , xna+2 , xna+3 ),

n=1

na = 3(n − 1)

(20)

and the quantum force is determined from its analytical gradient. Due to the limited form of the basis defined by Eq. (19), the exchange of energy between different quantum and classical particles is mediated only by the classical potential V . Such a description, however, is reasonable for P3HT given that the quantum particles are not bonded to each other and that there is no bond-breaking in the system. To relax this limitation, the first step is to include the nearest-neighbor linear basis functions into the basis for each quantum particle or do a two-step – global linear/local non-linear – fitting procedure as described in Ref. 27 and employed to describe the ZPE of solid He.

2.3

Quadratic basis quantum force for -CH bond in methane

In this section the accuracy of the quadratic basis approximation to the quantum force described in Section 2.2 has been numerically tested on a CH bond of a methane molecule, CH4 . Specifically, the expectation energies associated with the localization of a selected (single) hydrogen and its motion in CH4 have been compared with the ZPE modes. For the current analysis purposes, the molecule is taken at the equilibrium geometry and is oriented in the three-dimensional Cartesian space so that one of the CH bonds is aligned with the z-axis. This ’top’ hydrogen (or proton) is treated as the only quantum nucleus described by three coordinates labeled here for clarity as (x, y, z). The remaining atoms are frozen at the equilibrium geometry. The mass of the top hydrogen is m = 1836 a.u. To freeze the positions of the remaining atoms in the course of dynamics, their masses are simply set to infinity,

12


Page 13 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


which effectively decouples motion of the selected hydrogen from remaining hydrogens. The motion of hydrogen along z-axis corresponds to -CH stretching while the motion in the xy plane corresponds to CH bending. The initial distribution of quantum trajectories for this quantum hydrogen at time t = 0 has been set according to a normalized three-dimensional Gaussian wavefunction, ψ(x, y, z, 0) = N exp(−αx (x − x0 )2 − αy (y − y0 )2 − αz (z − z0 )2 ),

(21)

centered at the equilibirum geometry (x0 , y0 , z0 ) = (0, 0, 2.09) a0 . Thus the resulting initial distribution of the hydrogen position corresponds to an oblate spheroid symmetry of wavefunction of the decoupled hydrogen with a width parameter equal to αx = αy = 4.72 −2 a−2 0 and αz = 13.00 a0 . The values of parameters αx , αy and αz are determined from the

second derivative of the potential energy curves with respect to the displacement of H from its equilibrium position along the x, y and z directions. To further validate accuracy of the quadratic basis in description of the CH ground state, we have performed dynamics of the proton in the Morse potential describing the CH stretch in methane along the bond axis while freezing all the other atoms. The initial wavefunction is a Gaussian close to (but not quite equal to) the vibrational ground state as described above for the z-coordinate. The expectation value of the quantum potential, averaged over 5000 a.u. of the exact QM propagation, is hU iQM = 0.00362 ± 0.00014 Eh . The same quantity obtained within the QQF approximation is close to the exact value, hU iQQF = 0.00352 ± 0.00011 Eh . The original LQF, which includes all the trajectories including those which move to the dissociation region due to imbalance between the classical force and the LQF-approximated quantum force, becomes zero after 1000 a.u. If the LQF is computed using just the trajectories in the well region (displaced by no more than 1 a0 from the equilibrium bond distance) then the average quantum potential decays to approximately half of its original value by t = 500 a.u. and stabilized around this value; the expectation value

13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 14 of 44

Page 15 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


trajectories to represent ψ(x, y, z, 0). The potential V for the hydrogen motion is computed on the fly within the QTES-DFTB approach. To analyze the effect of the non-linear quantum force correction, the approximate quantum potential, computed within the linear and quadratic bases, is included into one QDOF at a time, while keeping the remaining DOFs frozen. The components of the expectation energy for the hydrogenic wavefunction are shown in Fig. 4. The results of the quantum correction in x (corresponds to CH bending) and z (corresponds to CH stretching) are displayed in panels (a) and (b) respectively. On each panel, the expectation values of the total energy E, the potential energy hV i and the quantum potential hU i calculated as average over all 3840 trajectories are plotted as a function of time. For an eigenstate, both hV i and hU i would be constant in time and hU + V i = E. In the QT formulation, the trajectories describing an eigenstate are stationary and their ˆ = h~p2 i/(2m) = 0. Our initial wavefunction is close expectation kinetic energy is zero, hKi to, but not an exact ground state of the system, thus, small oscillations of hV i and hU i and a decrease in hU i at short times are seen. In the x and y directions the classical potential V is well approximated by parabolic functions, therefore, the linear quantum force is fairly accurate during the time-evolution and the nonlinear correction has negligible effect. The expectation value of the quantum potential remains close to 0.001 Eh , which shows that the wavefunction remains localized and close to a Gaussian shape. We note that the total energy of the system is rigorously conserved. Once we consider dynamics along the z-axis (CH stretching), the potential V shows appreciable anharmonicity, and the linear quantum force is no longer accurate. Here the linear quantum force does not compensate for the classical force correctly, which means that the fringe QTs experience an unphysically large total force sending them to an unbound region of the potential. To circumvent this issue we artificially stop such QTs if their displacement distance from the average position exceeds a given cutoff ∆ (here ∆ was set to 1.0 a0 ) which effectively removes some small amount of kinetic energy from the system. This is the reason for the loss of total energy in the LQF

15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 44

calculation in panel (b). In the same calculation hU i drops with time because some of it was incorrectly transferred to the kinetic energy of the QTs. The nonlinear correction of the quadratic basis quantum force (QQF) – the quadratic basis is used when approximating the quantum force – mitigates this effect: the quantum energy remains constant around 0.003 Eh and the wavefunction stays localized as expected. The total energy is conserved as there are no QTs with ’frozen’ positions reaching the cutoff value. Overall, the use of the cutoff on the trajectory positions does not qualitatively change the dynamics of the system, because there is no bond-breaking at this energy in CH4 . One can expect qualitative differences in cases for CH bond breaking QT dynamics with quadratic vs linear basis used to define the quantum potential due to the significant changes in wavefunction distribution in time. In this case the linear or even quadratic bases may not be sufficiently flexible to describe splitting and branching of the wavefunction. This is, however, beyond the scope of our current studies. To summarize, the LQF approximation is sufficient for describing the motion that corresponds to bending modes of methane. This result is due to the fact that the harmonic approximation sufficiently represents the corresponding potential energy and the LQF is exact for such cases. For the CH bond stretching there is significant anharmonicity of the potential energy profile. In this case a QQF works well while a LQF is insufficient and causes fringe trajectories. We did not observe fringe trajectories when a QQF approximation was used. Overall, the nonlinear quantum force generated with the quadratic basis is sufficiently accurate to describe the wavefunction localization and associated ZPE of the CH bond. Therefore, we will use a QQF to treat the quantum protons in P3HT with the CH interactions similar to those in CH4 . Table 1: The average classical (V ) and quantum (U ) energies in the CH bond of CH4 obtained from the QTES-DFTB dynamics with a quadratic basis quantum force. The results (except the initial values) are averaged over the time interval t = [2000, 10000] a.u. The deviation of the total energy ∆E, averaged over the entire time-interval is listed in the last column. atom H D

V0 0.217 0.160

Vt 0.258 0.183

V change 18.9% 14.4%

U0 0.238 0.162 16

U 0.154 0.110


U change 36% 32%

∆E 2.6% 1.9%

(a)

0.008

0.004

E LQF LQF LQF

E QQF QQF QQF

0

Energy [Eh]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Energy [Eh]

Page 17 of 44

0.008

(b)

0.004

0 0

1000

2000

3000

4000

5000

Time [a.u.] Figure 4: Energy of the proton in CH4 as a function of time for a wavefunction evolved for (a) motion of H on the xy -plane (this corresponds to bending of CH bond) and for (b) motion of H along z-coordinate (which corresponds to stretching of the CH bond). The average total energy E, potential energy hV i and the quantum potential hU i obtained with the quadratic basis quantum force (QQF on the legend) are shown as a solid line, dash and dot-dash lines, respectively. The corresponding results obtained with the linearized quantum force (LQF on the legend) are marked with circles, squares and triangles, respectively.

17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

3

The QTES dynamics of P3HT-PCBM

The computational model of the P3HT/PCBM system consists of five thiophene monomer rings with each thiophene monomer functionalized with a hexyl group, referenced here as a ’side-chain arm’. Three out of five hexyl arms in our P3HT/PCBM model are wrapped around the fullerene of PCBM (see Fig. 1). The molecular geometry for P3HT/PCBM was obtained from DFTB optimization that included dispersion corrections. The side-chain H/D-atoms in these three arms nearest to the fullerene are treated as quantum protons or deuterons, respectively. The remaining atoms are described as classical particles which, in general, move according to the MQCB formalism. The simulations were first performed with all atoms allowed to move. We found, however, that for our time of interest, which is shorter than 500 fs, the motion of heavy atoms was negligible. Therefore, the dynamics discussed below pertains to all heavy atoms being frozen in space and time. With all heavy atoms in the P3HT/PCBM model being fixed during QTES-DFTB simulations, it is easy to distinguish the effects of the motion of quantum H/D nuclei on the optoelectronic properties from the effects of the remaining nuclei motion. The wavefunction for the 39 quantum protons (and deuterons) of the three arms closest to the fullerene (see Fig. 1) is initialized as a direct product of three-dimensional Gaussian wavefunctions. Each Gaussian describing H or D atoms is in the form of Eq. (21) and is centered at the H position from a DFTB+dispersion optimization. The local z-axis of each Gaussian is aligned to the equilibrium position of the CH bond. The values of α parameters in Eq. (21) describing the initial nuclear wavefunction of H/D in 3 hexyl arms are set to αx = αy = 4.72 and αz = 13.00 a−2 for H and to 0 αx = αy = 6.43 and αz = 17.74 a−2 0 for D. These values correspond, respectively, to a typical single-bond CH in methane CH4 (see Sec. 2.3) and in deuterated methane CD4 . Each threedimensional wavefunction is sampled by the initial QT positions in the respective local set of coordinates, defined by the equilibirum position of the H-atoms, according to the normal distribution of |ψ|2 , given by Eq. (21). The QT positions are transformed back into the global Cartesian coordinates to generate a representation of the full-dimensional wavefunction (the 18


Page 18 of 44

Page 19 of 44

number of QDOFs is 3 × 39 H/D) in terms of the QT ensemble. The initial wavefunction is real and the initial momenta of the QTs are zeroes. In order to compare the average geometries of the P3HT and its H/D substituted version SD-P3HT, which could lead to some differences in the electronic energy level structure of the two species, the dynamics simulation is performed for 10000-20000 a.u. (0.25-0.5 ps) until the average quantum potential exhibits stable behavior. A time-step of dt = 10 a.u. was used and the number of trajectories was Ntraj = 3840. Fig. 5 shows the expectation value of the quantum potential as a function of time for the P3HT and SD-P3HT in the top and bottom panels, respectively. A stable regime has been achieved within 1000 a.u. (50 fs). Fig. 6 shows the time dependence of H/D postion from a single, selected QTES trajectory.

[Eh]

0.25

Hydrogen

0.2 0.15 0.1 0.25

[Eh]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.2

linear diagonal quadratic

Deuterium

0.15 0.1 0

2000

4000

6000

8000

10000

Time [a.u.] Figure 5: The expectation value of the quantum potential hU i as a function of time for the P3HT/PCBM system obtained with three different fitting bases: linear, with just the diagonal terms, and quadratic with all terms. The upper and lower panels show the results for the H/D isotopologues. To maintain the accuracy of the approximate quantum potential over time, freezing of the fringe QTs as described in Section 2.3 was modified in the following way: if the distance 19



Displacement [Å]

0.8

(a)

middle chain H methyl group H middle chain D methyl group D

0.6

0.4

0.2

0 0

0.8

Displacement [Å]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 44

2000

4000

6000

8000

10000

2000

4000

6000

8000

10000

Time [a.u.]

(b)

0.6

0.4

0.2

0 0

Time [a.u.]

Figure 6: Distance of H to its initial position for (a) the ensemble-average and (b) for the single-trajectory geometries. The results for the middle-of-chain and methyl-group hydrogen atoms are shown with filled circles and squares respectively. The same results for the deuterium atoms are shown with open circles and squares, as stated in the legend.

20


Page 21 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of a QT to the average position of a certain nucleus deviates by more than ∆ = 0.5 ˚ A, then the corresponding momentum of that QT is set to the average value of the QT momenta for that nucleus. To study the effect of the non-linear quantum corrections, dynamics was performed using three basis sets: (i) the linear basis generating Linear Quantum Force f~l = {x, y, z, 1}; (ii) the ’diagonal’ quadratic basis with functions f~d = {x, y, z, x2 , y 2 , z 2 , 1}; (iii) the full quadratic basis f~q = {x, y, z, x2 , y 2 , z 2 , xy, xy, yz, 1}. The last basis yields the largest expectation value of the quantum potential at long t > 2000 a.u. times and these are the results we discuss in more detail. In the course of dynamics, the average quantum potential, hU i, drops by about 35% compared to its initial value as the initial wavefunction delocalizes, and the effect is similar for both isotopologues. The delocalization is accompanied by the increase in the average potential energy by 18.9% and 14.4% for P3HT and SDP3HT, respectively. The total energy in the system is not rigorously conserved because of the momentum change introduced for the fringe QTs, but the deviation of the total energy for the two systems, 2.6% and 1.9% is significantly smaller than the changes in the quantum potential energies. The average quantum potential, averaged over time interval √ t = [2000, 8000], is hU iH = 0.15 and hU iD = 0.11 Eh and their ratio of 1.36 is close to 2 expected from the Harmonic oscillator model. To rationalize the difference in the configuration of the two isotopologues averaged over their respective wavefunctions we examined the average CH bond-distance for the three ’arms’ with quantum H/D shown in Fig. 7. We observe distinguishable behavior in the dynamics of various hydrogens (or deuteriums) depending on their location within the hexyl groups of P3HT. Specifically, there is noticeable difference in the dynamics of three ’terminal’ hydrogens (deuteriums) located in the terminal methyl group (end of chain) of the each hexyl arm from the ten ’internal’ hydrogens located within -CH2 - groups of hexyl arm. The top panel in Fig. 7 shows the bond-distance for the -CH2 - internal H/D (located on -CH2 -), while the bottom panels shows the bond-distance for the terminal H/D (located on -CH3 ends). For the internal H/D the oscillation pattern is clearly visible with the period ratio of

21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

D to H close to

√

2. The average bond distances are 1.127 and 1.120 ˚ A, respectively, which

is consistent with the wavefunctions of the H atoms reaching further into the anharmonic region of the potential describing the stretched CH/CD bond. The difference in the average CH/CD bond distance is 0.6% and its variation is about twice as large for D compared to that for the H isotopologues. For the terminal H/D of the methyl groups (end of chains), the oscillations are less pronounced, and at the end of the propagation time, the stable regime of dynamics had not quite been reached. The CH/CD distances for the terminal H/D are somewhat shorter, 1.072 and 1.080 ˚ A, respectively, and their variation is in the same range for CH and CD. One possible explanation for these trends is stronger interaction of the methyl group H/D with C60 compared to that of the internal H/D in -CH2 - groups. However, given the limited accuracy of the dynamics and of the DFTB energies, the results should be interpreted with caution. We assume that given that the entire charge transfer process is ultrafast, i. e., on the order of 100 fs 38 – the difference in the isotopologue properties is traceable to the difference of the ground state nuclear wavefunctions, and thus we used more sophisticated ab initio methods to examine the electronic structure of P3HT/PCBM. Table 2: Average hydrogen-carbon bond distances obtained from the QTES-DFTB dynamics using a non-linear quantum correction in P3HT/PCBM with 39 quantum H/D. The terminal H/D are located at the methyl ends of the hexyl arms. The internal H/D are located within (CH2 )5 - section of the hexyl groups. The results are averaged over the time interval indicated in the last column. Quantum atom location internal H internal D terminal H terminal D

Bond distance [a0 ]

Period [a.u.]

Time interval [a.u.]

1.127 ± 0.0020 1.120 ± 0.0047 1.072 ± 0.0020 1.080 ± 0.0016

520 718 513 691

[500, 1000] [500, 1000] [800, 1000] [800, 1000]

22


Page 22 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Bond distance [Å] Bond distance [Å]

Page 23 of 44

middle chain H/D 1.15

1.1

methyl group H/D

H D

1.1

1.05

0

2000

4000

6000

Time [a.u.]

8000

10000

Figure 7: Average hydrogen-carbon bond length for the hexyl arms in P3HT/PCBM in the course of the QTES-DFTB dynamics as a function of time. Bottom panel: terminal H/D, located within methyl ends of hexyl groups. Top panel: internal H/D, located within -(CH2 )5 - of hexyl groups. The non-linear quantum force correction was used.

23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

4

Page 24 of 44

NQE on the electronic structure and charge transfer

The geometries obtained from the QTES-DFTB dynamics are used to analyze the energy level structure of P3HT/PCBM and SD-P3HT/PCBM using Density Functional Theory (DFT) 39,40 and time-dependent DFT (TD-DFT). 41 In general, for a given number of excitons (electron-hole pairs) generated by absorption of light in P3HT the power output of the P3HT/PCBM solar cell depends on three factors: (1) the rate constants for the tunneling of charge cariers (electrons and holes) across the heterojunction that leads to interfacial charge transfer states and similar rates for recombination, (2) the energy gap Eg between interfacial charge transfer state and the ground state, and (3) the mobilities of the charge carriers that determine the diffusion rates of dissociated holes and electrons away from the heterojunction to their respective electrodes. The charge mobilities (described through series resistance) mostly affect the output current from the solar cell and can be neglected in an open circuit situation (infinitely large load resistance and zero output current). It is mainly the first two factors that contribute to open circuit voltage (Voc ), namely, the energy gap and charge tunneling rates. The experimentally measured open-circuit voltage VOC of the protonic and deuteronic samples is 0.58 V and 0.34 V respectively, while the short-circuit current JSC remained essentially unchanged upon the SD H/D-substitution. Based on experimental data for a variety of polymer-PCBM blends, the following empirical relation of VOC to the energies of the frontier orbitals has been suggested: 5,42

VOC = Eg /e − ∆Vloss ,

A D Eg = ELU M O − EHOM O .

(22)

Here the quasiparticle energy gap Eg of P3HT/PCBM is approximated by Kohn-Sham gap between energy of HOMO localized on electron donor (P3HT on the Fig. 2(a)) and energy of LUMO localized on electron acceptor (PCBM on the (Fig. 2(c)). The HOMO-LUMO gap provides the simplest approximation of Eg . Other ways to estimate Eg are also possible, for 24


Page 25 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


example constrained DFT or TD-DFT can be used to directly estimate Eg as a difference between charge transfer and ground states. We stress that there are several interfacial charge transfer (CT) states that can be found from TD-DFT which are energetically below the first optically active valence-excited neutral π − π ∗ excitation of the P3HT/PCBM as can be seen in Figure 8. The lowest interfacial CT state from a B3LYP/6-31G(d) calculation is at ∼1.2 eV whereas the optical band gap (the lowest neutral excited state with large oscilator strength) is at ∼2.15 eV. However, charge transfer excited states are often underestimated in TD-DFT. In contrast the energies of neutral excited states are typically well described in TDDFT/B3LYP. To verify that, we estimated the lowest interfacial CT state and optical band gap from long range corrected CAM-B3LYP (LC-ωPBEh) functionals. The lowest interfacial CT state from CAM-B3LYP (LC-wPBEh) is 2.23 eV (2.35 eV) which is about 1eV above the B3LYP value. The optical band gap which is a neutral state from CAM-B3LYP (LCwPBEh) is 2.5 eV (2.6eV), ∼ 0.4eV above the B3LYP value. The better agrement between B3LYP and the long range corrected functionals for the optical band gap is expected because it is a neutral state. The ∆Vloss in the Eq. (22) is related to the specific details of charge transfer and other recombination mechanisms. But the empirical Eq. (22) was obtained for polymer/PCBM blends and should be corrected for the non-ideality of the donor-acceptor junction and temperature. 9 Note, that the exciton on the electron donor (HOMO and LUMO on the donor) are localized on the main chain of P3HT and are spatially separated from the H/D substitution of the side-chains at the same geometry. Therefore, the corresponding HOMO and LUMO energies (π and π ∗ orbitals in thiophene ring) are essentially independent of the H/D substitution of the side chains, even taking the ZPE of nuclear vibrations into account. Our model of a heterojunction consists of just a single PCBM and a short P3HT polymer, thus a qualitative, rather than quantitative description of NQE on VOC is all that can be expected from the electronic structure analysis. Also, to distill the effect resulting purely from quantum nuclei, the molecular structures used were obtained from QTES-DFTB with

25



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 26 of 44

Page 27 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


all heavy nuclei frozen. The sets of structures representing the evolving nuclear wavefunction that were used for electronic structure analysis differ only in the position of quantum nuclei. This approximation seems to be justified/supported by experimental measurements. 5 Atomic force microscope phase images and energy filtered transmission electron microscope images show very similar crystallinity between pristine P3HT and SD-P3HT. Both, pristine P3HT and SD-P3HT reveal complex interpenetrating domains of similar structure and size (see Fig. 3 in Ref 5 ), which is clearly different from the main thiophene chain deuterated P3HT (denoted as MD-P3HT in Ref 5 ). Table 3: The energies of the HOMO and LUMO and the energy gap obtained with DFT for the P3HT/PCBM and SD-P3HT/PCBM. The calculations are performed for the geometries averaged over the nuclear wavefunction (represented by an ensemble of 4320 quantum trajectories) and for a single quantum trajectory. Time [a.u.] B3LYP 0 10000

Ensemble average P3HT EHOM O [Eh ] ELU M O [Eh ] Eg [eV]

Ensemble average SD-P3HT EHOM O [Eh ] ELU M O [Eh ] Eg [eV]

-0.172 -0.113 1.605 -0.171 -0.113 1.578 Single trajectory P3HT

-0.172 -0.113 1.605 -0.171 -0.113 1.578 Single trajectory SD-P3HT

-0.171 -0.170 -0.171 -0.171 -0.170

-0.114 -0.114 -0.114 -0.114 -0.113

1.551 1.524 1.551 1.551 1.524

-0.171 -0.171 -0.171 -0.171 -0.171

-0.114 -0.114 -0.113 -0.114 -0.114

1.551 1.551 1.524 1.551 1.551

-0.215 -0.214 -0.214 -0.214 -0.214

-0.080 -0.080 -0.081 -0.080 -0.080

3.663 3.629 3.632 3.643 3.641

-0.215 -0.214 -0.214 -0.215 -0.214

-0.080 -0.080 -0.080 -0.080 -0.080

3.666 3.654 3.654 3.661 3.635

-0.238 -0.237 -0.238 -0.237 -0.237

-0.067 -0.068 -0.068 -0.067 -0.067

4.638 4.602 4.604 4.615 4.615

-0.238 -0.237 -0.237 -0.238 -0.237

-0.067 -0.067 -0.067 -0.067 -0.068

4.641 4.628 4.628 4.632 4.608

B3LYP 0 7000 8000 9000 10000 CAM-B3LYP 0 7000 8000 9000 10000 LC-ωPBEh 0 7000 8000 9000 10000

Quasiparticle energy gap. We examined several models of the quasiparticle energy 27



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

gap Eg for three different functionals (B3LYP, CAM-B3LYP, LC-ωPBEh) while utilizing the information about the quantum nuclei wavefunction from quantum trajectory simulations. In the first model Eg was calculated as difference between DFT energies of the HOMO of the donor and the LUMO of the acceptor for the averaged geometry over the QT ensemble of 4320 trajectories. The DFT calculations are performed every 1000 a.u. using Q-Chem4.2. 43 Eg computed with B3LYP/6-31G(d) at the average geometry for the two isotopologues is essentially independent of H/D substitution, as seen from Table 3, and is equal to 1.58 eV. In the second model, Eg was calculated along a single trajectory as a difference of DFT energies (B3LYP, CAM-B3LYP, LC-ωPBEh) of the HOMO of donor and the LUMO of acceptor and averaged over various structures along a single trajector. Even when taken along a single trajectory, the variation of Eg for the two species does not exceed 0.03 eV which is below 2% of Eg . In the third model Eg was calculated from TD-DFT (at the B3LYP, CAM-B3LYP, LC-ωPBEh theory and with 6-31G(d) basis set) along a single trajectory and averaged over structures along a single trajectory. Specifically, Eg was calculated as a difference between S1 and S0 singlet states of the P3HT/PCBM and then averaged over time. S0 is the ground state (neutral) and the S1 is the state that corresponds to the interfacial charge transfer state with the leading term involving a transition between the HOMO of donor and the LUMO of acceptor. The transition between S0 and S1 states is optically inactive (very low oscillator strength) which is a consequence of large spatial separation. Note, PCBM has three unoccupied orbitals (LUMO, LUMO+1 and LUMO+2) below the LUMO of P3HT. These three LUMOs on PCBM are good electron acceptors and appear as the three lowest interfacial CT states below the optical band gap (∼2.15 eV at the B3LYP/631G(d) level of theory in Figure 8). The TD-DFT model of Eg gives similar results to first two models based on the HOMO-LUMO energy gap: the H/D difference in Eg is on the order of 0.01 eV. Specifically, Eg at B3LYP theory is 1.16 eV and 1.17 for hydrogenated and deuterated species. Similarly, Eg at CAM-B3LYP (LC-ωPBEh) theory is 2.22 eV (2.346 eV) and 2.23 eV (2.353 eV) for hydrogenated and deuterated P3HT. Thus, the variation in Eg

28


Page 28 of 44

Page 29 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


does not explain the experimental observation for VOC lowering by 0.24 eV upon side-chain deuteration. This conclusion is consistent with the character of the HOMO orbital, localized on the thiophene main chain shown in Fig. 2(a), while the difference in the ground state nuclear wavefunction is limited to the hexyl side chains. The fourth model of Eg was similar to TD-DFT but the energy of the charge transfer state was estimated from the open-shell singlet constrained DFT calculations (CDFT) along a single trajectory. Our final approach based on constrained density functional theory led to a similar observation that the difference in VOC cannot be explained based on the difference between energies of the charge transferred P3HT(+)/PCBM(-) state and the neutral ground P3HT/PCBM state. 44 Table 4: Average values of charge transfer integrals W and energy difference ∆E for the neutral excited state (dominated by LUMO+3) and closest charge transfer state (dominated by LUMO+2) from various DFT fucntionals with 6-31G(d) basis set. The averaging was performed over QTES-DFTB structures for hydrogenated and deuterated P3HT. ∆ET D−DF T corresponds to the energy difference betnween excitaed states calculated from linear response TD-DFT. ∆EM O corresponds to energy difference calculated between adiabatic states (LUMO+3 and LUMO+2). ∆EDiabats corresponds to energy difference between diabatic states from four state model. All energies are in meV. Hydrogen Energy [meV] RMS [meV]

Deuterium Energy [meV] RMS [meV]

B3LYP ∆ET D−DF T ∆EM O ∆EDiabats W

163.1 769.7 841.7 72.2

67.8 7.6 12.5 0.19

137.6 759.4 824.9 73.7

31.9 6.5 9.6 0.20

∆ET D−DF T ∆EM O ∆EDiabats W

74.9 836.2 879.0 74.8

5.7 10.0 12.5 0.21

78.2 822.1 858.7 78.2

4.4 8.3 9.8 0.22

∆EM O ∆EDiabats W

893.3 879.0 33.8

8.3 49.3 0.05

881.3 890.7 56.2

7.6 38.3 0.12

CAM-B3LYP

LC-ωPBEh

Charge tunneling rates. After elimination of the effect of nuclear wavefuction on Eg we move our attention to the ∆Vloss term in Eq. (22). This term descibes rates for 29



interfacial tunneling between neutral excited states of P3HT/PCBM and interfacial charge transfer states (as shown in Figure 8) as well as recombination rates. We explored further the NQE within a more detailed picture of the electron transfer in the P3HT-PCBM model. In this case, the electron was first optically excited into a higher energy state localized on the donor (neutral state in Figure 8), such as that qualitatively depicted by a LUMO+3 in Fig. 2(b) for P3HT/PCBM heterojunction (or equivalently the LUMO on the donor), leaving a positive hole in the HOMO. We will refer to this excited state as the ’neutral state’ and assume that the electronic charge is transferred to the energetically nearest lower electronic state localized on the acceptor (such as that qualitatively depicted by a LUMO+2 in Fig. 2(d) ), which will be referred to as the ’charge transfer’ state. Fig. 8 shows three such CT states (black) coupled to neutral excited states (green) near 2 – 2.5eV energy range.

2.2

Side-chain deuterated 2.2

2.1

2.1 Energy [eV]

Protonated

Energy [eV]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 44

2

2

1.9

1.9

1.8

1.8 0

2000

4000 6000 Time [a.u.]

8000

0

10000

2000

4000 6000 Time [a.u.]

8000

10000

Figure 9: Time evolution of several excited states from TD-DFT at B3LYP/6-31G(d) level of theory for pristine P3HT/PCBM (left) and for side-chain deuterated SD-P3HT/PCBM (right). According to a recent first-principles theory, the coupling of the neutral and interfacial CT states 19 or tunneling from excitonic states to CT states 18 in a material or device, may determine the charge transfer process. According to the theory of Bittner and Silva, 18 noise may promote fast charge transfer via a neutral–CT states tunneling mechanism. Within the two-state model defined by the Hamiltonian,  

H= 

E



W 

W −E

, 

30

E = ∆E /2,


(23)

Page 31 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the off-diagonal element, W , couples the optically excited neutral state (localized on the donor) with the charge transfer state (localized on acceptor). This coupling promotes resonant tunneling from the excited neutral to charge transfer state in the limit of E ≪ W and suppresses such transitions in the opposite limit, E ≫ W , due to increased charge recombination on the donor. For a given (fixed) values of ∆E and W of a system that is initially in the neutral state (localized on donor) the probability of finding system in the acceptor state oscillates in time according to Rabi formula with probability amplitude P M AX =

W2 W 2 +E 2

and with frequency equal to half of difference between eigenvalues of the coupled two-level system: PCT (t) =

W2 t · sin2 |E+ − E− | 2 2 W +E 2¯h

(24)

After averaging over all times one can obtain a ”time-independent” probability as

P CT =

1 1 W2 = P M AX · 2 2 2 W +E 2

(25)

However, our model consists an ensemble of structures from QTES that collectively represent a nuclear wavefunction for H/D. Specifically, the nuclear wavefunction is described through the fluctuation of coordinates for quantum mechanically treated atoms H and D. Taking into account the fluctuations of H/D positions leads to the expectation value (nuclear wavefunction averaged) of charge transfer probabillity:

hP CT i = N

−1

Z N 0

P CT (ξ)dξ

(26)

where N is a normalization factor and ξ numbers the structures in the nuclear wavefunction ensemble. To analyze charge tunneling probabilities we estimated the energy difference ∆E and coupling W between excited neutral state of donor and charge transfer states of acceptor for a set of structures collected from QTES simulations for hydrogenated and deuterated

31



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

P3HT/PCBM. We used a modified version of the Valeev’s 46 procedure for determination of charge transfer integrals from the Löwdin orthogonalized 2x2 Hamiltonian expressed in the basis of the donor and acceptor diabats. Our modification accounted for the fact that there are three available ’target’ states on the acceptor (LUMO, LUMO+1, LUMO+2) below the optically active neutral excited state of the donor (LUMO+3). Thus, instead of 2×2 problem we deal with 4×4 problem (one donor state and three acceptor states). These four diabats are obtained from the ground state Hamiltonian of the full (donor+acceptor) system as the eigenstates of the diagonal donor-donor and acceptor-acceptor blocks. The reduced size (4×4) Hamiltonian and overlap matrix for the charge trasnfer problem are then computed in the diabatic basis set from the off-diagonal donor-acceptor block of the full Hamiltonian and full overlap matrix. Finally, the 4x4 Hamiltonian matrix was orthogonalized using the Löwdin scheme. The non-diagonal elements of the 4x4 Löwdin orthogonalized diabatic Hamiltonian correspond to the coupling terms W in the Eq. (23) while the difference between diagonal terms correspond to ∆E . Similarly to the Eg analysis in previous section, the analysis of 4x4 coupling Hamiltonian has been performed for B3LYP, CAM-B3LYP and LC-ωPBEh for several structures from QTES simulations. Our analysis of 4x4 charge transfer Hamiltonian shows that: (a) Couplings between donor and acceptor states are usually much smaller than the ∆E for the corresponding states (calculated as the difference between diagonal terms). For example, the ∆E from CAM-B3LYP is between the donor neutral state and the acceptor state is ∼0.8 eV, whereas the coupling (charge transfer integral), W, is usually ∼70 meV. (b) All three acceptor states have a significant W coupling, and the correct model should likely consist of all three acceptor states rather than a single state. For example, the charge transfer coupling elements from CAM-B3LYP/6-31G(d) for a selected (initial) hydrogenated structure is 80, 61 and 9 meV for the coupling of the donor state (LUMO+3) to, respectively, LUMO+2, LUMO+1 and LUMO of the acceptor. (c) The range of fluctuations of ∆E is much larger than the range of variation of W . (d) Due to the perturbation from the donor, the three acceptor states (LUMO,

32


Page 32 of 44

Page 33 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


LUMO+1 and LUMO+2) are coupled to each other. The magnitude of this couplings is small. It amounts to ∼0.1-1 meV and vanishes in the absence of the donor. The non-zero couplings between acceptor states in the presence of donor states suggests a possibility of charge relaxation cascade between acceptor states that follows the initial charge tunneling from donor to acceptor. Table 4 compares couplings W and energy difference ∆E = Eneutral − ECT calculated as an average over QTES structures for hydrogenated and deuterated species. We used three different approximations for caluclation of ∆E and three functionals (B3LYP, CAM-B3LYP and LC-ωPBEh). In the first model, denoted as ∆EM O in the Table 4, the excited neutral and charge transfer state are approximated as LUMO+3 and LUMO+2 (see Figure 2). In the second model, denoted as ∆Ediabats in the Table 4, the excited states were approximated from the states of the 4x4 charge transfer diabatic Hamiltonian. In both cases ∆E is calculated as a differenece between LUMO+3 neutral and LUMO+2. (LUMO+3 corresponds to excited neutral donor state; LUMO+2 corresponds to energetically closest charge transfer state) In the third model based on TD-DFT (denoted as ∆ET D−DF T ) ∆E is calculated as a difference between lowest optically active neutral excited state and energetically closest charge transfer excited state. The models based on frontier molecular orbitals ( ∆EM O ) and on the diabatic charge transfer Hamiltonian (∆Ediabats ) give values ∆E which are much larger than the value of coupling W , regardless of the functional used (B3LYP, CAM-B3LYP, LC-ωPBEh). TD-DFT analysis gives ∆E for the H/D species that is comparable with W . One can expect that for a given functional the true ∆E is likely somewhere in between frontier orbitals/diabatic and TD-DFT models. That is LUMOs are rather crude approximation of excited states charge transfer and they overestimate the excitation energy as can be seen in Table 4. Similarly TD-DFT states are only approximation of diabatic states (neutral excited state has usually a small contribution from acceptor while charge transfer state has usually small contribution from donor) and hence they understimate the energy of diabatic states.

33



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 44

We will now focus on the effect of fluctuations on the charge transfer probabilities in Eq. (25). As indicated in the Table 4 by the RMS values, the fluctuations of W are negligible compared to the fluctuations of ∆E . Next, the fluctuations of the ∆E for hydrogenated species (P3HT/PCBM) are larger (for all models) than for deuterated ones (SD-P3HT/PCBM). This observation is consistent with a larger span of nuclear wavefucntion for hydrogens than for deuteriums. That is, calculations of ∆E (see Figure 9) along the trajectories yield small variations between the protonic and deuteronic species. These isotope-specific fluctuations are essentially the isotope-specific noise, which is small compared to Eg , but is appreciable compared to ∆E . 45 The side-chain deuterium substitution, while not appreciably changing the energies of the electronic states, reduces carrier density in SD-P3HT due to reduced probability of donor-acceptor charge transfer and increased probability of charge recombination on the donor. Within Eq. (22) the drop in VOC comes from the increase in Vloss for SD-P3HT compared to P3HT. Thus, our calculations suggest that for the hydrogenated species the decreased value of E = ∆E /2 due to their larger fluctuations (RMS values in the Table 4) leads to increased charge transfer probabilities. The effect of the ∆E variance caused by the hydrogenic/deuteronics nuclear wavefunction can be included into Eq. 25 through the root-mean square of ∆E dependence of E: P CT (ξ) =

1 W2 1 W2 · 2 = · 2 W + E2 2 W 2 + ∆E +α cos(ξ) 2 2

where ’α cos(ξ)’ term describes variance of ∆E with the amplitude α =

(27) √

2 · RM S and ξ is a

parameter ranging from 0 to 2π. The charge transfer probabilities averaged over diffeerent configurations of P3HT/PCBM can then be evaluated as a definite integral

hP CT i = (2π)−1

Z 2π 0

P CT (ξ)dξ

(28)

Equation (27) can be interpreted as charge transfer probabilities for P3HT/PCBM at the influence of the H/D isotope sepcific RMS noise. Then Equation (28) describes charge 34


Page 35 of 44

transfer probabilities averaged over isotope specific nuclear wavefunction. Figure 10 shows P CT (ξ) and hP CT i for CAM-B3LYP based values of W and ∆E . 0.06

0.14

W=75 meV W=60 meV W=40 meV W=20 meV W=10 meV

Probability, 〈PCT〉ξ

0.12

Probability, PCT(ξ)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.04

0.02

0.1 0.08 0.06 0.04 0.02

RMS=0.3 eV RMS=0.2 eV

0 0

π/2

2π

3π/2

0 0 2π

Nuclear fluctuation, ξ

0.2 0.4 0.6 0.8 Amplitude of ∆E fluctuation (eV)

1

Figure 10: (a) Fluctuation of charge transfer probabilities for W =75 meV, ∆E =0.84 eV and RMS of ∆E equal to 0.2 and 0.3 eV from Eq. (27). Fluctuation of ∆E modulates charge transfer probabilities. (b) Increasing fluctuations of ∆E for a given, fixed values of W and ∆E increases ξ averaged probability for charge transfer hPCT i until it reaches its maximum. At the apex (near 0.8eV) the amplitude of fluctuations becomes equal to ∆E which results in vanishing gap. Here ∆E =0.84 eV (CAM-B3LYP). As can be seen on Figure 10 the fluctuations ∆E can modulate the the charge transfer probabilities even when W and ∆E are on average unchanged. Our simulations suggest that stronger fluctuations of ∆E for H than for D are responsible for increased charge transfer probabilities and lead to larger Voc. In this section, by employing a similar to Bittner two-state model of charge transfer in conjuction with Rabi’s analytical expressions for dynamical behavior of the two state model we have demonstrated that seemingly unrelated fluctuation can cause an observable change in the performance of P3HT/PCBM photovoltaic material. Next, by applying a 4 state model (a single donor state + 3 acceptor states), we have also demonstrated that more detailed studies of charge transfer in P3HT/PCBM studies should not be limited to coupling between a single donor state (neutral excited state) with a single acceptor state (energetically closest charge transter state). Our calculations show that the magnitude of 35



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

coupling between all three acceptor states (LUMO,LUMO+1, LUMO+2) with the excited donor state (LUMO+3) is similar. Ideally, all three acceptor states as well as the HOMO (to account for recombination rates) should be included in more detailed studies. The authors believe, that viable path forward is to combine the QTES approach 20 with real-time electron dynamics. 38,47

5

Summary

The goals of this study were threefold: (i) to develop a quantum dynamics methodology for multiple quantum nuclei suitable for simulations of the isotope effect in conducting polymers, (ii) to expand the QTES methodology, particularly the treatment of quantum forces, by including higher order terms, (iii) to model an unexpected NQE caused by side-chain H/D substitution in a P3HT/PCBM blend on the open circuit voltage. Our molecular model of P3HT/PCBM, involving a single PCBM and short P3HT polymer is a minimalistic representation of a heterojunction. However, the NQE itself should manifest already in a single-molecule model and thus yield valuable qualitative information for understanding the previously experimentally observed effects. Since the electronic structure is isotopeindependent we examined variations in the nuclear wavefunctions of the two isotopologues. In summary, the work presented in this paper established the following: 1. To investigate the effect of H/D isotope substitution on the performance of conducting polymers we proposed using a combination of quantum trajectories with two electronic structure methods. A cheap electronic structure method such as DFTB is used within QTES method to account for a quantum nuclear effect of hydrogen or deuterium atoms. Snapshots of molecular structure from QTES dynamics collectively represent the nuclear wavefunction and can be used with DFT or other electronic structure methods to study optoelectronic or magnetic properties. 2. The quadratic basis quantum force correction in QTES is an adequate extension of the 36


Page 36 of 44

Page 37 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


linear quantum force for quantum mechanical simulations of multiple quantum nuclei. Specifically, in order to describe multiple anharmonic CH bonds, we have developed the non-linear quantum corrections for the QT dynamics, where an ensemble of QTs represents the nuclear wavefunction evolving in time. The corrections were introduced into the local modes of the protons allowing for the quantum mode-coupling for each proton, but not between different protons. 3. The QTES-DFTB dynamics involving 39 H/D atoms of three hexyl side chains for up to 500 fs was used to assess the nuclear wavefunction localization. Given the ultrafast nature of the CT process (on the order of 100 fs), the differences in the ground state nuclear wavefunctions determine the differences corresponding to the electronic structure of the two isotopologues, leading to the difference in VOC of 0.24 eV. 4. We analyzed the electronic energy levels and charge transfer couplings along the quantum trajectories describing the H/D nuclear wavefunction and for the geometries averaged over those trajectories, with the latter showing more variations between the H/D species. The fluctuations due to nuclear wavefunctions associated with the substituted atoms, of the electronic energy levels in the optical transition range, are small compared to the HOMO-LUMO gap. Fluctuations of energy gap are larger for H than for D regardless of theory used. Fluctuations of charge coupling integrals are much smaller than the fluctuations of ∆E and can be neglected. 5. Our simulations strongly suggest that the isotope-dependent fluctuations may account for the observed NQE by promoting CT via tunneling in P3HT compared to SD-P3HT and increasing carrier loss in SD-P3HT. Our two-state model of charge transfer in conjuction with Rabi’s analytical expressions for charge transfer probabilities demonstrate that seemingly unrelated fluctuation of H/D modulates the energy of diabatic states involved in charge transfer and this in turn can cause an observable change in charge transfer probabilities that affect the performance of P3HT/PCBM photovoltaic cell. 37



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Experiments on the temperature dependence of the NQE in P3HT/PCBM would provide further insight into the mechanism of NQE in these type of systems. 6. We applied a modified Valeev’s model 46 to estimate charge transfer couplings between the optically active excited neutral state of the donor (dominated by LUMO+3) and three charge transfer acceptor states (dominated by LUMO, LUMO+1, LUMO+2). The analysis of this 4 state diabatic Hamiltonian model suggests: (a) that magnitude of coupling between all three acceptor states (LUMO,LUMO+1, LUMO+2) and the excited donor state (LUMO+3) is similar, and (b) more detailed studies of charge transfer in P3HT/PCBM should include several states of the donor and acceptor. One promising path forward would be to combine the QTES approach with real-time electron dynamics.

6

Acknowledgements

This material is based upon work partially supported by the National Science Foundation under Grants No. CHE-1056188 (S.G.). The work regading P3HT was conducted at the Center for Nanophase Materials Sciences, a U.S. Department of Energy Office of Science User Facility. The XSEDE allocation TG-DMR110037 and use of the USC HPC cluster funded by the National Science Foundation under Grant No. CHE-1048629 are also acknowledged. The authors declare no competing financial interest.

References (1) Nalwa, H. S. Ed. handbook of Organic Conductive Molecules and Polymer ; John Wiley, 1997. (2) Li, X. G.; Huang, M. R.; Duan, W.; Yang, Y. L. Chem. Rev. 2002, 102, 2925.

38


Page 38 of 44

Page 39 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(3) Yang, H.; Shibayama, M.; Stein, R. S.; Shimizu, N.; Hashimoto, T. Macromolecules 1986, 19, 1667–2447. (4) White, R. P.; Lipson, J. E. G.; Higgins, J. S. Macromolecules 2010, 43, 4287–4293. (5) Shao, M.; Keum, J.; Chen, J.; He, Y.; Chen, W.; Browning, J. F.; Jakowski, J.; Sumpter, B. G.; Ivanov, I. N.; Ma, Y.-Z.; Rouleau, C. M.; Smith, S. C.; Geohegan, D. B.; Hong, K.; Xiao, K. Nature Commun. 2014, 5, 3180. (6) Asadi, K.; Kronemeijer, A. J.; Cramer, T.; Koster, L. J. A.; Blom, P. W.; de Leeuw, D. M. Nature Commun. 2013, 4, 1710. (7) Jiang, Y.; Geng, H.; Shi, W.; Peng, Q.; Zheng, X.; Shuai, Z. J. Phys. Chem. Lett. 2014, 5, 2267–2273. (8) Tsuji, H.; Mitsui, C.; Nakamura, E. Chem. Commun. 2014, 50, 14870–14872. (9) W. J. Potscavage Jr.,; Sharma, A.; Kippelen, B. Acc. Chem. Res. 2011, 42 . (10) Perez, M. D.; Borek, C.; Forrest, S. R.; Thompson, M. E. J. Am. Chem. Soc. 2009, 131, 9281–9286. (11) Tong, C. C.; Hwang, K. C. J. Phys. Chem. C 2007, 111, 3490–3494. (12) Nguyen, T.; Hukic-Markosian, G.; Wang, F.; Wojcik, L.; Li, X.-G.; Ehrenfreund, E.; Vardeny, Z. V. Nature Materials. 2010, 9, 345. (13) Nietfeld, J. P.; Schwiderski, R. L.; Gonnella, T. P.; Rasmussen, S. C. J. Org. Chem. 2011, 76, 6383. (14) Bourass, M.; Fitri, A.; Benjelloun, A. T.; Benzakour, M.; Mcharfi, M.; Hamidi, M.; Serein-Spirau, F.; Jarrosson, T.; L` ere-Porte, J. P.; Sotiropoulos, J. M.; Bouachrine, M. Der Pharma Chemica 2013, 5, 144–153.

39



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(15) Kroon, R.; Lenes, M.; Hummelen, J. C.; Blom, P. W. M.; Boer, B. D. Polymer Reviews 2008, 48, 531–582. (16) Tong, C. C.; Hwang, K. C. J. Phys. Chem. C 2007, 111, 3490–3494. (17) Nguyen, T. D.; Markosian, G. H.; Wang, F. J.; Wojcik, L.; Li, X. G.; Ehrenfreund, E.; Vardeny, Z. V. Nat. Mater. 2010, 9, 345–352. (18) Bittner, E. R.; Silva, C. Nat. Commun. 2014, 5 . (19) Tamura, H.; Burghardt, I. J. Am. Chem. Soc. 2013, 135, 16364–16367. (20) Garashchuk, S.; Jakowski, J.; Wang, L.; Sumpter, B. G. J. Chem. Theory and Comp. 2013, 9, 5221–5235. (21) Porezag, D.; Fraunheim, T.; Kohler, T.; Serfert, G.; Kaschner, R. Phys. Rev. B 1995, 51, 12947–12957. (22) Seifert, G.; Porezag, D.; Frauenheim, T. Int. J. Quantum Chem. 1996, 58, 185–192, International Workshop on Electronic Structure Methods for Truly Large Systems Moving the Frontiers in Quantum Chemistry, Braungale, Germany, AUG 01-07, 1994. (23) Frauenheim, T.; Seifert, G.; Elstner, M.; Hajnal, Z.; Jungnickel, G.; Porezag, D.; Suhai, S.; Scholz, R. Phys. Status Solidi B 2000, 217, 41–62. (24) Zheng, G.; Lundberg, M.; Jakowski, J.; Vreven, T.; Frisch, M. J.; Morokuma, K. Int. J. Quantum Chem. 2009, 109, 1841–1854. (25) Garashchuk, S.; Volkov, M. V. Molecular Physics 2012, 110, 985–993. (26) Mazzuca, J. W.; Garashchuk, S.; Jakowski, J. Chem. Phys. Lett. 2014, 613, 104–109. (27) Gu, B.; Hinde, R. J.; Rassolov, V. A.; Garashchuk, S. J. Chem. Theory Comput. 2015, 11, 2891–2899.

40


Page 40 of 44

Page 41 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(28) Habershon, S.; Manolopoulos, D. E. J. Chem. Phys. 2009, 131, 244302. (29) Czako, G.; Kaledin, A. L.; Bowman, J. M. J. Chem. Phys. 2010, 132 . (30) Frauenheim, T.; Seifert, G.; Elstner, M.; Niehaus, T.; Kohler, C.; Amkreutz, M.; Sternberg, M.; Hajnal, Z.; Di Carlo, A.; Suhai, S. J. Phys.: Condens. Matter 2002, 14, 3015. (31) Seifert, G.; Joswig, J.-O. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2012, 2, 456–465. (32) Gaus, M.; Cui, Q.; Elstner, M. Interdiscip. Rev.: Comput. Mol. Sci. 2014, 4, 49–61. (33) Cui, Q.; Elstner, M. Phys. Chem. Chem. Phys. 2014, 16, 14368–14377. (34) Wang, L.; Jakowski, J.; Garashchuk, S. J. Phys. Chem. C 2014, 118, 16175–16187. (35) Madelung, E. Z. Phys. 1927, 40, 322–326. (36) Bohm, D. Phys. Rev. 1952, 85, 166–193. (37) Garashchuk, S.; Rassolov, V. A. J. Chem. Phys. 2004, 120, 1181–1190. (38) Jakowski, J.; Irle, S.; Sumpter, B. G.; Morokuma, K. The Journal of Physical Chemistry Letters 2012, 3, 1536–1542, PMID: 26285634. (39) Becke, Axel D., J. Chem. Phys. 1993, 98, 5648–5652. (40) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B 1988, 37, 785–789. (41) Burke, K.; Werschnik, J.; Gross, E. K. U. The Journal of Chemical Physics 2005, 123, 062206. (42) Scharber, M.; M¨ uhlbacher, D.; Koppe, M.; Denk, P.; Waldauf, C.; Heeger, A.; Brabec, C. Advanced Materials 2006, 18, 789–794. (43) Shao, Y.; Gan, Z.; Epifanovsky, E.; Gilbert, A. T. B.; Wormit, M.; Kussmann, J.; Lange, A. W.; Behn, A.; Deng, J.; Feng, X.; Ghosh, D.; Goldey, M.; Horn, P. R.; Jacobson, L. D.; Kaliman, I.; Khaliullin, R. Z.; Kus, T.; Landau, A.; Liu, J.; Proynov, E. I.; 41



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Rhee, Y. M.; Richard, R. M.; Rohrdanz, M. A.; Steele, R. P.; Sundstrom, E. J.; Woodcock, H. L., III; Zimmerman, P. M.; Zuev, D.; Albrecht, B.; Alguire, E.; Austin, B.; Beran, G. J. O.; Bernard, Y. A.; Berquist, E.; Brandhorst, K.; Bravaya, K. B.; Brown, S. T.; Casanova, D.; Chang, C.-M.; Chen, Y.; Chien, S. H.; Closser, K. D.; Crittenden, D. L.; Diedenhofen, M.; DiStasio, R. A., Jr.; Do, H.; Dutoi, A. D.; Edgar, R. G.; Fatehi, S.; Fusti-Molnar, L.; Ghysels, A.; Golubeva-Zadorozhnaya, A.; Gomes, J.; Hanson-Heine, M. W. D.; Harbach, P. H. P.; Hauser, A. W.; Hohenstein, E. G.; Holden, Z. C.; Jagau, T.-C.; Ji, H.; Kaduk, B.; Khistyaev, K.; Kim, J.; Kim, J.; King, R. A.; Klunzinger, P.; Kosenkov, D.; Kowalczyk, T.; Krauter, C. M.; Lao, K. U.; Laurent, A. D.; Lawler, K. V.; Levchenko, S. V.; Lin, C. Y.; Liu, F.; Livshits, E.; Lochan, R. C.; Luenser, A.; Manohar, P.; Manzer, S. F.; Mao, S.-P.; Mardirossian, N.; Marenich, A. V.; Maurer, S. A.; Mayhall, N. J.; Neuscamman, E.; Oana, C. M.; Olivares-Amaya, R.; O’Neill, D. P.; Parkhill, J. A.; Perrine, T. M.; Peverati, R.; Prociuk, A.; Rehn, D. R.; Rosta, E.; Russ, N. J.; Sharada, S. M.; Sharma, S.; Small, D. W.; Sodt, A.; Stein, T.; Stueck, D.; Su, Y.-C.; Thom, A. J. W.; Tsuchimochi, T.; Vanovschi, V.; Vogt, L.; Vydrov, O.; Wang, T.; Watson, M. A.; Wenzel, J.; White, A.; Williams, C. F.; Yang, J.; Yeganeh, S.; Yost, S. R.; You, Z.-Q.; Zhang, I. Y.; Zhang, X.; Zhao, Y.; Brooks, B. R.; Chan, G. K. L.; Chipman, D. M.; Cramer, C. J.; Goddard, W. A., III; Gordon, M. S.; Hehre, W. J.; Klamt, A.; Schaefer, H. F., III; Schmidt, M. W.; Sherrill, C. D.; Truhlar, D. G.; Warshel, A.; Xu, X.; Aspuru-Guzik, A.; Baer, R.; Bell, A. T.; Besley, N. A.; Chai, J.-D.; Dreuw, A.; Dunietz, B. D.; Furlani, T. R.; Gwaltney, S. R.; Hsu, C.-P.; Jung, Y.; Kong, J.; Lambrecht, D. S.; Liang, W.; Ochsenfeld, C.; Rassolov, V. A.; Slipchenko, L. V.; Subotnik, J. E.; Van Voorhis, T.; Herbert, J. M.; Krylov, A. I.; Gill, P. M. W.; HeadGordon, M. Molecular Physics 2015, 113, 184–215. (44) Kaduk, B.; Kowalczyk, T.; Van Voorhis, T. Chemical Reviews 2012, 112, 321–370.

42


Page 42 of 44

Page 43 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(45) Valiev, M.; Bylaska, E.; Govind, N.; Kowalski, K.; Straatsma, T.; Dam, H. V.; Wang, D.; Nieplocha, J.; Apra, E.; Windus, T.; de Jong, W. Computer Physics Communications 2010, 181, 1477 – 1489. (46) Valeev, E.; Coropceanu, V.; da Silva Filho, D. A.; Salman, S.; Bredas, J.-L. J. Am. Chem. Soc. 2006, 128, 9882-9886 . (47) Jakowski, J; Morokuma, K. J. Chem. Phys. 2009, 130, 224106 .

43


EDEUTERATED-P3HT

Electron donor Acceptor Journal of Chemical Theory and Page Computation 44 of 44 P3HT

(c)

1 2 3 4 5 6 7

(a)

EPCBM(d)

PCBM

CT probability Time

EHYDROGENATED-P3HT (b)

ACS Paragon Plus Environment EPCBM CT probability

PCBM: A

Recommend Documents