From Model Hamiltonians to ab Initio Hamiltonians and Back Again

Article pubs.acs.org/JCTC

From Model Hamiltonians to ab Initio Hamiltonians and Back Again: Using Single Excitation Quantum Chemistry Methods To Find Multiexciton States in Singlet Fission Materials Nicholas J. Mayhall* Department of Chemistry, Virginia Tech, Blacksburg, Virginia 24060, United States S Supporting Information *

ABSTRACT: Due to the promise of significantly enhanced photovoltaic efficiencies, significant effort has been directed toward understanding and controlling the singlet fission mechanism. Although accurate quantum chemical calculations would provide a detail-rich view of the singlet fission mechanism, this is complicated by the multiexcitonic nature of one of the key intermediates, the 1(TT) state. Being described as two simultaneous and singlet-coupled triplet excitations on a pair of nearest neighbor monomers, the 1(TT) state is inherently a multielectronic excitation. This fact renders most singlereference ab initio quantum chemical methods incapable of providing accurate results. This paper serves two purposes: (1) to demonstrate that the multiexciton states in singlet fission materials can be described using a spin-only Hamiltonian and with each monomer treated as a biradical and (2) to propose a very simple procedure for extracting the values for this Hamiltonian from single-reference calculations. Numerical examples are included for a number of different systems, including dimers, trimers, tetramers, and a cluster comprised of seven chromophores.

I. INTRODUCTION Improving device efficiency remains one of the most important challenges for the development of more economically competitive organic photovoltaic solar cells. By creating multiple excited states out of a single photon, singlet fission stands as one possible avenue for significantly improving the efficiency of solar cells. Singlet fission is a photophysical mechanism by which a system undergoes an internal conversion from a single exciton electronic state (a singlet electron−hole pair) to a biexciton state (two coupled electron−hole pairs). Although originally discovered many years ago,1−4 singlet fission has recently enjoyed a dramatic growth of interest due to the realization that the mechanism can provide a path for exceeding the conventional theoretical efficiency limitations for solar cells.5,6 IA. Mechanism of 1(TT) Formation. The basic mechanism for singlet fission is typically written as4,7,8 S0 → S1 → 1(TT ) → T + T

Figure 1. Basic mechanism of singlet fission. (a) Ground state. Both monomers in respective ground states. (b) Bright singly excited state. Although depicted as a purely local state for simplicity, Frenkel exciton delocalizes over several monomers. (c) Biexciton state. Both monomers in lowest energy triplet state, simultaneous excitations on both monomers.

(1)

as direct,9,10 indirect,11 mediated,12−21 coherent,22,23 or nonadiabatic:10,24−30where V2e denotes a transfer mechanism driven

and has been given a simplistic depiction in Figure 1. Following an initial excitation to the bright S1 state, systems that exhibit singlet fission rapidly transition to a biexciton state, 1(TT), which can be thought of as two spin-coupled triplet states residing on neighboring chromophores. It is then usually hypothesized that from the 1(TT) state, two uncorrelated triplet excitons emerge and diffusively migrate to the interface.8 Most efforts to understand singlet fission have focused primarily on defining the mechanism by which the 1(TT) state is formed, the S1 → 1(TT) step. From these studies, multiple mechanisms have been put forth, which can be roughly described © XXXX American Chemical Society

hν

V2e

Direct: S0 → S1 → 1(TT )

Received: May 26, 2016

A

DOI: 10.1021/acs.jctc.6b00545 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation hν

V1e

one as it correctly describes the dominant character of the respective adiabatic state. IC. Multiexcitons in Quantum Chemistry. Because the intermediate 1(TT) state is dark, experiments can typically (with some exceptions)31 only make indirect observations. This would ordinarily be an ideal scenario for theoretical chemistry to provide the needed analysis. However, being a two-electron transition from the S0 ground state, most quantum chemistry methods are not able to be used since they can typically only describe single-electron transitions. For instance, the most popular computational approach for computing excited states is time-dependent density functional theory (TDDFT). However, TDDFT is not capable of computing multiple electronic excitations. This is not simply a matter of the desired calculations being too computationally demanding; linear response within the adiabatic approximation formally does not permit multiple excitations.32 Thus, the most commonly used computational method for excited states cannot be used for processes involving multiexcitons. Multiple electronic excitations are, however, included in equation-of-motion coupled cluster methods,33−38 such as EOM-CCSD. Unfortunately, not only is this method is too expensive for application to the large molecules which exhibit singlet fission, but the accuracy for doubly excited states is often very poor. For electronic states with dominant double excitation character, EOM-CCSD can have errors of 1−2 eV.39 This is much larger than the errors typically expected for single excitations (typically around 0.1−0.2 eV). The source of these errors are understood and come from the neglect of triple excitations which provide the extra degrees of freedom necessary to describe the simultaneous orbital relaxation and double excitation. This results in an unbalanced treatment of electron correlation between the doubly excited state and the ground state. Triple excitations can sometimes be included, but the incredible computational cost of such methods precludes their application to all but the smallest systems. Active-space methods, such as complete active-space selfconsistent field (CASSCF), are the only tools available that can treat single and double excitations on a similar footing.40−43 However, with each chromophore added to the model system, the cost of the calculation increases dramatically. The task of modeling excited states that are delocalized over multiple chromophores quickly becomes intractable. Further, because the active-space excitations only provide the minimal number of configurations to qualitatively describe the system, additional correlation effects must then be included to obtain quantitative accuracy. This is usually done either with perturbation theory or multireference CI.44−50 Although, for the systems studied here, expansive multireference CI treatments are not tractable, and perturbation theory can be wrought with difficulties.7,21 An alternative approach, which has recently seen success in studying singlet fission mechanisms,14,24−26,51−55 is the restricted active-space spin-flip method (RAS-SF).56−58 Unlike most other spin-flip methods that start from a high-spin triplet reference and use α → β excitations to generate the desired low-spin states,59−62 RAS-SF has the ability to flip multiple spins to treat systems with several unpaired electrons. Because the 1(TT) intermediate can be characterized as a singlet state comprised of four unpaired electrons, a double spin-flip RAS-SF calculation (RAS-2SF) has the ability to efficiently describe the qualitative features of this state. By starting from the lowest energy quintet reference state, the double α → β excitations in the RAS-2SF method conveniently generate the six completely open-shell

V1e

Indirect: S0 → S1 → CTx → 1(TT ) hν

V1e

Mediated: S0 → S1 + CTx → 1(TT ) hν

VN

Nonadiabatic: S0 → S1 → 1(TT ) hν

Coherent: S0 → S1 + 1(TT )

