Modulating electronic coupling at the quantum dot/molecule interface by wavefunction engineering

Interfacial charge transfer (CT) is a fundamental process in photocatalysis,¹ biocatalysis,^2–5 conversion of solar energy to chemical and electrical energy,^6,7 and many other fields involving charge carrier dynamics. Efficiency of these processes depends on the charge generation in the utilized Quantum Dots (QDs), semiconductors, biocatalysts, and photocatalysts,⁸ as well as on the interfacial charge transfer dynamics and the rate of charge separation between these materials and adsorbed small molecules. Therefore, controlling the rate of interfacial charge (either electron or hole) transfer is critical for designing efficient catalytic processes and functional devices^9–11 (see Scheme 1).

Marcus theory¹² has long been a reliable tool for studying electron transfer (ET) in molecular and biomolecular dyads consisting of electron donor and acceptor moieties.¹³ The three primary factors controlling the ET rate, at a given temperature, are the electronic coupling strength, the change in the free energy (i.e., the driving force), and the solvent reorganization energy. While the origins of the driving force and solvent effects can be straightforwardly analyzed by classical thermodynamics and electrostatics theories,^14–17 scrutinizing electronic coupling, which is purely quantum mechanical in nature,^12,18–20 requires a deep knowledge of the multi-electron wavefunctions of the involved systems. A common empirical approach to control the nature of electronic coupling is to utilize the intrinsic quantum confinement effects in finite size nanocrystals,²¹ also called wavefunction engineering, by varying valence and conduction band (VB and CB) offsets via mixing different materials, changing the size and facilitating electron/hole delocalization effects.^22–25 Conversely, calculations based on first principles can shed light on the nature of electronic coupling. Noting that the wavefunctions of bulk semiconductors, photocatalysts, and molecules can be readily obtained by modern computational approaches [such as atomic-centered basis and plane-wave density functional theory (DFT)], the generation of QD wavefunctions requires special attention. It is expected that revealing the nature of electronic coupling at the QD/molecule interface by wavefunction engineering of both the QD and the molecule with accurate electronic structure calculations will be instrumental in controlling charge carrier dynamics in these systems.

Previously, we have studied ET between anthraquinone (AQ) and a family of CdSe/ZnS core-shell QDs using a hybrid molecular orbital approach.²⁶ In that study, we used QDs with a 3 nm diameter core and up to 2 nm of shell thickness, treated the core-shell QD and AQ with the effective mass approximation (EMA)²⁷ and a density functional theory, respectively. To overcome the known limitations of the approximate EMA theory (the envelope function is not well defined outside the lattice, e.g., in vacuum), recently, we developed an atomistic, Bulk-Adjusted Linear Combination of Atomic Orbitals (BA-LCAO)²⁸ approach that allows the generation of QD orbitals by targeted parameter fitting to the experiment-corrected DFT-quality bulk band structure. In the present paper, we explore the impact of the molecular orbital spatial modification, in other words wavefunction engineering, on the QD/molecule electronic coupling in spherical wurtzite-CdSe/AQ dyads. With the above, we stress that our primary goal is to characterize this impact rather than to calculate explicit CT rates in either the electron or hole channels. Thus, to achieve our goals, we integrate the BA-LCAO generated wavefunction of spherical wurtzite-CdSe QDs with that of the AQ molecule obtained at a DFT level of theory (see below) and calculate the coupling at the interface of CdSe quantum dot and anthraquinone-2,3-dicarboxylic acid. In order to alter the spatial probability distribution of acceptor (i.e., AQ) orbitals, we modify its highest occupied (HOMO) and lowest unoccupied (LUMO) molecular orbitals by introducing various substituents. Finally, to model the experimental conditions in which QD/molecule complexes are typically prepared, we carry out rigorous molecular dynamics (MD) simulations at room temperature in weak and strong solvent regimes by using a specifically designed non-bonding force field that has previously been tested by us in similar applications.²⁸ The results of these calculations are presented and discussed in comparison with available experimental measurements.

