Absorption Spectra for Disordered Aggregates of Chromophores

Similar agreement is obtained for a hybrid functional (B3LYP), a long-range ..... by quantum mechanics is often necessary for computing accurate prope...
0 downloads 0 Views 3MB Size
Article pubs.acs.org/JCTC

Absorption Spectra for Disordered Aggregates of Chromophores Using the Exciton Model Aleksey A. Kocherzhenko, Xochitl A. Sosa Vazquez, Joel M. Milanese, and Christine M. Isborn* Department of Chemistry and Chemical Biology, University of California, Merced, California 95343, United States S Supporting Information *

ABSTRACT: Optimizing the optical properties of large chromophore aggregates and molecular solids for applications in photovoltaics and nonlinear optics is an outstanding challenge. It requires efficient and reliable computational models that must be validated against accurate theoretical methods. We show that linear absorption spectra calculated using the molecular exciton model agree well with spectra calculated using time-dependent density functional theory and configuration interaction singles for aggregates of strongly polar chromophores. Similar agreement is obtained for a hybrid functional (B3LYP), a long-range corrected hybrid functional (ωB97X), and configuration interaction singles. Accounting for the electrostatic environment of individual chromophores in the parametrization of the exciton model with the inclusion of atomic point charges significantly improves the agreement of the resulting spectra with those calculated using all-electron methods; different charge definitions (Mulliken and ChelpG) yield similar results. We find that there is a size-dependent error in the exciton model compared with all-electron methods, but for aggregates with more than six chromophores, the errors change slowly with the number of chromophores in the aggregate. Our results validate the use of the molecular exciton model for predicting the absorption spectra of bulk molecular solids; its formalism also allows straightforward extension to calculations of nonlinear optical response.

1. INTRODUCTION Modeling of the linear and nonlinear optical response of molecular materials offers insights into the optimization of the structure and properties of these materials for photovoltaic,1 photosynthetic,2 and sensing applications.3 However, accurately predicting the optical response of disordered chromophore assemblies is challenging: neither ab initio cluster calculations nor band theory is appropriate. Ab initio all-electron quantumchemical calculations are useful for predicting the optical properties of individual chromophores, but they are too computationally expensive to be practical for aggregates consisting of more than a few chromophores. Band theory only predicts the properties of idealized periodic crystals and provides no insight into the effects of energetic and structural disorder in a material (unless the disorder is weak enough to be treated perturbatively).4 An alternative method is necessary for modeling large numbers of disordered chromophores. Numerical simulations based on the molecular exciton model5,6 are often successfully used to calculate the linear and nonlinear optical response of chromophore assemblies in complex molecular environments, e.g., in photosynthetic complexes.7−12 A multiscale modeling approach is employed, with each chromophore in the system represented by one or a few Frenkel exciton states that are often parametrized using ab initio calculations. The number of excitonic states per chromophore that need to be included in the calculation is much smaller than the number of atomic basis functions © 2017 American Chemical Society

needed to describe a chromophore. Therefore, unlike ab initio all-electron methods, calculations based on the exciton model are efficient enough to use on a large assembly consisting of thousands of chromophores. The underlying approximations of the molecular exciton model are valid only for weak chromophore−chromophore interactions.6 Many implementations of the exciton model consider only excitonic transitions and neglect charge-transfer transitions between chromophores. However, including intraor intermolecular charge-transfer states in the exciton model13−15 can be essential for reproducing experimentally observed optical and charge-transfer properties.16,17 Even in cases where including only excitonic states is sufficient, it may be necessary to account for interchromophore interactions in the parametrization of the excitonic Hamiltonian of an aggregate.18 Relatively weak interchromophore interactions in photosynthetic complexes ensure the applicability of the molecular exciton model to such systems. However, this theory may fail to reliably predict the properties of materials that have stronger interchromophore interactions. Nonlinear optical materials find applications in modulation of optical signals, sensing, and imaging.19−21 Organic chromophores exhibiting a strong nonlinear optical response have large hyperpolarizabilities. Such chromophores usually Received: May 8, 2017 Published: July 7, 2017 3787

DOI: 10.1021/acs.jctc.7b00477 J. Chem. Theory Comput. 2017, 13, 3787−3801

Article

Journal of Chemical Theory and Computation feature extended π-conjugated systems that allow ultrafast electron displacement induced by an applied electric field. They also often include electron-donating and electron-accepting groups and have significant ground- and excited-state dipole moments.22−26 A prototypical charge-transfer chromophore with a large ground-state dipole moment (∼30 D) and a large molecular hyperpolarizability is the CLD-type27 optical chromophore 2-[4-[(E,3E)-3-[3-[(E)-2-[4-[bis[2-[tert-butyl(dimethyl)silyl]oxyethyl]amino]phenyl]ethenyl]-5,5-dimethylcyclohex-2-en-1-ylidene]prop-1-enyl]-3-cyano-5-phenyl-5(trifluoromethyl)furan-2-ylidene]propanedinitrile (YLD124),28 shown in Figure 1A.

materials, thus serves as a challenging test for the molecular exciton model. To test the applicability of the molecular exciton model to calculations on chromophore aggregates, predictions of this theory must be benchmarked against predictions of ab initio calculations that account for interchromophore polarization as well as inter- and intramolecular charge transfer in the ground and excited states. Using time-dependent density functional theory (TDDFT), it is computationally feasible to compute excited-state properties for small chromophore aggregates, but care must be taken to avoid large energetic errors when considering charge-transfer transitions.33−35 The energies of charge-transfer transitions can be modeled accurately via longrange-corrected hybrid density functionals36,37 that include the correct electron−hole Coulombic interaction at long range.38−40 In contrast, pure DFT functionals lack this interaction, and TDDFT underestimates the energies of charge-transfer transitions compared with excitonic transitions. This error is remediated to some extent in global hybrid functionals that contain some fraction of exact exchange.33,34 Configuration interaction singles (CIS) calculations have full exact exchange and no correlation and are biased in the opposite direction, predicting charge-transfer transitions to be too high in energy compared with excitonic transitions.41 Because of their differing treatments of charge transfer compared with excitonic transitions, we might expect a hybrid TDDFT method, a long-range-corrected TDDFT method, and the CIS method to have varying agreement with molecular exciton theory that neglects charge-transfer transitions. Because of the large size of nonlinear optical chromophores such as YLD124, most electronic structure programs using hybrid TDDFT or CIS can only handle calculations on small aggregates (up to five chromophores). However, the TeraChem program42−44 implements algorithms optimized for massively parallel GPUs that enable large-scale ab initio calculations, allowing the computation of the absorption spectra of larger aggregates. The exciton model has previously shown good spectral agreement with all-electron TDDFT calculations using TeraChem for aggregates of six bacteriochlorophyll a chromophores.45 However, the exciton model is valid only in the limit of weak interchromophore interactions.5,6 Therefore, before the model is applied to aggregates of more polar chromophores with smaller interchromophore separations, it is necessary to verify that the good agreement with all-electron methods holds for such systems. In order to confidently predict the spectra of bulk molecular systems, it is also important to study how this agreement varies with the number of chromophores in the aggregate. In this paper, we investigate the applicability of the exciton model to calculations of linear absorption spectra for aggregates of the strongly polar chromophore YLD12428 selected from a single snapshot of an equilibrated coarse-grained Monte Carlo simulation46 of a bulk molecular solid. Examples of such aggregates are shown in Figure 1B,C. The aggregate spectra calculated using the exciton model are benchmarked against the aggregate spectra computed with TeraChem using CIS and TDDFT with a global hybrid functional (B3LYP) and a longrange corrected hybrid functional (wB97X). We find good agreement between the spectra calculated using the exciton model and using ab initio methods for aggregates with all numbers of chromophores considered, although we observe size-dependent errors in the exciton model. We also study how the agreement can be further improved by accounting for the

Figure 1. (A) Chemical structure of the YLD124 chromophore. (B, C) Sample chromophore aggregates: (B) dimer; (C) hexamer.

Generally, optoelectronic devices are built using thin films or bulk molecular materials, which consist of many individual chromophores. The electro-optical response is enhanced by increasing the chromophore alignment and the number density of the chromophores.29−31 However, chromophore−chromophore electrostatic interactions make the alignment of chromophores with large dipole moments very energetically unfavorable at high number densities.32 Therefore, bulk molecular materials constructed from polar molecules tend to be disordered and cannot be adequately described using band theory. Furthermore, polarization of the electron density of individual chromophores by the surrounding chromophores with large dipole moments is likely to be substantial. The electrostatic environment of the chromophores affects the optical properties of the sample and thus needs to be adequately taken into account by any computational model. The molecular exciton model has the potential to provide guidance for developing new materials that consist of polar chromophores. However, this model may not be applicable to materials where high number densities of polar chromophores could lead to strong electrostatic interactions and strong interchromophore coupling as well as to interchromophore charge transfer.6,17 The ability to predict absorption spectra for aggregates of chromophores with large dipole moments at high number densities, such as those used in organic electro-optic 3788

DOI: 10.1021/acs.jctc.7b00477 J. Chem. Theory Comput. 2017, 13, 3787−3801

Article

Journal of Chemical Theory and Computation M

electrostatic environment of chromophores in the parametrization of excitonic states.

Ĥ =

N

∑ ∑ εiv aiv̂ † aiv̂ i

i

i

vi = 1 i = 1 M