through a two-electron Hamiltonian matrix element, V1e a oneelectron Hamiltonian matrix element, and VN a nuclearelectronic nonadiabatic coupling element. The differences among the above mechanisms are mostly due to differences among methods of investigation (e.g., model Hamiltonians vs ab initio theoretical studies), with some differences being attributed to different model systems (e.g., dimer models vs larger cluster models, which has a large impact on the position and mixing of the charge transfer states).19 While significant progress has been made to work out the singlet fission mechanism, much less work has been dedicated to characterizing the structure of the multiexciton intermediate. The purpose of this paper is to def ine a simple computational strategy to aid the calculation and analysis of the 1(TT) multiexciton states in singlet f ission materials. IB. Structure of the 1(TT) State. So far, we have described the 1(TT) state as being a singlet-coupled pair of triplet states. This language implies that the Ms = 0 biexciton states can be written in terms of diabatic states coupled by Clebsch−Gordan coefficients: 1

(TT ) ≈ |S1 = 1, S2 = 1; S = 0, Ms = 0⟩ =

1 (| − 1, 1⟩|1, 1⟩ + |1, 1⟩| − 1, 1⟩ 3 − |0, 1⟩|0, 1⟩)

3

(TT ) ≈ |S1 = 1, S2 = 1; S = 1, Ms = 0⟩ =

1 (| − 1, 1⟩|1, 1⟩ − |1, 1⟩| − 1, 1⟩) 2

5

(TT ) ≈ |S1 = 1, S2 = 1; S = 2, Ms = 0⟩ =

1 (| − 1, 1⟩|1, 1⟩ + |1, 1⟩| − 1, 1⟩) 6

+

2 |0, 1⟩|0, 1⟩ 3

If this were a complete description of the states, then we would always find the energies of the biexciton states determined solely by direct exchange interactions following the order: 5(TT) ≤ 3 (TT) ≤ 1(TT), which we know is usually incorrect. Note that the underlying diabatic states are only perfectly degenerate in a nonrelativistic framework and void of external magnetic fields. Due to configuration mixing between 1(TT) and states of differing character, the singlet combination usually becomes the lowest energy biexciton state. However, despite the change in spin-state ordering, qualitatively the spin coupling model is typically a good approximation the final configurationn interaction (CI) eigenstates, and thus, the label 1(TT) is a useful B


Article

Journal of Chemical Theory and Computation configurations needed to properly describe the angular momentum coupling, in addition to the remaining 30 excitations which include partially and fully closed shell configurations. Further, unlike more expensive methods such as EOM-2SFCCSD,63 RAS-SF uses an orbital active space (typically comprised of the ROHF singly occupied orbitals) to reduce computational complexity, resulting in a number of configurations which only grows linearly with system size, assuming the number of spin-flips is kept constant. However, although RAS-SF is very efficient for systems with small active spaces, the calculations quickly become intractable for larger collections of chromophores. The reason is that each additional chromophore necessarily introduces more orbitals and electrons into the active space, making the cost of the calculation scale exponentially, and thus limiting that application of RAS-SF to systems with many chromophores. Furthermore, RAS-SF cannot provide quantitatively accurate results due to the neglect of dynamical correlation, and correction factors are often used to shift the computed energies.24,25,64 In light of the challenges in computing multiexciton states in molecular materials, this paper presents a simple strategy to obtain the multiexciton energies and wave functions. Herein, simple single excitation spin-flip calculations are projectively mapped onto a spin-Hamiltonian, which can then be solved by exact diagonalization, yielding the full multiexciton manifold. The overall methodology is outlined in the following section, with numerical examples in the following section.

1. Orbital optimization: Converge the reference high spin state (ROHF, DFT, CCSD, etc.), which has multiplicity 2n + 1 2. 1SF-CI: Calculate the 2n lowest energy eigenpairs of the ab initio Hamiltonian in the basis of Ms = S − 1 configurations.93 This is simply the conventional single spin-flip method performed on the (potentially much) higher spin reference state. Methods such as SF-TDDFT, SF-CIS, SF-CIS(D), RAS-1SF, EOM-SF-CCSD, etc. could all potentially be used here. 3. Eigenstate projection: Project the 1SF-CI eigenvectors onto the neutral determinant basis. This is done by first transforming the eigenvector into the localized-corresponding orbital basis, then deleting components coming from ionic contributions (those involving excitation from one site to another). The corresponding orbitals are singular vectors of the αβ overlap matrix:

AB

(3)

=Uis ΣsVa ̅ s

(4)

where Σ is the diagonal matrix of singular values, and Uis and Vas̅ are the left and right singular vectors, which act as projectors of the occupied α and virtual β spaces onto distinct closed-shell and open-shell subspaces. Because the orbital localization occurs only among the corresponding orbitals (these are identical to the singly occupied orbitals in an ROHF calculation), the notoriously difficult virtual orbital localization is not needed. In this paper, the Boys localization scheme is used.70−72 However, both Pipek−Mezey73 and Edmiston−Ruedenberg orbitals74 were found to give essentially identical results. Neither the corresponding orbital transformation nor the following orbital localization involve any mixed occupied and virtual orbital rotations, a fact necessary to maintaining a well-defined high-spin reference state. Using s and s′ to now indicate two possibly distinct, localized corresponding orbitals, the spin-flip amplitude vectors and then be transformed,

II. METHODS As described above, the description of the 1(TT) state as a pair of singlet-coupled triplet states, suggests that we should be able to find a localized basis in which a spin-Heisenberg−Dirac-Van Vleck Hamiltonian (HDvV) provides an accurate model of the system,65−67 HDvV Ĥ = −2 ∑ JAB SÂ SB̂

Pia ̅ = CμiSμνCνa ̅

(2)

where JAB is the effective exchange coupling constant between chromophores A and B, and SÂ (SB̂ ) is the local spin operator acting on site A (B). The sign of JAB determines the nature of the coupling; JAB < 0 indicates that the biexciton states are ordered from low-spin to high-spin, with the opposite true for JAB > 0. However, this single site model does not allow for the interaction of the states with configurations containing the local ground states. Furthermore, this would not make it possible to study the biexciton spectrum in systems with more than two chromophores because each additional chromophore would be required to increase the exciton number rank. To make these degrees of freedom available, we have chosen rather to adopt a two-site model, where each chromophore contains two radical sites. Thus, in a dimer system, there will be four sources of angular momentum, and the Ĥ HDvV Hamiltonian will have six distinct couplings. This is, of course, the natural model for this system, as it allows a one-to-one mapping between the lattice sites of the model to the localized singly occupied orbitals of the 5(TT) reference state. The process of mapping a high-spin reference single spin-flip wave function onto a spin-Hamiltonian to study ground state problems in magnetic systems has been described previously in detail.68,69 Below, a modified process is described for extracting an effective Hamiltonian which describes the multiexciton states in a system of n singlet-fission chromophores.

tss ′ =

∑ Uistia̅Va̅s ′

(5)

ia

To complete the projection of the spin-flip eigenstates, the above amplitudes are now multiplied by a Kronecker delta δs,s′ to project out the ionic components (those in which the spin-flip occurs between two separate orbitals). Thus, only the neutral tss excitations survive, and all ionic tss′ excitations are eliminated: tss ′ → ts