Briefly, the BA-LCAO approach for periodic structure calculations combines the features of environment-modified total-energy-based tight-binding (TETB)^29,30 and the standard orbital-based density-functional tight-binding (DFTB)^31,32 methods. It applies a non-local, environment dependent scaling function to the two-center interaction potentials (derived independently using DFT data) strictly according to atomic valences acquired naturally in a bulk structure, similar to that used by others.^29,30 To advance practical application of this method, the BA-LCAO Hamiltonian must be parameterized by fitting the band structure to a high-level DFT bulk reference, corrected for experimentally measured band edges. In the present work, the reference bulk band structure is derived from a projector-augmented-wave generalized-gradient-approximation pseudopotential calculation with the Perdew-Burke-Ernzerhof functional, as reported in our earlier work.²⁸ 

In BA-LCAO, the molecular orbitals of the QD are expressed as linear combinations of normalized atom-centered orbitals (AOs),

where j runs over all basis functions on atomic center a. The expansion coefficients are found by solving the secular equation

where H is the Hamiltonian matrix, S is the AO overlap matrix, c_I is the solution vector, and α_I is the corresponding eigenvalue, which is interpreted as the energy of ψ_I. The off-diagonal elements of the overlap matrix S are integrals involving two AOs at positions R_a and R_b,

Following the 2-center LCAO model,³³ the matrix representation of the Hamiltonian in an AO basis is

where a and b are the atomic sites with respective AO orbitals i and j; ε_ai is the ith energy level of free atom a, and $Vijab$ is the off-diagonal element representing the degree of AO mixing in the interacting cluster. The latter depends on the properties of the pair of atoms, their internuclear distance, and valences (see below). These potentials are computed for the three atom-pair sets considered in this work: Cd–Cd, Cd–Se, and Se–Se, where the atomic basis sets {φ_bj} are the valence Cd(4d5s5p) and Se(4s4p) functions. In Table I, we give the atomic energies and Slater exponents used in the calculations.

The method takes a semi-empirical approach to deriving the elements in Eq. (3b) guided by the goals of (i) properly describing the long-range tails of molecular orbitals of nanoparticles in the vacuum space and (ii) recovering the qualitative properties of these MOs at high and low energies, i.e., deep regions of the CBs and VBs, respectively. To this end, the diagonal elements ε_ai are taken from experimental values of IPs of neutral atoms, and the AO basis is hydrogen-like functions with the Slater exponents taken from non-relativistic all-electron SCF wavefunctions.^34,35 In other words, the first IP formally corresponds to the HOMO, for closed shell atoms, or to the singly occupied (SO) MO, for open-shell atoms, energy level. The off-diagonal matrix elements are derived from DFT calculations of diatomics. We note that this approach is different from that used in DFTB where the off-diagonal matrix elements are calculated directly by evaluating one-electron integrals between the two orbital centers.³¹ 

Environment-induced effects^29,30 are incorporated into Eq. (3b) by using a valence function for each atom in the nanoparticle,

where η is a distance-dependent bond function, bound on [0, 1], between a pair of atoms a and c, with the summation running over all atoms.²⁸ If a diatomic a-b is surrounded by other atoms in a cluster, it will have a high valence number N_a + N_b, and, therefore, its diatomic identity will be modified. To address this, we proposed a simple valence-dependent modifying function²⁸ 

where $Na′R=NaR−η(Rab)$ is a reduced valence of atom a (in the presence of atom b) and $γlalbλab$ is an adjustable parameter. The integers l_a and l_b are the AO angular momenta on atoms a and b, respectively, and λ is the corresponding sum-projection on the bond’s axis.³³ The modifying functions are used as scaling factors of the corresponding off-diagonal elements in Eq. (3b). This is illustrated on the specific example of a p_z–p_z potential matrix element for an arbitrary Cd–Se pair at distance r apart,

where r_x, r_y, r_z are the laboratory frame Cartesian components of the bond vector r; R is the collective position vector of all the atoms in the cluster, $υzzσCdSe$ and $υzzπCdSe$ are the σ-bond and π-bond potentials, and $γzzσCdSe$ and $γzzπCdSe$ are parameters. The above expression is equivalent to the original Slater-Koster formula (used also in the standard DFTB) if $κzzσCdSe=κzzπCdSe=1$⁠, that is, if the matrix element V depends only on the relative position of the two atoms irrespective of the positions of the other atoms in the cluster. Thus, in BA-LCAO, the non-local effects are introduced by scaling the individual AO interactions.

