Charge Transport in Nanostructured Materials: Implementation and

Apr 20, 2017 - The in silico design of novel complex materials for energy conversion requires accurate, ab initio simulation of charge transport. In t...
0 downloads 0 Views 2MB Size
Article pubs.acs.org/JCTC

Charge Transport in Nanostructured Materials: Implementation and Verification of Constrained Density Functional Theory Matthew B. Goldey,*,† Nicholas P. Brawand,† Márton Vörös,‡ and Giulia Galli†,‡ †

Institute for Molecule Engineering, University of Chicago, Chicago, Illinois 60637, United States Argonne National Laboratory, Lemont, Illinois 60439, United States

Downloaded via UNIV OF SOUTH DAKOTA on September 13, 2018 at 13:09:12 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.



ABSTRACT: The in silico design of novel complex materials for energy conversion requires accurate, ab initio simulation of charge transport. In this work, we present an implementation of constrained density functional theory (CDFT) for the calculation of parameters for charge transport in the hopping regime. We verify our implementation against literature results for molecular systems, and we discuss the dependence of results on numerical parameters and the choice of localization potentials. In addition, we compare CDFT results with those of other commonly used methods for simulating charge transport between nanoscale building blocks. We show that some of these methods give unphysical results for thermally disordered configurations, while CDFT proves to be a viable and robust approach.

1. INTRODUCTION Functional nanostructured systems have long been recognized as a promising platform for the development of materials with target properties, and they are an active area of theoretical and experimental research.1−8 Much progress has been made in recent years to predict structural and electronic properties of functional materials for energy conversion applications;9−17 however, robust, ab initio methods for the description and optimization of electronic transport properties of nanostructured materials are not yet available, although several useful codes18−27 for electronic transport have been developed, based on various models and approximations. Much effort is still required, anchored within first-principles theories, to describe electron and hole transport in materials for energy applications in order to define descriptors for the design of semiconductor nanostructures, e.g., with optimal solar and thermal energy conversion properties. Here we focus on hopping transport and the efficient calculation of first-principles parameters required either to compute transfer rates using Marcus theory28 or to evaluate currents from Green’s function approaches, e.g., using the Keldysh formalism.29,30 In the framework of Marcus theory, charge transport rates (k) are expressed in terms of free energy differences (ΔG) between sites A and B involved in charge hopping, the electronic coupling between these sites (HAB), and the reorganization energy (λ): kAB =

⎡ (ΔG + λ)2 /4λ ⎤ 2π 1 |HAB|2 exp⎢ − ⎥ ℏ 4kBTπλ kBT ⎦ ⎣

In eq 1, sites A and B may represent two nanoparticles of an ensemble (e.g., a nanoparticle-based inorganic solar cell), or two different cages in a clathrate solid or a nanoparticle and a region of space in the surrounding matrix. Depending on the specific nanostructured system and problem and on the strength of the HAB coupling, one may consider using transfer rate expressions derived from Landau−Zener theory33 which include information about the phonon density of states of the nanoparticles; e.g., ⎡ ⎛ − 2π 3/2H 2 ⎞⎤ ω ⎛ ⎛λ ⎞ ⎞ AB ⎟⎥ k = ⎢1 − exp⎜⎜ ⎟⎥ 2π exp⎜⎝− ⎜⎝ 4 − |HAB|⎟⎠ /kBT ⎟⎠ ⎢⎣ ⎝ ℏω λkBT ⎠⎦ (2)

where ω is a representative molecular or nanoparticle phonon frequency. In the cases of Marcus theory and the Landau−Zener formula, the electronic coupling can be obtained from first principles and in this work we focus on constrained density functional theory (CDFT) to obtain this quantity. The use of CDFT to obtain couplings to evaluate rates within Marcus theory has been reviewed, e.g., by Kaduk et al.36 In addition to CDFT, many methods have been proposed for predicting the electronic coupling, such as the generalized Mulliken−Hush (GMH) method,37 excited state wave function methods,38 tight-binding approaches,32,39,40 anticrossing methods,41 and empirical quasiclassical approximations.42 Many of these methods neglect polarization effects, and some of them do not contain any atomistic detail or orbital relaxation beyond the single particle picture,43 which may be important for materials with high densities of states near the conduction band

(1)

where kB is the Boltzmann constant and T is the temperature. For the details of Marcus theory, we refer readers to existing reviews.28,31−35 © 2017 American Chemical Society

Received: January 26, 2017 Published: April 20, 2017 2581

DOI: 10.1021/acs.jctc.7b00088 J. Chem. Theory Comput. 2017, 13, 2581−2590

Article

Journal of Chemical Theory and Computation

are determined self-consistently). Perhaps more importantly, CDFT involves ground state total energies (which DFT describes much better than single particle virtual states) and describes charge transfer between ground states of a given system in the presence of a nonlinear perturbation (the CDFT potential). We implemented CDFT in Quantum-ESPRESSO,73 using Hirshfeld partitioning74 of the electronic density to determine the real-space shape of the constraining potential. We chose Hirshfeld partitioning in order to compare directly with existing plane wave codes and, most importantly, because such partitioning has been shown to be more reliable in defining fragments (e.g., in the presence of possibly delocalized anionic states) than schemes not based on the density such as Löwdin’s and Mulliken’s. The weighting function used in Hirshfeld partitioning is defined by the selection of atoms donating (D) and accepting (A) electrons,

minimum or valence band maximum. Furthermore, none of these methods have been extended to include effects present in real devices arising from image charges within the surrounding electrodes and substrates. CDFT, on the other hand, includes the effects of polarization and orbital relaxation through the evaluation of many-body wave functions describing localized charges. These wave functions may be obtained within density functional theory (DFT) by modifying the Kohn−Sham equations to include a spatially dependent potential, whose strength is varied selfconsistently to satisfy the constraint of the desired number of electrons localized on a given fragment or building block.36,44−62 The use of constrained DFT within the quantum chemistry, and, since recently, the condensed matter physics, community is entering its fourth decade, since the first formulation in 1984 by Dederichs et al.,44 and it is often employed within the formulation proposed by Wu and Van Voorhis.46 In addition to total energies for localized states and electronic couplings, CDFT may be used to obtain charge transfer energies and Heisenberg spin couplings.36,63 This method has been applied to molecules,61,64−66 molecular solids,67,68 oxides,69 organic photovoltaics,70 solvated ions in liquid water,50,51 and defects in nanoparticle arrays,15,71 as well as to provide corrections to results where hybrid DFT fails to compute total energy differences and geometries accurately.62 The remainder of this work is organized as follows: we first present the details of our implementation of constrained DFT for periodic systems using plane wave basis sets. Second, we verify our implementation on molecular systems by comparing our results to existing literature data, and we include tests of the sensitivity of our implementation on numerical parameters. Finally, we present results for silicon quantum dots comparing constrained DFT against other commonly used methodologies for the calculation of electronic couplings, finding a robust performance for CDFT when other approaches lead to unphysical results for thermally disordered configurations.

D

w(r) =

∫Ω w(r) ρ(r) d3r − N0)

all

∑i ρi (r)

(4)