(6) 75,76

4. Bloch effective Hamiltonian: Following Bloch, and then Malrieu and co-workers,77 we use the 2n states with maximum projection to construct the following effective Hamiltonian: 2nStates

Hseff, t =

∑ m=0

:= − Jst

t sm̃ Emttm̃

(7) (8)

where tms̃ is the projected 1SF-CI eigenvectors, with the tilde indicating that the amplitude vectors have been C


Article

Journal of Chemical Theory and Computation

enough degrees of freedom to allow extension to higher number of chromophores since each additional chromophore would increase the multiexciton rank (number of simultaneously excited units). In contrast, the current approach treats an n chromophore cluster as a 2n-site lattice of S = 1/2 spins. This allows exciton quenching to occur, gaining access to biexcitons, triexcitons, etc. Effective exchange integrals have also been invoked to describe the biradical character present on a single chromophore, which was used as a design principle for maximizing singlet fission.80 b. Interaction with Charge-Transfer States. Although the projection onto the neutral determinant basis removes all explicit coupling to the charge-transfer (or charge-resonance) states, the energetic contributions of these couplings are still present, but are treated implicitly. For a concrete, yet simple, depiction of how the charge transfer states are folded into the effective Hamiltonian, we follow ref 77 and consider the effective Hamiltonian construction for a two-site model. The extension to four-site models (like the ones most relevant to the systems in this paper) is conceptually straightforward. In a localized basis, the ab initio two-orbital Hamiltonian can be directly mapped to a two-site Hubbard Hamiltonian, with the following matrix form:

symmetrically orthogonalized after the projection described above, and Em is the 1SF-CI energy for state m. Notice that in eq 7, variables s and t were used for the site labels, whereas in eq 2, variables A and B list the spin sites. This is highlights the essence of the current approach: Each localized corresponding orbital, s, is treated as a single site of a quantum spin−lattice. For the acene systems, the corresponding orbitals localize primarily along the long edges. An example is given in Figure 2 for a tetracene dimer cluster.

Figure 2. Four singly occupied orbitals of the quintet state for the tetracene dimer. Orbitals were computed at the unrestricted ωB97X-D/ 6-31G level of theory, and then the singly occupied orbitals were identified using a corresponding orbital transformation, followed by a Boys localization among the singly occupied orbitals. Because the localization only mixes the singly occupied orbitals, π − σ mixing is not a problem.

5. Diagonalize Ĥ HDvV: Obtain the energy spectrum of Ĥ HDvV by exact diagonalization. The size of the systems studied in this paper are limited due to the cost of the underlying SFCIS or SF-TDDFT calculation, with the largest cluster containing seven chromophores. This translates into diagonalizing a 14 × 14 spin−lattice, which only has 3432 configurations, and can easily be directly diagonalized without using subspace techniques. For this paper, a simple Octave code was written to build and diagonalize the Heisenberg Hamiltonian.78 For much larger clusters, different diagonalization techniques will be needed. Like full CI, the exact diagonalization of the Heisenberg Hamiltonian scales exponentially. However, because triexcitons and higher excitations are not expected to be important, truncated excitation methods should be effective for describing the biexcitonic manifold of states. a. Comparison with Other Approaches. With the plethora of computational approaches to studying singlet fission, it is important to clearly define the unique characteristics of this approach. Unlike most computational strategies which have employed model Hamiltonians,12,15−17,19,20,79 this approach uses a spin-only Hamiltonian and extracts the Hamiltonian from the results of ab initio spin-flip CI calculations through a one-to-one correspondence of an effective Hamiltonian and the matrix representation of the model Hamiltonian in the same space. This straightforward and well-defined approach makes it simple to apply routinely. Further, because the Bloch effective Hamiltonian theory is used, our model Hamiltonian is guaranteed to exactly reproduce the ab initio results in the Ms = S − 1 space. The approximation, therefore, lies in our assumption that this same Hamiltonian can be transferred to the Ms = S − n space. The use of a spin Hamiltonian to describe the biexciton states is not entirely new.8 Smith and Michl make arguments about the relative ordering of the 1(TT), 3(TT), and 5(TT) states by assuming an underlying Heisenberg Hamiltonian acting on a two-site lattice of S = 1 sites. An S = 1 model does not contain

Hubbard Ĥ

⎛ |12̅ ⟩ |2 1̅ ⟩ |1 1̅ ⟩ |22̅ ⟩⎞ ⎜ ⎟ ⎜ ⟨12̅ | 0 K t t ⎟ ⎜ ⎟ = ⎜⟨2 1̅ | K 0 t t ⎟ ⎜ ⎟ t U K ⎟ ⎜⟨1 1̅ | t ⎜ ⎟ t K U ⎠ ⎝ ⟨22̅ | t

(9)

In the above Hamiltonian, the values for direct exchange, K, electron transfer (hopping), t, and coulomb repulsion, U, can be taken directly from the ab initio integrals, i.e., K = ⟨12̅||21̅⟩. For weakly coupled systems, the lowest two eigenstates of the above Hamiltonian can be written as ⎛ 1−α 1⎞ ⎜ ⎟ − − α 1 1⎟ 1 ⎜ ⎜ ⎟ T= 2 0⎟ 2 ⎜ 1 − (1 − α) ⎜⎜ ⎟⎟ 2 0⎠ ⎝ 1 − (1 − α)

(10)

where columns one and two are the biradical singlet state triplet state, respectively, and α measures the amount of charge transfer present in the eigenstate. Note that the triplet state does not acquire any charge-transfer character because both chargetransfer configurations are of singlet spin symmetry and thus cannot mix with the triplet state. Projecting onto the neutral determinant basis and reorthogonalizing yields ⎛1 1 ⎞ ⎜ ⎟ 1 ⎜ 1 −1⎟ T̃ = 2 ⎜0 0 ⎟ ⎜ ⎟ ⎝0 0 ⎠

(11)

To form the effective Hamiltonian, we first project the Hamiltonian onto the ab initio eigenbasis (T) and then transform into the model space (T̃ ): eff ̃ †Ĥ HubbardTT̃ † Ĥ = TT

D

(12) DOI: 10.1021/acs.jctc.6b00545 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation ⎛ −4t(α − 1) −(α − 2)α − (α − 2)αU 2K − 4t(α − 1) −(α − 2)α − (α − 2)αU ⎜ ⎜ 1 −4t(α − 1) −(α − 2)α − (α − 2)αU = ⎜ 2K − 4t(α − 1) −(α − 2)α − (α − 2)αU 2⎜ 0 0 ⎜ ⎝ 0 0

From eq 12, it is clear that although the Hamiltonian exists only in the neutral determinant basis, the interaction with charge transfer configurations is implicitly accounted for by the introduction of the charge-hopping (t) and on-site repulsion integrals (U). We can also now write the Ĥ HDvV effective exchange coupling constant (eq 2) in terms of the Hubbard parameters, which clearly indicate the extent of charge transfer: JAB = K + 2t(1 − α) (2 − α)α +