It needs to be pointed out that the present treatment neglects particle correlation, such as electron-hole Coulomb and exchange interactions which are present in a single exciton. The interaction energies have been shown to be on the order of a few hundred meV for single excitons in CdSe and GaAs quantum dots by Zunger and co-workers.^36,37 These interactions can be incorporated into BA-LCAO at the level of singles CI whereby the leading few electron and hole orbitals, a small set of the unperturbed ψ_I wavefunctions from Eq. (1), are mixed to yield correlated one-electron excited states, i.e. single excitons, their energies and wavefunctions. Mixing excited states to 1S_e/1S_h orbitals will affect their behavior at the interface region (see discussion below), which may impact the overlaps and charge transfer probabilities. These corrections are out of scope of the present paper but are being considered for implementation in our future studies.

We note, in passing, the use of other approaches to correlation in excited states, such as those based on Green’s function formalism^38,39 and linear-response theories such as time-dependent HF and DFT,⁴⁰ can yield more accurate electronic spectra than a small-scale singles CI. However, application of these theories to large systems is problematic due to the self-consistent procedure for orbital calculation. More importantly, they do not provide the wavefunctions required for the calculation of coupling elements and, consequently, for charge transfer rates and interfacial charge transfer dynamics in the junction of nanostructures and molecules.

The electronic structure of the molecular electron/hole acceptors is described by a density functional theory. Here, we use the B3LYP functional with the 6-31++G(d,p) basis sets, with the diffuse functions added to provide a better description of the LUMO orbitals. All geometries were optimized from the starting structure generated by GaussView to the nearest minimum, where population analysis and orbital information, in the form of cube files, was performed and stored for electron transfer calculations. These calculations were performed with Gaussian 09 software.⁴¹ 

To demonstrate the effect of QD’s orbital extent on coupling with the molecule, we considered CdSe quantum dots of two sizes: (CdSe)₂₀₁ and (CdSe)₆₉₃, approximately corresponding to 2.75 nm and ∼4.08 nm diameters, respectively. These QDs were cut out of bulk wurtzite phase CdSe and subsequently refined to remove any single coordinated surface atoms. The remaining 2- and 3-coordinated atoms were “ligated” by applying shifting potentials, as described in our earlier work,²⁸ to remove any surface states entering the bandgap region. Using the valence basis sets, as described above, the sizes of the respective secular equations were 2613 and 9009.

In Fig. 1, we show spatial analysis of two frontier orbitals in each band, LUMO/LUMO + 1 in the CB and HOMO/HOMO-1 in the VB, in the two quantum confinement regimes. The energies are summarized in Table S1. The LUMOs and HOMOs correspond to the nodeless 1S and 1S_h type orbitals, while the LUMO + 1 and HOMO-1 correlate with 1P and 1P_h types of the single-nodal-plane, respectively. Visual inspection allows us to discern the slightly greater spatial extent of the 1P-type QD orbitals than of the 1S-type. Furthermore, consistent with spatial confinement theories, the electron density of the smaller dot has a significantly larger value at the surface and penetrates farther into the vacuum than the corresponding density of the larger dot, demonstrating that BA-LCAO properly captures this well-known phenomenon. Thus, the size effect allows tuning, in a very simple way, of electron’s wavefunction amplitude outside the dot. Another important feature observed in Fig. 1 is the difference in the spatial extent between the two bands, i.e., the electron (LUMO/LUMO + 1) extends slightly further into the vacuum than does the hole (HOMO/HOMO-1). From this, one expects a stronger coupling of the molecule with the QD’s conduction band than the valence band, which implies a more efficient electron transfer than hole transfer (HT). This behavior has been described previously^22,26 and is fully predictable on the basis of the effective masses of the electron and the hole.