2. METHODS 2.1. Choice of Chromophore Aggregates. We compared linear optical absorption spectra calculated using the exciton model to spectra calculated using TDDFT and CIS for aggregates consisting of N = 2−10 YLD124 chromophores.28 The aggregates were selected from a single snapshot (see Figure S1 in the Supporting Information) of a coarse-grained Monte Carlo simulation46 of a bulk molecular solid (refractive index n = 1.7; dielectric constant ϵ ≈ 3) performed at 1 atm and 300 K with 108 YLD124 chromophores in a cubic simulation cell and periodic boundary conditions. Aggregates with various numbers of chromophores were obtained by selecting from the snapshot N chromophores with centers of mass closest to 27 grid points, (0.25Lnx , 0.25Lny, 0.25Lnz), with periodic boundary conditions taken into account, where (0, 0, 0) are the coordinates of the center, L = 60.448 Å is the length of the side of the cubic cell, and nx, ny, nz = 0, ±1. The two grid points with the largest distance between chromophore centers for N = 2 were discarded in order to retain only aggregates with significant interchromophore interactions. For any two of the aggregates with N = 10 chromophores, on average only 4.8% of chromophores belonged to both aggregates, with a standard deviation of 6.2%. Thus, all of the aggregates can be considered independent. 2.2. Calculation of Linear Optical Absorption Spectra Using the Exciton Model. We calculated the optical absorption spectra for disordered aggregates using the standard exciton model, which is valid in the limit of weak interchromophore interactions.5 In the exciton model for a system of N chromophores with M excited states on each chromophore a basis can be constructed from the ground and excited states of individual chromophores as follows:

M

N

i−1

∑ ∑ ∑ ∑ bijvv (aiv̂ † + aiv̂ )(ajv†̂ i j

+

i

i

j

vi = 1 vj = 1 i = 1 j = 1

+ ajv̂ j) (2)

εvi i

where is the vith excitation energy for the ith chromophore, bvijivj is the excitonic coupling between the vith excitation on the ith chromophore and the v j th excitation on the jth chromophore, and â†ivi and âivi are creation and annihilation operators for the vith excitation on the ith chromophore: â†ivi|0⟩i = |vi⟩i , âivi|0⟩i = 0. The Hamiltonian given by eq 2 neglects charge transfer between chromophores in both the ground and excited states; in other words, for all of its eigenstates, the total charge of each chromophore is equal to its ionization state. The dimensions of the Hamiltonian given by eq 2 are (M + 1)N × (M + 1)N, so its diagonalization quickly becomes computationally intractable as N and M increase. To enable calculation of the linear absorption spectrum, we simplify this Hamiltonian by making a number of approximations. First, we assumed that the excitation energies are considerably greater than absolute values of the interchromophore excitonic couplings (εvi i, εvj j ≫ bvijivj ∀ i, j, vi , vj), which is typically the case in organic materials.48,49 Second, because under typical experimental conditions the spatial density of excitations is relatively low, we assumed that interaction between excitons is negligible. Third, we assumed that the lowest excited bright state for each YLD124 chromophore is well-separated in energy from its higher excited states. Consequently, the lowest excited states of all of the chromophores in the aggregate form a distinct excitonic band in the spectrum in the limit of a bulk solid. We restricted our calculations to excitations within this band. Under these conditions, the manifold formed by the lowest excitations on individual chromophores can be considered separately, and the Hamiltonian given by eq 2 reduces to5,6,48 N

Ĥ =

N

i−1

∑ εiaî †aî + ∑ ∑ bij(aî †aĵ + aî aj†̂ ) i=1

(3)

i=1 j=1

where we have discarded the indices vi = vj = 1 (e.g., âj = âj1). The excitation energies for this tight-binding Hamiltonian are size-intensive by construction. Equation 3 is the standard tight-binding Hamiltonian in the Heitler−London approximation:48,50 its ground-state wave function is assumed to be |Ø⟩ (the product of the groundstate wave functions of single chromophores), and its excitedstate wave functions are of the form N

|Ψ⟩ i =

∑ cji|Φ⟩j (4)

j=1

⊗k N= 1|δjk⟩k

where i = 1, ..., N; |Φj⟩ = (in which the jth chromophore is in its lowest excited state, |1⟩j , and the rest of the chromophores are in their ground states, |0⟩k≠j , j, k = 1, ..., N); and the expansion coefficients cij for the ith eigenstate can be selected to be real. The oscillator strength for the transition between the ground state |Ø⟩ and the state |Ψi⟩ is defined as51

where |0⟩i and |vi⟩i are the ground and vith excited state of the ith chromophore, respectively, and δij is the Kronecker delta. Thus, in the state |Φvi i⟩, the ith chromophore is in its vith excited state, |vi⟩i , and all of the other chromophores are in their ground states, |0⟩j≠i; in the state |Φvijivj⟩, the ith and jth chromophores are in their vith and vjth excited states, |vi⟩i and |vj⟩j , respectively, and all of the other chromophores are in their ground states; etc. We constructed the excitonic Hamiltonian as follows:5,6,47,48

fØ → i = 3789

2m Ei ⟨Ψ|i μ̂ |Ø⟩ 2 3e 2ℏ2

(5) DOI: 10.1021/acs.jctc.7b00477 J. Chem. Theory Comput. 2017, 13, 3787−3801

Article

Journal of Chemical Theory and Computation where m is the electron mass, ℏ is the reduced Planck’s constant, e is the elementary charge, Ei is the eigenvalue that corresponds to the eigenvector |Ψi⟩, and μ̂ is the transition dipole operator that couples these states. Since interaction with an external electric field can produce an excitation on any of the chromophores in the system, the transition dipole operator can be written as

tridiagonal form using Householder reduction; its eigenvalues and eigenvectors are found using the QL algorithm with implicit shifts.53 The oscillator strengths for transitions from the ground state to each of the excited states are then calculated using eq 8. 2.3. Parametrization of the Exciton Model Hamiltonian Based on TDDFT and CIS Calculations. Optical absorption spectra of YLD124 chromophore aggregates calculated using the molecular exciton model were compared to spectra of the same aggregates obtained from ab initio calculations. Three ab initio methods were used: TDDFT with the B3LYP global hybrid functional (default exact exchange value of cx = 0.2),54 TDDFT with the long-range-corrected ωB97X functional (default parameters, range-separation parameter ω = 0.3 bohr−1, short-range exact exchange cx = 0.157706),55 and CIS.56 These three methods were selected because of their differing treatments of charge-transfer excitations. Because of the lack of full exact exchange at long range, TD-B3LYP underestimates charge-transfer excitation energies, whereas TD-ωB97X predicts charge-transfer excitation energies more accurately. CIS overestimates chargetransfer excitation energies relative to excitonic excitation energies.41 Consequently, spectra calculated using the molecular exciton model that neglects charge-transfer transitions could be expected to agree better with spectra calculated using CIS than with those calculated using TD-ωB97X or TDB3LYP. The TeraChem program was needed for calculations on the larger chromophore aggregates, but parametrization of the molecular exciton Hamiltonian only requires calculations on single chromophores and can be performed using any quantum-chemistry program. We performed this parametrization using the popular Gaussian 09 package.57 The excitation energies εi and the transition dipole moments μi were calculated using CIS and the Tamm−Dancoff approximation to TDDFT.58 The 6-31G* basis set was used for all of the calculations, except as noted in the text.59,60 The “Ultrafine” integration grid was used for all of the Gaussian TDDFT calculations. This ensured that the lowest excitation energies and corresponding oscillator strengths for single chromophores from the Gaussian calculations agreed with the values from TeraChem calculations with the wave function convergence threshold set to 10−8, the two-electron integral threshold set to 10−15, and the finest DFT integration grid (the “5” setting) up to the hundredth place. It should be noted that transition dipole moments μi are defined up to an arbitrary phase factor. Because Gaussian selects μi to be real, transition dipole moments for calculations on individual chromophores are defined up to their sign. To obtain correct interchromophore interactions in an aggregate, the signs of μi for all of the chromophores must be selected consistently. We ensured that this condition was satisfied by requiring that the orientation of the transition dipole moment of a chromophore relative to its ground-state dipole moment be similar for all chromophores in the aggregate. The excitonic couplings bij were calculated using the transition density cube method.61,62 This method captures only the Coulomb contribution to interchromophore couplings, neglecting the effects of electron exchange and correlation. However, at small interchromophore distances where these effects could become pronounced, they are typically overshadowed by the effect of interchromophore charge transfer.63 Thus, the approximation of purely Coulomb

N

∑ (exμkx + eyμky + ezμkz )(ak̂ † + ak̂ )

μ̂ =

k=1

(6)

μαk

where (α = x, y, z) is the component of the kth chromophore’s transition dipole moment in the α direction, â†k (âk) are creation (annihilation) operators for an excitation on the kth chromophore, and eα is the unit vector in the α direction. Introducing eqs 4 and 6 into eq 5 and taking into account that â†k |Ø⟩ = |Φk⟩ and âk|Ø⟩ = 0, we obtain the following expression for the oscillator strength: fØ → i =

2m Ei 3e 2ℏ2 ⎡N ⎤2 N ⎢∑ (e μ x + e μ y + e μ z ) ∑ c i⟨Φ|Φ ⟩⎥ x k y k z k j j k ⎢⎣ ⎥⎦ k=1 j=1

(7)

According to our definition of the states |Φj⟩ and |Φk⟩, ⟨Φj|Φk⟩ = δjk , so eq 7 becomes fØ → i

⎡ N ⎞2 ⎛ N ⎞2 2m ⎢⎛ i x⎟ i y⎟ ⎜ ⎜ = 2 2 Ei ⎜∑ ckμk ⎟ + ⎜∑ ckμk ⎟ 3e ℏ ⎢⎣⎝ k = 1 ⎠ ⎝k=1 ⎠ ⎛N ⎞2 ⎤ i z + ⎜⎜∑ ckμk ⎟⎟ ⎥ ⎝k=1 ⎠ ⎥⎦

(8)

For a single chromophore with transition dipole moment (μx , μy, μz) and excitation energy ε, the expression for the oscillator strength given by eq 8 reduces to f|0⟩→|1⟩ =

2m ε(μ 2 + μy 2 + μz 2 ) 3e 2ℏ2 x

(9)

