Effects of Dispersion Forces on Structure and Photoinduced Charge

Subscriber access provided by University of Texas Libraries

Article

Effects of Dispersion Forces on Structure and Photoinduced Charge Separation in Organic Photovoltaics Juan Pablo Martinez, Daniel Eduardo Trujillo-González, Andreas W. Götz, Fray de Landa Castillo-Alvarado, and Juan Ignacio Rodriguez Hernandez J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.7b05107 • Publication Date (Web): 25 Aug 2017 Downloaded from http://pubs.acs.org on August 28, 2017

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.

The Journal of Physical Chemistry C 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 24

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 Journal of Physical Chemistry

Effects of Dispersion Forces on Structure and Photoinduced Charge Separation in Organic Photovoltaics Juan Pablo Martínez,∗,† Daniel Eduardo Trujillo-González,† Andreas W. Götz,∗,‡ Fray L. Castillo-Alvarado,† and Juan I. Rodríguez∗,† †Escuela Superior de Física y Matemáticas, Instituto Politécnico Nacional, Edificio 9, U.P. A.L.M., Col. San Pedro Zacatenco, C.P. 07738, Ciudad de México, México ‡San Diego Supercomputer Center, University of California San Diego, 9500 Gilman Drive, La Jolla, CA 92093-0505, USA E-mail: [email protected]; [email protected]; [email protected]

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

Abstract We present a theoretical study on the role of van der Waals (vdW) interactions on the structure and, as a consequence, the photoinduced charge separation (CS) of a series of dimer complexes formed by the polymer P3HT and the fullerene derivative PCBM. CS rate constants for P3HT/PCBM dimer structures in which vdW interactions are taken into account agree well with experimental data. Without proper treatment of vdW interactions during geometry optimizations, the predicted CS rates can be too low by up to three orders of magnitude. These variations in computed CS rates are not due to changes in the Gibbs energy for CS. Instead, the electronic coupling increases by up to two orders of magnitude for structures obtained with dispersion-corrected density functionals that lead to deformations in the P3HT oligomer with pronounced π-π stacking interactions with PCBM.

Introduction The use of π-conjugated polymers as electron donors and fullerene derivatives as electron acceptors in bulk heterojunction (BHJ) organic photovoltaic devices (OPVs) represents an attractive alternative to silicon-based solar cells due to their potential low cost, flexibility, negligible toxicity, and short energy payback time. 1–7 Fullerene C60 and its derivatives perform well in OPVs due to their excellent electron acceptor capability, 8,9 and also their small internal reorganization energy. 10,11 Due to its high power conversion efficiency (PCE≈5%) 12,13 and internal quantum efficiency (≈0.8), 14 experimental studies based on transient absorption spectroscopy have been focused on the OPV blends made of poly(3-hexylthiophene) (P3HT) and [6,6]-phenyl-C61 -butyric acid methyl ester (PCBM), 15–17 reporting rate constants of charge separation (CS) kCS > 4 × 1012 s−1 . 18,19 Morphology of the BHJ active layer mainly at the donor-acceptor interface, the intermolecular distance (D) between donor and acceptor molecules, and exciton formation/dissociation/recombination are key factors that determine the CS process. 20–23 2


Page 2 of 24

Page 3 of 24

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


Based on Marcus theory and its extensions, 24 density functional theory (DFT) based theoretical studies have been focused on the calculation of the CS and charge recombination (CR) rate constants of the P3HT/PCBM system, modeling P3HT with oligomers of several sizes. 22,25,26 The common-denominator conclusions of these studies are: 1) there is an uncertainty in the theoretically predicted structure at the P3HT/PCBM interface; 2) the rate constants (kCS and kCR ) can vary orders of magnitude depending on a particular structure of the P3HT/PCBM interface; 3) the theoretically computed values for kCS and kCR are usually 1 to 2 orders of magnitude below the experimental values. 25,26 Although there is strong experimental and theoretical evidence that van der Waals interactions (vdW) are indispensable to reach the so-called chemical accuracy in the structure and interaction energies of supramolecular systems formed by fullerenes and aromatic molecules, 27–31 surprisingly none of the DFT studies that reports CS and CR rates considered explicitly vdW interactions for the P3HT/PCBM interface. In this letter, the influence of vdW interactions on the electron transfer parameters that directly affect the photoinduced CS in the supramolecular dimer formed by an 8-unit oligomer of P3HT and PCBM are quantum chemically analyzed in detail for a series of five P3HT/PCBM dimers (see Figure 1).

Methodology and computational details The rate constant kCS was computed via Marcus theory for five P3HT/PCBM dimers (Figure 1), in two of which (I1 and I2) the P3HT adopts a U-shaped structure surrounding to PCBM due to vdW interactions as recently reported by us 32,33 (see Figure 1 (top) and discussion below). Isomer dimers I1 and I2 are local minima of the PES. Geometries were optimized in the gas phase at the DFT level using the Perdew-Burke-Ernzerhof (PBE) 34,35 exchange-correlation (xc) potential as implemented in the Amsterdam Density Functional (ADF2016) package. 36–38 The vdW interactions were taken into account via Grimme’s dis-

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