As a molecular electron/hole acceptor, we take anthraquinone-2,3-dicarboxylic acid, C₁₄H₆O₂(CO₂H)₂, which is abbreviated as AQ throughout the paper (see Fig. 2). In order to study the impact of the frontier orbitals (i.e., electronic properties) of the acceptor on the QD/Acceptor interfacial charge transfer rate, we introduce two types of modifications to the original AQ molecule, namely, substituting (a) 7-H by electron withdrawing groups OH and OF that result in 7-ol-AQ and 7-hypofluorite-AQ or 7-OH-AQ and 7-OF-AQ, respectively, and (b) 10-oxo by electron donating groups BH and CH₂, resulting in 10-boraneylidene-AQ and 10-methylene-AQ or 10-BH-AQ and 10-CH₂-AQ, respectively, as shown in Fig. 2. For the sake of simplicity, these species are abbreviated as X-AQ, where X = 7-ol, 7-hypofluorite, 10-BH, and 10-CH₂.

Visual inspection of the MOs reveals that the HOMO of AQ is a σ-type orbital, symmetric about reflection in the plane, and with all its nodes orthogonal to the molecule’s plane (see Fig. 2). On the other hand, the HOMOs of all the substituted X-AQ are π-type orbitals, antisymmetric about reflection in the plane, i.e., with one of the nodal planes passing through the molecular plane. As we show below using numerical analysis, this creates a major effect on the extension of AQ’s HOMO in the direction perpendicular to its plane. The LUMO of AQ is a π-type orbital, and the 7-H substitution by OH and OF only slightly affects its nature. By contrast, the 10-oxo substitutions by BH and CH₂ lead to the substantial change of LUMO. In the case of 10-BH-AQ, the LUMO is almost entirely localized on the B-center as an in-plane 2p orbital, making the LUMO a σ-type molecular orbital. Similarly, but to a lesser extent, the shape of the LUMO of 10-CH₂-AQ is asymmetric toward CH₂ and retains its π-character.

To better characterize the changes to the HOMO/LUMO reference pair, we perform envelope orbital analysis by fitting an elliptic, 1S-type exponential to each of these MOs,

where ζ_k are the three principal exponents, determined as the square roots of the three eigenvalues of the inverse variance matrix σ⁻¹ and η_k are the three principal axes as the corresponding eigenvectors of σ_MO. A 3D ellipse, as is Eq. (7), is able to closer reflect the shape of a bulky polyatomic molecule. The “laboratory frame” Cartesian elements α = (x, y, z) and β = (x, y, z) of the variance matrix are defined as averages over the exact MO density,

with r₀ being the center of the MO density distribution, taken from the DFT calculation (see above). It can be shown that with the choice of ζ_k as above, the variance matrix σ is the same for the original MO and its fit, Eq. (7). In other words, the Cartesian variances of the MOs have the physical meaning that they correspond to Slater-type exponents of the orbital envelope. In Table II, we summarize the results of the fit for the MOs shown in Fig. 2.

For a clearer picture, the exponents are arranged according to the spatial orientation of the three principal axes (as seen in Fig. 2) in the same order for all MOs, that is, η₁ is best described as the out-of-plane direction, η₂ as the vertical direction, and η₃ as the horizontal direction. We note immediately that the out-of-plane tail of the AQ’s HOMO is dramatically extended by all of the substituents. This effect was alluded to previously from a simple visual inspection of the MOs. In the vertical direction, the HOMO is contracted, most noticeably by the OH substituent, which is also consistent with a visual analysis. And in the horizontal direction, the HOMO is moderately extended by the substitution. In the LUMO, the effect is subtle. The electron withdrawing substituents (OH and OF) do not seem to affect the LUMO’s tail to a significant extent. The electron donating substituents (BH and CH₂), however, do affect the LUMO, albeit in different ways. BH contracts the LUMO in the out-of-plane direction by making it a σ-type orbital, while CH₂ expands it by increasing its π-character (see above). In the in-plane directions, CH₂ has almost no effect, and BH elongates the tail in the vertical direction while contracting it in the horizontal. Figure 3 provides a more detailed representation of MO variation with substituent in the out-of-plane direction.