Calculating the optical absorption spectrum for an aggregate of chromophores requires three types of parameters: (1) the excitation energies εi for all of the chromophores that constitute the aggregate, (2) the excitonic couplings bij between all pairs of chromophores, and (3) the transition dipole moments μi = (μxi , μyi , μzi ) for all of the chromophores. The excitation energies and excitonic couplings are used to construct the Hamiltonian given by eq 3. Diagonalizing this Hamiltonian produces the excitation energies of the aggregate, Ei , and the eigenstate expansion coefficients, cij , from eq 4. These quantities, together with the transition dipole moments μi for all of the chromophores, are then used to calculate the oscillator strengths of all excitations in the aggregate according to eq 8. For two coupled chromophores, we solve the timeindependent Schrödinger equation analytically. This solution can be found in the Supporting Information, with examples of dimers in several orientations shown in Figure S2. Although general analytical solutions also exist for clusters of chromophores with N ≤ 4,52 for calculations involving more than two chromophores we diagonalize the Hamiltonian numerically. The Hamiltonian matrix is transformed to 3790

DOI: 10.1021/acs.jctc.7b00477 J. Chem. Theory Comput. 2017, 13, 3787−3801

Article

Journal of Chemical Theory and Computation couplings is appropriate within the applicability range of the standard exciton model. Transition density matrices from Gaussian calculations were converted to average transition densities for small cubic spatial regions that spanned the rectangular volume encompassing each chromophore using the Multiwfn program64 with a coarse grid (the “1” setting). The couplings were then calculated according to61 bij =

1 2

∑∑ k

l

MkiMlj 4π ϵ0rklij

(10)

where Mki and Mlj are the average transition densities for the kth cubic volume for the ith chromophore and the lth cubic volume for the jth chromophore, respectively, rijkl is the distance between the centers of these cubic volumes, the summations run over all of the cubic volumes into which the spatial regions of both chromophores are fragmented, ϵ0 is the vacuum permittivity, and the factor of 1/2 is necessary to prevent double counting of interactions. In section 3.1 we compare absorption spectra calculated using all-electron methods to those calculated using the exciton model parametrized by all-electron calculations on individual chromophores in vacuum. In section 3.2 the electrostatic effects of surrounding polar chromophores are included in the parametrization of the exciton-model Hamiltonian by electrostatic embedding of partial atomic point charges in a quantummechanical/molecular-mechanical (QM/MM) calculation.65 The chromophore for which parameters are calculated is described quantum-mechanically, with the rest of the chromophores in the aggregate represented by either Mulliken66,67 or ChelpG68 atomic point charges calculated for the ground states of individual chromophores. 2.4. Calculations of Linear Optical Absorption Spectra Using All-Electron TDDFT and CIS. For small clusters of up to 10 chromophores, we compared the spectra calculated using the exciton model to spectra found from calculations using the Tamm−Dancoff approximation to TDDFT and CIS. These calculations were performed using the TeraChem package, which enables the analysis of the electronic structure of large molecular systems by exploiting algorithms optimized for massively parallel GPUs.42−44 Calculations on the aggregates employed the same all-electron methods (TD-B3LYP, TDωB97X, and CIS) and the same basis set (6-31G*) as used to parametrize the model Hamiltonian. 2.5. Metrics for Spectral Comparisons. The absorption spectra were computed as a sum over Gaussian functions centered at the excitation energies Ei and weighted by the corresponding oscillator strengths fØ→i: N

f (E ) =

∑ i=1

⎡ (E − Ei) ⎤ exp⎢ − ⎥ ⎣ σ 2 2σ 2 ⎦

Figure 2. Sample hexamer optical absorption spectra calculated at the TD-ωB97X/6-31G* level (magenta) and with the exciton model parametrized using ChelpG charges to represent the environment of each chromophore (blue). The area under the spectrum is its integrated intensity. The dashed lines indicate the mean absorption energy; ΔE is the relative shift in the mean absorption energy between two spectra. The Pearson correlation coefficient for two spectra, r, is equal to 1 if all of the absorption maxima and minima are at the same energy values for both spectra. The Pearson correlation coefficient for the derivatives of two spectra, r′, is equal to 1 if all of the shoulders are within the same energy range for both spectra.

We compared the exciton model and TDDFT or CIS spectra (Figure 2) using the shift in the mean energy, ΔE = Eexc − Eall‑electron, and the integrated intensity ratio, Iexc/Iall‑electron. To quantify the similarities in the shapes of the exciton-model and all-electron spectra, we also calculated the Pearson productmoment correlation coefficients:69 r=

r′ =

(11)

(12)

+∞

∫−∞

f (E ) d E

dfall‐electron (E) ⎤ dE

⎦⎥

⎡ dfexc (E − ΔE) ⎤

⎡ dfall‐electron (E) ⎤

dE

dE

⎦⎥σ ⎣⎢

⎦⎥

(15)

3. RESULTS AND DISCUSSION 3.1. Validity of the Exciton Model for Calculating Optical Absorption Spectra of Chromophore Aggregates. We first establish the applicability of the exciton model to calculations of the optical absorption spectra of aggregates of

To simplify notation, hereafter we omit the overline from the mean energy (E̅ → E). The integrated intensity is I=

(14)

where the numerators are the covariances and the denominators are the products of standard deviations of the exciton model and all-electron spectra (eq 14) or of their first derivatives with respect to energy (eq 15). The Pearson correlation coefficients r and r′ were calculated for the excitonmodel spectrum shifted to have the same mean energy as the all-electron spectrum. A value of 1 for the Pearson correlation coefficient indicates that two functions are perfectly correlated, a value of 0 indicates no correlation, and a value of −1 indicates perfect anticorrelation. For the spectra shown in Figure 2, the ratio of the integrated intensities of the spectra calculated using the exciton model and TDDFT is Iexc/ITDDFT = 1.124, and the difference in the mean absorption energies is Eexc − ETDDFT = 0.094 eV. The Pearson correlation coefficient for the spectra is rexc,TDDFT = 0.982, and the Pearson correlation coefficient for the derivatives of the spectra is rexc,TDDFT ′ = 0.932.

fØ → i

Ef (E) dE

⎡ df (E − ΔE) cov⎣⎢ exc dE , σ ⎣⎢

+∞

∫−∞

σ[fexc (E − ΔE)]σ[fall‐electron (E)]

and

for E ∈ [1.5, 3.0) eV with the broadening parameter σ = 0.08 eV and an energy discretization step of 0.001 eV. Sample spectra are shown in Figure 2. The mean absorption energy is

E̅ =

cov[fexc (E − ΔE), fall‐electron (E)]

(13) 3791

DOI: 10.1021/acs.jctc.7b00477 J. Chem. Theory Comput. 2017, 13, 3787−3801

Article

Journal of Chemical Theory and Computation

predict charge-transfer states with significantly higher energies than those of Frenkel excitons. Thus, charge-transfer states are not expected to significantly affect the aggregate spectra for these methods. However, calculations using TD-B3LYP predict large numbers of dark states that have charge-transfer character and occur at unphysically low energies.34 For 25 aggregates consisting of N = 2 chromophores (dimers) and 25 aggregates consisting of N = 4 chromophores (tetramers), there are 5−7 and 21−27 states with energies below 2.5 eV, respectively. Because the number of low-energy charge-transfer states predicted by TD-B3LYP calculations increases rapidly with the number of chromophores in the aggregate, calculations using this method were performed only for N ≤ 4. If the mixing of charge-transfer and Frenkel exciton states significantly alters the spectrum, one would expect considerably worse agreement of TD-B3LYP calculations with the exciton model compared with that seen for TD-ωB97X and CIS. For 25 dimers and 25 tetramers of YLD124 chromophores, we compared spectra predicted by the exciton model parametrized using TD-B3LYP, TD-ωB97X, and CIS to spectra calculated using the same methods on the full dimers or tetramers. The parameters of the molecular exciton Hamiltonian, eq 3, for dimers calculated using all three methods can be found in Figure S3. The four similarity metrics described in Methods are shown for dimer spectra in Figure 4. These metrics for tetramers exhibit features similar to those seen for dimers (Figure S4). From Figure 4A,B it is evident that most of the values of the Pearson correlation coefficients are close to 1, indicating good agreement in spectral shape for all of the methods. The Pearson correlation coefficients for spectra calculated using the exciton model and TDDFT/CIS on the dimer are greater than 0.99 for 18, 16, and 11 dimers for CIS, TD-ωB97X, and TD-B3LYP calculations, respectively. They are greater than 0.90 for 25 (CIS), 24 (TD-ωB97X), and 21 (TD-B3LYP) dimers. The Pearson correlation coefficients for the derivatives of spectra are greater than 0.99 for 14 (CIS), 13 (TD-ωB97X), and 7 (TDB3LYP) dimers and greater than 0.90 for 19 (CIS and TDωB97X) and 17 (TD-B3LYP) dimers. This agreement can be improved further by using a larger basis set that includes diffuse functions, such as 6-31+G* (Figure S5A,B). For both dimers and tetramers, the agreement in spectral shape between calculations using the exciton model and the all-electron methods is on average slightly better for methods with more exact exchange (Table 1). This may be the case because the Frenkel exciton model is less able to reproduce the more delocalized excitonic states predicted by methods with less exact exchange. However, these differences should not be considered significant because the average Pearson correlation coefficients both for the spectra and their derivatives calculated using all three methods are well within a standard deviation of one another (Table 1). One might also note that the Pearson correlation coefficients are not consistently larger for any of the three methods: for instance, although for most dimers the exciton model performs better with CIS and TD-ωB97X than with TD-B3LYP, the opposite is true for dimers #9 and #11 (Figure 4A,B). On average, the exciton model parametrized with any of the three methods consistently overestimates the integrated intensity of the spectra compared with calculations on dimers and tetramers performed with the same method (Table 1). This overestimation remains virtually the same when the basis set is increased to 6-31+G* (Figure S5C). It is greater for tetramers