1 (2 − α)αU 2

(13)

methods used in this paper include, SF-CAS, SF-CIS, and SF-TDDFT, and in principle, more sophisticated methods such as SF-CIS(D) and EOM-CCSD could also be used, cost permitting. 2. Phenomenological requirement: The physics of the multiexciton states must be describable by the model Hamiltonian, Ĥ HDvV. While it is quite obvious that this requirement must be met, it is not so obvious that this requirement will be met. The ground state of a cluster of closed shell chromophores certainly would not be well described using a spin-Hamiltonian. However, the biexciton states are most significantly comprised of open-shell components which, in a localized orbital basis, can be considered as interacting spins. In order to determine if the model Hamiltonian is accurately capturing the physics, we carry out calculations for which the full ab initio calculation can be compared against. All spin-flip calculations have been carried out with a development version of QChem.81 IIIA. Dimers. The assumption that the 1(TT) state will be correctly described by the same spin Hamiltonian which is extracted from the one spin-flip calculations is not necessarily an obvious one. To demonstrate the legitimacy of this claim, the biexciton binding energy

(14)

The eigenvalues of this effective Hamiltonian are identical to the corresponding eigenstates of the original Hamiltonian, E(Triplet) = −K

0 0⎞ ⎟ ⎟ 0 0⎟ 0 0⎟ ⎟ 0 0⎠

(15)

E(Singlet) = K + 4t (1 − α) (2 − α)α + (2 − α)αU (16)

The energy gap, E(Singlet) − E(Triplet), is therefore related to the extent of charge transfer present in the wave function. As α → 0, the gap tends to 2K, which is the result one would obtain if the bare Hamiltonian was diagonalized in the basis of neutral determinants. As α → 1, the gap tends to 2K + U, indicating complete charge transfer. c. Limitations. Although efficient, the Hamiltonian in eq 2 is far too simple to answer all questions about the multiexciton states. Most importantly, this approach cannot directly compute rates of singlet fission or 1(TT) formation or describe the relevant singlet and charge transfer states. The reason is that, by design, this model essentially integrates out the interaction with all bright states, such that this interaction is folded into the effective integrals, JAB. Thus, the spectrum of Ĥ HDvV does not contain the bright S1 state. However, as stated above, despite the significant recent effort to understand the mechanisms which generate 1 (TT), much less has been achieved in terms of understanding the nature and structure of the manifold of 1(TT) states. The approach in this paper is designed solely to address these questions, with the target quantities being the energy gaps between the spin states of the biexciton states Eb = E1(TT) − E5(TT) and the character and relative ordering of biexciton states in larger clusters.

Eb = E 1(TT ) − E 5(TT )

(17)

is computed for a series of chromophore dimers, which are shown in Figure 3. In recent studies,51 Eb was found to be the critical quantity in determining the rate of the 1(TT) → T + T process, which was defined as k 2[1(TT ) → T + T ] ≈ e αβEb

(18)

where α is a parameter (0.5) which comes from a linear free energy analysis, and β is 1/kBT.94 However, Eb had to be determined using 2SF calculations. We consider a broad range of chromophore coupling strengths, ranging from the very weakly coupled, nonbonded chromophore dimers, to covalently linked chromophores26,83 to the strong coupling limit observed in βcarotene undergoing intramolecular triplet pair formation. In Figure 4, the values of Eb are plotted for the conventional 2SF approach and the proposed 1SF+Ĥ HDvV approach. In order to make this comparison, a spin-flip method which is capable of performing both 1SF and 2SF calculations is needed. Although this prevents the use of SF-TDDFT for the evaluation of the spin model, other methods such as SF-CAS can be used. If our proposed methodology works perfectly, then the 1SF+Ĥ HDvV results should provide an identical Eb to the more expensive direct 2SF-CAS results. As shown in Figure 4, the proposed method provides incredibly accurate values of Eb, across a broad range of values reaching nearly 4 orders of magnitude. Because the proposed approach is most applicable in the weakcoupling limit, it is perhaps most surprising to see the accuracy obtained for β-carotene. Unlike the other systems studied, all

III. NUMERICAL RESULTS Having described all the necessary components of the method, we are now in a position to discuss performance. There are two independent requirements that must be satisfied for this approach to be of useful accuracy. 1. Ab initio requirement: The one spin-flip method used should be reliably accurate. In principle, one should perform a full CI calculation to obtain the Ms = n − 1 eigenstates (where n is again equal to the number of chromophores in the cluster). Clearly, this is generally not possible, and approximations must be made here. Spin-flip methods are perfectly suited to describing these particular states, and all of our results make use of this. The spin-flip E


Article


IIIB. Dimers: SF-TDDFT. Achieving a balanced description of both dynamical and static correlation is an outstanding problem in quantum chemistry, and several of the calculations of the 1 (TT) states in singlet fission have added in dynamical correlation corrections a posteriori. Because the current approach attempts to fold much of the physics up into the effective Hamiltonian, it is possible that dynamical correlation effects can also be part of this. As mentioned in the previous section, any one spin-flip method can be used to construct the effective Hamiltonian, this includes SF-TD-DFT. In Figure 5, we compare five different density functionals that are used to compute the biexciton binding energy for the same

Figure 3. Structures of the dimer systems used in evaluating the validity of the spin Hamiltonian. The 1,3 diphenylisobenzofuran (DPIBF) structures are taken from ref 82, with the three structures obtained by translating once along each lattice vector. β-Carotene is not actually a dimer, as the 1(TT) state exists on a monomer. The five structures, BETX-excimer, BET-X1, BET-X2, BET-B, BET-B-excimer are taken from ref 51.

Figure 5. Biexciton binding energies for dimers using various density functionals. Bottom axis: Eb in eV. Top axis: Reduction in norm after projection of all four one spin-flip eigenstates (0 indicates perfect projection, 4 indicates the eigenstate is orthogonal to the model space.) 6-31G* basis set used. All SF-TD-DFT calculations have been performed using the noncollinear kernel of Ziegler and co-workers (NC-SF-TDDFT).85−87

set of dimers shown in Figure 3. The functionals are ordered according the percentage of exact exchange included, with the range-separated hybrids coming first. The SF-CAS results are also included as a comparison. Of course, unlike with the 1SFCAS results in Figure 4, there exists no analogous two spin-flip calculation to compare against. As an admittedly weak measure of performance, the deviation in the norms of the projected eigenstates are also included. This corresponds to the value α in eq 10 The most obvious trend seen in Figure 5 is that upon increasing exact exchange the biexciton binding energy increases. In addition to the raw percentage of HF exchange, the manner in which it is introduced also matters. Both range-separated functionals (CAM-B3LYP and ωB97X-D) have relatively smaller binding energies. These results are rather easily explained by considering two aspects of the calculations. 1. Increasing the exact exchange preferentially stabilizes the high spin-states. Because we are computing the binding energy as the energy gap between the singlet and quintet biexciton states, increasing the exact exchange stabilizes the quintet state, thus decreasing the binding energy. 2. Range-separated functionals shift charge transfer excitations to higher energies. It is the coupling of the neutral configurations to charge transfer excitations which ultimately drives electron delocalization and antiferromagnetic coupling. By ramping up the percentage of exact