To calculate quantum coupling between QD and X-AQ, we introduce two approximations. First, we treat QD/(X-AQ) dyad as a one-particle system which implies 1 electron in the conduction band/LUMO space and 1 hole in the valence band/HOMO space. In the present work, these are non-interacting singly excited states arising from a non-correlated electron-hole pair. Second, we treat it as a weakly interacting dyad, i.e., we neglect any additional interactions between QD and X-AQ that can affect the zero-order orbitals. This approach rests on the assumption that the QD and the molecule do not form a covalent bond, but rather interact via non-covalent interactions. To validate this assumption, we performed electronic structure calculations using a smaller, computationally manageable model of QD, Cd₃₃Se₃₃. We calculated geometry, complexation energy and examined the nature of the frontier orbitals of the [Cd₃₃Se₃₃]/AQ dyad. Briefly, we have found that the effect of complexation on the fragments frontier orbitals, e.g., the LUMO, is very small (see Sec. S2 of the supplementary material). The Hamiltonian of such a weakly interacting system can be written as²⁶ 

in terms of the zero-order states, |n⟩, with corresponding energies E_n being the one-electron molecular orbitals of QD and X-AQ. The matrix element coupling any two orbitals across the QD/(X-AQ) interface, at a single geometric configuration, has a formal extended Huckel notation as follows:

where |m^QD⟩ and |l^X−AQ⟩ are the two electron orbitals of the QD and X-AQ, respectively. Under normal experimental conditions, a QD/(X-AQ) system exists in a thermal equilibrium, with the ligands attached to the QD as well as with the solvent molecules. Accordingly, the rate constant, as well as the electron transfer probability, must be evaluated as an ensemble average. Recalling the relationship of the coupling element to the charge transfer rate,¹³ $kETT∝HQDX-AQT2$⁠, we emphasize the thermal statistical nature of the coupling strength factor through its dependence on temperature. Utilizing the expression for the quantum coupling matrix element, Eq. (10), we define the coupling strength as an ensemble average over a Boltzmann distribution,

In this formulation, the fragment geometries (i.e., geometries of QD and X-AQ) are frozen in their respective bulk and global minimum structures, and thus, the orbital energies are constant and can be taken out of the ensemble average. The only variable is the orbital overlap integral, which depends on the distance between the fragments and their orientation.

To approach the task of evaluating the squared overlap integral appearing in Eq. (11) as a Boltzmann average, we perform molecular dynamics simulations of QD/(X-AQ) at 300 K. The force field describing the QD–(X-AQ) interactions is a previously developed fragmented molecular mechanics (FraMM) model²⁸ that considers the two main components retained in the weak interaction (non-bonding) limit used in the present work, (i) the partial charge interactions and (ii) the long-range van der Waals forces acting on the two, or more, as we discuss below, rigid fragments. We performed the additional calculations, as stated above, (see also Sec. S2 of the supplementary material), to validate the rigid fragment assumption and the use of the FraMM force field. Briefly, we have found that the nature of the interaction between QD and AQ molecule is non-covalent and the effect of complexation on fragments’ geometries is very small, less than 1% relative r.m.s. difference.

For the simplest case of two rigid fragments A = QD and B = X-AQ, the non-bonding interaction potential is

where R_AB is a six-dimensional vector describing relative orientation of the two fragments, N_A and N_B are the number of atoms in respective fragments, r^A and r^B are the atomic Cartesian coordinates, and ε is the solvent dielectric. The partial q^A and q^B charges are derived from Mulliken population analysis from the BA-LCAO calculated valence electron density for (CdSe)_n, and from the DFT calculations of X-AQ. The Lennard-Jones binding energies υ^A, υ^B and radii σ^A, σ^B are taken from the AMBER-99⁴² database for the various F, O, C, H atom types and from the Universal Force Field⁴³ database for Cd, Se, and B atoms (see Table III). It is imperative to point out that the rigid fragment model neglects intra-fragment dynamical effects, such as orbital energy and electron density, i.e., the atomic charges, response to internal bond fluctuations, as investigated by others.^44–46 Nevertheless, these fluctuations occur on a much faster time scale than inter-fragment overlap oscillations and are expected to be averaged out over the long simulation time (1 ns) and thus to have the same effect for all the AQ-X fragments. Other specific effects caused by structural fluctuations in real time, such as electron/hole trapping, do play a role in ET, but describing or understanding these phenomena is out of scope of the present work.