YLD124 chromophores. The Frenkel exciton Hamiltonian with one excitation included per chromophore, eq 3, can be used within the following three approximations: (a) the first and higher excited states of individual chromophores are wellseparated in energy; (b) the couplings between chromophores are much smaller than the excitation energies; and (c) the presence of charge-transfer states in the system can be neglected. Let us now verify that these conditions hold for the system under consideration. In this section we neglect polarization of the electron density of chromophores by their electrostatic environment. The calculated lowest excitation energy of the YLD124 chromophore depends on the quantum-chemical method used and increases with the amount of exact exchange. For 87 YLD124 chromophores from which the aggregates were built, we calculated average excitation energies of 2.083 ± 0.012 eV (TD-B3LYP), 2.535 ± 0.022 eV (TD-ωB97X), and 2.833 ± 0.025 eV (CIS). The energy separation between the two lowest excited states of YLD124 was found to be 0.676 ± 0.027 eV (TD-B3LYP), 1.366 ± 0.028 eV (TD-ωB97X), and 1.541 ± 0.020 eV (CIS). The excitonic couplings calculated using the three methods have similar distributions (Figure 3). The mean

Figure 3. Probability density distributions for interchromophore couplings in 25 tetramers calculated using TD-B3LYP (magenta), TDωB97X (black), and CIS (orange).

coupling value for all three methods is 0.021 eV, with standard deviations of 0.021 eV for TD-B3LYP and 0.024 eV for TDωB97X and CIS. Despite the large transition dipole moments and the close proximity of YLD124 chromophores, the magnitude of the coupling between them is similar to that typically found in photosynthetic systems, where the exciton model is commonly used.7−12 For all three methods, the energy difference between the first and second excitations of the chromophore significantly exceeds the interchromophore couplings. Therefore, it should be reasonable to consider the band formed by the lowest excitations on single chromophores separately, neglecting any effects of higher excitonic states. Similarly to photosynthetic systems,8,9 the excitonic couplings are 2 orders of magnitude smaller than the excitation energies, which is within the limits of applicability of the molecular exciton model in the Heitler− London approximation.48 Therefore, conditions (a) and (b) for the applicability of the Frenkel exciton Hamiltonian, eq 3, are satisfied. To verify condition (c), we must study how the neglect of charge-transfer states in the exciton model affects the agreement of its predictions with those of all-electron quantum-chemical methods. Calculations on aggregates of YLD124 chromophores using TD-ωB97X and CIS correctly 3792

DOI: 10.1021/acs.jctc.7b00477 J. Chem. Theory Comput. 2017, 13, 3787−3801

Article

Journal of Chemical Theory and Computation

Figure 4. Similarity metrics for spectra of 25 chromophore dimers calculated with the exciton model parametrized using TD-B3LYP (circles, magenta), TD-ωB97X (squares, black), or CIS (triangles, orange) calculations with no point charges relative to aggregate spectra calculated with the corresponding quantum-chemical method. (A) Pearson correlation coefficients for spectra. (B) Pearson correlation coefficients for the first derivatives of spectra. (C) Ratios of integrated intensities of the exciton-model and all-electron spectra. (D) Shifts of the spectral means of the exciton-model spectra relative to the all-electron spectra.

Table 1. Mean Values and Standard Deviations of Pearson Correlation Coefficients of Spectra (r) and Their First Derivatives (r′), Ratios of Integrated Intensities (Iexc/Iall‑electron), and Absolute Mean Energy Shifts (|ΔE|) for Spectra of 25 Chromophore Dimers and 25 Chromophore Tetramers Calculated Using the Exciton Model Relative to Spectra of the Same Chromophore Aggregates Calculated Using TDDFT or CIS

method dimers

tetramers

TD-B3LYP TD-ωB97X CIS TD-B3LYP TD-ωB97X CIS

r ± σr 0.9636 0.9760 0.9868 0.9672 0.9874 0.9902

± ± ± ± ± ±

0.0536 0.0501 0.0232 0.0304 0.0137 0.0174

r′ ± σr′ 0.9166 0.9157 0.9547 0.9051 0.9386 0.9492

than for dimers (Table 1 and Figure S4) and increases further for larger chromophore aggregates. It suggests that dark chargetransfer states, neglected in the exciton model, may weakly contribute to Hamiltonian eigenstates of dimers and tetramers calculated using TDDFT or CIS. The integrated intensities calculated using CIS and TD-ωB97X are very similar (Figure 4C). The integrated intensities calculated using TD-B3LYP vary more than those calculated with the other two methods. This could be caused by a larger contribution of charge-transfer states, some higher and some lower in energy than the Frenkel exciton states, to the spectrum calculated with TD-B3LYP. Compared with all-electron aggregate calculations, the exciton model tends to slightly overestimate the mean excitation energy for dimers and tetramers (Table 1). This overestimation also does not significantly change when the basis set is increased to 6-31+G* (Figure S5D). We will see that this error can be reduced by accounting for the electrostatic environment of the chromophores. Thus, the error is partially due to the neglect of stabilization of the

± ± ± ± ± ±

0.1206 0.1538 0.0749 0.0724 0.0741 0.0851

⎛ Iexc ⎞ ⎟ ± σI ⎜ ⎝ Iall‐electron ⎠ 1.0402 1.0633 1.0600 1.0732 1.1074 1.1012

± ± ± ± ± ±

0.0645 0.0435 0.0417 0.0868 0.0449 0.0419

|ΔE| ± σΔE (eV) 0.0259 0.0341 0.0293 0.0352 0.0553 0.0517

± ± ± ± ± ±

0.0240 0.0267 0.0232 0.0224 0.0404 0.0358

electronically excited states by the environment. Like the integrated intensities, the mean energies are overestimated more for larger aggregates (Table 1 and Figure S4). The origin of this size dependence is likely similar in these two cases and may be related to the neglect of charge-transfer states in the exciton model. For all four metrics considered, we find good agreement between the spectra of YLD124 dimers and tetramers calculated using TDDFT or CIS and those calculated using the molecular exciton model. This agreement is similar regardless of the ab initio method used; differences between methods are not significant. Thus, we conclude that the neglect of charge-transfer states results only in small changes to the calculated spectra of chromophore aggregates for all of the methods considered. The use of the molecular exciton model is therefore justified for the system studied and for similar systems. In the next section we will discuss how the agreement between spectra calculated using the exciton model and allelectron methods can be improved by accounting for the 3793

DOI: 10.1021/acs.jctc.7b00477 J. Chem. Theory Comput. 2017, 13, 3787−3801

Article

Journal of Chemical Theory and Computation

Figure 5. Similarity metrics for spectra of 25 chromophore dimers calculated using the exciton model with no point charges (squares, black), with Mulliken point charges (triangles, green), or with ChelpG charges (circles, red) relative to spectra calculated using TD-ωB97X. (A) Pearson correlation coefficients for spectra. (B) Pearson correlation coefficients for the first derivatives of spectra. (C) Ratios of integrated intensities of the exciton-model and all-electron spectra. (D) Shifts of the spectral means of the exciton-model spectra relative to the all-electron spectra.

In this section, we compare the agreement of spectra calculated using the exciton model parametrized with no point charges, with Mulliken charges, and with ChelpG charges with spectra calculated using TD-ωB97X on the aggregate. In section 3.3 we will see that the effects of point charges are similar for B3LYP and CIS calculations. Atomic point charges were obtained from ground-state electronic structure calculations on individual chromophores. There are considerable differences between the values of the Mulliken charges (derived from Mulliken’s population analysis) and ChelpG charges (derived from fitting the electrostatic potential in the vicinity of the chromophore) for most atoms in a chromophore. However, the dipole moments that mostly determine the electrostatic potential far away from a chromophore are similar for the two sets of charges. Figure 5 shows the four similarity metrics for the spectra of chromophore dimers calculated using the exciton model relative to the spectra calculated using TD-ωB97X. For 14 of the 25 chromophore dimers considered, the spectral shapes calculated using the exciton model without the inclusion of point charges are very similar to those calculated using TDωB97X: the values of the Pearson correlation coefficients for both the spectra, r, and their first derivatives, r′, exceed 0.99. For two chromophore dimers, #2 and #16, inclusion of point charges leads to greatly improved agreement between the shapes of spectra calculated using the exciton model and spectra calculated using TD-ωB97X (Figure 5A,B). For chromophore dimer #2, inclusion of point charges increases r from 0.7699 to 0.9993 (Mulliken charges) or 0.9998 (ChelpG charges) and r′ from 0.3783 to 0.9974 (Mulliken charges) or 0.9993 (ChelpG charges). For chromophore dimer #16, inclusion of point charges increases r and r′ from 0.9191 and 0.7983, respectively, to over 0.9999 for both Mulliken and ChelpG charges. For four more chromophore dimers (#7, #15, #17, and #25) where r or r′ or both coefficients are below

electrostatic environment of individual chromophores in the parametrization of the exciton-model Hamiltonian. 3.2. Effects of the Electrostatic Environment on Spectra Calculated Using the Exciton Model. Including a quantum-mechanical description of the chromophore environment takes into account the polarization of the electron density and charge transfer between molecules; treatment of a large region of the environment by quantum mechanics is often necessary for computing accurate properties.65,70 All-electron TDDFT and CIS calculations on chromophore aggregates account for the polarization of the electronic density by the other chromophores in the aggregate as well as for charge transfer between chromophores. In calculations using the molecular exciton model we have so far neglected both polarization and ground-state charge transfer. Polarization can be approximately included in the exciton-model Hamiltonian, eq 3, by calculating the excitation energies and transition densities of individual chromophores in the presence of the electrostatic environment of all other chromophores via an electrostatic embedding hybrid QM/MM approach.71−74 The electrostatic environment of a chromophore can be approximated by the partial atomic charges of other chromophores in the aggregate. Partial atomic charges are not observables and can be defined in a variety of ways. Common definitions are based on population analysis of wave functions (e.g., Mulliken charges66,67), on partitioning of the electron density distribution, on fitting of the electrostatic potential (e.g., ChelpG charges68) or the electrostatic multipole moments corresponding to the electron density distribution, or on fitting spectroscopic or other experimental data.75 The accuracy of the computed spectra depends on the accuracy of the chromophore polarization. Therefore, it is important to determine how strongly the choice of a specific definition for partial atomic charges affects the calculated optical absorption spectrum. 3794