persion correction scheme, termed DFT-D3, 39 in which the dispersion coefficient C6 for each pair of elements are obtained from ab-initio calculations using the Casimir-Polder formula considering the hydrates of each element, and the C8 coefficient is obtained recursively from C6 . All other parameters (cutoff radii, damping function parameter, etc.) are obtained empirically and only two of them are xc potential dependent (PBE-D3). 39 An all-electron uncontracted set of Slater-type orbitals (STOs) of triple-ζ quality (TZP) containing one set of polarization and diffuse functions was used to expand molecular orbitals for the geometry optimization procedure. Via time-dependent DFT (TDDFT) calculations, electron transfer parameters were determined in the gas phase at the PBE-D3/TZP optimized geometries for isomers I1 and I2. To better describe excitations with charge transfer (CT) character, the range-separated xc functional CAMY-B3LYP 40,41 was chosen to evaluate the energy terms involved in the photoinduced CS between P3HT and PCBM. 42 The frozen core approximation (FCA) was included to decrease computer demand, which freezes the core orbitals (1s for carbon and oxygen, and up to 2p for sulfur). 36 It is known that FCA usually has a negligible effect in the electronic transitions of interest since the main contributions come from frontier molecular orbitals (see below for our benchmark calculations on the adenine dimer electronic coupling). It has been reported that the P3HT alkyl side chains have a significant effect on the structure of isolated P3HT 43 and P3HT in the P3HT/PCBM dimer. 32,33 However, alkyl groups have hardly any effect on the optical properties and CS parameters. 26,44,45 Thus, in order to save computational resources, each hexyl group is replaced by a methyl group in the P3HT chain for the TDDFT calculations reported in this work. For dimers I3, I4 and I5, the CS parameters were computed via TDDFT single point calculations. It is worth mentioning that we tried to optimize I3 using the same xc potential PBE but without Grimme’s corrections to properly quantify how the lack of the vdW interactions affects the value of the charge separation rate constant kCS . However, these caculations did not lead to any binding between P3HT and PCBM. So we decided to consider these dimers (I3-I5) with ”idealistic” structures to have a reference point that is representative of approaches

4


Page 4 of 24

Page 5 of 24


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


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 6 of 24

the electronic coupling, respectively; T is the temperature, kB and h ¯ are the Boltzmann and reduced Planck constants. The parameter λ is usually divided into two contributions: internal and external reorganization energies. The former, λint , is the energy required to rearrange all the nuclei of the system due to CS from a neutral to a charged state. 24 λint can be calculated by considering isolated donor and acceptor fragments as a separated contribution to the internal reorganization, which is a well known approximation, 46

λint = λD + λA

λD(A) =

1 ′ ′ En − En + Eion − Eion , 2

(2)

(3)

where λD(A) is the reorganization energy of either the donor (D) or the acceptor (A). En(ion) is the energy of the neutral (ionic) state computed in the equilibrium geometry of the neutral (ionic) state; and En(ion) is the energy of the neutral (ionic) state computed in the ′

equilibrium geometry of the ionic (neutral) state. The ionic states are +1e− cation and -1e− anion states for the donor and acceptor fragments, respectively. The optimized geometries for both the neutral and ion for each fragment, donor (P3HT) and acceptor (PCBM), were obtained at the DFT level in order to apply Eq (3). The external reorganization energy corresponds to changes in the environment. However, since we are mainly interested in the influence of the structural conformation of the P3HT/PCBM dimer in the CS process, environmental effects are excluded from the current work. It should be noted that in some cases the environment can stabilize CT states such that the CS process enters the so called far inverted regime ( |∆G| ≫ λ) for which Marcus theory (Eq. (1)) can not be aplied. 47,48 Instead other more proper theories must be used. 49,50 Electron transfer coupling Vij was assessed using the orbital approximation within the energy splitting in dimer (ESID) method. 51–53 An all-electron DZP basis set was used to obtain Vij . As a benchmark calculation, the CS electronic coupling Vij was computed for the adenine dimer obtaining a value of 0.057 eV, which

6


Page 7 of 24

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


is in relatively good agreement with the reported value of 0.068 eV obtained via the fragment charge difference (FCD) method (the same value of 0.057 eV was obtained considering the FCA basis set approximation). It is worth mentioning that more accurate values are expected from a FCD approach since adiabatic states are computed using multi-configuration excited states as generated by TDDFT; therefore, it does not use orbital approximations to represent locally excited and CS states. 54,55 The good performance is attributed to the use of the range-separated CAMY-B3LYP/DZP, since ESID via PBE/DZP results in a Vij of 0.036 eV; that is, half of the expected value. The CS reaction under study is understood as follows:

P3HT/PCBM + hν −→ P3HT∗ /PCBM P3HT∗ /PCBM −→ P3HT+ /PCBM−

(4) (5)

wherein a molecular exciton localized at the P3HT fragment is initially formed by photon absorption (Equation 4), then the CS occurs after the dissociation of such a Frenkel exciton (Equation 5). Assuming a negligible entropic component, which may be valid under the Franck-Condon approximation, changes in electronic energy, ∆ECS , for CS reactions nearly correspond to ∆GCS . With the aim of assigning the excited states in Equations 4 and 5, the 30 lowest-in-energy electronic transitions were estimated under the Tamm-Dancoff approximation. 56