At present, no explicit solvent molecules or capping ligands are introduced. The solvent effects are measured by performing the calculations in two distinct regimes, (a) vacuum, with a dielectric constant of ε = 1.0 and (b) strong solvent regime (water) with a dielectric constant of ε = 78.3. To facilitate the sampling of X-AQ adsorption sites on the non-uniform and rugged (CdSe)_n surfaces, we simultaneously let 12 and 26 (roughly proportional to the surface area) X-AQ molecules (evenly distributed on the each surface) interact with (CdSe)₂₀₁ and (CdSe)₆₉₃, respectively. For all of these complexes and both solvent regimes, we ran a 1 ns trajectory at a fixed temperature (300 K), having first equilibrated the system for 100 ps (also at 300 K) and recorded the coordinates at every 0.1 ps, thus collecting 10 000 structures, per (CdSe)_n/(X-AQ)_J super-complex with J = 12 for n = 201 and J = 26 for n = 693, for statistical analysis and evaluating the average squared overlap appearing in Eq. (11). The MD simulations were done with TINKER-6.1 suite of codes.⁴⁷ The final value of the overlap is calculated as the time average along the equilibrated trajectory and then averaged over all X-AQ instances,

where τ = 1 ns. At every time point, the overlap integral is calculated by merging two MO cube files and integrating on the grid of the cube with the coarser mesh. The missing MO values at non-grid points of the other cube are found by a local 3D Lagrange interpolation on the finer mesh (refer to Scheme S4.1 in the supplementary material).

In Fig. 4, we show a representative snapshot from the MD trajectory for one of the complexes, (CdSe)₂₀₁/(AQ)₁₂, taken in both solvent environment regimes. Similar snapshots for all the other complexes are reported in the supplementary material. The common adsorption profile in both cases is for an AQ molecule to lay near flat on the QD’s surface at an average distance of ∼3 Å. There is an occasional edge-on orientation at the surface in vacuum, as can be seen in Fig. 4(a), where partial charges play a more significant role in the interaction and make such a configuration energetically feasible.

The results of the MD simulations recovered by Eq. (13) are compiled in Fig. 5. For a general analysis, we found it most revealing to represent the data as a function of the out-of-plane Slater exponent (ζ₁) of the HOMO and LUMO of X-AQ as it serves as the best measurement of the molecular size effect on the overlap. As seen in Fig. 5, there is a well-defined correlation between the QD and X-AQ orbital overlap and the extent of the MO’s “tail” going out of the X-AQ molecular plane: the farther the MO extends, the bigger its overlap with the QD. This correlation is present for both the electron transfer process (via 1S_e/LUMO and 1P_e/LUMO couplings) and the hole transfer process (via 1S_h/HOMO and 1P_h/HOMO couplings) regardless of the solvent regime. The linear correlation illustrated in the plots should not be interpreted in a physically meaningful way, however. The relative values of the coupling strength in the two processes are consistent with our earlier predictions based on the QD orbital properties (see Fig. 1); that is, the ET process is about two orders of magnitude more efficient than HT, assuming the driving forces, ΔG, are equal. Furthermore, the 1P excited states of the QD couple to the molecule consistently stronger than the 1S states, supporting the previously expressed prediction based on QD radial probability densities (Fig. 1). In the same fashion, the smaller of the two QDs exhibits a much stronger coupling to the molecule than does the larger QD, as can be seen by a direct comparison of Figs. 5(a) and 5(b).