DOI: 10.1021/acs.jctc.7b00477 J. Chem. Theory Comput. 2017, 13, 3787−3801

Article

Journal of Chemical Theory and Computation

Table 2. Mean Values and Standard Deviations for Pearson Correlation Coefficients of Spectra (r) and Their First Derivatives (r′), Ratios of Integrated Intensities (Iexc/Iall‑electron), and Absolute Mean Energy Shifts (|ΔE|) for Spectra of 25 Chromophore Dimers, Tetramers, Hexamers, Octamers, and Decamers Calculated Using the Exciton Model Relative to Spectra of the Same Chromophore Aggregates Calculated Using TD-ωB97X

charges dimers

tetramers hexamers octamers decamers

none Mulliken ChelpG none ChelpG none ChelpG none ChelpG none ChelpG

r ± σr 0.9760 0.9916 0.9911 0.9874 0.9926 0.9776 0.9919 0.9695 0.9930 0.9653 0.9901

± ± ± ± ± ± ± ± ± ± ±

0.0501 0.0192 0.0203 0.0137 0.0077 0.0200 0.0090 0.0362 0.0085 0.0474 0.0128

r′ ± σr′ 0.9157 0.9611 0.9587 0.9386 0.9665 0.8806 0.9577 0.8428 0.9619 0.8188 0.9407

± ± ± ± ± ± ± ± ± ± ±

0.1538 0.0874 0.0932 0.0741 0.0354 0.1168 0.0448 0.1703 0.0455 0.1850 0.0650

⎛ Iexc ⎞ ⎟ ± σI ⎜ ⎝ Iall‐electron ⎠ 1.0633 1.0633 1.0635 1.1074 1.1058 1.1278 1.1258 1.1326 1.1314 1.1395 1.1384

± ± ± ± ± ± ± ± ± ± ±

0.0435 0.0460 0.0460 0.0449 0.0429 0.0489 0.0481 0.0463 0.0467 0.0392 0.0398

|ΔE| ± σΔE (eV) 0.0341 0.0227 0.0236 0.0553 0.0448 0.0717 0.0574 0.0750 0.0618 0.0797 0.0662

± ± ± ± ± ± ± ± ± ± ±

0.0267 0.0126 0.0129 0.0404 0.0191 0.0417 0.0173 0.0346 0.0180 0.0391 0.0185

Figure 6. Similarity metrics for spectra of 25 chromophore octamers calculated using the exciton model with no point charges (squares, black) and with ChelpG charges (circles, red) relative to spectra calculated using TD-ωB97X. (A) Pearson correlation coefficients for spectra. (B) Pearson correlation coefficients for the first derivatives of spectra. (C) Ratios of integrated intensities of the exciton-model and all-electron spectra. (D) Shifts of the spectral means of the exciton-model spectra relative to the all-electron spectra.

the chromophore dimers, except for dimer #10, where the ratio of these intensities is 0.9949. Inclusion of point charges in the parametrization of the exciton-model Hamiltonian does not significantly change the ratio of integrated intensities of the exciton model and TD-ωB97X spectra for any of the chromophore dimers (Figure 5C and Table 2). Therefore, the differences in the integrated intensities calculated using the exciton model and TD-ωB97X are likely due to the neglect of charge-transfer states in the exciton model rather than to the polarization of the electron density by the environment. Figure 5D shows that when no charges are included in the parametrization of the exciton-model Hamiltonian, the mean absorption energy for 17 out of 25 chromophore dimers is above that found from TD-ωB97X calculations by an average of 0.0434 eV. For the remaining eight chromophore dimers, the

0.9900 when no point charges are included, the values of both Pearson correlation coefficients exceed 0.9900 upon inclusion of either Mulliken or ChelpG point charges. For other chromophore dimers, inclusion of Mulliken or ChelpG point charges leads to either slight improvements or to no significant changes in the spectral shape. For instance, for dimer #9, a slight improvement is observed: r increases from 0.9107 to 0.9249 (Mulliken charges) or 0.9209 (ChelpG charges), and r′ increases from 0.6231 to 0.6866 (Mulliken charges) or 0.6670 (ChelpG charges). On average, inclusion of point charges, whether Mulliken or ChelpG, increases the values of r and r′ and decreases their standard deviations (see Table 2). Figure 5C shows that the integrated intensities of spectra calculated using the exciton model are larger than the integrated intensities calculated using TD-ωB97X for all of 3795

DOI: 10.1021/acs.jctc.7b00477 J. Chem. Theory Comput. 2017, 13, 3787−3801

Article

Journal of Chemical Theory and Computation

estimates the mean excitation energy for octamers relative to TD-ωB97X, although this is mitigated to a certain extent when atomic point charges are included in the calculation (Figure 6D). Thus, we find that accounting for individual chromophore polarization by other chromophores in the environment via inclusion of atomic point charges in the parametrization of the exciton-model Hamiltonian improves the agreement with the all-electron results. We next examine whether allowing charge transfer between chromophores in the environment further improves this agreement. In our approach, charges are calculated separately for each chromophore in the aggregate, and therefore, individual chromophores have no net charge. However, charge delocalization among chromophores in their ground states could change the absorption spectra of chromophore aggregates. For 25 YLD124 hexamers, we calculated the excitation energies, transition densities, and transition dipole moments for each of the monomers that constitute the hexamer with the electrostatic environment represented by Mulliken point charges obtained from a TD-ωB97X calculation on the other five monomers. In these calculations, the monomer for which parameter calculations were performed was neutral, but charge delocalization across the other five monomers was allowed. The spectra of Hamiltonians parametrized using this procedure show little variation compared to the spectra of Hamiltonians parametrized using point charges calculated for individual chromophores to represent the electrostatic environment (see Figure S8). The negligible electrostatic effects of interchromophore charge delocalization in the ground state on the absorption spectra may be attributed to the small extent of this delocalization. Despite the close proximity of chromophores in the aggregate, the mean absolute chromophore charge for the chromophores comprising 25 hexamers is only 0.034e per chromophore, with a standard deviation of 0.027e. We have seen that accounting for the electrostatic environment of chromophores in the parametrization of the excitonmodel Hamiltonian of chromophore aggregates improves the agreement of the resulting absorption spectra with those obtained from TD-ωB97X calculations on entire aggregates. Electrostatic polarization of the chromophores by their environment significantly affects the shape of the absorption spectrum. It also changes the mean absorption energy but has little effect on the absorption intensity. With the methods used here, the choice of definition for atomic point charges that represent the electrostatic environment of the chromophore makes little difference for the spectra calculated using the exciton model. Different definitions of point charges result in very different electrostatic potentials in the direct vicinity of a molecule but in similar values of the ground-state molecular dipole moment. Therefore, we conclude that chromophores are sufficiently separated to neglect the effects of the higher terms in the multipole expansions of their electrostatic potentials on the aggregate spectra. In the following section, for calculations based on the exciton model we will only present results with no point charges and with ChelpG point charges. 3.3. Size-Dependent Errors in the Exciton Model. The molecular exciton model can be used to calculate the optoelectronic properties of very large aggregates: up to several thousand chromophores. The optical response for aggregates with these numbers of chromophores is expected to converge to that of the bulk molecular material. However, before the molecular exciton model is applied to very large systems, one must verify that the good agreement between spectra calculated

mean absorption energy is below that found from TD-ωB97X calculations by an average of 0.0144 eV. Interestingly, when point charges are included, the mean absorption energy is overestimated for all 25 chromophore dimers by an average of 0.0227 eV (Mulliken charges) or 0.0236 eV (ChelpG charges). However, the absolute value of the difference between the mean values of the spectra calculated using TD-ωB97X and the exciton model, |ΔE|, decreases for 17 out of 25 chromophore dimers when point charges are included, by an average of 0.0229 eV for Mulliken charges and 0.0207 eV for ChelpG charges (Figure 5D). For the other eight chromophore dimers, including point charges increases |ΔE|, but only by an average of 0.015 eV for Mulliken charges and 0.013 eV for ChelpG charges. The mean value of |ΔE| and its standard deviation both decrease when point charges are included in the calculation of model Hamiltonian parameters (Table 2). The improved agreement between the mean excitation energies calculated using the exciton model and TD-ωB97X suggests the importance of accounting for the polarization of the electron density on the chromophores by their environment. Although the polarization of a chromophore can be described by the inclusion of the atomic point charges of surrounding chromophores in the parametrization of the exciton model, in this procedure the polarization of the chromophore and of its environment is not self-consistent. Insight into the significance of such self-consistency can be gained from calculations based on the exciton model parametrized using monomer calculations with a polarizable continuum model (PCM).76 Using a PCM does not improve the agreement with TDDFT in the spectral shape compared with the exciton model parametrized with no pont charges (Figure S6A,B). However, it does lower both the integrated intensity and the mean excitation energy predicted by the exciton model (Figure S6C,D). Therefore, one may conclude that the spectral shape is primarily determined by local electrostatic interactions between chromophores in the aggregate, whereas the excitation energy and intensity depend considerably on the mutual polarization between the chromophore and its environment. It should be noted that in contrast to the case of bulk solids, using a PCM to describe polarization in small chromophore aggregates is not physically reasonable. Any discrepancies between TDDFT and the exciton model that remain upon accurately accounting for the effects of polarization would likely be due to the neglect of chargetransfer states in the exciton model. The qualitative trends observed for all numbers of chromophores in the aggregates are similar to those observed for dimers. Figure 6 shows the four similarity metrics for the spectra of chromophore octamers calculated using the exciton model relative to the spectra calculated using TD-ωB97X. The same metrics for tetramers, hexamers, and decamers can be found in Figures S7−S9, respectively. From a comparison of Figure 5A,B to Figure 6A,B it is evident that the spectra calculated using TD-ωB97X and the exciton model without atomic point charges have less similar shapes for octamers than they do for dimers. However, the inclusion of point charges significantly improves the agreement between the spectral shapes, as evidenced by increases in the Pearson correlation coefficients for both the spectra, r, and their derivatives, r′. As in the case of dimers, the exciton model systematically overestimates the integrated intensity of octamer spectra relative to TD-ωB97X (Figure 6C). The extent to which this occurs is not affected by the inclusion of atomic point charges. The exciton model also systematically over3796