Figure 4. Biexciton binding energies, Eb, for the series of chromophore dimers, in Figure 3. Red bars: current approach, 1SF+Ĥ HDvV. Black bars: direct ab initio approach. Agreement between the two approaches provides evidence of the validity of the spin Hamiltonian, Ĥ HDvV. The “∗” next to the DPIBF-112 indicates a positive exciton binding energy but is depicted as a negative number on the log plot. The ab initio methods used are 1SF-CAS and 2SF-CAS, respectively, using the 631G* basis set.

four spin-sites on beta-carotene reside on a single molecule. This increases the molecular orbital overlap and, thus, the effective superexchange interaction compared to the other dimers studied. Using the direct ab initio calculation (2SF-CAS), the biexciton binding energy was 0.348 eV, compared to 0.343 eV using the proposed 1SF+Ĥ HDvV approach. On the basis of these results, it appears as though the Ĥ HDvV contains the necessary physics to describe the 1(TT) states. However, this conclusion is made with the assumption that different ab initio spin-flip methods will behave similarly. As is shown below, this is not always the case, particularly concerning the use of DFT methods using functionals that do not incorporate long-range corrections. For a comparison with methods incorporating external excitations such as RAS-SF and SF-CAS(h,p),56,84 see the Supporting Information. F


Article


Figure 6. Biexciton binding energies for trimers and tetramers of DPIBF. Binding energies relative to lowest quintet 5(TT) state: Ebn = 1(TT)n - min (5(TT)n) in meV. Red bars: Current approach, 1SF+Ĥ HDvV. Black bars: Direct ab initio approach. The ab initio methods used are 1SF-CAS, 3SF-CAS, and 4SF-CAS, with the 6-31G*basis set.

mixing of the charge-transfer and triplet biexciton state creating an erroneously delocalized biexciton state. This is ultimately the cause for the rather large binding energy computed with ωB97XD. We note that no intruder states where found for non-DFT methods based on ROHF orbitals. In lieu of a concrete assessment of the accuracy of different functionals (a study which would require the existence of 2SFTDDFT methods), it seems to be the case that range-separated functionals are better behaved, yield better projected norms, have fewer intruder states, and provide more consistent results. IIIC. DPIBF Trimers and Tetramers. Because the 1(TT) states exist on pairs of chromophores, the density of states for these biexcitons increases quadratically with the size of the cluster. For the currently proposed method to be useful for studying larger clusters, it is necessary that the distribution of the numerous 1(TT) states be accurately described. In Figure 6, the energies of the 1(TT) manifold of states, referenced to the lowest energy quintet 5(TT) state, are depicted for a series of trimer and tetramer clusters of DPIBF. The clusters are generated and named by extending a monomer along different combinations of lattice vectors. Due to the very small coupling of chromophores, the density of states that lie above the lowest quintet state increases with number of chromophores, with the largest binding energy of −2.5 meV occurring for the 221 cluster. For each trimer (or tetramer), 3SF-CAS (or 4SF-CAS) calculations are compared to 1SF-CAS + Ĥ HDvV results. In all cases, both the qualitative trends are reproduced well, in addition to a quantitative agreement with lowest energy 1(TT) states. The 122 cluster which has the largest deviation between the 1SF+Ĥ HDvV and 4SF-CAS results. Here, one can see that the second and third biexciton states are noticeably underestimated compared to the 4SF-CAS energies. Because the norms of the 1SF-CAS eigenstates for these clusters were all near 1 (>0.99), any deviations from the exact ab initio results should be interpreted as deficiencies in the model Hamiltonian in eq 2. Future work will explore the impact of introducing more terms into the model Hamiltonian, such as a biquadratic interaction.88−90 D. Tetracene 7-mer. As a final demonstration, a comparatively large cluster comprised of seven tetracene monomers is considered. Following the direct approach using spin-flip wave functions, one would need a seven spin flip

exchange for larger electron−electron distances, chargetransfer excitations are preferentially destabilized relative to the neutral configurations, thus diminishing the coupling responsible for antiferromagnetic coupling. This preferentially destabilizes the 1(TT) state compared to 5(TT). In addition to the general trends, there are a couple of notable failures that need mentioning. For both BET-B and BET-Xexcimer, both PBE0 and B3LYP predict that the 1(TT) state is unbound with respect to the 5(TT) state. This is neither intuitive nor consistent with the results using the other functionals and SF-CAS. By inspecting the projection norms of BET-X-excimer, it becomes apparent that the PBE0 and B3LYP eigenstates have significantly less overlap with the model space neutral determinants. This large decrease in the norm after projection is a good indication that the results are not as reliable. Unfortunately, however, the projection norms for BET-B are not obviously any worse for PBE0 and B3LYP than for the other functionals, and a simple inspection of the norms would not indicate a failure of the method. To better explain the failures of some of the DFT results, one must analyze the energies of the target states, in addition to the magnitude of the norms as described above. Typically, one would expect that the states of maximal projection onto the model space should fall toward the bottom of the 1SF spectrum. While this is usually the case, intruder states (states with near zero projection on the model space that fall below the energy of states with large model space projections) can arise from artificially low-energy charge transfer states. This is intimately related to electron selfinteraction problems in DFT. For both B3LYP and PBE0, every single calculation in Figure 5 suffered from the presence of intruder states. However, when a greater amount of exact exchange is added (BHHLYP), or when range-separated functionals are used (ωB97X-D, CAM-B3LYP), all of the intruder states disappear, with the exception of the tetracene dimer. The mere presence of intruder states is not necessarily ruinous. Often, the intruder states fall well within the relatively large window (≈1.5 eV) that separates the biexciton and single exciton 1SF states, with only modest coupling to the target states. However, for tetracene, the charge-transfer state was closer in energy to the biexciton target states, resulting in a non-negligible G


Article


cluster, whereas the fourth and fifth AA states are linear combinations of dimer states comprised of vertical combinations of A and A′. The final AA state is clearly the single A′A′ state, as depicted by the diagonal correlation lines in the graph. A similar analysis holds for each of the remaining states, which have been labeled accordingly. Several of these states involve the long-range spin-correlation described by Scholes.92 2. Relative Energies of 1(TT) States. Despite the small interchromophore coupling, a relatively broad energy window is observed for the 1(TT) states. On the basis of the above analysis of the correlation functions, the dimer states retain their monomer identities, and thus, the broadening of 1(TT) states comes from differences along the diagonal of the Hamiltonian and not due to significant coupling between the biexciton states. Thus, the state broadening seen here would be described as “static disorder” such that the individual triplet states are seen to have differing energies. This suggests that the 1(TT) states do not delocalize much beyond a dimer.