Results and discussion The five P3HT/PCBM dimer structures under study are schematized in Figure 1. For Isomer 1 (I1) and Isomer 2 (I2) (see Figure 1, top), vdW interactions were taken into account via the PBE-D3 scheme in a local geometry optimization procedure. As reported previously by us (see reference 32 for a detailed description of these isomers), when considering the vdW correction, the initially planar P3HT chain gets bent towards PCBM. In the final stable 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

structures of I1 and I2, the P3HT chain is in a U-shaped conformation surrounding PCBM (the alkyl side chains get also bent towards PCBM). In isomers I1 and I2, the adjacent thiophene rings in the P3HT backbone chain show an alternate torsion (≈ 15 deg), which leads not only to a loss of planarity but also to a loss of the local translational symmetry that is related to the head-tail arrangement between neighboring thiophene rings. 32 The other three dimers considered in this letter were constructed by combining the optimized geometries of PCBM and planar P3HT, from I3 to I5. The intermolecular distance between PCBM and P3HT for the two U-shaped isomers and the planar dimer I3 is the typical π − π stacking distance (≈ 3.5 Å). We intentionally locate the planar P3HT at 4.5 and 5.5 Å from PCBM for dimers I4 and I5, respectively (see Figure 1, bottom). So, for P3HT/PCBM dimers I3 to I5, the vdW interactions were not taken into account. Table 1 shows the frontier orbital localization according to the fragments P3HT and PCBM for isomer I1. Interestingly, this distribution is nearly identical for all other dimers. 57 It can be observed that, in all dimers, molecular orbitals from HOMO-2 to HOMO are fully located at the P3HT fragment. Lower energy HOMOs are located at the PCBM moiety. Besides, frontier LUMOs are localized at the PCBM fragment; except for LUMO+4 (or LUMO+5 in the other dimers). In view of that, a low-in-energy electronic transition leading to the formation of an exciton fully situated at P3HT can take place from HOMO (or HOMO1, or HOMO-2) to LUMO+4. The resulting state can be represented as P3HT∗ /PCBM (see Equation 4). On the other hand, an electronic transition characterized by strong CS takes place from one of the frontier HOMOs located at P3HT to a LUMO located at PCBM, thus giving rise to a CS state represented as P3HT+ /PCBM− . Table 2 summarizes the ten lowest-energy electronic transitions in isomer I1. The formation of P3HT∗ /PCBM is given by the electronic transition with the highest oscillator strength value (f = 1.48), resulting in an excitation energy of 3.02 eV mainly originated from a HOMO-to-LUMO+4 transition; which also resembles the experimental maximum of 2.8 eV for pure P3HT films. 58 According to Tables 1 and 2, for isomer I1, there are two CS states,

8


Page 8 of 24

Page 9 of 24

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


Table 1: Orbital localization according to fragments in isomer I1 of the P3HT/PCBM dimer. Orbital HOMO-6 HOMO-5 HOMO-4 HOMO-3 HOMO-2 HOMO-1 HOMO LUMO LUMO+1 LUMO+2 LUMO+3 LUMO+4 LUMO+5 LUMO+6

Fragment PCBM PCBM PCBM PCBM P3HT P3HT P3HT PCBM PCBM PCBM PCBM P3HT PCBM PCBM

P3HT+ /PCBM− . The first one (2.75 eV), CS1, is formed via the first electronic transition because of the main HOMO-to-LUMO contribution. The second one, CS2, is formed via the HOMO-to-LUMO+1 electronic transition characterized by an excitation energy of 2.93 eV and also a vanishing oscillator strength. Thus, there are two channels for the dissociation of the exciton P3HT∗ /PCBM into P3HT+ /PCBM− , via the photoinduced states CS1 and CS2 with a favorable driving force ∆GCS1 = −0.27 eV and ∆GCS2 = −0.09 eV, respectively. This second lowest-energy state CS2 has been reported previously for the P3HT/PCBM dimer and a set of other donor-acceptor pairs. 59 However, we found these two CS states with favorable driving force only for isomer I1. For isomer I2 and dimer I3 only one CS state with favorable driving force was found considering the ten lowest-energy electronic transitions. 57 For dimers I4 and I5 the Gibbs energy is positive. In Table 3, the transitions leading to the charge-localized exciton P3HT∗ /PCBM and P3HT+ /PCBM− are also reported. In the former, the reduced value of the oscillator strength in the U-shaped isomers I1 and I2 is due to a loss of symmetry resulting in a diminuition of the transition dipole moment of symmetry-allowed transitions. Hence the intensity of the main absorption peaks in these isomers is decreased as compared to the respective peaks 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

Table 2: The ten lowest-energy electronic transitions for isomer I1. Excitation energies Eexc in eV, HOMO (H) to LUMO (L) orbital contribution along with its weight, and oscillator strength f . Contribution weights lower than 0.10 are ignored. Eexc 2.75 2.80 2.81 2.84 2.89 2.93 3.02 3.07 3.10 3.14

Transition H→L (0.34) H-3→L+1 (0.20) H-4→L (0.12) H-3→L (0.59) H→L (0.15) H-3→L+1 (0.39) H→L (0.18) H-3→L (0.16) H-5→L (0.49) H-4→L (0.31) H-4→L (0.31) H-5→L (0.28) H→L+1 (0.17) H→L+1 (0.46) H→L+2 (0.15) H→L+4 (0.34) H→L+5 (0.17) H-6→L+1 (0.45) H-4→L+1 (0.13) H→L+2 (0.13) H→L+2 (0.23) H-4→L+1 (0.22) H→L+1 (0.12) H-5→L+1 (0.32) H→L+2 (0.16) H-3→L+2 (0.14)