DOI: 10.1021/acs.jctc.7b00477 J. Chem. Theory Comput. 2017, 13, 3787−3801

Article

Journal of Chemical Theory and Computation

Figure 7. Average similarity metrics for spectra of chromophore aggregates calculated using the exciton model with no point charges (solid lines) and with ChelpG charges (dotted lines) relative to spectra calculated using CIS (triangles, orange), TD-ωB97X (squares, black), and B3LYP (circles, magenta) as functions of the number of chromophores in the aggregate. (A) Pearson correlation coefficients for spectra. (B) Pearson correlation coefficients for the first derivatives of spectra. (C) Ratios of integrated intensities of the exciton-model and all-electron spectra. (D) Shifts of the spectral means of the exciton-model spectra relative to the all-electron spectra. The averaging is over 25 chromophore aggregates for each number of chromophores in the aggregate.

using this model and using all-electron methods is not restricted to aggregates consisting of a small number of chromophores. TeraChem’s algorithms for ab initio calculations that are optimized for massively parallel GPUs make it possible to benchmark the performance of the molecular exciton model against all-electron methods for relatively large chromophore aggregates. Using TeraChem on our available GPU machines, we performed CIS and TD-ωB97X calculations on aggregates consisting of up to 10 YLD124 chromophores (1250 atoms, 10 640 basis functions). As discussed in section 3.1, calculations using TD-B3LYP were performed only for aggregates with N ≤ 4 chromophores because of the large number of excited states that had to be included as a result of the presence of chargetransfer states with unphysically low energies. In this section, we consider how increasing the number of chromophores in the aggregate affects the agreement between the linear absorption spectra calculated using CIS, TD-ωB97X, or TDB3LYP and using the molecular exciton model. For each comparison, the molecular exciton model was parametrized with the equivalent all-electron method. The mean values and standard deviations for the four similarity metrics between spectra calculated using TD-wB97X and the exciton model (without partial atomic charges and with ChelpG charges included in the calculations) on aggregates of up to N = 10 chromophores are given in Table 2. These metrics for calculations using CIS and TD-B3LYP can be found in Tables S1 and S2, respectively. Figure 7 shows how the mean values of the similarity metrics for all three methods depend on the number of chromophores in the aggregate. The agreement between the exciton model and different all-electron methods varies slightly because of their varying treatments of exchange and correlation. The agreement in spectral shape seems to be slightly better for methods with more exact exchange (Figure

7A,B). Interestingly, for the integrated intensity (Figure 7C) and the mean excitation energy (Figure 7D), the agreement for TD-ωB97X is slightly worse than for either CIS or TD-B3LYP. However, the differences between methods for all of the similarity metrics are within a standard deviation and should not be considered significant. Figure 7A,B suggests that as the number of chromophores in the aggregate, N, increases, the agreement between the shapes of the spectra calculated using the exciton model with no atomic point charges and the spectra calculated using allelectron methods gets worse (except for N = 4). While r changes relatively little and remains close to 1, r′ for CIS decreases from 0.9547 for dimers to 0.8738 for decamers; the decrease is even larger for TD-ωB97X. However, the inclusion of atomic point charges in the calculation using the exciton model improves the agreement with the all-electron methods considerably. With point charges, for CIS and TD-ωB97X both r and r′ change little with the number of chromophores in the aggregate and remain greater than 0.99 and 0.95, respectively, for all numbers of chromophores in the aggregate. Thus, when point charges are included to account for interchromophore polarization, the agreement between the spectral shapes remains good as the number of chromophores in the aggregate is increased. The differences between the spectra predicted by the exciton model and by all-electron methods, as indicated by both the ratio of integrated absorption intensity, Iexc/Iall‐electron (Figure 7C), and the mean absolute energy, |ΔE| (Figure 7D), increase with the number of chromophores in the aggregate. The presence of atomic point charges makes no difference for Iexc/Iall‐electron and somewhat mitigates the discrepancy in |ΔE|. 3797

DOI: 10.1021/acs.jctc.7b00477 J. Chem. Theory Comput. 2017, 13, 3787−3801

Article

Journal of Chemical Theory and Computation

Figure 8. Average similarity metrics for spectra of chromophore aggregates with varying numbers of chromophores relative to spectra for aggregates with N = 20 chromophores. All of the spectra were calculated using the exciton model with no point charges parametrized with TD-ωB97X. (A) Pearson correlation coefficients for spectra. (B) Pearson correlation coefficients for the first derivatives of spectra. (C) Ratios of integrated intensities (per chromophore). (D) Shifts of the spectral means. The averaging was performed over 25 chromophore aggregates for each number of chromophores in the aggregate.

of the spectra for this figure were calculated using the exciton model with no point charges parametrized with TD-ωB97X. It is evident that the Pearson correlation coefficients for the spectra and their derivatives (Figure 8A,B) as well as the mean excitation energy (Figure 8D) change considerably as the number of chromophores in the aggregate increases, at least up to N = 20. Thus, the optical absorption spectra have not yet converged to that of the bulk material for these numbers of chromophores. For YLD124 chromophore aggregates with numbers of chromophores N ≤ 10, we find good agreement between optical absorption spectra calculated using all-electron methods and those calculated using the exciton model parametrized with point charges representing the chromophore environment. The agreement for larger aggregates is considerably worse when point charges are not included. The discrepancies between the spectra calculated using the exciton model and those calculated using all-electron methods become somewhat larger as the number of chromophores in the aggregate increases. Exploiting TeraChem’s algorithms optimized for massively parallel GPUs allowed us to calculate all-electron spectra for aggregates with up to N = 10 chromophores. We found that the increase in discrepancies between spectra calculated using the exciton model and those calculated using all-electron methods becomes slower for aggregates with N > 6. Therefore, we expect that the exciton model can predict the spectroscopic properties of large chromophore aggregates, thin films, and bulk molecular solids that are too large to study using all-electron methods.

For N > 6 chromophores, Iexc/Iall‐electron and |ΔE| change relatively slowly. Because the approximations made in the exciton model concern interchromophore interactions, they are most significant at small interchromophore distances. Therefore, the rate of change of all four spectral similarity metrics decreases as the number of chromophores in the aggregate increases. The variations in excitation energies among chromophores (∼0.01 eV) could lead to Anderson-localized excitonic states.77 However, these variations in the excitation energies are similar to the magnitudes of the interchromophore couplings, so excitonic states in the aggregates could be expected to be rather delocalized. Indeed, we find that the eigenstates of the molecular exciton Hamiltonian, eq 3, have contributions from excited states on all chromophores in the aggregate for N ∼ 10. For linear chains of helium atoms, convergence of errors in the excitation energy with the size of the system has been observed in calculations using an ab initio realization of the molecular exciton model for numbers of atoms N ∼ 10.78 However, one should not assume that the low rates of change of the similarity metrics shown in Figure 7 indicate that for YLD124 chromophore aggregates with N ∼ 10 the optoelectronic properties are approaching those of a bulk material. For a three-dimensional system with strong interchromophore interactions, such quick convergence is very unlikely. While ab initio calculations of the optical absorption spectra for aggregates with N > 10 YLD124 chromophores are difficult even using TeraChem, such calculations can be easily performed for up to several thousand chromophores using the exciton model. Figure 8 shows the average similarity metrics for spectra of 25 YLD124 aggregates with each number of chromophores up to N = 20 relative to the largest number of chromophores in the aggregate, N = 20. All

4. CONCLUSIONS In this paper we have examined the validity of the molecular exciton model for calculating linear absorption spectra of aggregates that consist of strongly polar YLD124 chromo3798

DOI: 10.1021/acs.jctc.7b00477 J. Chem. Theory Comput. 2017, 13, 3787−3801

Article

Journal of Chemical Theory and Computation phores. We have found that the excitation energies and excitonic couplings for chromophores in these aggregates are similar to those typical for photosynthetic complexes, where the exciton model is commonly used. We have also found that there is little ground-state charge transfer between chromophores in an aggregate. Therefore, for these strongly polar chromophores at high number density, the approximations of the exciton model are valid. To test the accuracy of the molecular exciton model for this system, we have systematically compared the linear absorption spectra of aggregates calculated using all-electron TD-B3LYP, TD-ωB97X, and CIS to those calculated using the exciton model parametrized with the corresponding methods. The spectral agreement between the exciton model and the allelectron methods does not significantly depend on the electronic structure method employed. Aggregate spectra calculated using TD-B3LYP contain many low-energy chargetransfer transitions, but these dark transitions do not significantly affect the agreement with the exciton-model spectra, despite the fact that the exciton model does not account for charge-transfer transitions. This suggests that little mixing between Frenkel exciton and charge-transfer states occurs, supporting the applicability of the exciton model to the systems studied. Accounting for the electrostatic environment of individual chromophores via electrostatic embedding of atomic point charges significantly improves the agreement between spectra calculated using the exciton model and those calculated using all-electron methods. We performed calculations using Mulliken and ChelpG definitions for atomic point charges and found that the specific definition had a negligible effect on the calculated spectrum. Thus, it is likely that chromophores in the aggregates are sufficiently well-separated to describe the electrostatic effects on their electron densities by the dipole moments of other chromophores. These dipole moments are similar for various definitions of point charges. The agreement between the integrated intensities and the mean excitation energies for spectra calculated using allelectron methods and using the exciton model is best for smaller chromophore aggregates. Using the massively parallel GPU-optimized algorithms implemented in the TeraChem program, we were able to perform all-electron calculations of spectra for aggregates with up to N = 10 chromophores. We found that for aggregates consisting of N > 6 chromophores, similarity metrics between spectra calculated using the exciton model and spectra calculated using all-electron methods change slowly with the number of chromophores in the aggregate. Thus, the exciton model should be applicable for calculations of spectral properties for large aggregates. As the number of chromophores increases, the spectra of very large aggregates are expected to converge to the spectra of bulk molecular solids that can be measured experimentally. The formalism of the molecular exciton model used in this paper allows straightforward extension to calculations of nonlinear optical response. In future work, we will assess the performance of the molecular exciton model in calculations of nonlinear optical properties of molecular aggregates.