calculation to obtain the 1(TT) states in this cluster. This is a rather large calculation, involving over 11 million excitations. In comparison, the current approach requires only a one spin-flip calculation (196 excitations), followed by a diagonalization of the Ĥ HDvV Hamiltonian in a basis of 3432 spin configurations. For a cluster of seven chromophores, there are 21 different 1 (TT) states, and all of these are obtained in the diagonalization of Ĥ HDvV. In addition to the state energies, the eigenvectors of the Ĥ HDvV are in a convenient basis for analysis. To characterize the eigenstates, the spin−spin correlation functions are computed52,53,91 Cijm = 4(⟨SîzSjẑ ⟩m − ⟨Sîz⟩m ⟨Sjẑ ⟩m )

(19)

and plotted alongside the excitation energies in Figure 7, where Cmij is the spin−spin correlation between radical sites i and j for

Figure 7. Graphical depictions of the spin−spin correlation functions for each the 21 singlet biexciton states (1(TT)) in the tetracene 7-mer cluster. The eigenstates and energies are obtained from diagonalizing the Heisenberg Hamiltonian built from a 1SF-CAS/6-31G*calculation. Each graph vertex (dots) corresponds to one localized open-shell orbital on the 7-mer cluster in the top right of the figure. The size and weight of graph edge (lines) indicates the relative spin−spin correlation Cmij between two orbitals. Heavy blue (orange) lines indicate high orbital pair (anti-) correlation. Figure 8. Comparison of the 21 biexciton state energies (red lines) with the summed energies of pairs of triplet states (black lines).

state m, Ŝzi is the spin projection operator for site i, and the subscript m denotes an expectation value for state m. In Figure 7, the structure of the cluster is shown, with the chromophores labeled, A (and A′), B, and C. Analysis of the spin−spin correlation functions for each state provides a clear characterization of each of the 1(TT) states. Here, each singly occupied orbital i is represented by a vertex, and the weight of the lines connecting the vertices i and j denotes the magnitude of the spin−spin correlation function, Cmij for state m, with orange (blue) lines indicating negative (positive) correlation. 1. Characterization of 1(TT) States. Unsurprisingly, the ground state displays strong negative correlations between the electrons within a chromophore and negligible correlation between chromophores. The correlation functions for the 1 (TT) states are much more interesting, and each can be described as having the “intra-chromophore” correlation traded for “inter-chromophore” correlation. This is consistent with the description of the 1(TT) states being dimer states. For some of the states, this characterization is not immediately obvious due to degeneracies between similar states resulting in linear combinations. For example, the second and third AA states are linear combinations of dimer states existing on the top an bottom of the

To make this more concrete, in Figure 8, the energies of the (TT) biexciton states are compared with the energies obtained by adding together pairs of independent triplet states. Each of the single exciton eigenstates of Ĥ HDvV can be seen to exist primarily on a single chromophore, and thus, one can label the lowest seven triplet eigenstates accordingly: 4A, 2B, and C. For example, in the weak coupling limit, the 1(TT) states described as “AB” should have identical energies to the sum of two triplet states in which one is described as “A” and one is “B”. As shown in Figure 8, the 1(TT) spectrum is almost entirely determined by the energies of the underlying triplet state energies. Thus, the large 1 (TT) bandwidth observed in these calculations (SF-CAS being significantly larger than ωB97X-D) should be interpreted as a truncation effect coming from the use of a cluster model. This is an important result which urges caution when using larger cluster sizes in studying singlet fission. Periodic boundary conditions or density embedding approaches might be used to minimize trunctation affects. 1

H


Article

Journal of Chemical Theory and Computation Let us now attempt to justify the 1SF+Ĥ HDvV approach for the 7-mer cluster. In the previous calculations on dimers, trimers, and tetramers, the systems were small enough that the direct ab initio calculation (for SF-CAS) was available to compare against. For the 7-mer cluster, this is no longer the case. Thus, in lieu of a direct comparison of the 1(TT) states, we choose to compare a different set of states that are common to the spectrum obtained from both the direct and 1SF+Ĥ HDvV methods: the single exciton triplet states. This comparison should be adequately challenging due to the fact that in the 1SF+Ĥ HDvV approach the single exciton triplet states are referenced to the high-spin (15-et) multiplicity reference state, whereas in the straightforward TDDFT approach the triplet states are referenced to the closed shell singlet reference state. This comparison is made for the ωB97X-D results and shown in Figure 9. Here, triplet state energies from

calculation using the Bloch effective Hamiltonian theory. The proposed approach allows insight into the relative energies and spatial distributions of multiexciton states in large clusters which are typically unaccessible from direct ab initio calculations. An analysis of the numerical examples in this paper yields the following conclusions: 1. The electronic structure of biexciton states in singlet fission materials are easily mapped onto a 2n spin−lattice, where n is the number of chromophores in a cluster. 2. For dimer clusters, the exciton binding energy (Eb computed as the energy difference between the 1(TT) and 5(TT) states) is very accurately approximated with the proposed method. This was established by directly comparing multiple spin-flip calculations to the results obtained with the proposed method. 3. When the norm of the projection onto neutral spin configurations deviates significantly from one, legitimacy of the model must be questioned, as this indicates that charge-transfer (or ionic) configurations are increasingly important. 4. For SF-TDDFT, the percentage of exact exchange can significantly affect the results. Functionals with high amounts of exchange (BHHLYP) or range-separated functionals (CAM-B3LYP, ωB97X-D) seem to behave most reliably. 5. Most reliable results are expected for weakly coupled systems, which also happen to be the most efficient singlet fission systems. 6. For larger clusters, truncation effects can lead to site energy differences which artificially broaden the biexciton energy spectrum. On the basis of these results, we anticipate that the current approach will prove complementary to ab intio calculations and helpful to efforts in screening candidate structures for possible singlet fission technologies.

Figure 9. Comparison of the lowest energy triplet states for the ω-B97XD/6-31G* results. Ab initio TDDFT triplet states (black X labels). SFTDDFT/HDvV triplet states (blue solid line, + labels). Shifted SFTDDFT/HDvV triplet states (red dashed line, + labels). Energies are shifted (red dashed line, + labels) by a constant to remove the average deviation. States 1−4 are labeled A. State 5 is labeled C. States 6 and 7 are labeled B.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b00545. All Cartesian coordinates for all structures used in this paper and a figure detailing the comparison of the 2SFCAS(h,p) and 1SF-CAS(h,p)+Ĥ HDvV results. (PDF)