10

f 0.05 0.03 0.03 0.00 0.04 0.01 1.48 0.18 0.10 0.02


Page 10 of 24

Page 11 of 24

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


in the planar conformation. The intense peaks at 2.85 eV in the planar dimers I3 to I5 are in better agreement with the experimental peak of 2.8 eV of a pure P3HT film; as expected in view of the lack of distortions in the chain of P3HT. CS states are not accessed by light absorption (f = 0) but by the migration of an electron to a phase boundary through spatial translation upon dissociation of the state P3HT∗ /PCBM. The most favored driving force, ∆GCS1 = −0.27 eV, is calculated for isomer I1; even though I2 and I3 deviate by no more than 0.16 eV. Disfavored driving forces correspond to the dimers I4 and I5 which is attributed to the lack of interactions between frontier molecular orbitals; thus, CS states should not be modeled via isolated fragments and vdW interactions must be considered. Table 3 summarizes all the electron transfer parameters that influence kCS . Notice that the fastest CS reaction, which can be computed as the sum of the contribution of the two channels, 60 kCS−I1 = kCS1 + kCS2 = 5.97 × 1013 + 2.48 × 1013 = 8.46 × 1013 s−1 is calculated for the U-shaped isomer I1. This is up to two and three orders of magnitude higher than kCS in I2 (3.0 × 1011 s−1 ) and I3 (7.0 × 1010 s−1 ), respectively. The slow CS in dimers I4 and I5 is evidently due to their poor Vij and endergonic ∆GCS . Thus, there is an increase of kCS of 3 (1) orders of magnitude for isomers I1 (I2) with respect to I3 for which dispersion forces were not taken into consideration, although I3 does have a similar intermolecular distance (D = 3.5 Å). There is a difference of 6 orders of magnitude between the extreme cases, I1 and I5. Actually, notice that kCS−I1 = 8.46 × 1013 s−1 is the only value that reaches the experimental range (> 4 × 1012 s−1 ), 18,19 which is for the most stable isomer I1 for which vdW forces were taken into account. This underlines the necessity to work with realistic dimer structures, which are predicted when vdW forces are properly taken into account. In order to further analyze the cause of such high increase of kCS , let us focus on the main parameters that influence its value (λint , ∆GCS , and Vij ). For the internal reorganization energy λint (Equations 2-3 ) we obtained a value of 0.245 eV, which resembles the value of 0.27 eV obtained by spectroscopic studies 61 and the typical small value that characterizes systems involving fullerenes. 10,11,21 The variation of

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 24

Table 3: Electronic couplings, Vij ; excitation energies, Eexc ; main HOMO (H) to LUMO (L) orbital contributions; oscillator strength values f ; and Gibbs energy for charge separation, ∆GCS , for each dimer under study. All energy terms in eV. The rate constant kCS is also provided in 1/s calculated with λint = 0.245 eV. Isomer I1

P3HT∗ /PCBM Eexc Transition f 0.0427 3.020 H→L+4 1.48

I2 I3 I4 I5

0.0031 0.0020 0.0004 0.0003

P3HT+ /PCBM− ∆GCS Eexc Transition f 2.745 H→L 0.05 -0.274 2.926 H→L+1 0.01 -0.094 3.062 H→L 0.00 -0.209 2.730 H→L+1 0.00 -0.116 2.858 H→L+1 0.00 0.008 2.977 H→L 0.00 0.126

Vij

3.271 2.846 2.850 2.852

H→L+5 H→L+5 H→L+5 H→L+5

0.75 3.00 3.13 3.17

kCS 5.97 × 1013 2.48 × 1013 3.02 × 1011 7.00 × 1010 3.68 × 108 1.30 × 107

λint over the five dimers was very small as expected due to the structural stability of C60 and its anion. For the P3HT chain, we found that the reorganization energy in the U-shaped and planar P3HT varies in the range of 0.159 to 0.163 eV, thus suggesting that λint remains constant through the five dimers, ∆λint = 4.1 meV. In view of that, considering a fixed value of λint = 0.245 eV for any dimer is a reasonable approximation, which in turn “eliminates” one independent variable for determining the variation of kCS . From Figure 2 (top), we can see an almost perfect linear correlation (R2 = 0.98) of ∆GCS for the five P3HT/PCBM dimers. Notice that this correlation is not just ∆GCS versus D, like previously reported, 21,62 since isomers I1, I2, and dimer I3 can be considered to approximately have the same average intermolecular distance (D ≈ 3.0 − 3.5 Å). What really makes the difference between these isomers (I1 and I2) is the drastic effect that vdW interactions produce on the dimer structures and, as a consequence, on the electronic and CS properties. Although there is a favorable driving force for I3 (∆GCS = −0.116 eV), which is not the case for the other two dimers with “idealistic” structures (I4 and I5), this driving force is increased by 100% (∆GCS−I1 = −0.274 eV, ∆GCS−I2 = −0.209 eV) for the isomers with ”realistic” structures (I1 and I2). In conclusion, there is not only an effect of the intermolecular distance on ∆GCS as previously remarked 62 but also an effect of the molecular orbital relaxation when the P3HT chain is bent in the equilibrium structures of isomers I1 and I2. An increase in the CS driving force directly results in a favorable change in kCS (see Equation 1). Our data show, 12