Snapshot from a Monte Carlo simulation from which molecular aggregates were selected; analytical calculation of absorption spectra for chromophore dimers; parameters of the exciton-model Hamiltonian and magnitudes of single-chromophore transition dipole moments for dimers; similarity metrics for spectra of chromophore tetramers calculated using the exciton model with no point charges parametrized using TD-B3LYP, TDωB97X, or CIS calculations relative to aggregate spectra calculated using the corresponding all-electron method; similarity metrics for spectra of 25 chromophore tetramers calculated using the exciton model with no point charges and with ChelpG point charges relative to spectra calculated using TD-ωB97X; similarity metrics for spectra of 25 chromophore hexamers calculated using the exciton model with no point charges, with ChelpG point charges, and with Mulliken point charges calculated for pentamers relative to spectra calculated using TDωB97X; mean similarity metrics and their standard deviations for spectra of 25 chromophore dimers, tetramers, hexamers, octamers, and decamers calculated using the exciton model relative to spectra of the same chromophore aggregates calculated using CIS; mean similarity metrics and their standard deviations for spectra of 25 chromophore dimers and tetramers calculated using the exciton model relative to spectra of the same chromophore aggregates calculated using B3LYP (PDF)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Christine M. Isborn: 0000-0002-4905-9113 Funding

This work was supported by the U.S. Department of Defense (Proposal 67310-CH-REP) under the Air Force Office of Scientific Research Organic Materials Division. This work used the XStream and MERCED computational resources supported by the NSF MRI Program (Grants ACI-1429830 and ACI1429783). Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The authors thank Drs. Bruce Robinson and Andreas Tillack for providing the Monte Carlo simulation snapshot from which the structures of the chromophore aggregates were selected.



REFERENCES

(1) De Angelis, F. Modeling Materials and Processes in Hybrid/ Organic Photovoltaics: From Dye-Sensitized to Perovskite Solar Cells. Acc. Chem. Res. 2014, 47, 3349−3360. (2) Mirkovic, T.; Ostroumov, E.; Anna, J.; van Grondelle, R.; Govindjee; Scholes, G. Light Absorption and Energy Transfer in the Antenna Complexes of Photosynthetic Organisms. Chem. Rev. 2017, 117, 249−293. (3) Brehmer, L. Ultrathin Films for Sensorics and Molecular Electronics. In Polymer Sensors and Actuators; Osada, Y., De Rossi, D. E., Eds.; Springer: Berlin, 2010; Chapter 2, pp 15−90. (4) Soukoulis, C.; Economou, E. Electronic localization in disordered systems. Waves Random Media 1999, 9, 255−269.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b00477. 3799

DOI: 10.1021/acs.jctc.7b00477 J. Chem. Theory Comput. 2017, 13, 3787−3801

Article

Journal of Chemical Theory and Computation (5) Davydov, A. S. Theory of Molecular Excitons; Plenum Press: New York, 1971. (6) Agranovich, V. M. Excitations in Organic Solids; International Series of Monographs on Physics, Vol. 142; Oxford University Press: Oxford, U.K., 2008. (7) Tretiak, S.; Middleton, C.; Chernyak, V.; Mukamel, S. Exciton Hamiltonian for the Bacteriochlorophyll System in the Lh2 Antenna Complex of Purple Bacteria. J. Phys. Chem. B 2000, 104, 4519−4528. (8) Frähmcke, J.; Walla, P. Coulombic couplings between pigments in the major light-harvesting complex LHC II calculated by the transition density cube method. Chem. Phys. Lett. 2006, 430, 397−403. (9) Ishizaki, A.; Fleming, G. Theoretical examination of quantum coherence in a photosynthetic system at physiological temperature. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 17255−17260. (10) Rivera, E.; Montemayor, D.; Masia, M.; Coker, D. F. Influence of Site-Dependent Pigment-Protein Interactions on Excitation Energy Transfer in Photosynthetic Light Harvesting. J. Phys. Chem. B 2013, 117, 5510−5521. (11) Duan, H.-G.; Stevens, A.; Nalbach, P.; Thorwart, M.; Prokhorenko, V.; Miller, R. Two-Dimensional Electronic Spectroscopy of Light-Harvesting Complex II at Ambient Temperature: A Joint Experimental and Theoretical Study. J. Phys. Chem. B 2015, 119, 12017−12027. (12) Curutchet, C.; Mennucci, B. Quantum Chemical Studies of Light Harvesting. Chem. Rev. 2017, 117, 294−343. (13) Li, X.; Parrish, R. M.; Liu, F.; Kokkila Schumacher, S. I. L.; Martínez, T. J. An Ab Initio Exciton Model Including Charge-Transfer Excited States. J. Chem. Theory Comput. 2017, DOI: 10.1021/ acs.jctc.7b00171. (14) Lee, D.; Forsuelo, M. A.; Kocherzhenko, A. A.; Whaley, K. B. Higher-Energy Charge Transfer States Facilitate Charge Separation in Donor-Acceptor Molecular Dyads. J. Phys. Chem. C 2017, 121, 13043−13051. (15) Kocherzhenko, A. A.; Lee, D.; Forsuelo, M. A.; Whaley, K. B. Coherent and Incoherent Contributions to Charge Separation in Multichromophore Systems. J. Phys. Chem. C 2015, 119, 7590−7603. (16) Painelli, A.; Terenziani, F. Optical Spectra of Push-Pull Chromophores in Solution: A Simple Model. J. Phys. Chem. A 2000, 104, 11041−11048. (17) Beljonne, D.; Yamagata, H.; Brédas, J.; Spano, F.; Olivier, Y. Charge-Transfer Excitations Steer the Davydov Splitting and Mediate Singlet Exciton Fission in Pentacene. Phys. Rev. Lett. 2013, 110, 226402. (18) Lee, D.; Greenman, L.; Sarovar, M.; Whaley, K. Ab Initio Calculation of Molecular Aggregation Effects: A Coumarin-343 Case Study. J. Phys. Chem. A 2013, 117, 11072−11085. (19) Kurtz, S. New nonlinear optical materials. IEEE J. Quantum Electron. 1968, 4, 578−584. (20) Bass, M.; Bua, D.; Mozzi, R.; Monchamp, R. Optical secondharmonic generation in crystals of organic dyes. Appl. Phys. Lett. 1969, 15, 393−396. (21) Marder, S. Organic nonlinear optical materials: where we have been and where we are going. Chem. Commun. 2006, 131−134. (22) Hermann, J.; Ducuing, J. Third-order polarizabilities of longchain molecules. J. Appl. Phys. 1974, 45, 5100−5102. (23) Oudar, J. Optical nonlinearities of conjugated molecules Stilbene derivatives and highly polar aromatic compounds. J. Chem. Phys. 1977, 67, 446−457. (24) Oudar, J.; Chemla, D. Hyperpolarizabilities of the nitroanilines and their relations to the excited state dipole moment. J. Chem. Phys. 1977, 66, 2664−2668. (25) Bailey, R. T.; Cruickshank, F. R.; Halfpenny, P. J.; Pugh, D.; Sherwood, J. N.; Simpson, G. S. Molecular Crystals. In Principles and Applications of Nonlinear Optical Materials; Munn, R. W., Ironside, C. N., Eds.; CRC Press: Boca Raton, FL, 1993; Chapter 7, pp 143−188. (26) Benight, S. J.; Robinson, B. H.; Dalton, L. R. Nano-Engineering of Molecular Interactions in Organic Electro-Optic Materials. In Molecular Interactions; Meghea, A., Ed.; InTech: Rijeka, Croatia, 2012; Chapter 8, pp 183−208.