the direct ab initio TDDFT calculation (black X) and the 1SF +Ĥ HDvV calculation (blue +) are compared. As is clear from this graph, there is a significant deviation between the two approaches, with the 1SF+Ĥ HDvV energies lying about 0.06 eV higher than the direct TDDFT. However, this deviation is very consistent from state to state; if one shifts the 1SF+Ĥ HDvV energies by a constant value to remove the average deviation (red dashed line), the results are nearly identical to the direct TDDFT. This is valuable in that while the average position of the triplet states contains some error, the relative energies of the triplet states is very accurately reproduced. Since we have already established that the relative energies of the 1(TT) states are primarily determined by the relative energies of the underlying triplet states, it is reasonable to expect that the 1SF+Ĥ HDvV relative energies for the 1(TT) states are good approximations to what would be obtained from a direct ab initio calculation with the same level of theory

■

AUTHOR INFORMATION

Corresponding Author

*Electronic address: [email protected]. Notes

The author declares no competing financial interest.

■

ACKNOWLEDGMENTS Support for this work was provided by the Department of Chemistry at Virginia Tech. The author acknowledges Advanced Research Computing at Virginia Tech for providing computational resources and technical support that have contributed to the results reported within this paper. URL: http://www.arc.vt. edu

IV. CONCLUSIONS In this paper, a new method is described for calculating biexciton excited states in singlet fission materials. In short, the biexciton states are obtained as eigenstates of a Heisenberg spin Hamiltonian which is extracted from a relatively simple ab initio

■

REFERENCES

(1) Geacintov, N.; Pope, M.; Vogel, F. Phys. Rev. Lett. 1969, 22, 593− 596.

I


Article

Journal of Chemical Theory and Computation (2) Swenberg, C.; Stacy, W. Chem. Phys. Lett. 1968, 2, 327−328. (3) Singh, S.; Jones, W. J.; Siebrand, W.; Stoicheff, B. P.; Schneider, W. G. J. Chem. Phys. 1965, 42, 330. (4) Merrifield, R.; Avakian, P.; Groff, R. Chem. Phys. Lett. 1969, 3, 386− 388. (5) Shockley, W.; Queisser, H. J. J. Appl. Phys. 1961, 32, 510. (6) Hanna, M. C.; Nozik, A. J. J. Appl. Phys. 2006, 100, 074510. (7) Zimmerman, P. M.; Zhang, Z.; Musgrave, C. B. Nat. Chem. 2010, 2, 648−52. (8) Smith, M. B.; Michl, J. Chem. Rev. 2010, 110, 6891−936. (9) Burdett, J. J.; Bardeen, C. J. J. Am. Chem. Soc. 2012, 134, 8597−607. (10) Yost, S. R.; et al. Nat. Chem. 2014, 6, 492−7. (11) Johnson, J. C.; Akdag, A.; Zamadar, M.; Chen, X.; Schwerin, A. F.; Paci, I.; Smith, M. B.; Havlas, Z.; Miller, J. R.; Ratner, M. A.; Nozik, A. J.; Michl, J. J. Phys. Chem. B 2013, 117, 4680−95. (12) Yamagata, H.; Norton, J.; Hontz, E.; Olivier, Y.; Beljonne, D.; Brédas, J. L.; Silbey, R. J.; Spano, F. C. J. Chem. Phys. 2011, 134, 204703. (13) Teichen, P. E.; Eaves, J. D. J. Phys. Chem. B 2012, 116, 11473−81. (14) Casanova, D. J. Chem. Theory Comput. 2014, 10, 324−334. (15) Beljonne, D.; Yamagata, H.; Brédas, J. L.; Spano, F. C.; Olivier, Y. Phys. Rev. Lett. 2013, 110, 226402. (16) Berkelbach, T. C.; Hybertsen, M. S.; Reichman, D. R. J. Chem. Phys. 2013, 138, 114102. (17) Berkelbach, T. C.; Hybertsen, M. S.; Reichman, D. R. J. Chem. Phys. 2013, 138, 114103. (18) Chan, W.-L.; Berkelbach, T. C.; Provorse, M. R.; Monahan, N. R.; Tritsch, J. R.; Hybertsen, M. S.; Reichman, D. R.; Gao, J.; Zhu, X.-Y. Acc. Chem. Res. 2013, 46, 1321−9. (19) Berkelbach, T. C.; Hybertsen, M. S.; Reichman, D. R. J. Chem. Phys. 2014, 141, 074705. (20) Parker, S. M.; Seideman, T.; Ratner, M. A.; Shiozaki, T. J. Phys. Chem. C 2014, 118, 12700−12705. (21) Zeng, T.; Hoffmann, R.; Ananth, N. J. Am. Chem. Soc. 2014, 136, 5755−64. (22) Greyson, E. C.; Vura-Weis, J.; Michl, J.; Ratner, M. A. J. Phys. Chem. B 2010, 114, 14168−77. (23) Chan, W.-L.; Ligges, M.; Jailaubekov, A.; Kaake, L.; Miaja-Avila, L.; Zhu, X.-Y. Science 2011, 334, 1541−5. (24) Zimmerman, P. M.; Bell, F.; Casanova, D.; Head-Gordon, M. J. Am. Chem. Soc. 2011, 133, 19944−52. (25) Feng, X.; Luzanov, A. V.; Krylov, A. I. J. Phys. Chem. Lett. 2013, 4, 3845−3852. (26) Kolomeisky, A. B.; Feng, X.; Krylov, A. I. J. Phys. Chem. C 2014, 118, 5188−5195. (27) Tao, G. J. Phys. Chem. C 2014, 118, 17299−17305. (28) Akimov, A. V.; Prezhdo, O. V. J. Am. Chem. Soc. 2014, 136, 1599− 608. (29) Tao, G. J. Chem. Theory Comput. 2015, 11, 28−36. (30) Alguire, E. C.; Subotnik, J. E.; Damrauer, N. H. J. Phys. Chem. A 2015, 119, 299−311. (31) Tayebjee, M. J. Y.; et al. Phys. Chem. Chem. Phys. 2013, 15, 14797. (32) Reviews in Computational Chemistry; Lipkowitz, K. B., Cundari, T. R., Eds.; John Wiley & Sons, Inc.: Hoboken, NJ, 2008; Vol. 26. (33) Paldus, J.; Shavitt, I.; Cizek, J. Phys. Rev. A: At., Mol., Opt. Phys. 1972, 5, 50. (34) Cizek, J. J. Chem. Phys. 1966, 45, 4256. (35) Coester, F.; Kümmel, H. Nucl. Phys. 1960, 17, 477−485. (36) Saeh, J. C.; Stanton, J. F. J. Chem. Phys. 1999, 111, 8275. (37) Koch, H.; Jensen, H. J. A.; Jorgensen, P.; Helgaker, T. J. Chem. Phys. 1990, 93, 3345. (38) Koch, H.; Jorgensen, P. J. Chem. Phys. 1990, 93, 3333. (39) Schreiber, M.; Silva-Junior, M. R.; Sauer, S. P. A.; Thiel, W. J. Chem. Phys. 2008, 128, 134110. (40) Roos, B. O. Int. J. Quantum Chem. 1980, 18, 175−189. (41) Schmidt, M. W.; Gordon, M. S. Annu. Rev. Phys. Chem. 1998, 49, 233−66. (42) Li Manni, G.; Aquilante, F.; Gagliardi, L. J. Chem. Phys. 2011, 134, 034114.