Page 13 of 24

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


however, that this favorable change is not the main factor that impacts the large increase in kCS . Changes in the electronic coupling have a significantly larger impact, while differences in ∆GCS result in an increase of kCS of less than one order of magnitude. The electronic coupling Vij = 0.002 eV (see Table 3) calculated for the planar dimer I3 is in qualitative agreement with the previously computed value of 0.013 eV; this latter was determined through a direct approach assuming that φLUMO in each fragment represents the nonadiabatic states, hφi |F |φj i, where F is the Kohn-Sham-Fock operator of the total system evaluated with the density matrix of the isolated fragments. 26 Increasing the distance between the fragments exponentially decreases the value of Vij , as described elsewhere 21 and confirmed by observing Vij from I3 to I5. However, Vij = 0.0427 eV for the U-shaped isomer I1 is one order of magnitude greater than Vij for I2 and two orders of magnitude greater than the one for dimers I3 to I5. After analyzing the quantities that determine Vij , 51–53,57 we found that this high increase in Vij is because the charge transfer matrix element for I1 is one order of magnitude greater than the one for I2 and I3, and two orders of magnitude greater than the one for I4 and I5. Even though I2 has a structure that includes effects of vdW interactions (see Figure 1, top), the electronic coupling is probably reduced in this isomer because of its local translational symmetry breaking related with the head-tail arrangement in the thiophene chain (isomer I1 does not show such a symmetry breaking). Notice that all the computed values of Vij (see Table 3) are sufficiently small so that the nonadiabatic approach, Equation 1, can be applied to assess kCS . Overall, there is an exponential increase of Vij with the five dimers. Notice again, like for ∆GCS , that this trend is not just Vij vs D because I1, I2, and I3 have on average a similar intermolecular distance D = 3.0 − 3.5 Å. Figure 2 (bottom) shows that there is a strong exponential correlation between ∆GCS and Vij (R2 = 0.84). This interesting correlation can be used for estimating Vij if the values of ∆GCS are known (or vice versa).

13



0.20

∆GCS (eV)

0.10 0.00 −0.10 −0.20 −0.30 −0.40 1

2 3 4 5 P3HT/PCBM dimer isomer