(27) Shi, Y.; Zhang, C.; Zhang, H.; Bechtel, J. H.; Dalton, L. R.; Robinson, B. H.; Steier, W. H. Low (Sub-1-V) Halfwave Voltage Polymeric Electro-optic Modulators Achieved by Controlling Chromophore Shape. Science 2000, 288, 119. (28) Jen, A.; Luo, J.; Kim, T.-D.; Chen, B.; Jang, S.-H.; Kang, J.-W.; Tucker, N. M.; Hau, S.; Tian, Y.; Ka, J.-W.; Haller, M.; Liao, Y.; Robinson, B.; Dalton, L.; Herman, W. Exceptional electro-optic properties through molecular design and controlled self-assembly. Proc. SPIE 2005, 5935, 593506. (29) Sekkat, Z.; Wood, J.; Aust, E. F.; Knoll, W.; Volksen, W.; Miller, R. D. Light-induced orientation in a high glass transition temperature polyimide with polar azo dyes in the side chain. J. Opt. Soc. Am. B 1996, 13, 1713−1724. (30) Dalton, L.; Sullivan, P.; Bale, D.; Olbricht, B. Theory-inspired nano-engineering of photonic and electronic materials: Noncentrosymmetric charge-transfer electro-optic materials. Solid-State Electron. 2007, 51, 1263−1277. (31) Olbricht, B.; Sullivan, P.; Wen, G.; Mistry, A.; Davies, J.; Ewy, T.; Eichinger, B.; Robinson, B.; Reid, P.; Dalton, L. Laser-assisted poling of binary chromophore materials. J. Phys. Chem. C 2008, 112, 7983−7988. (32) Dalton, L. R.; Harper, A. W.; Robinson, B. H. The role of London forces in defining noncentrosymmetric order of high dipole momenthigh hyperpolarizability chromophores in electrically poled polymeric thin films. Proc. Natl. Acad. Sci. U. S. A. 1997, 94, 4842− 4847. (33) Dreuw, A.; Weisman, J. L.; Head-Gordon, M. Long-range charge-transfer excited states in time-dependent density functional theory require non-local exchange. J. Chem. Phys. 2003, 119, 2943− 2946. (34) Dreuw, A.; Head-Gordon, M. Single-Reference ab Initio Methods for the Calculation of Excited States of Large Molecules. Chem. Rev. 2005, 105, 4009−4037. (35) Maitra, N. T. Perspective: Fundamental aspects of timedependent density functional theory. J. Chem. Phys. 2016, 144, 220901. (36) Tawada, Y.; Tsuneda, T.; Yanagisawa, S.; Yanai, T.; Hirao, K. A long-range-corrected time-dependent density functional theory. J. Chem. Phys. 2004, 120, 8425−8433. (37) Vydrov, O. A.; Scuseria, G. E. Assessment of a long-range corrected hybrid functional. J. Chem. Phys. 2006, 125, 234109. (38) Rohrdanz, M. A.; Herbert, J. M. Simultaneous benchmarking of ground- and excited-state properties with long-range-corrected density functional theory. J. Chem. Phys. 2008, 129, 034107. (39) Autschbach, J. Charge-Transfer Excitations and Time-Dependent Density Functional Theory: Problems and Some Proposed Solutions. ChemPhysChem 2009, 10, 1757−1760. (40) Stein, T.; Kronik, L.; Baer, R. Reliable Prediction of Charge Transfer Excitations in Molecular Complexes Using Time-Dependent Density Functional Theory. J. Am. Chem. Soc. 2009, 131, 2818−2820. (41) Subotnik, J. E. Communication: Configuration interaction singles has a large systematic bias against charge-transfer states. J. Chem. Phys. 2011, 135, 071104. (42) Ufimtsev, I.; Martínez, T. Quantum Chemistry on Graphical Processing Units. 3. Analytical Energy Gradients and First Principles Molecular Dynamics. J. Chem. Theory Comput. 2009, 5, 2619. (43) Titov, A.; Ufimtsev, I.; Luehr, N.; Martínez, T. Generating Efficient Quantum Chemistry Codes for Novel Architectures. J. Chem. Theory Comput. 2013, 9, 213. (44) Isborn, C.; Luehr, N.; Ufimtsev, I.; Martínez, T. Excited-State Electronic Structure with Configuration Interaction Singles and Tamm-Dancoff Time-Dependent Density Functional Theory on Graphical Processing Units. J. Chem. Theory Comput. 2011, 7, 1814−1823. (45) Sisto, A.; Glowacki, D. R.; Martinez, T. J. Ab Initio Nonadiabatic Dynamics of Multichromophore Complexes: A Scalable GraphicalProcessing-Unit-Accelerated Exciton Framework. Acc. Chem. Res. 2014, 47, 2857−2866. (46) Tillack, A.; Johnson, L.; Eichinger, B.; Robinson, B. Systematic Generation of Anisotropic Coarse-Grained Lennard-Jones Potentials 3800

DOI: 10.1021/acs.jctc.7b00477 J. Chem. Theory Comput. 2017, 13, 3787−3801

Article

Journal of Chemical Theory and Computation and Their Application to Ordered Soft Matter. J. Chem. Theory Comput. 2016, 12, 4362−4374. (47) Krugler, J.; Montgomery, C.; McConnell, H. Collective Electronic States in Molecular Crystals. J. Chem. Phys. 1964, 41, 2421. (48) Bakalis, L.; Knoester, J. Optical properties of one-dimensional exciton systems: Beyond the Heitler-London approximation. J. Chem. Phys. 1997, 106, 6964. (49) Kocherzhenko, A. A.; Dawlaty, J.; Abolins, B. P.; Herrera, F.; Abraham, D. B.; Whaley, K. B. Collective Effects in Linear Spectroscopy of Dipole-Coupled Molecular Arrays. Phys. Rev. A: At., Mol., Opt. Phys. 2014, 90, 062502. (50) Agranovich, V.; Basko, D. Frenkel excitons beyond the HeitlerLondon approximation. J. Chem. Phys. 2000, 112, 8156. (51) Craig, D. P.; Thirunamachandran, T. Molecular Quantum Electrodynamics: An Introduction to Radiation Molecule Interactions; Dover Publications: Mineola, NY, 1998; p 94. (52) Abel, N. H. Démonstration de l’impossibilité de la résolution algébrique des équations générales qui passent le quatrième degré. In Œuvres Complètes de Niels Henrik Abel, reissue ed.; Sylow, L., Lie, S., Eds.; Cambridge University Press: Cambridge, U.K., 2012; Vol. 1, Chapter 7, pp 66−94. (53) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in C: The Art of Scientific Computing, 2nd ed.; Cambridge University Press: Cambridge, U.K., 1992; Chapters 11.2 and 11.3, pp 469−480. (54) Becke, A. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648−5652. (55) Chai, J.-D.; Head-Gordon, M. Systematic optimization of longrange corrected hybrid density functionals. J. Chem. Phys. 2008, 128, 084106. (56) Hirata, S.; Head-Gordon, M.; Bartlett, R. Configuration interaction singles, time-dependent Hartree-Fock, and time-dependent density functional theory for the electronic excited states of extended systems. J. Chem. Phys. 1999, 111, 10774−10786. (57) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, revision E.01; Gaussian, Inc.: Wallingford, CT, 2009. (58) Hirata, S.; Head-Gordon, M. Time-dependent density functional theory within the Tamm-Dancoff approximation. Chem. Phys. Lett. 1999, 314, 291−299. (59) Hehre, W.; Ditchfield, R.; Pople, J. Self-Consistent Molecular Orbital Methods. XII. Further Extensions of Gaussian-Type Basis Sets for Use in Molecular-Orbital Studies of Organic Molecules. J. Chem. Phys. 1972, 56, 2257. (60) Hariharan, P.; Pople, J. Influence of polarization functions on molecular-orbital hydrogenation energies. Theor. Chem. Acc. 1973, 28, 213−222. (61) Krueger, B.; Scholes, G.; Fleming, G. Calculation of Couplings and Energy-Transfer Pathways between the Pigments of LH2 by the ab Initio Transition Density Cube Method. J. Phys. Chem. B 1998, 102, 5378−5386. (62) Czader, A.; Bittner, E. R. Calculations of the exciton coupling elements between the DNA bases using the transition density cube method. J. Chem. Phys. 2008, 128, 035101. (63) Harcourt, R.; Scholes, G.; Ghiggino, K. Rate expressions for excitation transfer. II. Electronic considerations of direct and through-

configuration exciton resonance interactions. J. Chem. Phys. 1994, 101, 10521. (64) Lu, T.; Chen, F. Multiwfn: A multifunctional wavefunction analyzer. J. Comput. Chem. 2012, 33, 580−592. (65) Isborn, C.; Götz, A.; Clark, M.; Walker, R.; Martínez, T. Electronic Absorption Spectra from MM and ab Initio QM/MM Molecular Dynamics: Environmental Effects on the Absorption Spectrum of Photoactive Yellow Protein. J. Chem. Theory Comput. 2012, 8, 5092−5106. (66) Mulliken, R. Electronic Population Analysis on LCAO-MO Molecular Wave Functions. I. J. Chem. Phys. 1955, 23, 1833. (67) Mulliken, R. Criteria for the Construction of Good SelfConsistent-Field Molecular Orbital Wave Functions, and the Significance of LCAO-MO Population Analysis. J. Chem. Phys. 1962, 36, 3428. (68) Breneman, C.; Wiberg, K. Determining atom-centered monopoles from molecular electrostatic potentials The need for high sampling density in formamide conformational analysis. J. Comput. Chem. 1990, 11, 361. (69) Randolph, K. A.; Myers, L. L. Basic Statistics in Multivariate Analysis; Oxford University Press, 2013; Chapter 2, pp 11−34. (70) Milanese, J. M.; Provorse, M. R.; Alameda, E., Jr.; Isborn, C. M. Convergence of Computed Aqueous Absorption Spectra with Explicit Quantum Mechanical Solvent. J. Chem. Theory Comput. 2017, 13, 2159. (71) Kulik, H. J.; Zhang, J.; Klinman, J. P.; Martínez, T. J. How Large Should the QM Region Be in QM/MM Calculations? The Case of Catechol O-Methyltransferase. J. Phys. Chem. B 2016, 120, 11381− 11394. (72) Karelina, M.; Kulik, H. J. Systematic Quantum Mechanical Region Determination in QM/MM Simulation. J. Chem. Theory Comput. 2017, 13, 563−576. (73) Zuehlsdorff, T. J.; Haynes, P. D.; Hanke, F.; Payne, M. C.; Hine, N. D. M. Solvent Effects on Electronic Excitations of an Organic Chromophore. J. Chem. Theory Comput. 2016, 12, 1853−1861. (74) Daday, C.; Curutchet, C.; Sinicropi, A.; Mennucci, B.; Filippi, C. Chromophore-Protein Coupling beyond Nonpolarizable Models: Understanding Absorption in Green Fluorescent Protein. J. Chem. Theory Comput. 2015, 11, 4825−4839. (75) Meister, J.; Schwarz, W. Principal Components of Ionicity. J. Phys. Chem. 1994, 98, 8245. (76) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, 105, 2999−3093. (77) Anderson, P. Absence of Diffusion in Certain Random Lattices. Phys. Rev. 1958, 109, 1492−1505. (78) Morrison, A. F.; You, Z.-Q.; Herbert, J. M. Ab Initio Implementation of the Frenkel-Davydov Exciton Model: A Naturally Parallelizable Approach to Computing Collective Excitations in Crystals and Aggregates. J. Chem. Theory Comput. 2014, 10, 5366− 5376.

3801

DOI: 10.1021/acs.jctc.7b00477 J. Chem. Theory Comput. 2017, 13, 3787−3801