(43) Li Manni, G.; Ma, D.; Aquilante, F.; Olsen, J.; Gagliardi, L. J. Chem. Theory Comput. 2013, 9, 3375−3384. (44) Hirao, K. Chem. Phys. Lett. 1993, 201, 59−66. (45) Hirao, K. Chem. Phys. Lett. 1992, 196, 397−403. (46) Hirao, K. Chem. Phys. Lett. 1992, 190, 374−380. (47) Roos, B. O.; Andersson, K.; Fülscher, M. P. Chem. Phys. Lett. 1992, 192, 5−13. (48) Angeli, C.; Cimiraglia, R.; Malrieu, J.-P. Chem. Phys. Lett. 2001, 350, 297−305. (49) Angeli, C.; Cimiraglia, R.; Evangelisti, S.; Leininger, T.; Malrieu, J.-P. J. Chem. Phys. 2001, 114, 10252. (50) Angeli, C.; Calzado, C. J. J. Chem. Phys. 2012, 137, 034104. (51) Korovina, N. V.; Das, S.; Nett, Z.; Feng, X.; Joy, J.; Haiges, R.; Krylov, A. I.; Bradforth, S. E.; Thompson, M. E. J. Am. Chem. Soc. 2016, 138, 617. (52) Luzanov, A. V.; Casanova, D.; Feng, X.; Krylov, A. I. J. Chem. Phys. 2015, 142, 224104. (53) Casanova, D.; Krylov, A. I. J. Chem. Phys. 2016, 144, 014102. (54) Chien, A. D.; Molina, A. R.; Abeyasinghe, N.; Varnavski, O. P.; Goodson, T.; Zimmerman, P. M. J. Phys. Chem. C 2015, 119, 28258− 28268. (55) Zimmerman, P. M.; Musgrave, C. B.; Head-Gordon, M. Acc. Chem. Res. 2013, 46, 1339−1347. (56) Casanova, D.; Head-Gordon, M. Phys. Chem. Chem. Phys. 2009, 11, 9779−90. (57) Zimmerman, P. M.; Bell, F.; Goldey, M.; Bell, A. T.; HeadGordon, M. J. Chem. Phys. 2012, 137, 164110. (58) Bell, F.; Zimmerman, P.; Casanova, D.; Goldey, M.; HeadGordon, M. Phys. Chem. Chem. Phys. 2013, 15, 358−366. (59) Shao, Y.; Head-Gordon, M.; Krylov, A. I. J. Chem. Phys. 2003, 118, 4807. (60) Krylov, A. Chem. Phys. Lett. 2001, 350, 522−530. (61) Krylov, A. I.; Sherrill, C. D. J. Chem. Phys. 2002, 116, 3194. (62) Krylov, A. I. Chem. Phys. Lett. 2001, 338, 375−384. (63) Casanova, D.; Slipchenko, L. V.; Krylov, A. I.; Head-Gordon, M. J. Chem. Phys. 2009, 130, 044103. (64) Casanova, D. J. Chem. Theory Comput. 2015, 11, 2642−50. (65) Dirac, P. A. M. Proc. R. Soc. London, Ser. A 1926, 112, 661. (66) Heisenberg, W. Eur. Phys. J. A 1928, 49, 619−636. (67) van Vleck, J. H. The Theory of Electric and Magnetic Susceptibilities; Clarendon Press: Oxford, 1932. (68) Mayhall, N. J.; Head-Gordon, M. J. Chem. Phys. 2014, 141, 134111. (69) Mayhall, N. J.; Head-Gordon, M. J. Phys. Chem. Lett. 2015, 6, 1982−1988. (70) Boys, S. Rev. Mod. Phys. 1960, 32, 296−299. (71) Subotnik, J. E.; Shao, Y.; Liang, W.; Head-Gordon, M. J. Chem. Phys. 2004, 121, 9220−9. (72) Subotnik, J. E.; Sodt, A.; Head-Gordon, M. Phys. Chem. Chem. Phys. 2007, 9, 5522−30. (73) Pipek, J.; Mezey, P. G. J. Chem. Phys. 1989, 90, 4916. (74) Edmiston, C.; Ruedenberg, K. Rev. Mod. Phys. 1963, 35, 457−464. (75) Bloch, C. Nucl. Phys. 1958, 6, 329−347. (76) des Cloizeaux, J. Nucl. Phys. 1960, 20, 321−346. (77) Malrieu, J. P.; Caballol, R.; Calzado, C. J.; de Graaf, C.; Guihéry, N. Chem. Rev. 2014, 114, 429−492. (78) Eaton, J. W.; Bateman, D.; Hauberg, S. GNU Octave, version 3.0.1 manual: A high-level interactive language for numerical computations; CreateSpace Independent Publishing Platform, 2009. (79) Aryanpour, K.; Shukla, A.; Mazumdar, S. J. Phys. Chem. C 2015, 119, 6966−6979. (80) Minami, T.; Ito, S.; Nakano, M. J. Phys. Chem. Lett. 2013, 4, 2133− 2137. (81) Shao, Y.; et al. Mol. Phys. 2015, 113, 184−215. (82) Schwerin, A. F.; et al. J. Phys. Chem. A 2010, 114, 1457−73. (83) Feng, X.; Krylov, A. I. Phys. Chem. Chem. Phys. 2016, 18, 7751−61. (84) Mayhall, N. J.; Goldey, M.; Head-Gordon, M. J. Chem. Theory Comput. 2014, 10, 589−599. (85) Wang, F.; Ziegler, T. J. Chem. Phys. 2005, 122, 074109. J


Article

Journal of Chemical Theory and Computation (86) Wang, F.; Ziegler, T. J. Chem. Phys. 2004, 121, 12191−6. (87) Bernard, Y. A.; Shao, Y.; Krylov, A. I. J. Chem. Phys. 2012, 136, 204103. (88) Huang, N.; Orbach, R. Phys. Rev. Lett. 1964, 12, 275−276. (89) Harris, E.; Owen, J. Phys. Rev. Lett. 1963, 11, 9−10. (90) Rodbell, D.; Jacobs, I.; Owen, J.; Harris, E. Phys. Rev. Lett. 1963, 11, 10−12. (91) Hachmann, J.; Dorando, J. J.; Avilés, M.; Chan, G. K.-L. J. Chem. Phys. 2007, 127, 134309. (92) Scholes, G. D. J. Phys. Chem. A 2015, 119, 12699−705. (93) Note that unlike the situation for the ground state magnetic systems69 we need twice as many eigenstates, 2n, to model the multiexciton manifold. (94) Note that the definition of Eb used here is the negative of that used in ref 51.

K


From Model Hamiltonians to ab Initio Hamiltonians and Back Again

Recommend Documents