10−1 Vij = 7 × 10−4 × e−11.34 ∆GCS Vij (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 14 of 24

R2 = 0.8435

10−2

10−3

10−4 −0.3

−0.2

−0.1 0 ∆GCS (eV)

0.1

0.2

Figure 2: (top) Gibbs energy for photoinduced charge separation, ∆GCS , as a function of the dimer number as defined in 1. (bottom) Correlation between the electronic coupling for charge separation Vij (logarithmic scale) and ∆GCS . Trend lines and their corresponding equations are shown in the inset.

14


Page 15 of 24

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


Conclusions In this letter, the effects of the vdW interactions on the structure and, as a consequence, on the CS rate constant of the P3HT/PCBM dimer have been quantified. For the P3HT/PCBM isomers with realistic structures obtained by properly modeling the vdW interactions, I1 and I2, the CS reaction is predicted to be up to 3 (6) orders of magnitude faster than for P3HT/PCBM dimer structures with intermolecular distances of 3.5 Å (5.5 Å) in which vdW forces are not explicitly accounted for (I3 to I5). The fastest CS reaction is predicted for the vdW modified isomer I1, which is the only isomer whose CS rate constant value (kCS−I1 = 8.46 × 1013 s−1 ) falls into the experimentally measured range (kCS > 4 × 1012 s−1 ) for P3HT/PCBM active OPV layers. Small variations in the Gibbs energy for CS due to vdW effects, ∆GCS ≈ 0.2 eV, in the U-shaped as compared to the planar conformation do not have a significant impact on kCS , not affecting its order of magnitude. The electronic coupling Vij in isomers I1 and I2 increases by up to two and one orders of magnitude, respectively, with respect to the corresponding values for I3 to I5. This larger value of the electronic coupling is the main factor that leads to an increase of kCS by up to 3 and 6 orders of magnitude, respectively, with respect to the corresponding values of dimers I3 and I5. Therefore, we conclude that faster photoinduced CS is attained by states that are properly coupled electronically, which is the case here for the U-shaped isomers I1 and I2. These conclusions give quantitative evidence that vdW interactions should be taken into account for the proper quantum physical modeling of CS reactions, and thus computational predictions of the efficiency of BHJ organic solar cells.

Acknowledgement We acknowledge financial support from UC-MEXUS-CONACYT (Project 15-1462). The excellent service by the Dutch National Supercomputer (SURFsara Services), as well as the San Diego Supercomputer Center, is gratefully acknowledged. This work used the Ex15



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

treme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575 (award TG-CHE130010 to AWG). JIR thanks support from SIP-IPN (project 20171617).

Supporting Information Available The following files are available free of charge. • SI1: Orbital localization. Excitation Energies. Electronic coupling integrals. This material is available free of charge via the Internet at http://pubs.acs.org/.

References (1) Jariwala, D.; Sangwan, V. K.; Lauhon, L. J.; Marks, T. J.; Hersam, M. C. Carbon nanomaterials for electronics, optoelectronics, photovoltaics, and sensing. Chem. Soc. Rev. 2013, 42, 2824–2860. (2) Yan, J.; Saunders, B. R. Third-generation solar cells: a review and comparison of polymer:fullerene, hybrid polymer and perovskite solar cells. RSC Adv. 2014, 4, 43286– 43314. (3) Chen, G.; Seo, J.; Yang, C.; Prasad, P. N. Nanochemistry and nanomaterials for photovoltaics. Chem. Soc. Rev. 2013, 42, 8304–8338. (4) Delgado, J. L.; Bouit, P.-A.; Filippone, S.; Herranz, M. A.; Martín, N. Organic photovoltaics: a chemical approach. ACS Chem. Comm. 2010, 46, 4853–4865. (5) Mauter, M. S.; Elimelech, M. Environmental Applications of Carbon-Based Nanomaterials. Environ. Sci. Technol. 2008, 42, 5843–5859.

16


Page 16 of 24

Page 17 of 24

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


(6) Kumar, P.; Chand, S. Recent progress and future aspects of organic solar cells. Prog. Photovolt.: Res. Appl. 2012, 20, 377–415. (7) Darling, S. B.; You, F. The case for organic photovoltaics. RSC Adv. 2013, 3, 17633– 17648. (8) Illescas, B. M.; Martín, N. [60]Fullerene Adducts with Improved Electron Acceptor Properties. J. Org. Chem. 2000, 65, 5986–5995. (9) Martín, N.;

Sánchez, L.;

Illescas, B.;

Pérez, I. C60 -Based Electroactive

Organofullerenes. Chem. Rev. 1998, 98, 2527–2548. (10) Kawashima, Y.; Ohkubo, K.; Fukuzumi, S. Small reorganization energies of photoinduced electron transfer between spherical fullerenes. J. Phys. Chem. A 2013, 117, 6737–6743. (11) Ihly, R.; Mistry, K. S.; Ferguson, A. J.; Clikeman, T. T.; Larson, B. W.; Reid, O.; Boltalina, O. V.; Strauss, S. H.; Rumbles, G.; Blackburn, J. L. Tuning the driving force for exciton dissociation in single-walled carbon nanotube heterojunctions. Nat. Chem. 2016, 8, 603–609. (12) Hauch, J. A.; Schilinsky, P.; Choulis, S. A.; Childers, R.; Biele, M.; Brabec, C. J. Flexible organic P3HT:PCBM bulk-heterojunction modules with more than 1 year outdoor lifetime. Sol. Energy Mat. Sol. Cells 2008, 92, 727–731. (13) Scharber, M. C.; Mühlbacher, D.; Koppe, M.; Denk, P.; Waldauf, C.; Heeger, A. J.; Brabec, C. J. Design Rules for Donors in Bulk-Heterojunction Solar Cells-Towards 10 % Energy-Conversion Efficiency. Adv. Mater. 2006, 18, 789–794. (14) Dennler, G.; Forberich, K.; Scharber, M. C.; Brabec, C. J.; Tomiš, I.; Hingerl, K.; Fromherz, T. Angle dependence of external and internal quantum efficiencies in bulkheterojunction organic solar cells. J. Appl. Phys. 2007, 102, 054516. 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

(15) Guo, J.; Ohkita, H.; Benten, H.; Ito, S. Charge generation and recombination dynamics in poly(3-hexylthiophene)/fullerene blend films with different regioregularities and morphologies. J. Am. Chem. Soc. 2010, 132, 6154–6164. (16) Ohkita, H.; Cook, S.; Astuti, Y.; Duffy, W.; Tierney, S.; Zhang, W.; Heeney, M.; McCulloch, I.; Nelson, J.; Bradley, D. D. C. et al. Charge carrier formation in polythiophene/fullerene blend films studied by transient absorption spectroscopy. J. Am. Chem. Soc. 2008, 130, 3030–3042. (17) Hwang, I.-W.; Moses, D.; Heeger, A. J. Photoinduced Carrier Generation in P3HT/PCBM Bulk Heterojunction Materials. J. Phys. Chem. C 2008, 112, 4350– 4354. (18) Choi, J. H.; Son, K.-I.; Kim, T.; Kim, K.; Ohkubo, K.; Fukuzumi, S. Thienylsubstituted methanofullerene derivatives for organic photovoltaic cells. J. Mater. Chem. 2010, 20, 475–482. (19) Piris, J.; Dykstra, T. E.; Bakulin, A. A.; van Loosdrecht, P. H.; Knulst, W.; Trinh, M. T.; Schins, J. M.; Siebbeles, L. D. Photogeneration and Ultrafast Dynamics of Excitons and Charges in P3HT/PCBM Blends. J. Phys. Chem. C 2009, 113, 14500–14506. (20) Yi, Y.; Coropceanu, V.; Brédas, J.-L. Exciton-Dissociation and Charge-Recombination Processes in Pentacene/C 60 Solar Cells: Theoretical Insight into the Impact of Interface Geometry. J. Am. Chem. Soc. 2009, 131, 15777–15783. (21) Rudolf, M.; Kirner, S. V.; Guldi, D. M. A multicomponent molecular approach to artificial photosynthesis - the role of fullerenes and endohedral metallofullerenes. Chem. Soc. Rev. 2016, 45, 612–630. (22) Pelzer, K. M.; Darling, S. B. Charge generation in organic photovoltaics: a review of theory and computation. Mol. Syst. Des. Eng. 2016, 1, 10–24. 18


Page 18 of 24

Page 19 of 24

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


(23) Vella, E.; Li, H.; Grégoire, P.; Tuladhar, S. M.; Vezie, M. S.; Few, S.; Bazán, C. M.; Nelson, J.; Silva-acuña, C.; Bittner, E. R. Ultrafast decoherence dynamics govern photocarrier generation efficiencies in polymer solar cells. Sci. Rep. 2016, 6, 1–12. (24) Marcus, R. A. Electron transfer reactions in chemistry. Theory and experiment. Rev. Mod. Phys. 1993, 65, 599–610. (25) Liu, T.; Cheung, D. L.; Troisi, A. Structural variability and dynamics of the P3HT/PCBM interface and its effects on the electronic structure and the chargetransfer rates in solar cells. Phys. Chem. Chem. Phys. 2011, 13, 21461. (26) Liu, T.; Troisi, A. Absolute Rate of Charge Separation and Recombination in a Molecular Model of the P3HT/PCBM Interface. J. Phys. Chem. C 2011, 115, 2406–2415. (27) Woods, L. M.; Dalvit, D. A. R.; Tkatchenko, A.; Rodriguez-Lopez, P.; Rodriguez, A. W.; Podgornik, R. Materials perspective on Casimir and van der Waals interactions. Rev. Mod. Phys. 2016, 88, 045003. (28) Hermann, J.; DiStasio, R. A.; Tkatchenko, A. First-Principles Models for van der Waals Interactions in Molecules and Materials: Concepts, Theory, and Applications. Chem. Rev. 2017, 117, 4714–4758. (29) DiStasio, R. A.; Gobre, V. V.; Tkatchenko, A. Many-body van der Waals interactions in molecules and condensed matter. J. Phys. Condens. Matter 2014, 26, 213202. (30) Tkatchenko, A. Current Understanding of Van der Waals Effects in Realistic Materials. Adv. Funct. Mater. 2015, 25, 2054–2061. (31) Grimme, S.; Antony, J.; Schwabe, T.; Mück-Lichtenfeld, C. Density functional theory with dispersion corrections for supramolecular structures, aggregates, and complexes of (bio)organic molecules. Org. Biomol. Chem. 2007, 5, 741–758.

19



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

(32) Gutiérrez-González, I.; Molina-Brito, B.; Götz, A. W.; Castillo-Alvarado, F.; Rodríguez, J. I. Structural and electronic properties of the P3HT-PCBM dimer: A theoretical Study. Chem. Phys. Lett. 2014, 612, 234–239. (33) Rodríguez, J. I.; Matta, C. F.; Uribe, E. A.; Götz, A. W.; Castillo-Alvarado, F.; MolinaBrito, B. A QTAIM topological analysis of the P3HT-PCBM dimer. Chem. Phys. Lett. 2016, 644, 157–162. (34) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. (35) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple [Phys. Rev. Lett. 77, 3865 (1996)]. Phys. Rev. Lett. 1997, 78, 1396–1396. (36) te Velde, G.; Bickelhaupt, F. M.; Baerends, E. J.; Fonseca Guerra, C.; van Gisbergen, S. J. A.; Snijders, J. G.; Ziegler, T. Chemistry with ADF. J. Comput. Chem. 2001, 22, 931–967. (37) van Gisbergen, S.; Snijders, J.; Baerends, E. Implementation of time-dependent density functional response equations. Comput. Phys. Commun. 1999, 118, 119–138. (38) Fonseca Guerra, C.; Snijders, J. G.; te Velde, G.; Baerends, E. J. Towards an order- N DFT method. Theoret. Chem. Acc. 1998, 99, 391–403. (39) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (40) Yanai, T.; Tew, D. P.; Handy, N. C. A new hybrid exchange-correlation functional using the Coulomb-attenuating method (CAM-B3LYP). Chem. Phys. Lett. 2004, 393, 51–57.

20


Page 20 of 24

Page 21 of 24

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


(41) Seth, M.; Ziegler, T. Range-Separated Exchange Functionals with Slater-Type Functions. J. Chem. Theory Comput. 2012, 8, 901–907. (42) Peach, M. J. G.; Helgaker, T.; Satek, P.; Keal, T. W.; Lutnaes, O. B.; Tozer, D. J.; Handy, N. C. Assessment of a Coulomb-attenuated exchange-correlation energy functional. Phys. Chem. Chem. Phys. 2006, 8, 558–562. (43) Darling, S. B.; Sternberg, M. Importance of Side Chains and Backbone Length in Defect Modeling of Poly ( 3-alkylthiophenes ). J. Phys. Chem. B 2009, 113, 6215–6218. (44) Goeb, S.; De Nicola, A.; Ziessel, R. Oligomeric Ligands Incorporating Multiple 5,5’Diethynyl-2,2’-bipyridine Moieties Bridged and End-Capped by 3,4-Dibutylthiophene Units. J. Org. Chem. 2005, 70, 1518–1529. (45) The length of the oligomer has already been discussed in a previous work by Liu and Troisi, wherein the 6-unit oligomer was concluded to have the suitable length in CT reactions between P3HT and PCBM (see Ref. 26). The 8-unit oligomer was discarded because of an inappropriate HOMO-LUMO electronic structure, as determined via UDFT in the triplet-spin state. However, from a previous work, both singlet and triplet excited CT states were indeed conveniently determined via U-DFT by altering the ground-state electronic structure (see Ref. 48); as a result, triplet excited CT states can be in principle found in the 8-unit oligomer P3HT/PCBM dimer. We believe that the better performance of CAMY-B3LYP over hybrid functionals (used by Liu and Troisi) to describe excited states with strong CS character may be the reason of their misleading conclusion. (46) Nelsen, S. F.; Silas, B. C.; Kim, Y. Estimation of Inner Shell Marcus Terms for Amino Nitrogen Compounds by Molecular Orbital Calculations . J. Am. Chem. Soc. 1987, 109, 677–682.

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

(47) Norton, J. E.; Brédas, J.-L. Polarization energies in oligoacene semiconductor crystals. J. Am. Chem. Soc. 2008, 130, 12377–12384. (48) McMahon, D. P.; Troisi, A. Evaluation of the External Reorganization Energy of Polyacenes. J. Phys. Chem. Lett. 2010, 1, 941–946. (49) Lee, M. H.; Dunietz, B. D.; Geva, E. Calculation from First Principles of Intramolecular Golden-Rule Rate Constants for Photo-Induced Electron Transfer in Molecular DonorAcceptor Systems. J. Phys. Chem. C 2013, 117, 23391–23401. (50) Lee, M. H.; Geva, E.; Dunietz, B. D. Calculation from First-Principles of Golden Rule Rate Constants for Photoinduced Subphthalocyanine/Fullerene Interfacial Charge Transfer and Recombination in Organic Photovoltaic Cells. J. Phys. Chem. C 2014, 118, 9780–9789. (51) Senthilkumar, K.; Grozema, F. C.; Bickelhaupt, F. M.; Siebbeles, L. D. A. Charge transport in columnar stacked triphenylenes: Effects of conformational fluctuations on charge transfer integrals and site energies. J. Chem. Phys. 2003, 119, 9809–9817. (52) Senthilkumar, K.; Grozema, F. C.; Fonseca Guerra, C.; Bickelhaupt, F. M.; Lewis, F. D.; Berlin, Y. A.; Ratner, M. A.; Siebbeles, L. D. A. Absolute Rates of Hole Transfer in DNA. J. Am. Chem. Soc. 2005, 127, 14894–14903. (53) Coropceanu, V.; Cornil, J.; da Silva Filho, D. A.; Olivier, Y.; Silbey, R.; Brédas, J.-L. Charge transport in organic semiconductors. Chem. Rev. 2007, 107, 926–952. (54) Martínez, J. P.; Solà, M.; Voityuk, A. A. Theoretical estimation of the rate of photoinduced charge transfer reactions in triphenylamine C60 donor-acceptor conjugate. J. Comp. Chem. 2016, 37, 1396–1405. (55) Voityuk, A. A. Electronic coupling for charge transfer in donor-bridge-acceptor systems.

22


Page 22 of 24

Page 23 of 24

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


Performance of the two-state FCD model. Phys. Chem. Chem. Phys. 2012, 14, 13789– 13793. (56) Hirata, S.; Head-Gordon, M. Time-dependent density functional theory within the Tamm-Dancoff approximation. Chem. Phys. Lett. 1999, 314, 291–299. (57) See Supporting Information. (58) Böckmann, M.; Schemme, T.; de Jong, D. H.; Denz, C.; Heuer, A.; Doltsinis, N. L. Structure of P3HT crystals, thin films, and solutions by UV/Vis spectral analysis. Phys. Chem. Chem. Phys. 2015, 17, 28616–28625. (59) Peet, J.; Kim, J. Y.; Coates, N. E.; Ma, W. L.; Moses, D.; Heeger, A. J.; Bazan, G. C. Efficiency enhancement in low-bandgap polymer solar cells by processing with alkane dithiols. Nat. Mater. 2007, 6, 497–500. (60) Liu, T.; Troisi, A. What Makes Fullerene Acceptors Special as Electron Acceptors in Organic Solar Cells and How to Replace Them. Adv. Mater. 2013, 25, 1038–1041. (61) Vandewal, K.; Tvingstedt, K.; Gadisa, A.; Inganäs, O.; Manca, J. V. Relating the opencircuit voltage to interface molecular properties of donor:acceptor bulk heterojunction solar cells. Phys. Rev. B 2010, 81, 125204. (62) Ma, H.; Troisi, A. Modulating the Exciton Dissociation Rate by up to More than Two Orders of Magnitude by Controlling the Alignment of LUMO + 1 in Organic Photovoltaics. J. Phys. Chem. C 2014, 118, 27272–27280.

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 ACS Paragon Plus Environment

Page 24 of 24

Effects of Dispersion Forces on Structure and Photoinduced Charge

Recommend Documents