and is generated from the square moduli (ρi) of pseudoatomic wave functions, consistent with the chosen pseudopotentials. The sums in eq 4 run over pseudoatomic densities belonging to donors (D), acceptors (A), or all atoms (all). These definitions are general and may be applied to both the total electronic density and the spin density. CDFT yields the total energy and density of diabatic states; the calculation of the electronic coupling requires further approximations since many-body wave functions are not readily available from DFT approaches. Following Wu and Van Voorhis,48 we computed electronic couplings from the overlap of Kohn−Sham single particle wave functions, thus effectively employing a single Slater determinant built from Kohn−Sham orbitals to approximate the many-body wave function. In general, the approximation for the calculation of couplings must be tested numerically, e.g., by comparing with the results of systematically improvable, wave function-based methodologies.65,66 After generating two diabatic states (e.g., corresponding to a charge (electron or hole) localized on two different quantum dots A and B), the coupling matrix H′ is evaluated:

2. METHODOLOGY Within constrained DFT, an external potential is added to the self-consistent potential entering the Kohn−Sham equations, and its strength is varied self-consistently in order to localize a desired integer number of charges N0 on a specified site.36,46 Specifically, one defines a free energy F[ρ(r)](N), a functional of the density ρ(r), and the number of localized electrons N, belonging to a volume Ω, which is minimized with respect to the total charge density of the system, under a constraint (expressed as a Lagrange multiplier) of charge localization: F[ρ(r)](N ) = E[ρ(r)] + V (

A

∑i ρi (r) − ∑i ρi (r)

H′AB =

FA + FB SAB − 2

ΦA

VÂ + VB̂ ΦB 2

(5)

Here F is the free energy defined in eq 3, SAB is the overlap between the two Slater determinants (⟨ΦA | and ⟨ΦB |) constructed from the Kohn−Sham single particle energies from diabatic states A and B, and V̂ = V w (r) (see eq 3). The electronic coupling (HAB) is obtained from eq 5 using Löwdin’s orthogonalization:75

(3)

Here E is the Kohn−Sham total energy, ρ is the charge density, and w(r) is a spatially dependent potential of strength V. Given a trial potential strength V, the solution of the Kohn− Sham equations yields a density ρ and a corresponding number of localized electrons N = ∫ Ω w(r) ρ(r) d3r. V is then updated and the procedure repeated until the error in localized charges ΔN = N − N0 is within a chosen, small threshold. As shown by Kaduk et al.,36 CDFT charge transfer states give the correct 1/R asymptotic distance dependence and the correct asymptotic limit. These states can be correctly described even with semilocal functionals for which, for example, timedependent DFT fails,72 since CDFT involves nonlinear response of the density (the CDFT potential and the density

H = S −1/2H′S −1/2

(6)

When investigating hopping transport in ensembles of quantum dots representing a given device, or in any other system of nanostructured building blocks, one is usually faced with the problem of including the effect of metallic leads and/ or substrates. The derivation of methods to calculate the renormalization of molecular electronic levels at metallic interfaces is an area of active research in condensed matter physics and physical chemistry.76−78 Methods based on classical point charges neglect orbital relaxation, dipole interactions, and orbital reordering and are applied as post hoc corrections to 2582

DOI: 10.1021/acs.jctc.7b00088 J. Chem. Theory Comput. 2017, 13, 2581−2590

Article

Journal of Chemical Theory and Computation orbital energies. Other methods require an atomistic simulation of the metallic surface or substrate and thus provide much greater detail; however, they are computationally demanding because they require that the surface be simulated along with the quantum dot or molecular ensemble being adsorbed.76−78 In addition, these methods often require many-body perturbation theory to accurately model the relative position of single particle energy levels, and hence their applicability has been limited to systems with a small number of atoms, for example a single benzene molecule adsorbed on graphite.78 In order to include self-consistently image charge effects in CDFT calculations, we added an image charge potential Vimg. to the Kohn−Sham potential: Vimg[ρ ; r] = −

qtot |r|



ptot ⃗ [ρ]·r r3

+ ...

Figure 1. Electronic coupling HAB (eV) as a function of distance R (Å) for hole transfer between two Zn atoms belonging to a Zn2+ cation, obtained with the generalized Mulliken−Hush method37 (denoted GMH), with constrained density functional theory using Becke weight partitioning (Wu and Van Voorhis),48,85 and in this work using Hirshfeld partitioning74 (see text). Calculated β decay rates are 2.55, 2.69, and 2.72 Å−1 for the three methods, respectively.

(7)

where r is centered on the site containing the constrained charge defined by CDFT, qtot is the total charge of the system, and p⃗tot is the dipole moment of the system. Vimg is a multipole expansion of the electrostatic potential of the image, which is added to the Kohn−Sham potential. In our procedure, the charge obtained at the end of each self-consistent Kohn−Sham iteration is used to update Vimg[ρ;r] until the image interaction energy, E img =

∫Ω Vimg[ρ; r] ρ(r) dr

(8)

is converged within a predefined accuracy.15 We now turn to presenting calculations carried out with the procedure outlined above. These calculations were performed with the Quantum-ESPRESSO package,73 using plane wave basis sets and optimized norm-conserving Vanderbilt (ONCV) pseudopotentials79,80 and Troullier−Martins pseudopotentials.81 We used a wave function energy cutoff of 80 Ry, except for molecular dynamics and CDFT calculations of systems containing only Si and H, where a cutoff of 30 Ry was used. Our calculations employed the generalized gradient approximation with the PBE exchange-correlation functional.82 The Martyna−Tuckerman correction83 was applied to all CDFT calculations of molecular systems to remove any remaining fictitious electrostatic interactions between periodic images. We first verified our results for simple molecules using literature reference data. We then compared our calculations on complex systems with those of other approaches commonly used at zero temperature, including quantum dots at finite temperature.

Figure 2. Dependence of the energy and potential strength on the number of charges transferred N = ∫ Ω w(r) ρ(r) d3r (top) and % error

(

in electronic coupling 100% ×

HAB(N ) − HAB(N0) HAB(N0)

) versus the error in

the number of charges transferred (ΔN = N − N0), relative to results obtained with N0 (bottom) for Zn2+ at 5 Å separation (the equilibrium distance of Zn+2 within PBE is 2.55 Å).

3. VERIFICATION As pointed out by Lejaeghere et al.,84 the consistent reproduction of benchmark values is necessary for confidence in any new code and is a sign of the technological maturity of computational methodologies. In this section, we present first the electronic coupling for hole transfer in a simple dimer of Zn, and then we investigate the molecular data sets presented by Kubas et al.65,66 Results for Zn2+ are shown in Figure 1. We found excellent agreement with the results of Wu and Van Voorhis48 in absolute magnitude of hole transfer couplings. Given the exponential decay of wave functions, the dependence of electronic couplings on distance is commonly parametrized to the form HAB ∝ e−βR /2 to provide a characteristic decay rate for a charge transfer process.

Figure 3. Electronic couplings (HAB; see eq 6) calculated using constrained density functional theory (CDFT) implemented in Quantum-ESPRESSO (QE) versus values obtained with the CPMD code86 (see Kubas et al.65,66) for the HAB11 and HAB7 sets of πconjugated molecules.

Figure 2 presents the total energy relative to the ground state energy of unmodified DFT and the applied potential strength as a function of the number of electrons transferred for Zn2+ at 5 Å of separation, as well as the error in HAB versus the error in 2583

DOI: 10.1021/acs.jctc.7b00088 J. Chem. Theory Comput. 2017, 13, 2581−2590

Article

Journal of Chemical Theory and Computation Table 1. Electronic Couplings (HAB/meV) Calculated Using CDFT for the HAB1165 and HAB766 Setsa system

separation/Å

QE/ONCV

QE/TM

CPMD

anthracene anthracene anthracene anthracene cyclopentadiene cyclopentadiene cyclopentadiene cyclopentadiene cyclopropene cyclopropene cyclopropene cyclopropene ethylene ethylene ethylene ethylene furane furane furane furane imidazole imidazole imidazole imidazole pentacene pentacene pentacene pentacene perfluoroanthracene perfluoroanthracene perfluoroanthracene perfluoroanthracene perylene perylene perylene perylene perylene diimide perylene diimide perylene diimide perylene diimide phenol phenol phenol phenol porphin porphin porphin porphin pyrrole pyrrole pyrrole pyrrole tetracene tetracene tetracene tetracene thiophene thiophene thiophene thiophene

3.50 4.00 4.50 5.00 3.50 4.00 4.50 5.00 3.50 4.00 4.50 5.00 3.50 4.00 4.50 5.00 3.50 4.00 4.50 5.00 3.50 4.00 4.50 5.00 3.50 4.00 4.50 5.00 3.50 4.00 4.50 5.00 3.50 4.00 4.50 5.00 3.50 4.00 4.50 5.00 3.50 4.00 4.50 5.00 3.50 4.00 4.50 5.00 3.50 4.00 4.50 5.00 3.50 4.00 4.50 5.00 3.50 4.00 4.50 5.00

651.07 329.15 171.49 91.60 699.45 335.82 159.32 75.03 810.11 358.90 157.57 69.42 611.67 306.29 153.01 N/A 594.68 285.25 136.33 65.36 586.23 278.88 132.22 63.01 645.97 309.60 156.82 80.19 479.46 219.90 99.96 45.01 654.05 332.83 178.72 100.26 559.38 266.88 132.64 67.30 557.29 264.94 123.52 56.87 590.25 288.54 146.85 76.17 629.02 309.13 151.59 74.56 644.67 318.79 163.00 84.36 676.13 328.24 158.10 75.61

648.78 327.96 170.80 90.98 696.57 334.03 158.06 74.04 806.87 357.16 156.58 68.87 607.08 302.77 150.12 82.97 590.53 281.40 132.82 62.23 581.57 274.77 128.69 59.99 N/A 308.54 156.12 79.74 478.35 219.41 99.76 44.76 651.42 331.43 177.86 99.59 557.01 265.55 131.94 66.83 552.43 260.46 119.44 53.55 588.03 287.50 146.65 75.71 623.59 304.20 147.11 70.59 642.28 317.65 162.24 83.84 672.15 323.50 153.35 71.45

637.00 324.50 169.40 87.90 702.10 346.40 167.40 80.90 816.10 369.60 164.80 73.90 621.50 314.80 158.30 78.40 598.20 292.50 141.00 67.00 590.10 286.30 136.80 64.30 618.30 303.00 154.00 77.70 479.60 227.20 107.90 49.60 633.70 324.50 174.00 94.40 541.80 261.00 131.20 65.10 557.30 271.30 129.40 58.90 577.70 285.00 146.60 74.70 629.80 314.80 155.50 75.70 628.80 313.10 160.80 81.50 669.90 332.30 162.40 77.40

2584

CP2K-BW

CP2K-BW+A

CP2K-FBB+A

653.07 304.77 138.78 65.31 772.80 334.70 146.94 68.03 598.65 304.77 152.38 73.47

647.63 299.33 136.06 65.31 770.08 331.98 146.94 68.03 601.37 304.77 152.38 73.47

745.59 386.40 190.48 92.52 900.70 440.82 209.53 92.52 702.05 361.91 174.15 78.91

541.51 247.62 117.01 54.42

541.51 247.62 117.01 54.42

636.75 318.37 152.38 68.03

487.08 217.69 95.24 43.54

487.08 214.97 95.24 43.54

566.00 280.28 130.61 59.87

574.16 272.11 125.17 57.14

574.16 269.39 125.17 57.14

663.96 340.14 165.99 76.19

609.54 280.28 127.89 59.87

612.26 277.56 127.89 59.87

699.33 353.75 165.99 73.47

DOI: 10.1021/acs.jctc.7b00088 J. Chem. Theory Comput. 2017, 13, 2581−2590

Article

Journal of Chemical Theory and Computation Table 1. continued a

QE, CPMD, and CP2K denote the Quantum-ESPRESSO,89 CPMD,86 and CP2K61,87 codes, respectively. CP2K results were only available for the HAB11 set, and the molecules of the HAB7 set were not available. TM and ONCV denote the Troullier−Martins81 and optimized norm-conserving pseudopotentials,79,80 respectively (see text). BW, BW+A, and FBB+A denote the Becke weight partitioning, Becke weight partitioning with atomic size adjustments, and fragment-based Becke weights with atomic size adjustments, respectively (see Holmberg and Laasonen61). Two calculations carried out with QE and marked N/A did not converge without fractional occupations.

Table 4. Mean Error (ME/Å−1), Mean Absolute Error (MAE/Å−1), Root-Mean-Square Deviation (RMSD/Å−1), and Mean Absolute Percent Error (MAPE/%) for Decay Rates (β/Å−1) for Electronic Couplings Calculated Using CDFT Implemented in Quantum-ESPRESSO (QE) with Troullier−Martins (QE/TM) and Optimized NormConserving Pseudopotentials (QE/ONCV) and Those Obtained from the CPMD Code86 (see Kubas et al.65,66) for the HAB11 and HAB7 Setsa

Table 2. Mean Error (ME/meV), Mean Absolute Error (MAE/meV), Root-Mean-Square Deviation (RMSD/meV), and Mean Absolute Percent Error (MAPE/%) for Electronic Couplings (HAB) Calculated Using CDFT Implemented in Quantum-ESPRESSO (QE) for Troullier−Martins (QE/ TM) and Optimized Norm-Conserving Pseudopotentials (QE/ONCV) and Those Obtained from the CPMD Code86 (See Kubas et al.65,66) for the HAB11 and HAB7 Sets comparison

ME

MAE

RMSD

MAPE

QE/TM vs CPMD QE/ONCV vs CPMD QE/ONCV vs QE/TM

−2.49 0.20 2.34

6.83 5.94 2.34

7.93 7.79 2.84

3.53 2.68 1.30

comparison

ME

MAE

RMSD

MAPE

QE/TM vs CPMD QE/ONCV vs CPMD QE/ONCV vs QE/TM

0.040 0.029 −0.011

0.067 0.038 0.034

0.079 0.052 0.048

2.303 1.292 1.191

β values were computed by fitting the results of Table 1 to the functional form HAB = A exp−βR/2 on a log scale.

a

the number of charges transferred ΔN. The slope of the potential strength changes sharply in correspondence to an amount of charge transferred equal to one electron. A 1% deviation in the number of constrained electrons results in up to a 50 meV change in total energies. The same 1% deviation produces up to a 15% error in the calculated coupling, suggesting that tight convergence is necessary to achieve accurate CDFT couplings. Convergence within 1 × 10−4 electrons yields errors in the estimated couplings within 1%. On the basis of these results, we recommend a precision of 5 × 10−5 electrons for CDFT self-consistency, which we employ for our following calculations. Recent results employing a precision of 1 × 10−2 electrons may need to be revisited. In Figure 3 and Table 1, we present results for literature data sets (HAB7 and HAB1165,66) of electronic couplings between identical molecular dimers. The sets consist of organic

molecular dimers (varying in size from 8 to 80 atoms) at different separations (3.5, 4.0, 4.5, and 5.0 Å). Molecular dimer geometries for the HAB11 and the HAB7 sets, composed of planar and π-conjugated molecules, were taken from the literature.65,66 Further details about these molecules are given by Kubas et al.65,66 Comparisons between our results for Troullier−Martins (TM) pseudopotentials81 and optimized norm-conserving (ONCV) pseudopotentials,79,80 denoted QE/ TM and QE/ONCV, are presented in Tables 2, 3, and 4 and compared with those obtained with the Car−Parrinello molecular dynamics code CPMD86 and TM pseudopotentials by Kubas et al.65,66 Tables 1 and 3 also contain calculated values for the HAB11 set from Holmberg and Laasonen,61 which were

Table 3. Decay Rates (β/Å−1) for Electronic Couplings Calculated Using CDFT for the HAB1165 and HAB766 Sets of Molecular Dimersa system

QE/ONCV

QE/TM

CPMD

anthracene cyclopentadiene cyclopropene ethylene furane imidazole pentacene perfluoroanthracene perylene perylene diimide phenol porphin pyrrole tetracene thiophene

2.618 2.989 3.283 2.669 3.001 3.029 2.706 3.158 2.503 2.824 3.112 2.729 2.905 2.712 2.988

2.614 2.977 3.278 2.771 2.945 2.975 2.776 3.154 2.499 2.821 3.044 2.727 2.844 2.709 2.921

2.637 2.884 3.205 2.759 2.919 2.955 2.760 3.021 2.534 2.818 2.993 2.721 2.824 2.718 2.876

CP2K-BW

CP2K-BW+A

CP2K-FBB+A

3.078 3.245 2.795

3.068 3.238 2.800

2.787 3.028 2.915

3.057

3.057

2.978

3.228

3.223

3.001

3.079

3.075

2.885

3.099

3.100

3.007

a β values were computed by fitting the results of Table 1 to the functional form HAB = A exp−βR/2 on a log scale. QE, CPMD, and CP2K denote the codes Quantum-ESPRESSO,89 CPMD,86 and CP2K,87 respectively. CP2K results were only available for the HAB11 set. TM and ONCV denote the Troullier−Martins81 and optimized norm-conserving pseudopotentials,79,80 respectively (see text). BW, BW+A, and FBB+A denote the Becke weight partitioning, Becke weight partitioning with atomic size adjustments, and fragment-based Becke weights with atomic size adjustments, respectively (see Holmberg and Laasonen61).

2585

DOI: 10.1021/acs.jctc.7b00088 J. Chem. Theory Comput. 2017, 13, 2581−2590

Article

Journal of Chemical Theory and Computation Table 5. Root-Mean-Square Deviation (meV) for Electronic Couplings (HAB) between Different Implementations of Constrained DFT for Molecules of the HAB11 Database66a QE/TM QE/ONCV CPMD CP2K-BW CP2K-BW+A CP2K-FBB+A

QE/TM

QE/ONCV

CPMD

CP2K-BW

CP2K-BW+A

CP2K-FBB+A

0.0 3.7 8.6 29.92 30.95 40.67

3.7 0.0 5.75 33.58 34.61 38.7

8.6 5.75 0.0 36.43 37.52 33.53

29.92 33.58 36.43 0.0 2.04 64.83

30.95 34.61 37.52 2.04 0.0 65.9

40.67 38.7 33.53 64.83 65.9 0.0

a QE, CPMD, and CP2K denote the codes Quantum-ESPRESSO,89 CPMD,86 and CP2K,61,87 respectively. TM and ONCV denote the Troullier− Martins81 and optimized norm-conserving pseudopotentials,79,80 respectively (see text). BW, BW+A, and FBB+A denote the Becke weight partitioning, Becke weight partitioning with atomic size adjustments, and fragment-based Becke weights with atomic size adjustments, respectively (see Holmberg and Laasonen61).

Table 6. Mean Absolute Percent Error (%) for Electronic Couplings (HAB) between Different Implementations of Constrained DFT for Molecules of the HAB11 Database66a QE/TM QE/ONCV CPMD CP2K-BW CP2K-BW+A CP2K-FBB+A

QE/TM

QE/ONCV

CPMD

CP2K-BW

CP2K-BW+A

CP2K-FBB+A

0.0 2.05 4.7 10.24 10.5 13.03

1.98 0.0 2.61 11.78 12.05 11.27

4.45 2.52 0.0 13.7 13.96 8.34

11.78 13.96 16.57 0.0 0.36 25.61

12.12 14.32 16.92 0.36 0.0 26.0

11.15 9.69 7.4 19.98 20.22 0.0

a QE, CPMD, and CP2K denote the codes Quantum-ESPRESSO,89 CPMD,86 and CP2K,61,87 respectively. TM and ONCV denote the Troullier− Martins81 and optimized norm-conserving pseudopotentials,79,80 respectively (see text). BW, BW+A, and FBB+A denote the Becke weight partitioning, Becke weight partitioning with atomic size adjustments, and fragment-based Becke weights with atomic size adjustments, respectively (see Holmberg and Laasonen61).

obtained using the CP2K code.87 For a few systems (acetylene, benzene, and cyclobutadiene), direct comparison could not be made with CPMD since the electronic couplings were calculated between different orbitals.65,88 We found excellent agreement between the CPMD and Quantum-ESPRESSO implementations; the mean absolute percent error (MAPE) in the couplings between the two codes is 2.7−3.5%, depending on the choice of pseudopotentials. We also found agreement for the computed β values (ME = 0.040 and 0.029 Å−1 for TM and ONCV results, correspondingly). The different implementations of Hirshfeld partitioning, including intrinsic thresholds and cutoff radii for the atomic densities, are likely responsible for the small differences in the results obtained with these two codes. In Tables 5 and 6, we present pairwise root-mean-square deviations (RMSDs) and mean absolute percent errors (MAPEs) for the molecules in the HAB11 database65 as computed with our implementation, that of CPMD,51,86 and the recent implementation in CP2K.61,87 Our and the CPMD results are within 10 meV in RMSD for this database, while differing by 30−40 meV from those of the CP2K results. We note that the choice of constraining potentials may produce noticeable variations in the computed couplings, with MAPEs of 10−20% and RMSDs up to 66 meV between different localization approaches. Having verified our methodology, we now turn to a case study of silicon quantum dots.

Figure 4. Semilog plot of the electronic coupling for holes (top) and electrons (bottom) versus distance R/Å−1 between two Si35H36 QDs for configurations taken from 300 K molecular dynamics trajectory. Electronic couplings are computed using the splitting in dimer (SID), anticrossing (AX), and constrained DFT (CDFT) methods (see text).

properties and few address the effects of thermal disorder on the electronic coupling between QDs. For example, at 0 K a perfectly symmetric Si35H36 QD has degenerate HOMO and LUMO energy levels. These degeneracies are lifted by disorder at finite T or by the interaction with nearby QDs. Realistic simulation of charge transport in QD solids requires the inclusion of these effects. In this section, we examine the effect of structural distortions induced by finite temperature on computed electronic couplings when using several different methods for estimating charge transfer parameters. Two methods often used for nanoscale materials are the splitting in dimer method and the

4. APPLICATION TO SILICON QUANTUM DOTS The electronic coupling as well as other electron transfer properties such as driving force and reorganization energy are paramount to understanding hopping transport in QD devices and yet few ab initio studies of QD solids focus on these 2586

DOI: 10.1021/acs.jctc.7b00088 J. Chem. Theory Comput. 2017, 13, 2581−2590

Article

Journal of Chemical Theory and Computation

separate sites. As a result, the difference in the energy levels is largely independent of the intersite distance. The example shown here indicates that CDFT is a more robust framework to compute electronic couplings, especially when thermally disordered configurations are involved.

5. CONCLUSIONS In summary, we described a robust method for studying charge transport at the nanoscale using periodic constrained DFT. We compared our results to those found with other codes and localization potentials, showing minor differences coming from the use of different pseudopotentials and plane wave codes. Substantial differences were instead found for different localization potentials and for different convergence thresholds. Finally, our results for nanodots demonstrated key limitations of methods such as the energy splitting in dimer or the anticrossing methods. For silicon quantum dot dimers, constrained DFT predicts physical (i.e., exponential) decay rates for the electronic coupling of disordered geometries. By contrast, the energy splitting in dimer and anticrossing methods fail for hole transfers, resulting in an unphysical distance dependence of the electronic coupling. Work is in progress to simulate more complex nanoscale materials using the methodology presented here.

Figure 5. Isosurfaces of the highest occupied single particle molecular orbital (HOMO) of a Si35H36 dimer, as obtained using DFT and the PBE functional, for a configuration relaxed at T = 0 (top) and one extracted from an MD simulation performed at 300 K (bottom). The separation between the nanoparticles is 3.75 Å. White and yellow spheres represent hydrogen and silicon, respectively.

anticrossing method. Within the splitting in dimer (SID) method, the electron (hole) transfer coupling is computed from the difference of the lowest unoccupied molecular orbital, or LUMO (highest occupied molecular orbital, or HOMO) and the LUMO+1 (HOMO−1) one-electron energy levels: 1 HAB = [E LUMO + 1(HOMO) − E LUMO(HOMO − 1)] (9) 2 41,90 The anticrossing (hereafter AX) method utilizes the response of one-electron wave functions to applied electric fields to estimate the electronic coupling. The applied electric field is varied to find the nearest approach of the HOMO (LUMO+1) and HOMO−1 (LUMO) energy levels for hole (electron) transfer processes; then the coupling is evaluated using the SID method. We refer interested readers to representative publications for more details.41,90 In Figure 4, we use the CDFT, SID, and AX methods to calculate the electron coupling between a dimer of Si35H36 QDs with a geometry taken from a 300 K molecular dynamics trajectory. First, we consider the electronic couplings for hole transfer, on the top. While CDFT electronic couplings for holes −1 decrease exponentially (βhole CDFT = 4.10 Å ) for the disordered geometry, the SID and AX couplings behave unphysically, −1 remaining either mostly constant (βhole SID = 0.15 Å ) or behaving erratically. For the electron transfer (see bottom of Figure 4) the decay rates are 1.76, 0.97, and 1.35 Å−1 for the CDFT, SID, and AX methods, possibly implying different physical mechanisms responsible for the decay in the three methods. Two factors contribute to the failure of the SID and AX methods for the hole transfer: polarization effects from intermolecular interactions, as pointed out in the literature for the SID method;35,91 and disorder in the QD geometry leading to symmetry breaking. The single particle wave functions obtained for the relaxed configuration at T = 0 and for a snapshot extracted from an MD system demonstrate this symmetry breaking, as shown in Figure 5. As expected, the HOMO wave function is delocalized across both sites in the relaxed configuration; however, the HOMO and HOMO−1 wave functions of the thermally disordered dot are localized on



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Matthew B. Goldey: 0000-0002-2390-9554 Nicholas P. Brawand: 0000-0002-0624-3666 Márton Vörös: 0000-0003-1321-9207 Funding

This work (M.B.G.) was supported by the Center for Hierarchical Materials Design (CHiMaD) from the U.S. Department of Commerce, National Institute of Standards and Technology Award 70NANB14H012. Work by N.P.B. and G.G. was supported by the Midwest Integrated Center for Computational Materials (MICCoM), as part of the Computational Materials Sciences Program funded by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division. Work by M.V. was supported by the Laboratory Directed Research and Development (LDRD) funding from Argonne National Laboratory, provided by the Director, Office of Science, of the U.S. DOE under Contract No. DE-AC02-06CH11357. Computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program; the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231 and the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE AC02 06CH11357. Additional resources and computer time from the University of Chicago Research Computing Center were used in part of this research. Notes

The authors declare no competing financial interest. 2587

DOI: 10.1021/acs.jctc.7b00088 J. Chem. Theory Comput. 2017, 13, 2581−2590

Article

Journal of Chemical Theory and Computation



(19) Pizzi, G.; Volja, D.; Kozinsky, B.; Fornari, M.; Marzari, N. BoltzWann: A code for the evaluation of thermoelectric and electronic transport properties with a maximally-localized Wannier functions basis. Comput. Phys. Commun. 2014, 185, 422−429. (20) Madsen, G. K.; Singh, D. J. BoltzTraP. A code for calculating band-structure dependent quantities. Comput. Phys. Commun. 2006, 175, 67−71. (21) Li, W. Electrical transport limited by electron-phonon coupling from Boltzmann transport equation: An ab initio study of Si, Al, and MoS2. Phys. Rev. B: Condens. Matter Mater. Phys. 2015, 92, 075405. (22) Poncé, S.; Margine, E.; Verdi, C.; Giustino, F. EPW: Electronphonon coupling, transport and superconducting properties using maximally localized Wannier functions. Comput. Phys. Commun. 2016, 209, 116−133. (23) Ferrer, J.; Lambert, C. J.; García-Suárez, V. M.; Manrique, D. Z.; Visontai, D.; Oroszlany, L.; Rodríguez-Ferradás, R.; Grace, I.; Bailey, S.; Gillemot, K.; Sadeghi, H.; Algharagholy, L. A. GOLLUM: a nextgeneration simulation tool for electron, thermal and spin transport. New J. Phys. 2014, 16, 093029. (24) Stokbro, K.; Taylor, J.; Brandbyge, M.; Ordejon, P. TranSIESTA: a spice for molecular electronics. Ann. N. Y. Acad. Sci. 2003, 1006, 212−226. (25) Bernardi, M. First-principles dynamics of electrons and phonons. Eur. Phys. J. B 2016, 89, 239. (26) Bässler, H. Charge Transport in Disordered Organic Photoconductors a Monte Carlo Simulation Study. Phys. Status Solidi B 1993, 175, 15−56. (27) Carbone, I.; Carter, S. A.; Zimanyi, G. T. Monte Carlo Modeling of Transport in PbSe Nanocrystal Films. J. Appl. Phys. 2013, 114, 193709. (28) Marcus, R. A. Electron Transfer Reactions in Chemistry. Theory and Experiment. Rev. Mod. Phys. 1993, 65, 599−610. (29) Van Dyke, J. S.; Morr, D. K. Controlling the flow of spin and charge in nanoscopic topological insulators. Phys. Rev. B: Condens. Matter Mater. Phys. 2016, 93, 081401. (30) Morr, D. K. Crossover from quantum to classical transport. Contemp. Phys. 2016, 57, 19−45. (31) Tessler, N.; Preezant, Y.; Rappaport, N.; Roichman, Y. Charge Transport in Disordered Organic Materials and Its Relevance to ThinFilm Devices: A Tutorial Review. Adv. Mater. 2009, 21, 2741−2761. (32) Ortmann, F.; Bechstedt, F.; Hannewald, K. Charge transport in organic crystals: Theory and modelling. Phys. Status Solidi B 2011, 248, 511−525. (33) Troisi, A. Charge transport in high mobility molecular semiconductors: classical models and new theories. Chem. Soc. Rev. 2011, 40, 2347−2358. (34) Bässler, H.; Köhler, A. Charge transport in organic semiconductors. Top. Curr. Chem. 2011, 312, 1−65. (35) 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. (36) Kaduk, B.; Kowalczyk, T.; Van Voorhis, T. Constrained Density Functional Theory. Chem. Rev. 2012, 112, 321−370. (37) Cave, R. J.; Newton, M. D. Calculation of Electronic Coupling Matrix Elements for Ground and Excited State Electron Transfer Reactions: Comparison of the Generalized Mulliken Hush and Block Diagonalization Methods. J. Chem. Phys. 1997, 106, 9213. (38) Subotnik, J. E.; Vura-Weis, J.; Sodt, A. J.; Ratner, M. A. Predicting Accurate Electronic Excitation Transfer Rates via Marcus Theory with Boys or Edmiston-Ruedenberg Localized Diabatization. J. Phys. Chem. A 2010, 114, 8665−75. (39) Ortmann, F.; Bechstedt, F.; Hannewald, K. Charge transport in organic crystals: interplay of band transport, hopping and electronphonon scattering. New J. Phys. 2010, 12, 023011. (40) Maiti, S. K.; Nitzan, A. Mobility Edge Phenomenon in a Hubbard Chain: A Mean Field Study. Phys. Lett. A 2013, 377, 1205− 1209.

ACKNOWLEDGMENTS We are grateful for discussions with R. L. McAvoy and A. P. Gaiduk.



REFERENCES

(1) Kovalenko, M. V.; Manna, L.; Cabot, A.; Hens, Z.; Talapin, D. V.; Kagan, C. R.; Klimov, V. I.; Rogach, A. L.; Reiss, P.; Milliron, D. J.; Guyot-Sionnnest, P.; Konstantatos, G.; Parak, W. J.; Hyeon, T.; Korgel, B. A.; Murray, C. B.; Heiss, W. Prospects of Nanoscience with Nanocrystals. ACS Nano 2015, 9, 1012−1057. (2) Gutsch, S.; Laube, J.; Hartel, A. M.; Hiller, D.; Zakharov, N.; Werner, P.; Zacharias, M. Charge Transport in Si Nanocrystal/SiO2 Superlattices. J. Appl. Phys. 2013, 113, 133703. (3) Stavrinadis, A.; Konstantatos, G. Strategies for the Controlled Electronic Doping of Colloidal Quantum Dot Solids. ChemPhysChem 2016, 17, 632−644. (4) Schnabel, M.; Canino, M.; Kühnhold-Pospischil, S.; LópezVidrier, J.; Klugermann, T.; Weiss, C.; López-Conesa, L.; ZschintzschDias, M.; Summonte, C.; Löper, P.; Janz, S.; Wilshaw, P. R. Charge Transport in Nanocrystalline SiC with and without Embedded Si Nanocrystals. Phys. Rev. B: Condens. Matter Mater. Phys. 2015, 91, 195317. (5) Kagan, C. R.; Murray, C. B. Charge Transport in Strongly Coupled Quantum Dot Solids. Nat. Nanotechnol. 2015, 10, 1013− 1026. (6) Wippermann, S.; He, Y.; Vörös, M.; Galli, G. Novel silicon phases and nanostructures for solar energy conversion. Appl. Phys. Rev. 2016, 3, 040807. (7) Cheng, K.-Y.; Anthony, R.; Kortshagen, U. R.; Holmes, R. J. Hybrid Silicon NanocrystalOrganic Light-Emitting Devices for Infrared Electroluminescence. Nano Lett. 2010, 10, 1154−1157. (8) Priolo, F.; Gregorkiewicz, T.; Galli, M.; Krauss, T. F. Silicon Nanostructures for Photonics and Photovoltaics. Nat. Nanotechnol. 2014, 9, 19−32. (9) Luppi, M.; Ossicini, S. Ab Initio Study on Oxidized Silicon Clusters and Silicon Nanocrystals Embedded in SiO2: Beyond the Quantum Confinement Effect. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 71, 035340. (10) Voznyy, O.; Zhitomirsky, D.; Stadler, P.; Ning, Z.; Hoogland, S.; Sargent, E. H. A Charge-Orbital Balance Picture of Doping in Colloidal Quantum Dot Solids. ACS Nano 2012, 6, 8448−8455. (11) Lin, Z.; Franceschetti, A.; Lusk, M. T. Size Dependence of the Multiple Exciton Generation Rate in CdSe Quantum Dots Size Dependence of the Multiple Exciton Generation Rate in CdSe Quantum Dots. ACS Nano 2011, 5, 2503−2511. (12) Li, H.; Wu, Z.; Lusk, M. T. Dangling Bond Defects: The Critical Roadblock to Efficient Photoconversion in Hybrid Quantum Dot Solar Cells. J. Phys. Chem. C 2014, 118, 46−53. (13) Rocca, D.; Vö r ö s , M.; Gali, A.; Galli, G. Ab Initio Optoelectronic Properties of Silicon Nanoparticles: Excitation Energies, Sum Rules, and TammDancoff Approximation. J. Chem. Theory Comput. 2014, 10, 3290−3298. (14) Govoni, M.; Galli, G. Large Scale GW Calculations. J. Chem. Theory Comput. 2015, 11, 2680−2696. (15) Brawand, N. P.; Goldey, M. B.; Vörös, M.; Galli, G. Defect States and Charge Transport in Quantum Dot Solids. Chem. Mater. 2017, 29, 1255−1262. (16) Brawand, N. P.; Vörös, M.; Govoni, M.; Galli, G. Generalization of Dielectric-Dependent Hybrid Functionals to Finite Systems. Phys. Rev. X 2016, 6, 041002. (17) Kilina, S. V.; Tamukong, P. K.; Kilin, D. S. Surface Chemistry of Semiconducting Quantum Dots: Theoretical Perspectives. Acc. Chem. Res. 2016, 49, 2127−2135. (18) Myöhänen, P.; Stan, A.; Stefanucci, G.; van Leeuwen, R. Kadanoff-Baym approach to quantum transport through interacting nanoscale systems: From the transient to the steady-state regime. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 80, 115107. 2588

DOI: 10.1021/acs.jctc.7b00088 J. Chem. Theory Comput. 2017, 13, 2581−2590

Article

Journal of Chemical Theory and Computation

(62) Melander, M.; Jȯnsson, E. Ö .; Mortensen, J. J.; Vegge, T.; García Lastra, J. M. Implementation of Constrained DFT for Computing Charge Transfer Rates within the Projector Augmented Wave Method. J. Chem. Theory Comput. 2016, 12, 5367−5378. (63) Lao, K. U.; Herbert, J. M. Energy Decomposition Analysis With a Stable Charge-Transfer Term for Interpreting Intermolecular Interactions. J. Chem. Theory Comput. 2016, 12, 2569−2582. (64) Gajdos, F.; Valner, S.; Hoffmann, F.; Spencer, J.; Breuer, M.; Kubas, A.; Dupuis, M.; Blumberger, J. Ultrafast Estimation of Electronic Couplings for Electron Transfer Between Π-Conjugated Organic Molecules. J. Chem. Theory Comput. 2014, 10, 4653−4660. (65) Kubas, A.; Hoffmann, F.; Heck, A.; Oberhofer, H.; Elstner, M.; Blumberger, J. Electronic Couplings for Molecular Charge Transfer: Benchmarking CDFT, FODFT, and FODFTB Against High-Level Ab Initio Calculations. J. Chem. Phys. 2014, 140, 104105. (66) Kubas, A.; Gajdos, F.; Heck, A.; Oberhofer, H.; Elstner, M.; Blumberger, J. Electronic Couplings for Molecular Charge Transfer: Benchmarking CDFT, FODFT, and FODFTB Against High-Level Ab Initio Calculations. II. Phys. Chem. Chem. Phys. 2015, 17, 14342− 14354. (67) Gajdos, F.; Oberhofer, H.; Dupuis, M.; Blumberger, J. On the Inapplicability of Electron-Hopping Models for the Organic Semiconductor Phenyl-C61-Butyric Acid Methyl Ester (PCBM). J. Phys. Chem. Lett. 2013, 4, 1012−1017. (68) Oberhofer, H.; Blumberger, J. Revisiting Electronic Couplings and Incoherent Hopping Models for Electron Transport in Crystalline C60 at Ambient Temperatures. Phys. Chem. Chem. Phys. 2012, 14, 13846. (69) Blumberger, J.; McKenna, K. P. Constrained Density Functional Theory Applied to Electron Tunnelling Between Defects in MgO. Phys. Chem. Chem. Phys. 2013, 15, 2184−96. (70) Goldey, M. B.; Reid, D.; de Pablo, J.; Galli, G. Planarity and multiple components promote organic photovoltaic efficiency by improving electronic transport. Phys. Chem. Chem. Phys. 2016, 18, 31388−31399. (71) Vörös, M.; Brawand, N. P.; Galli, G. Hydrogen Treatment as a Detergent of Electronic Trap States in Lead Chalcogenide Nanoparticles. Chem. Mater. 2017, 29, 2485−2493. (72) Dreuw, A.; Head-Gordon, M. Failure of Time-Dependent Density Functional Theory for Long-Range Charge-Transfer Excited States: The Zincbacteriochlorin-Bacterlochlorin and Bacteriochlorophyll-Spheroidene Complexes. J. Am. Chem. Soc. 2004, 126, 4007− 4016. (73) Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; Dal Corso, A.; de Gironcoli, S.; Fabris, S.; Fratesi, G.; Gebauer, R.; Gerstmann, U.; Gougoussis, C.; Kokalj, A.; Lazzeri, M.; Martin-Samos, L.; Marzari, N.; Mauri, F.; Mazzarello, R.; Paolini, S.; Pasquarello, A.; Paulatto, L.; Sbraccia, C.; Scandolo, S.; Sclauzero, G.; Seitsonen, A. P.; Smogunov, A.; Umari, P.; Wentzcovitch, R. M. QUANTUM ESPRESSO: A Modular and Open-Source Software Project for Quantum Simulations of Materials. J. Phys.: Condens. Matter 2009, 21, 395502. (74) Hirshfeld, F. L. Bonded-Atom Fragments for Describing Molecular Charge Densities. Theor. Chem. Accounts Theory, Comput. Model. (Theoretica Chim. Acta) 1977, 44, 129−138. (75) Löwdin, P. On the Non-Orthogonality Problem Connected With the Use of Atomic Wave Functions in the Theory of Molecules and Crystals. J. Chem. Phys. 1950, 18, 365−375. (76) Roychoudhury, S.; Motta, C.; Sanvito, S. Charge Transfer Energies of Benzene Physisorbed on a Graphene Sheet From Constrained Density Functional Theory. Phys. Rev. B: Condens. Matter Mater. Phys. 2016, 93, 045130. (77) Souza, A. M.; Rungger, I.; Pemmaraju, C. D.; Schwingenschloegl, U.; Sanvito, S. Constrained-DFT Method for Accurate Energy-Level Alignment of Metal/Molecule Interfaces. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 88, 165112.

(41) Shan, W.; Walukiewicz, W.; Ager, J.; Haller, E.; Geisz, J.; Friedman, D.; Olson, J.; Kurtz, S. Band Anticrossing in GaInNAs Alloys. Phys. Rev. Lett. 1999, 82, 1221−1224. (42) Chen, T.; Reich, K. V.; Kramer, N. J.; Fu, H.; Kortshagen, U. R.; Shklovskii, B. I. Metal-Insulator Transition in Films of Doped Semiconductor Nanocrystals. Nat. Mater. 2016, 15, 299−303. (43) Janak, J. Proof That dE/dn_i = E_i in Density-Functional Theory. Phys. Rev. B: Condens. Matter Mater. Phys. 1978, 18, 7165− 7168. (44) Dederichs, P. H.; Blügel, S.; Zeller, R.; Akai, H. Ground States of Constrained Systems: Application to Cerium Impurities. Phys. Rev. Lett. 1984, 53, 2512−2515. (45) Prezhdo, O. V.; Kindt, J. T.; Tully, J. C. Perturbed ground state method for electron transfer. J. Chem. Phys. 1999, 111, 7818−7827. (46) Wu, Q.; Van Voorhis, T. Direct Optimization Method to Study Constrained Systems within Density-Functional Theory. Phys. Rev. A: At., Mol., Opt. Phys. 2005, 72, 024502. (47) Wu, Q.; Van Voorhis, T. Constrained Density Functional Theory and Its Application in Long-Range Electron Transfer. J. Chem. Theory Comput. 2006, 2, 765−774. (48) Wu, Q.; Van Voorhis, T. Extracting Electron Transfer Coupling Elements from Constrained Density Functional Theory. J. Chem. Phys. 2006, 125, 164105. (49) Wu, Q.; Van Voorhis, T. Direct Calculation of Electron Transfer Parameters Through Constrained Density Functional Theory. J. Phys. Chem. A 2006, 110, 9212−8. (50) Blumberger, J.; Tavernelli, I.; Klein, M. L.; Sprik, M. Diabatic Free Energy Curves and Coordination Fluctuations for the Aqueous Ag+/Ag2+ Redox Couple: A Biased Born-Oppenheimer Molecular Dynamics Investigation. J. Chem. Phys. 2006, 124, 064507. (51) Oberhofer, H.; Blumberger, J. Charge Constrained Density Functional Molecular Dynamics for Simulation of Condensed Phase Electron Transfer Reactions. J. Chem. Phys. 2009, 131, 064101. (52) Oberhofer, H.; Blumberger, J. Electronic Coupling Matrix Elements From Charge Constrained Density Functional Theory Calculations Using a Plane Wave Basis Set. J. Chem. Phys. 2010, 133, 244105. (53) Van Voorhis, T.; Kowalczyk, T.; Kaduk, B.; Wang, L.-P.; Cheng, C.-L.; Wu, Q. The Diabatic Picture of Electron Transfer, Reaction Barriers, and Molecular Dynamics. Annu. Rev. Phys. Chem. 2010, 61, 149−70. (54) Ghosh, P.; Gebauer, R. Computational approaches to charge transfer excitations in a zinc tetraphenylporphyrin and C70 complex. J. Chem. Phys. 2010, 132, 104102. (55) Sena, A. M. P.; Miyazaki, T.; Bowler, D. R. Linear Scaling Constrained Density Functional Theory in CONQUEST. J. Chem. Theory Comput. 2011, 7, 884−889. (56) Ř ezác,̌ J.; Lévy, B.; Demachy, I.; de la Lande, A. Robust and Efficient Constrained DFT Molecular Dynamics Approach for Biochemical Modeling. J. Chem. Theory Comput. 2012, 8, 418−427. (57) Souza, A. M.; Rungger, I.; Pemmaraju, C. D.; Schwingenschloegl, U.; Sanvito, S. Constrained-DFT method for accurate energy-level alignment of metal/molecule interfaces. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 88, 165112. (58) Ratcliff, L. E.; Grisanti, L.; Genovese, L.; Deutsch, T.; Neumann, T.; Danilov, D.; Wenzel, W.; Beljonne, D.; Cornil, J. Toward Fast and Accurate Evaluation of Charge On-Site Energies and Transfer Integrals in Supramolecular Architectures Using Linear Constrained Density Functional Theory (CDFT)-Based Methods. J. Chem. Theory Comput. 2015, 11, 2077−2086. (59) Ma, P.-W.; Dudarev, S. L. Constrained density functional for noncollinear magnetism. Phys. Rev. B: Condens. Matter Mater. Phys. 2015, 91, 054420. (60) Ramos, P.; Pavanello, M. Constrained subsystem density functional theory. Phys. Chem. Chem. Phys. 2016, 18, 21172−21178. (61) Holmberg, N.; Laasonen, K. Efficient Constrained Density Functional Theory Implementation for Simulation of Condensed Phase Electron Transfer Reactions. J. Chem. Theory Comput. 2017, 13, 587−601. 2589

DOI: 10.1021/acs.jctc.7b00088 J. Chem. Theory Comput. 2017, 13, 2581−2590

Article

Journal of Chemical Theory and Computation (78) Neaton, J. B.; Hybertsen, M. S.; Louie, S. G. Renormalization of Molecular Electronic Levels at Metal-Molecule Interfaces. Phys. Rev. Lett. 2006, 97, 216405. (79) Hamann, D. R. Optimized Norm-Conserving Vanderbilt Pseudopotentials. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 88, 085117. (80) Schlipf, M.; Gygi, F. Optimization Algorithm for the Generation of ONCV Pseudopotentials. Comput. Phys. Commun. 2015, 196, 36− 44. (81) Troullier, N.; Martins, J. Efficient pseudopotentials for planewave calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 1991, 43, 1993. (82) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (83) Martyna, G. J.; Tuckerman, M. E. A Reciprocal Space Based Method for Treating Long Range Interactions in Ab Initio and ForceField-Based Calculations in Clusters. J. Chem. Phys. 1999, 110, 2810. (84) Lejaeghere, K.; Bihlmayer, G.; Björkman, T.; Blaha, P.; Blügel, S.; Blum, V.; Caliste, D.; Castelli, I. E.; Clark, S. J.; Dal Corso, A.; de Gironcoli, S.; Deutsch, T.; Dewhurst, J. K.; Di Marco, I.; Draxl, C.; Dułak, M.; Eriksson, O.; Flores-Livas, J. A.; Garrity, K. F.; Genovese, L.; Giannozzi, P.; Giantomassi, M.; Goedecker, S.; Gonze, X.; Grånäs, O.; Gross, E. K. U.; Gulans, A.; Gygi, F.; Hamann, D. R.; Hasnip, P. J.; Holzwarth, N. A. W.; Iuşan, D.; Jochym, D. B.; Jollet, F.; Jones, D.; Kresse, G.; Koepernik, K.; Kücu̧ ̈kbenli, E.; Kvashnin, Y. O.; Locht, I. L. M.; Lubeck, S.; Marsman, M.; Marzari, N.; Nitzsche, U.; Nordström, L.; Ozaki, T.; Paulatto, L.; Pickard, C. J.; Poelmans, W.; Probert, M. I. J.; Refson, K.; Richter, M.; Rignanese, G.-M.; Saha, S.; Scheffler, M.; Schlipf, M.; Schwarz, K.; Sharma, S.; Tavazza, F.; Thunström, P.; Tkatchenko, A.; Torrent, M.; Vanderbilt, D.; van Setten, M. J.; Van Speybroeck, V.; Wills, J. M.; Yates, J. R.; Zhang, G.-X.; Cottenier, S. Reproducibility in Density Functional Theory Calculations of Solids. Science 2016, 351, aad3000. (85) Becke, A. D. A Multicenter Numerical Integration Scheme for Polyatomic Molecules. J. Chem. Phys. 1988, 88, 2547. (86) Andreoni, W.; Curioni, A. New Advances in Chemistry and Material Science With CPMD and Parallel Computing. Parallel Comput. 2000, 26, 819−842. (87) Lippert, G.; Hutter, J.; Parrinello, M. A hybrid Gaussian and plane wave density functional scheme. Mol. Phys. 1997, 92, 477−488. (88) Due to degeneracies in highly symmetric molecules, the chargelocalized state can occupy one of several energetically equivalent orbitals, which have varying overlaps and electronic couplings with charge-localized states on neighboring molecules. This problem was addressed by Kubas et al.,65 who chose the specific orbitals a priori in these cases. Since reading and reordering fragment-derived oneelectron orbitals is not a routine procedure for Quantum ESPRESSO, we neglect these systems. As a result, these data points were excluded from Tables 2, 3, and 4 and Figure 3. (89) Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; Dal Corso, A.; de Gironcoli, S.; Fabris, S.; Fratesi, G.; Gebauer, R.; Gerstmann, U.; Gougoussis, C.; Kokalj, A.; Lazzeri, M.; Martin-Samos, L.; Marzari, N.; Mauri, F.; Mazzarello, R.; Paolini, S.; Pasquarello, A.; Paulatto, L.; Sbraccia, C.; Scandolo, S.; Sclauzero, G.; Seitsonen, A. P.; Smogunov, A.; Umari, P.; Wentzcovitch, R. M. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. J. Phys.: Condens. Matter 2009, 21, 395502. (90) Chu, I.-H.; Radulaski, M.; Vukmirovic, N.; Cheng, H.-P.; Wang, L.-W. Charge Transport in a Quantum Dot Supercrystal. J. Phys. Chem. C 2011, 115, 21409−21415. (91) Valeev, E. F.; Coropceanu, V.; da Silva Filho, D. A.; Salman, S.; Brédas, J.-L. Effect of Electronic Polarization on Charge-Transport Parameters in Molecular Organic Semiconductors. J. Am. Chem. Soc. 2006, 128, 9882−9886.

2590

DOI: 10.1021/acs.jctc.7b00088 J. Chem. Theory Comput. 2017, 13, 2581−2590