Turning to group-specific analysis of the overlap data, we note that all of the substituents used in this work for AQ have a dramatic effect on the QD/(X-AQ) HT overlaps, compared with the near zero (∼10⁻⁹) 1S_h (1P_h)/HOMO overlaps in the QD/AQ complex. This follows precisely from the Slater orbital analysis of X-AQ HOMOs summarized in Table II and Fig. 3. For the QD/(X-AQ) electron transfer probability, the orbital effects are more pronounced only in the cases of X = BH (ζ₁ = 16.12 nm⁻¹) and X = CH₂ (ζ₁ = 11.78 nm⁻¹). Their effects on the overlap are, however, rather different: CH₂ enhancing electron transfer with the more diffuse LUMO, while BH slightly impeding it with the more compact LUMO, relative to the LUMO of the unmodified AQ. The electron donating groups OH and OF have little effect on the ET process, consistent with their unchanged ζ₁ exponents relative to AQ.

Further inspection of the data in Fig. 5 suggests that the solvent effect on orbital overlaps is more subtle than the orbital effect (presented above). One notable difference is in the ET curves, which show that “water” tends to minimize the 1P_e/LUMO overlap. All the other curves, on average, show that the effect of water environment on the coupling strength is likely spurious. In other words, solvent appears to play a minor role in the electronic part of interfacial charge transfer although its effect on the driving force and, more likely, on the reorganization energy may be far more dramatic, as demonstrated by Hyun et al.¹⁶ for PbS nanocrystals.

In the present work, we investigated the impact of molecular orbital spatial modification of an electron/hole acceptor molecule, as well as the size variation of a quantum dot on charge transfer dynamics at the QD/acceptor interface. Here, we applied our recently developed BA-LCAO approach to nanostructures and conventional DFT to the acceptor molecules. The presented electronic structure calculations, combined with extensive molecular dynamics simulations, revealed intimate details of charge transfer at the QD/acceptor interface. We used spherical wurtzite-(CdSe)₂₀₁ and (CdSe)₆₉₃ nanostructures (with ∼2.8 and ∼4.1 nm diameters) as model QDs, and anthraquinone-2,3-dicarboxylic acid (AQ) and its derivatives X-AQ with X = 7-OH, 7-OF, 10-BH, and 10-CH₂ substituents as model molecular acceptors. Briefly, we found the following:

1. Introducing both the electron donating (10-BH and 10-CH₂) and withdrawing (7-OH and 7-OF) substituents greatly enhances, roughly by a few orders of magnitude, HT by means of extending the HOMO of the acceptor molecule (i.e., AQ) in the dimension orthogonal to the molecular plane.

2. Electron transfer is affected only by the electron donating groups; in the case of CH_2, the ET probability is only moderately enhanced (∼1-2 orders of magnitude), while in the case of BH, it is slightly reduced. Analyses show that both cases correlate with the extended and contracted nature of LUMOs, respectively, in the dimension orthogonal to the AQ plane.

3. Solvent effects are largely negligible for the orbital overlaps and the coupling strength although MD simulations with explicit solvent molecules would help elucidate this question in a more consistent manner.

4. Consistent with spatial confinement theories, the electron density of the smaller dot has a significantly larger value at the surface and penetrates farther into the vacuum than the corresponding density of the larger dot. Thus, by altering the QD size, one can tune the electron’s wavefunction amplitude outside the dot.

These findings suggest that (a) one can effectively control details of charge transfer in the QD/molecule interface either by changing the size of the QD or by tuning diffuseness of frontier orbitals of the acceptor molecule; and (b) the combination of recently developed BA-LCAO approach, for nanostructures, with a DFT, for the acceptor molecules, and the fragmented molecular mechanics method, for the extensive molecular dynamics simulations, provide qualitatively accurate description of charge transfer dynamics at the QD/acceptor interface.

See supplementary material for XYZ structures of the X-AQ, representative MD trajectory snapshots, overlap as functions of time, and numerical details of calculations.

The FORTRAN codes used for performing these calculations can be downloaded from www.emerson.emory.edu/research/staff.html.

This work was funded by the U.S. Department of Energy, Office of Basic Energy Sciences, Solar Photochemistry Program (Grant No. DE-FG02-07ER-15906). We also gratefully acknowledge NSF MRI-R2 Grant No. (CHE-0958205) and the use of the resources of the Cherry Emerson Center for Scientific Computation.