Nonlinear response time-dependent density functional theory combined with the effective fragment potential method Federico Zahariev and Mark S. Gordon Citation: The Journal of Chemical Physics 140, 18A523 (2014); doi: 10.1063/1.4867271 View online: http://dx.doi.org/10.1063/1.4867271 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/18?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Time-dependent density functional theory excited state nonadiabatic dynamics combined with quantum mechanical/molecular mechanical approach: Photodynamics of indole in water J. Chem. Phys. 135, 054105 (2011); 10.1063/1.3622563 Implementation of the analytic energy gradient for the combined time-dependent density functional theory/effective fragment potential method: Application to excited-state molecular dynamics simulations J. Chem. Phys. 134, 054111 (2011); 10.1063/1.3523578 Solvent effects on optical properties of molecules: A combined time-dependent density functional theory/effective fragment potential approach J. Chem. Phys. 129, 144112 (2008); 10.1063/1.2992049 Density functional self-consistent quantum mechanics/molecular mechanics theory for linear and nonlinear molecular properties: Applications to solvated water and formaldehyde J. Chem. Phys. 126, 154112 (2007); 10.1063/1.2711182 Two-photon absorption in solution by means of time-dependent density-functional theory and the polarizable continuum model J. Chem. Phys. 122, 244104 (2005); 10.1063/1.1944727
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.120 On: Tue, 20 May 2014 20:29:47
THE JOURNAL OF CHEMICAL PHYSICS 140, 18A523 (2014)
Nonlinear response time-dependent density functional theory combined with the effective fragment potential method Federico Zahariev and Mark S. Gordona) Department of Chemistry, Iowa State University, Ames, Iowa 50011, USA
(Received 31 December 2013; accepted 18 February 2014; published online 17 March 2014) This work presents an extension of the linear response TDDFT/EFP method to the nonlinear-response regime together with the implementation of nonlinear-response TDDFT/EFP in the quantumchemistry computer package GAMESS. Included in the new method is the ability to calculate the two-photon absorption cross section and to incorporate solvent effects via the EFP method. The nonlinear-response TDDFT/EFP method is able to make correct qualitative predictions for both gas phase values and aqueous solvent shifts of several important nonlinear properties. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4867271] I. INTRODUCTION
While the polarization of a molecular system under the influence of a relatively weak oscillating electric field can be well described by a linear correction to the Hamiltonian within time-dependent perturbation theory, stronger fields typical in present-day laser experiments necessitate higherorder corrections characterizing a multitude of nonlinear optical effects. The frequency-dependent coefficients of the linear, quadratic, etc., corrections are, respectively, known as the dynamic polarizability, dynamic first hyperpolarizability, etc.1 The nonlinear optical effects of molecular systems are of practical interest in many areas of science, technology, and medicine. The second harmonic generation (SHG), electrooptical Pockels effect (EOPE), and optical rectification (OR) are second-order effects characterized by the dynamic first hyperpolarizability. The dynamic first hyperpolarizability has two frequencies as independent variables and the aforementioned second-order effects are each characterized by a specific combination of these frequencies. The version of the dynamic first hyperpolarizability corresponding to the SHG, i.e., the frequency-doubling effect, can also describe the process of a concerted two-photon absorption (TPA) by an atom or molecule. TPA was theoretically predicted by Goeppert-Mayer2 and first observed 30 years later by the generation of a blue fluorescent light upon laser excitation of CaF2 :Eu2+ crystals with a red light.3, 4 The individual photons participating in the TPA process are not necessarily resonant with any of the transitions in the atom or molecule; it is only their sum that is resonant. A simple conceptual model can be employed to explain TPA in terms of a short-lived virtual excited state intermediate between the ground and final excited states.5 A nonnegligible probability of such a process (measured by the TPA cross-section) occurs only when a high-intensity light source is used to assure a quick succession of photons. It is therefore not surprising that TPA was experimentally observed only after the advent of lasers. a) Electronic mail: [email protected]
0021-9606/2014/140(18)/18A523/10/$30.00
In addition to being of fundamental interest, TPA has numerous important applications, including 3D microscopy, 3D micro- and nano-fabrication, 3D optical storage, optical limiting, and pumped up-conversion lasing. An interesting application of TPA is in photodynamic cancer therapy,6 a technique that makes use of the high reactivity of singlet oxygen produced by one- or two-photon absorption. The advantage of using TPA in this context is that singlet oxygen formation by TPA occurs only when the incident photon influx is of high intensity, a constraint that allows for precise spatial and temporal control in destroying the cancer-infected tissue as well as an improved possibility of depth penetration in the body. There are several methods for the computation of linear and nonlinear molecular response properties. Among these, time-dependent density functional theory (TDDFT) presents a balance between accuracy and efficiency. High-level computational methods, such as coupled-cluster (CC),7–10 configuration-interaction (CI),11 and multi-reference (MR)12 methods provide high accuracy, but they generally have very steep computational scaling with the size of the molecular system. On the other hand, there are simplified approaches such as the few-state models for TPA13–15 using just a limited number of states in the sum-over-states (SOS) approach16, 17 that originates from the perturbative treatment of the effect18–22 and the finite-field approach.23–25 The finitefield approach is valid only for static (hyper-)polarizabilities. The SOS approach and its few-states truncation are sometimes necessary due to the excessive computational cost of the high-level electronic structure methods applied to relatively large molecules. Unfortunately, the simplified methods typically have very low accuracy. Modern density functionals provide good quality results for many properties and processes, with modest computational cost. However, there are several limitations caused by (1) the TDDFT adiabatic approximation, i.e., the neglect of explicit time-dependence (or frequency-dependence in the Fourier representation) in the TDDFT exchange-correlation term resulting from the use of a time-independent, ground state density functional approximation;26 (2) the well-known deficiencies of most present-day density functionals, such as
140, 18A523-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.120 On: Tue, 20 May 2014 20:29:47
18A523-2
F. Zahariev and M. S. Gordon
the self-interaction error,27, 28 the functional derivative discontinuity of the exchange-correlation potential with the particle number29, 30 and the incorrect long-range asymptotic behavior of the exchange-correlation potential.31 These limitations affect the accuracy of the computed excitations, especially those having charge-transfer, doubly excited or Rydberg character.32 Although nonlinear-response TDDFT has not yet reached the theoretical maturity and widespread practical use that linear-response TDDFT currently enjoys,33 there has recently been a growing interest in nonlinear response. An additional complication in comparing predicted nonlinear molecular response properties with experimental values is that these properties are typically measured in solution, not in the gas phase, while most calculations predict gas phase values. Two types of computational models are commonly used to account for solvent effects on electronic properties: implicit models and models that explicitly take into account the presence of solvent molecules. Explicit methods are usually employed by treating the solute with quantum mechanics (QM) and the solvent with model potentials (e.g., molecular mechanics (MM) force fields). The polarizable continuum model (PCM)34, 35 is a popular continuum solvent method. In this approach, the solvent effects are directly incorporated into the QM Hamiltonian. Although the TDDFT/PCM method36 has achieved some success in reproducing absorption and fluorescence spectra in solution, it has had difficulty describing the effect of microscopic interactions. Optical properties of solutes in a solvent environment are sensitive to the nature of the solvent-solvent and solute-solvent interactions, especially those that involve hydrogen-bonding interactions. The effective-fragment potential (EFP) method37, 38 is a sophisticated explicit solvent approach. The EFP method is based in quantum mechanics and provides a semi-classical description of the solvent-solvent and solute-solvent interactions including self-consistent induction, Coulombic interactions expressed as distributed multipoles, and exchange repulsion + charge transfer. In the initial implementation of EFP (EFP1), the latter term is fitted to either the Hartree-Fock (EFP1/HF) or DFT (EFP1/DFT) water dimer potential. The EFP1/DFT implementation is employed in the present work. EFP can be combined with DFT to incorporate solvent effects on QM ground-state properties39, 40 and with TDDFT to describe solvent effects on QM excited-state properties.41–44 A nonlinear-response TDDFT formalism, with the inclusion of EFP contributions for solvent effects, is presented here as recently implemented in the popular quantum-chemistry computer package GAMESS (General Atomic and Molecular Electronic Structure System).45, 46 The organization of this paper is as follows. A description of the density matrix formulation for nonlinear response TDDFT is followed by a brief overview of the EFP method and the manner in which the EFP method is interfaced with the nonlinear TDDFT method, in order to describe solvent effects on nonlinear properties. This is followed by a numerical analysis of the computational scheme that is developed and presented here. This is followed by a summary and conclusions.
J. Chem. Phys. 140, 18A523 (2014)
II. DENSITY-MATRIX BASED LINEARAND NONLINEAR-RESPONSE TDDFT
The present work employs the density matrix based approach to nonlinear-response TDDFT.47, 48 The method of Hirata and Head-Gordon47 is extended to the second-order response in the applied electric field. A density-matrix approach to second-order-response TDDFT was previously developed by using slightly different formalisms.48, 49 Although the linear-order response within this formalism recovers the well-known TDDFT matrix equations,47 it is presented here for consistency of the overall presentation. The system is initially in the ground state and its reduced one-electron density matrix ρ(r, r ) can be expanded in the Kohn-Sham spin-orbitals φp (r) as ρ(r, r ) =
Ppq φp (r)φq∗ (r ).
(1)
pq
The expansion coefficients Ppr (the discrete-index representation of the density matrix) satisfy the idempotence relation
Ppq Pqr = Ppr .
(2)
q
A simple derivation of this form of the idempotence relation for closed-shell systems is given in Appendix A. The electron density ρ(r) can easily be recovered from the density matrix ρ(r, r ): ρ(r) = ρ(r, r) =
Ppq φp (r)φq∗ (r).
(3)
pq
In addition to the idempotence relation, Ppr obeys the Heisenberg equation-of-motion (see p. 440 of Ref. 50): ∂Ppr . {Fpq Pqr − Ppq Fqr } = i ∂t q
(4)
For a non-hybrid density-functional, the Fock-like matrix F is Fpq =
−ZA 1 − ∇2 + |r − RA | 2 A ρ(r ) δE XC [ρ] φq (r)d r. + + d r |r − r| δρ(r) φp∗ (r)
(5)
The right-hand side of Eq. (5) contains the electronic kinetic, nuclear attraction, Coulomb repulsion, and exchangecorrelation energy components. For a hybrid densityfunctional, F is
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.120 On: Tue, 20 May 2014 20:29:47
18A523-3
F. Zahariev and M. S. Gordon
J. Chem. Phys. 140, 18A523 (2014)
⎞ 1 2 −ZA − ∇ + ⎟ ⎜ 2 |r − RA | ⎟ ⎜ A ⎟ φq (r)d r, = φp∗ (r) ⎜ N ⎜ ∗ XC φi (r )(1 − cH P (r, r ))φi (r ) δE [ρ] ⎟ ⎠ ⎝+ d r + (1 − cH ) |r − r| δρ(r) i=1 ⎛
Fpq
where P (r, r ) is the permutation operator that creates the exact Hatree-Fock (HF) exchange, and cH is a coefficient that mixes in HF exchange. The cH = 0 limit recovers the pure DFT expression of Eq. (5). The cH = 1 limit recovers the pure HF Fock matrix. Now, suppose there is an external perturbation, due to a time-dependent electric field, with matrix elements
1 gpq = λ hpq e−iωt + h∗qp eiωt , 2 where λ controls the strength of the perturbation and 3 Eξ μξ,pq hpq = φp∗ (r) (E · μ) φq (r)d r =
(7)
(8)
is a matrix element that couples the dipole-moment operator μˆ ξ = −ˆrξ (r1 = x, r2 = y, r3 = z), to the external electric field Eξ . For simplicity of notation, only one frequency ω is used but it is not difficult to generalize to the case of multiple frequencies. Both the density matrix and the Fock matrix can be expanded in powers of λ, thereby controlling the strength of the external perturbation (Eqs. (9) and (10)):
Fpq =
(0) Fpq
+
(1) λFpq
+λ
2
(2) Fpq
+λ
3
(3) Fpq
+ ....
(10)
q
is trivially satisfied due to (12)
(the indices i, j, . . . are used for occupied orbitals, a, b, . . . for unoccupied orbitals, and p, q, . . . for all orbitals), and Pij(0) = δij .
(13)
The right-hand side of Eq. (13) is the Kroneker delta. The ground state Kohn-Sham equation, Eq. (14), is a special case of the equation-of-motion, Eq. (4), such that (0) (0) (0) (0) Fpq (14) Pqr − Ppq Fqr = 0, q
where (0) Fpq = δpq εp .
q
(2) (0) (1) (1) (0) (2) (2) Ppq Pqr + Ppq Pqr + Ppq Pqr = Ppr ,
(17)
q
(3) (0) (2) (1) (1) (2) (0) (3) (3) Ppq Pqr + Ppq Pqr + Ppq Pqr + Ppq Pqr = Ppr , (18) and so on. The higher orders of the equations-of-motion are (0) (1) (1) (0) (0) (1) (1) (0) Fpq Fpq Pqr − Ppq Pqr − Ppq Fqr + Fqr q
(15)
q
=i
(1) ∂Ppr
∂t
(19)
,
(1) (1) (2) (0) (0) (2) (1) (1) Fpq Fpq Pqr − Ppq Pqr − Ppq Fqr + Fqr q
q
+
(9)
In the following, λ = 1 is assumed for simplicity. The expansions in Eqs. (9) and (10) result in a sequence of idempotence and equation-of-motion relations. In the zeroth order, the idempotence relation (0) (0) (0) Ppq Pqr = Ppr (11)
(0) Pia(0) = Pai(0) = Pab =0
In Eq. (15) εp are the Kohn-Sham orbital energies. The higher orders of the idempotence relation can also be derived from Eqs. (2) and (8): (1) (0) (0) (1) (1) Ppq Pqr − Ppq Pqr = Ppr , (16)
q
ξ =1
(0) (1) (2) (3) + λPpq + λ2 Ppq + λ3 Ppq + ..., Ppq = Ppq
(6)
(2) ∂Ppr (0) (2) (2) (0) , Fpq Pqr − Ppq Fqr = i ∂t q
(20)
(0) (3) (1) (2) (2) (1) (0) (3) Fpq Fpq Pqr − Ppq Pqr − Ppq Fqr + Fqr q
q
=i
(3) ∂Ppr
∂t
(21)
,
and so on. In order to relate the above density-matrix formalism to the molecular (hyper-) polarizabilities, consider the timedependent induced polarization pξ (t) [for simplicity, assume ∗ , i.e., the external perturbation is of the form gpq (t) fpq = fqp = hpq (ω)cos(ωt)], where the index ξ ∈ {x, y, z} specifies the Cartesian component of the vector: pξ (t) = μ(0) ξ + pξ (0) + pξ (ω) cos (ωt) + pξ (2ω) cos (2ωt) + . . . ,
(22)
where μ(0) ξ is the permanent ground-state dipole-moment, and the first three Fourier amplitudes of the time-dependent part are 3 1 βξ ζ η (0; −ω, ω) Eζ (−ω)Eη (ω) + . . . , (23) pξ (0) = 4 ζ,η=1
pξ (ω) =
3
αξ ζ (−ω; ω) Eζ (ω) + . . . ,
(24)
ζ =1
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.120 On: Tue, 20 May 2014 20:29:47
18A523-4
F. Zahariev and M. S. Gordon
pξ (2ω) =
J. Chem. Phys. 140, 18A523 (2014)
3 1 βξ ζ η (−2ω; ω, ω) Eζ (ω)Eη (ω) + . . . , 4 ζ,η=1
(25) α ξ ζ is the polarizability tensor and β ξ ζ η is the first hyperpolarizability tensor. The polarizability tensor α ξ ζ (−ω; ω) has two variables. The frequency ω to the right of the semicolon is the frequency of the incoming electric field Eζ (ω). The outgoing frequency −ω to the left of the semicolon is the frequency of the induced polarization pξ (ω). The sum of all of the variables is taken to be zero. Similarly, the hyperpolarizability tensor β ξ ζ η (−2ω; ω, ω) has two frequencies ω to the right of the semicolon that correspond to the two incoming fields Eζ (ω) and Eη (ω). There is one “outgoing” frequency −2ω to the left of the semicolon that corresponds to the induced polarization pξ (2ω). Once again the sum of all of the variables inside the parentheses is assumed to be zero. β ξ ζ η (−2ω; ω, ω) describes the SHG process and β ξ ζ η (0; −ω, ω) describes the OR process. Note that pξ (0),pξ (ω), pξ (2ω), etc. are the Fourier amplitudes of (1) (2) μξ,pq Ppq + Ppq + ... . (26) pξ (t) − μ(0) ξ =
further requires that dij = 0 and dab = 0, i.e., the Fourier component of the density matrix elements are non-zero only for the occupied-unoccupied and unoccupied-occupied pairs of indices. The first-order response of the Fock matrix is due to the original external perturbation, as well as to the dependence of the Fock matrix on the density (1) = gpq + Fpq
∂Fpq r,s
∂Prs
Prs(1) .
The response of the Fock matrix due to the density matrix response is δFpq =
∂Fpq r,s
+
∂Prs
Prs +
1 ∂ 2 Fpq Prs Pr s 2! r,s ∂Prs ∂Pr s r ,s
∂ 3 Fpq 1 Prs Pr s Pr s + . . . , 3! r,s ∂Prs ∂Pr s ∂Pr s r ,s r ,s
p,q
(33)
Based on Eqs. (23)–(25), the components of the (hyper-) polarizability tensors can be obtained by differentiation ∂pξ (ω) αξ ζ (−ω; ω) = , (27) ∂Eζ E=0 ∂ 2 pξ (2ω) βξ ζ η (−2ω; ω, ω) = 4 , (28) ∂Eζ ∂Eη E=0 ∂ 2 pξ (0) βξ ζ η (0; −ω, ω) = 4 . (29) ∂Eζ ∂Eη E=0 (1) Only Ppq from the perturbative expansion on the righthand side of Eq. (26) contributes to α ξ ζ (−ω; ω), as both are (2) contributes to β ξ ζ η (−2ω; ω, ω) and linear terms and only Ppq β ξ ζ η (0; −ω, ω) as these quantities are all quadratic terms, for example, (1) ∂ μξ,pq Ppq p,q αξ ζ (−ω; ω) = ∂Eζ E=0
=
(32)
μξ,pq
p,q
(1) ∂Ppq ∂Eζ
.
(30)
E=0
III. LINEAR-RESPONSE TDDFT
where Ppq is the full density matrix from the left-hand side of Eq. (9) and Fpq is the full Fock matrix from the left-hand side of Eq. (10). In order to obtain the Fock matrix response to a given order, the density matrix expansion from the right-hand side of Eq. (9) must be inserted into Eq. (33), and the result terminated at the desired order. Next, the explicit forms of the frequency-dependent density matrix and Fock matrix up to first order, given in Eqs. (31) and (32), respectively, are inserted into the firstorder equation-of-motion given in Eq. (19), and all of the terms in front of e−iωt on both sides of the resulting equation are collected. These two steps give ∂Fpq (0) (0) (0) dst Pqr Fpq dqr − dpq Fqr + hpq + ∂P st q st ∂Fqr (0) − Ppq dst (34) hqr + = ωdpr . ∂Pst st A similar procedure for eiωt will produce the complex conjugate of Eq. (34). Denoting Xai = dai and Yai = dia , the following linear-response matrix equation is obtained A B 1 0 X (ω) h −ω = . (35) B A 0 −1 Y (ω) h The solutions X and Y are vectors in terms of the index pairs “ai”, and the matrix elements of A and B are given by
The first-order (linear) response of the density-matrix is expressed as 1 ∗ iωt dpq e−iωt + dqp . (31) e 2 This form follows from standard time-dependent perturbation theory arguments. The first-order idempotence relation
Aia,j b = (εa − εi ) δij δab + (ia|j b) − cH (ab|j i) XC + (1 − cH )fia,j b
(1) = Ppq
(36)
and XC Bia,j b = (ai|bj ) − cH (aj |bi) + (1 − cH )fia,j b.
(37)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.120 On: Tue, 20 May 2014 20:29:47
18A523-5
F. Zahariev and M. S. Gordon
J. Chem. Phys. 140, 18A523 (2014)
The matrix elements of the exchange-correlation kernel fXC are 2 XC δ E [ρ] ∗ XC φ (r)φi (r)φb∗ (r )φj (r )d rd r . fai,bj = δρ(r)δρ(r ) a (38) The matrix elements of the two-electron integral are ∗ φp (r)φq (r)φs∗ (r )φr (r ) (pq|sr) = d rd r. (39) |r − r | The εp are KS orbital energies. The external perturbation 3 Eξ μξ,ai (see Eq. (8)) is also a vector in terms of the hai = ξ =1
index pairs “ai”. The linear response of the density matrix is Pai(1) = Xai + Yai .
(40)
The matrix elements of the solution to Eq. (35) are used on the right-hand side of Eq. (40). To obtain the polarizabil ∂P (1) ity tensor α ξ ζ (−ω; ω) using Eq. (30), one needs ∂Eaiξ E=0
(1) instead of Pai itself. In order to obtain a response equation (1) ∂Pai ∂ for ∂Eξ the derivative ∂Eξ is taken on both sides of E=0 E=0 Eq. (35). Referring to Eq. (8), one can see that after the differentiation is carried out, the dipole moment μξ becomes the perturbation on the right-hand side of Eq. (35). Since the matrices A and B do not depend on the electric field, the differentiation on the left-hand side of Eq. (35) acts only on the X and Y tensors that depend implicitly on the electric field. Hence, if a component of the dipole moment μξ is used as a perturbation in the linear-response equation, then A B 1 0 Xξ (ω) μξ =− . −ω Yξ (ω) B A 0 −1 μξ
(41) In Eq. (41), the dependence of the solution on the frequency ω and on the ξ -component of the dipole moment μ is explicitly noted. The polarizability tensor can be obtained from the solutions Xξ , ai (ω) and Yξ , ai (ω): αξ ζ (−ω; ω) = μξ,ai (Xζ,ai (ω) + Yζ,ai (ω)). (42)
The eigenvalue ωn corresponds to the nth excitation energy. The eigenvalues of Eq. (44) give the excitation spectrum and the eigenvectors form an orthonormal basis. The solutions of the inhomogenous matrix equation (44) can be expressed in terms of the eigenvectors and eigenvalues of the homogenous matrix equation by employing the spectral resolution method (see Appendix C). If the spectral resolution method is used, the averaged polarizability becomes 3
α¯ (−ω; ω) =
ξ =1
n
f0n , ωn − ω
(45)
IV. QUADRATIC-RESPONSE TDDFT
The second-order response of the density-matrix is 1 (2,−2) −i2ωt (2) (2,2) i2ωt (2,0) = + dpq e . (46) + dpq Ppq dpq e 2 The first superscript index of the frequency-independent (2,−2) (2,2) density matrix elements (Fourier coefficients) dpq , dpq (2,0) and dpq refers to the order of the response (2 means second order). The second superscript index refers to the frequency (2,−2) is the second order coin the exponent; for example, dpq i(−2)ωt (2,0) efficient multiplying e and dpq is the coefficient muli0ωt (2,−2) (2,2)∗ (2,0) (2,0)∗ tiplying e = 1. Also, dqp = dpq and dqp = dpq . (2,0) Since dpq contributes only to the zero-frequency component of the dipole moment and hence to β ξ ζ η (0; −ω, ω), it will not be considered here. The second-order idempotence relation determines the exact form of the occupied-occupied and unoccupied-unoccupied blocks of the second order density matrix directly in terms of the first-order density matrix components: (1) (1) Pia Paj (47) Pij(2) = − a
and (2) Pab =
Pai(1) Pib(1) .
(48)
i
The second-order response of the Fock-matrix is (2) = Fpq
∂Fpq r,s
∂Prs
Prs(2) +
where ∂ 2 Fpq = ∂Prs ∂Pr s
1 ∂ 2 Fpq P (1) P (1) , 2! r,s ∂Prs ∂Pr s rs r s
(49)
r ,s
E=0
The polarizability tensor α ξ ζ (−ω; ω) diverges at resonances. It is therefore assumed that when ω approaches an excitation frequency (resonance) ωn , the external perturbation fai can be adiabatically removed due to the resonant nature of the excitation. Then, the right-hand side of Eq. (41) is zero, and Eq. (41) becomes 1 0 A B Xn Xn = ωn . (44) Yn 0 −1 B A Yn
f0n is an oscillator strength for the 0 → n transition. Equation (45) reveals the divergent nature of the averaged polarizability near the excitation frequencies ωn . In practice, an eigenvalue matrix equation of half the dimension of Eq. (44) is solved (see Appendix B).
ai
Equation (42) is a direct consequence of Eqs. (30) and (40). Note that (1) ∂Ppq = Xξ,ai (ω) + Yξ,ai (ω) . (43) ∂Eξ
αξ ξ (−ω; ω) =
δ 3 EXC [ρ] φ ∗ (r)φq (r) δρ(r)δρ(r )δρ (r ) p
× φr∗ (r)φs (r)φr∗ (r)φs (r)d rd r d r .
(50)
Equation (50) follows from Eqs. (5) and (6). (2) Using the explicit form of Fpq from Eq. (49) and inserting it into Eq. (20) results in an explicit form of the secondorder equation-of-motion:
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.120 On: Tue, 20 May 2014 20:29:47
18A523-6
F. Zahariev and M. S. Gordon
J. Chem. Phys. 140, 18A523 (2014)
⎡
∂Fpq
∂Fqr
⎤
− + hpq + dst dqr − dpq hqr + dst ⎥ ∂Pst ∂Pst ⎥ ⎢ st st ⎥ ⎢ ⎛ ⎞ ⎛ ⎞ ⎥ ⎢ ⎢ ⎥ ⎢ ⎜ ⎟ ⎜ ⎟⎥ 2 2 ⎥ ⎢ ⎜ ∂Fpq (2) ⎟ ⎜ ⎟ ∂ F ∂F ∂ F 1 1 pq qr (2) qr q ⎢ (0) (0) ⎜ ⎥ Pqr dst + dst ds t ⎟ − Ppq dst + dst ds t ⎟ ⎥ ⎢+⎜ ⎜ ⎟ ⎜ ⎟ t 2! ∂Pst ∂Ps t ∂P 2! ∂P ∂P st st s ⎣ ⎝ s,t ∂Pst ⎠ ⎝ s,t ⎠⎦ (0) (2) ⎢ Fpq dqr
(2) (0) dpq Fqr
s,t s ,t
(51)
s,t s ,t
(2) = 2ωdpr .
(2) Recall that only the occupied-unoccupied dai and (2) unoccupied-occupied dia components of the solution are of interest. Using the notation
h(2) ai
=
(2) (2) Xai = dai ,
(52)
(2) Yai(2) = dia ,
(53)
hab +
b
∂Fab ck
∂Fab xck + yck ∂Pck ∂Pkc
(0) (0) and the diagonal form of Fpq and Ppq (see Eqs. (12), (13), and (15)). Equation (51) can be cast into a matrix equation form reminiscent of Eq. (35) (2) (2) A B h 1 0 X (ω) . =− − 2ω h(2) Y (2) (ω) B A 0 −1 (54) The components of h(2) are
xbi −
xaj hj i +
j
∂Fj i ck
∂ 2 Fai ∂ 2 Fai ∂ 2 Fai 1 xbj xb j + 2 xbj yb j + + 2! ∂Pbj ∂Pb j ∂Pbj ∂Pb j ∂Pbj ∂Pb j
∂Fj i xck + yck ∂Pck ∂Pkc ybj yb j ,
(55)
b,j b ,j
where x and y in Eq. (55) are the solutions to the nonhomogeneous linear-response Eq. (35). In Sec. III, the linear response solution was used to compute the polarizability based on Eq. (27) for the polarizability α ξ ζ (−ω; ω). In a similar manner, the second-order response solution gives rise to the SHG hyperpoliarizability β ξ ζ η (−2ω; ω, ω) based on Eq. (28) and to the OR hyperpolarizability β ξ ζ η (0; −ω, ω) based on Eq. (29). The most general form of the hyperpolarzability β ξ ζ η (ωξ ; ωζ , ωη ) can be obtained using an appropriate generalization of Eqs. (28) and (29). The dependence of the hyperpolarizability on the second-order response solution can be further transformed into a dependence on just the linear response solution by using the “(2n+1) theorem”:1 βξ ζ η (ωξ ; ωζ , ωη ) = A + B + C,
(56)
A=
perm.ξ,ζ,η
aij
Xξ,ai (ωξ ) −μζ,ij +
! +Yζ,ck (ωζ )) Yη,aj (ωη ) ,
C=
perm.ξ,ζ,η
iab
perm.ξ,ζ,η
H XC Xξ,ai (ωξ ) −μζ,ab + fab,ck (Xζ,ck (ωζ ) ck
(58)
XC gai,bj,ck (Xξ,ai (ωξ )+Yξ,ai (ωξ ))(Xζ,bj (ωζ )
aibj ck
+Yζ,bj (ωζ ))(Xη,ck (ωη ) + Yη,ck (ωη )) .
(59)
The “permutational sum” (perm. in Eqs. (57)–(59)) of a given quantity Qξ , ζ , η is
Qξ,ζ,η = Qξ,ζ,η + Qξ,η,ζ + Qζ,η,ξ + Qζ,ξ,η
perm.ξ,ζ,η
+Qη,ξ,ζ + Qη,ζ,ξ ,
H XC fij,ck (Xζ,ck (ωζ )
ck
! +Yζ,ck (ωζ )) Yη,bi (ωη ) ,
where any of the indices ξ , ζ , η can be x, y, or z and A, B, and C in Eq. (56) are defined as
B =−
(60)
and (57)
H XC fab,cd = (ia|j b) − cH (ab|j i) + (1 − cH )f XC ia,j b
(61)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.120 On: Tue, 20 May 2014 20:29:47
18A523-7
F. Zahariev and M. S. Gordon
J. Chem. Phys. 140, 18A523 (2014)
are the matrix elements of the combined Hartree-exchangecorrelation kernel and δ 3 E XC [ρ] XC gai,bj,ck = φ ∗ (r)φi (r)φb∗ (r ) δρ(r)δρ(r )δρ(r ) a × φj (r )φc∗ (r )φk (r )d rd r d r
lim (2ω − ωn ) βξ ζ η (−2ω; ω, ω) = σξ ζ (ωn )0|rη |n. (63)
2ω→ωn
Hence, the TPA cross-section itself is σξ ζ (ωn ) = A + B + C.
(64)
The three terms on the right-hand side of Eq. (64) are A=− Xξ,ai (ωξ ) − μζ,ij (1 − δζ n ) +
aij
H XC fij,ck (Xζ,ck (ωζ )
! + Yζ,ck (ωζ )) Yn,aj , (65)
ck
"
B=
perm.ξ,ζ,n
+
Xξ,ai (ωξ ) −ηζ,ab (1 − δζ n )
iab
H XC fabck (Xζ,ck (ωζ )
#
(66) C=−
perm.ξ,ζ,n
XC gai,bj,ck (Xξ,ai (ωξ ) + Yξ,ai (ωξ ))
aibj ck
×(Xζ,bj (ωζ ) + Yζ,bj (ωζ ))(Xn,ck
+ Yn,ck ) .
(67)
The value of the TPA cross section matrix element in Eq. (64) depends on the choice of coordinate system. However, orientationally averaged TPA transition moment quantities can be constructed from the matrix elements of the TPA. For a linear polarization of the light such an orientationally averaged TPA transition moment is 2Df + 4Dg ,
In the ground-state DFT/EFP1 method,39 the interaction between the DFT solute and the EFP water molecules39, 40 is represented by an effective one-body potential that has a Coulomb, polarization, and remainder component. The remainder term, which is fitted to a DFT water dimer potential, contains exchange repulsion, charge transfer, and short-range correlation contributions. The effective potential v EF P 1 (r) is added to the usual Kohn-Sham potential, but since only the EFP polarization term v P OL [ρ](r) depends on the electron density,39 only the EFP1 polarization term is updated during the SCF cycles. A key ingredient in the TDDFT linear response equation XC , is the exchange-correlation kernel f XC (r, r ) = δv δρ(r[ρ](r) ) which is a functional derivative of the exchange-correlation potential v XC [ρ] (r) with respect to the density. Since only the polarization component of the EFP1 effective potential v P OL [ρ] (r) has a non-zero functional derivative, the EFP1 correction to the linear-response TDDFT matrix equation41–44 is given by the following replacement: (73)
The last term in Eq. (73) is the EFP polarization kerP OL . The EFP polarization kernel nel f P OL (r, r ) = δv δρ(r[ρ](r) ) f P OL (r, r ), the derivative of the polarization potential with respect to the density, is expressed through the EFP polarizability tensor.41, 42 For applications of the combined TDDFT/EFP method, see Refs. 43 and 44. In order to extend the TDDFT/EFP1 scheme to nonlinear-response TDDFT, the replacement as given in Eq. (73) is performed in the expressions for the first hyperpolarizabilty and TPA cross-section. Since the EFP polarization potential v P OL [ρ](r) has only a linear dependence on the density, the EFP polarization kernel f P OL (r, r ) has no effective density dependence and hence g P OL (r, r , r ) P OL [ρ](r,r) = δf δρ(r = 0. Therefore, there is no need to correct ) g XC (r, r , r ) =
δf XC [ρ](r,r ) δρ(r )
in Eqs. (59) and (67).
(68) VI. COMPUTATIONAL DETAILS
where Df and Dg are defined as Df =
1 σξ ξ σζ ζ 30 ξ ζ
(69)
Dg =
1 σξ ζ σξ ζ , 30 ξ ζ
(70)
and
while the orientationally averaged TPA transition moment for a circular polarization is −2Df + 6Dg .
(72)
V. THE EFP METHOD
f XC (r, r ) → f XC (r, r ) + f P OL (r, r ).
+ Yζ,ck (ωζ )) Yn,bi ,
ck
(−Df + 3Dg ) . (Df + 2Dg )
(62)
are the matrix elements of the third functional derivative of the exchange-correlation energy density functional. The TPA absorption cross-section is proportional to the following residue:
perm.ξ,ζ,n
The polarization ration is
(71)
As an illustration, the difficult to model hyperpolarizability and two-photon absorption of bulk water will be computed by the new nonlinear-response TDDFT/EFP method as implemented in GAMESS. A system of 400 EFP molecules was equilibrated using a molecular dynamics (MD) simulation in the NVT ensemble. Upon equilibration, a smaller spherical sub-system of 150 EFP waters with a radius of about 9 A was extracted and one water molecule approximately in the middle of the sub-system was selected to be described by B3LYP/DH(d,p). A harmonic restraint potential with a force field of
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.120 On: Tue, 20 May 2014 20:29:47
18A523-8
F. Zahariev and M. S. Gordon
J. Chem. Phys. 140, 18A523 (2014)
The computed average TPA transition moment to the first excited state by two photons of equal energy is 9.9 a.u. with a solvent shift of 4.9 a.u. with respect to the gas phase value of 5.0 a.u. In Ref. 56, the same quantities were computed with a combined coupled cluster/molecular mechanics method to be, respectively, 7.2 a.u., 2.2 a.u., and 5.0 a.u. Reference 56 discusses the favorable qualitative comparison of the computed liquid phase TPA values with experiment57, 58 but there is no experimental basis to validate the computed TPA solvent shifts. VII. CONCLUSION
FIG. 1. The first excited state energy, second harmonic generation (SHG) hyper-polarizability (at ω = 0.043 a.u.) and the two-photon absorption transition moment and polarization ratio corresponding to the first excited state of water as a function of the simulation time (blue). The time averages are in green.
10 kcal/mol/A2 , effectively keeping the system in the NPT ensemble, was applied at the spherical boundary in order to prevent the evaporation of the water droplet.51, 52 The combined DFT/EFP system was equilibrated again for 55 ps and the next 45 ps production run was used to collect and average the optical properties of interest. All of the MD runs were done at 298 K using the velocity-Verlet algorithm with a time step of 0.5 fs. TDDFT (B3LYP)/EFP linear and nonlinear optical properties of water were computed on snapshots from the production run of the ground state DFT/EFP MD simulation. The d-aug-cc-pVDZ basis set was used for a proper description of the NLO properties. The snapshots from the MD trajectory were chosen with a separation of 250 fs to allow enough time for de-correlation between the snapshots to occur. Figure 1 shows both the running values (blue) and their time-averages (green) for several properties to illustrate the type and magnitude of statistical fluctuations during the MD simulation. All of the gas-phase optical values were computed with TDDFT using B3LYP/d-aug-cc-pVDZ. The average first excited state (11 B1 ) excitation energy is computed to be 7.72 eV with a solvent shift of 0.5 eV with respect to the gas phase first excited state energy of 7.22 eV. The corresponding experimental values are 8.2 eV for the liquid phase with a solvent shift of 0.8 eV relative to the 7.4 eV gas phase value.53 According to the experimental data, the sign of the SHG hyperpolarizability of liquid water changes relative to that of gas phase water. The average computed SHG hyperpolarizability of liquid water (11.8 a.u.) at the frequency ω = 0.043 a.u. compares reasonably well with the experimental value of 15.5 a.u. at the same frequency.54 The gas phase TDDFT/B3LYP/d-aug-cc-pVDZ SHG hyperpolarizability of water at ω = 0.0656 is predicted to be −5.9 a.u., again in reasonable agreement with the experimental value of −10.9 a.u. at the same frequency.55 Importantly, the calculations presented here predict the observed sign change.
A nonlinear response TDDFT computational scheme has been described that can successfully be combined with the EFP method to account for solvent effects. The combined method is able to correctly predict both gas phase values and aqueous solvent shifts for several important nonlinear properties. More details on the prediction of solvent effects (including non-aqueous solvents) within the presented TDDFT/EFP scheme will be described in subsequent publications. ACKNOWLEDGMENTS
This research is supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences, through the Ames Laboratory Chemical Physics program. The Ames Laboratory is operated for the U.S. Department of Energy by Iowa State University under Contract No. DE-AC02–07CH11358. APPENDIX A: THE IDEMPOTENCE RELATION OF THE DENSITY MATRIX
Expand the density-matrix in the basis of natural spinorbitals ρ(r, r ) =
N
φi (r)φi∗ (r ).
(A1)
i=1
Using this expansion it is easy to see that ρ r, r = ρ(r, r )ρ(r , r )d r
(A2)
holds. This is the idempotence relation of the density matrix in the coordinate representation. In order to get the form used in the text, insert the expansion of the density matrix in terms of the Kohn-Sham spinorbitals in the above equation. APPENDIX B: THE MODIFIED LINEAR-RESPONSE TDDFT EQUATION AND THE BIORTHOGONALITY CONDITION
The matrix equation can be simplified by a unitary transformation, which is equivalent to forming new block rows by adding and subtracting the original two block rows and thus forming two matrix equations of half-size, i.e., ( A + B) (X + Y ) = ω (X − Y )
(B1)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.120 On: Tue, 20 May 2014 20:29:47
18A523-9
F. Zahariev and M. S. Gordon
J. Chem. Phys. 140, 18A523 (2014)
and ( A − B) (X − Y ) = ω (X + Y ) .
(B2)
Multiplying both sides of the first matrix equation on the left by ( A − B) and using the second matrix equation results in a matrix eigenvalue equation ( A − B) ( A + B) (X + Y ) = ω2 (X + Y ) .
(B3)
Equation (B3) is half the size of the original matrix eigenvalue equation. This is still a non-Hermitian problem though. (X + Y )i and (X − Y )i are the right- and left-eigenvectors of the non-Hermitian matrix ( A − B) ( A + B)and the usual biorthonormality condition (X + Y )i · (X − Y )j = δij
(B4)
is imposed upon them, which condition could be also expressed as X i · X j − Y i · Y j = δij .
(B5)
The matrix equation can be modified further by multiplying both sides on the left by ( A − B)−1/2 and redefining the eigenvectors ( A − B)1/2 ( A + B) ( A − B)1/2 X˜ + Y˜ = ω2 X˜ + Y˜ , (B6) where X˜ + Y˜ = ( A − B)−1/2 (X + Y ). This is a Hermitian eigenvalue problem. The actual linear-response TDDFT implementation in GAMESS uses the Davidson algorithm based on both the Hermitian and non-Hermitian eiqenproblem formulations (see Ref. 59 for details). APPENDIX C: THE SPECTRAL RESOLUTION METHOD
The eigenvectors φi (r) of a Hermitian operator L, i.e., Lφi (r) = λi φi (r)
(C1)
form an orthogonal basis, so any function, say u(r), can be expanded in it, i.e., u(r) =
∞
αi φi (r).
(C2)
i=1
The action of L on u(r) becomes Lu(r) =
∞
αi λi φi (r).
(C3)
i=1
This representation is convenient for obtaining a function of L, i.e., f (L) u(r) =
∞
αi f (λi ) φi (r)
(C4)
i=1
and in the particular case of f (L) = (L − λ)−1 , it becomes −1
(L − λ)
u(r) =
∞ αi φi (r) i=1
λi − λ
.
(C5)
For additional information see Ref. 60. The solutions of the non-homogenous matrix equation in Eq. (41) can be expressed in terms of the eigenvectors and
eigenvalues of the homogenous matrix equation in Eq. (44) by employing the spectral resolution method: Yn Xn + (C6) [μξ · (Xn + Yn )] Xξ (ω) = ωn − ω ωn + ω n and Yξ (ω) =
n
[μξ · (Xn + Yn )]
Yn Xn + . (C7) ωn − ω ωn + ω
The scalar products on the right-hand side of Eqs. (46) and (47) are expressed in terms of the ai-matrix elements. The above expansions have a formal similarity with the sum-overstate (SOS) method, but the basis of the usual SOS expansion is different. 1 D.
Pugh, “Electric multipoles, polarizabilities and hyperpolarizabilities,” in Chemical Modeling: Applications and Theory, edited by A. Hinchliffe (RSC Publishing, 2000), Vol. 1. 2 M. Göppert-Mayer, Ann. Phys. 401, 273 (1931). 3 W. Kaiser and C. G. B. Garret, Phys. Rev. Lett. 7, 229 (1961). 4 I. D. Abella, Phys. Rev. Lett. 9, 453 (1962). 5 F. A. Moscatelli, Am. J. Phys. 54, 52 (1986). 6 J. D. Bhawalkar, G. S. He, and P. N. Prasad, Rep. Prog. Phys. 59, 1041 (1996); C. R. Shea, Y. Hefetz, R. Gillies, J. Wimberly, G. Dalickas, and T. Hasan, J. Biol. Chem. 265, 5977 (1990); P. Lenz, Photochem. Photobiol. 62, 333 (1995); J. D. Bhawalkar, N. D. Kumar, C. F. Zhao, and P. N. Prasad, J. Clin. Laser Med. Surg. 15, 201 (1997); W. G. Fisher, W. P. Partridge, C. Dees, and E. A. Wachter, Photochem. Photobiol. 66, 141 (1997). 7 D. Tunega and J. Noga, Theor. Chem. Acc. 100, 78 (1998). 8 A. B. Kumar, N. Vaval, and S. Pal, Chem. Phys. Lett. 295, 189 (1998). 9 O. Christiansen, A. Halkier, H. Koch, P. Jorgensen, and T. Helgaker, J. Chem. Phys. 108, 2801 (1998). 10 O. Chrisitiansen, J. Gauss, and J. F. Stanton, Chem. Phys. Lett. 292, 437 (1998). 11 P. B. Rozyczko and R. J. Bartlett, J. Chem. Phys. 109, 9201 (1998). 12 B. Fernandez, C. Hattig, H. Koch, and A. Rizzo, Chem. Phys. Lett. 288, 677 (1998). 13 P. Cronstrand, Y. Luo, and H. Agren, Chem. Phys. Lett. 352, 262 (2002). 14 P. Cronstrand, Y. Luo, and H. Agren, J. Chem. Phys. 117, 11102 (2002). 15 M. Drobizhev, Y. Stepanenko, Y. Dzenis, A. Karotki, A. Rebane, P. N. Taylor, and H. L. Anderson, J. Am. Chem. Soc. 126, 15352 (2004). 16 S. Delysse, P. Raimond, and J. M. Nunzi, Chem. Phys. 219, 341 (1997). 17 B. J. Orr and J. F. Ward, Mol. Phys. 20, 513 (1971). 18 C. Lambert, E. Schmalzlin, K. Meerholz, and C. Brauchle, Chem. – Eur. J. 4, 512 (1998). 19 U. S. Choi, T. W. Kim, S. W. Jung, and C. J. Kim, Bull. Korean Chem. Soc. 19, 299 (1998). 20 F. W. Vance, L. Karki, J. K. Reigle, J. T. Hupp, and M. A. Ratner, J. Phys. Chem. A 102, 8320 (1998). 21 Y. K. Lee, S. J. Jeon, and M. H. Cho, J. Am. Chem. Soc. 120, 10921 (1998). 22 P. Wang, C. G. Wang, P. W. Zhu, and C. Ye, Sci. China, Ser. B: Chem. 41, 156 (1998). 23 F. Grozema, R. Telesca, H. Jonkman, L. Siebbeles, and J. Snijders, J. Chem. Phys. 115, 10014 (2001). 24 F. Grozema, R. Telesca, J. Snijders, and L. Siebbeles, J. Chem. Phys. 118, 9441 (2003). 25 F. C. Grozema and P. T. van Duijnen, J. Phys. Chem. A 102, 7984 (1998). 26 A. Zangwill and P. Soven, Phys. Rev. A 21, 1561 (1980). 27 A. Dreuw, J. L. Weisman, and M. Head-Gordon, J. Chem. Phys. 119, 2943 (2003). 28 A. Dreuw and M. Head-Gordon, J. Am. Chem. Soc. 126, 4007 (2004). 29 D. J. Tozer, J. Chem. Phys. 119, 12697 (2003). 30 O. Gritsenko and E. J. Baerends, J. Chem. Phys. 121, 655 (2004). 31 M. E. Casida, C. Jamorsky, K. C. Casida, and D. R. Salahub, J. Chem. Phys. 113, 8918 (2000). 32 A. Dreuw and M. Head-Gordon, Chem. Rev. 105, 4009 (2005). 33 P. Elliott, F. Furche, and K. Burke, “Excited states from time-dependent density functional theory,” in Reviews in Computational Chemistry, edited by K. B. Lipkovitz and T. R. Cundari (Wiley, Hoboken, NJ, 2009), p. 91.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.120 On: Tue, 20 May 2014 20:29:47
18A523-10 34 J.
F. Zahariev and M. S. Gordon
Tomasi and M. Persico, Chem. Rev. 94, 2027 (1994). Tomasi, B. Mennucci, and R. Cammi, Chem. Rev. 105, 2999 (2005). 36 G. Scalmani, M. J. Frisch, B. Mennucci, J. Tomasi, R. Cammi, and V. Barone, J. Chem. Phys. 124, 094107 (2006). 37 M. S. Gordon, L. Slipchenko, L. Hui, and J. H. Jensen, Annu. Rep. Comput. Chem. 3, 177 (2007). 38 M. S. Gordon, M. A. Freitag, P. Bandyopadhyay, V. Kairys, J. H. Jensen, and W. J. Stevens, J. Phys. Chem. A 105, 293 (2001). 39 I. Adamovic, M. A. Freitag, and M. S. Gordon, J. Chem. Phys. 118, 6725 (2003); H. Li, H. M. Netzloff, and M. S. Gordon, ibid. 125, 194103 (2006); P. N. Day, J. H. Jensen, M. S. Gordon, S. P. Webb, W. J. Stevens, M. Krauss, D. Garmer, H. Basch, and D. Cohen, ibid. 105, 1968 (1996). 40 J. H. Jensen and M. S. Gordon, J. Chem. Phys. 108, 4772 (1998). 41 S. Yoo, F. Zahariev, S. Sok, and M. S. Gordon, J. Chem. Phys. 129, 144112 (2008). 42 N. Minezawa, N. De Silva, F. Zahariev, and M. S. Gordon, J. Chem. Phys. 134, 054111 (2011). 43 S. Sok, S. Y. Willow, F. Zahariev, and M. S. Gordon, J. Chem. Phys. A 115, 9801 (2011). 44 A. DeFusco, N. Minezawa, L. V. Slipchenko, and M. S. Gordon, J. Phys. Chem. Lett. 2, 2184 (2011). 45 M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. J. Su, T. L. Windus, M. Dupuis, J. A. Montgomery, J. Comput. Chem. 14, 1347 (1993). 35 J.
J. Chem. Phys. 140, 18A523 (2014) 46 M.
S. Gordon, M. W. Schmidt, in Theory and Applications of Computational Chemistry: The First Forty Years, edited by C. E. Dykstra, G. Frenking, K. S. Kim, and G. E. Scuseria (Elsevier, Amsterdam, 2005), Chap. 41, pp. 1167–1189. 47 S. Hirata and M. Head-Gordon, Chem. Phys. Lett. 314, 291 (1999). 48 F. Furche, J. Chem. Phys. 114, 5982 (2001). 49 F. Wang, C. Y. Yam, and G. Chen, J. Chem. Phys. 126, 244102 (2007). 50 R. McWeeny, Methods of Molecular Quantum Mechanics (Academic Press, 1992). 51 D. Beglov and B. Roux, J. Chem. Phys. 100, 9050 (1994). 52 C. H. Choi, S. Re, M. Feig, and Y. Sugita, Chem. Phys. Lett. 539–540, 218 (2012). 53 W. F. Chan, G. Cooper, and C. E. Brion, Chem. Phys. 178, 387 (1993). 54 B. F. Levine and C. G. Bethea, J. Chem. Phys. 65, 2429 (1976). 55 J. W. Ward and K. Miller, Phys. Rev. A 19, 826 (1979). 56 M. J. Paterson, J. Kongsted, O. Christiansen, K. V. Mikkelsen, and C. B. Nielsen, J. Chem. Phys. 125, 184501 (2006). 57 D. N. Nikogosyan, A. A. Oraevsky, and V. I. Rapusov, Chem. Phys. 77, 131 (1983). 58 C. L. Thomsen, D. Madsen, S. R. Keiding, J. Thogersen, and O. Christiansen, J. Chem. Phys. 110, 3453 (1999). 59 R. E. Stratmann and G. E. Scuseria, J. Chem. Phys. 109, 8218 (1998). 60 J. P. Keener, Principles of Applied Mathematics: Transformation and Approximation (Perseus Books, Cambridge, MA, 2000).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.120 On: Tue, 20 May 2014 20:29:47
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.120 On: Tue, 20 May 2014 20:29:47
THE JOURNAL OF CHEMICAL PHYSICS 140, 18A523 (2014)
Nonlinear response time-dependent density functional theory combined with the effective fragment potential method Federico Zahariev and Mark S. Gordona) Department of Chemistry, Iowa State University, Ames, Iowa 50011, USA
(Received 31 December 2013; accepted 18 February 2014; published online 17 March 2014) This work presents an extension of the linear response TDDFT/EFP method to the nonlinear-response regime together with the implementation of nonlinear-response TDDFT/EFP in the quantumchemistry computer package GAMESS. Included in the new method is the ability to calculate the two-photon absorption cross section and to incorporate solvent effects via the EFP method. The nonlinear-response TDDFT/EFP method is able to make correct qualitative predictions for both gas phase values and aqueous solvent shifts of several important nonlinear properties. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4867271] I. INTRODUCTION
While the polarization of a molecular system under the influence of a relatively weak oscillating electric field can be well described by a linear correction to the Hamiltonian within time-dependent perturbation theory, stronger fields typical in present-day laser experiments necessitate higherorder corrections characterizing a multitude of nonlinear optical effects. The frequency-dependent coefficients of the linear, quadratic, etc., corrections are, respectively, known as the dynamic polarizability, dynamic first hyperpolarizability, etc.1 The nonlinear optical effects of molecular systems are of practical interest in many areas of science, technology, and medicine. The second harmonic generation (SHG), electrooptical Pockels effect (EOPE), and optical rectification (OR) are second-order effects characterized by the dynamic first hyperpolarizability. The dynamic first hyperpolarizability has two frequencies as independent variables and the aforementioned second-order effects are each characterized by a specific combination of these frequencies. The version of the dynamic first hyperpolarizability corresponding to the SHG, i.e., the frequency-doubling effect, can also describe the process of a concerted two-photon absorption (TPA) by an atom or molecule. TPA was theoretically predicted by Goeppert-Mayer2 and first observed 30 years later by the generation of a blue fluorescent light upon laser excitation of CaF2 :Eu2+ crystals with a red light.3, 4 The individual photons participating in the TPA process are not necessarily resonant with any of the transitions in the atom or molecule; it is only their sum that is resonant. A simple conceptual model can be employed to explain TPA in terms of a short-lived virtual excited state intermediate between the ground and final excited states.5 A nonnegligible probability of such a process (measured by the TPA cross-section) occurs only when a high-intensity light source is used to assure a quick succession of photons. It is therefore not surprising that TPA was experimentally observed only after the advent of lasers. a) Electronic mail: [email protected]
0021-9606/2014/140(18)/18A523/10/$30.00
In addition to being of fundamental interest, TPA has numerous important applications, including 3D microscopy, 3D micro- and nano-fabrication, 3D optical storage, optical limiting, and pumped up-conversion lasing. An interesting application of TPA is in photodynamic cancer therapy,6 a technique that makes use of the high reactivity of singlet oxygen produced by one- or two-photon absorption. The advantage of using TPA in this context is that singlet oxygen formation by TPA occurs only when the incident photon influx is of high intensity, a constraint that allows for precise spatial and temporal control in destroying the cancer-infected tissue as well as an improved possibility of depth penetration in the body. There are several methods for the computation of linear and nonlinear molecular response properties. Among these, time-dependent density functional theory (TDDFT) presents a balance between accuracy and efficiency. High-level computational methods, such as coupled-cluster (CC),7–10 configuration-interaction (CI),11 and multi-reference (MR)12 methods provide high accuracy, but they generally have very steep computational scaling with the size of the molecular system. On the other hand, there are simplified approaches such as the few-state models for TPA13–15 using just a limited number of states in the sum-over-states (SOS) approach16, 17 that originates from the perturbative treatment of the effect18–22 and the finite-field approach.23–25 The finitefield approach is valid only for static (hyper-)polarizabilities. The SOS approach and its few-states truncation are sometimes necessary due to the excessive computational cost of the high-level electronic structure methods applied to relatively large molecules. Unfortunately, the simplified methods typically have very low accuracy. Modern density functionals provide good quality results for many properties and processes, with modest computational cost. However, there are several limitations caused by (1) the TDDFT adiabatic approximation, i.e., the neglect of explicit time-dependence (or frequency-dependence in the Fourier representation) in the TDDFT exchange-correlation term resulting from the use of a time-independent, ground state density functional approximation;26 (2) the well-known deficiencies of most present-day density functionals, such as
140, 18A523-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.120 On: Tue, 20 May 2014 20:29:47
18A523-2
F. Zahariev and M. S. Gordon
the self-interaction error,27, 28 the functional derivative discontinuity of the exchange-correlation potential with the particle number29, 30 and the incorrect long-range asymptotic behavior of the exchange-correlation potential.31 These limitations affect the accuracy of the computed excitations, especially those having charge-transfer, doubly excited or Rydberg character.32 Although nonlinear-response TDDFT has not yet reached the theoretical maturity and widespread practical use that linear-response TDDFT currently enjoys,33 there has recently been a growing interest in nonlinear response. An additional complication in comparing predicted nonlinear molecular response properties with experimental values is that these properties are typically measured in solution, not in the gas phase, while most calculations predict gas phase values. Two types of computational models are commonly used to account for solvent effects on electronic properties: implicit models and models that explicitly take into account the presence of solvent molecules. Explicit methods are usually employed by treating the solute with quantum mechanics (QM) and the solvent with model potentials (e.g., molecular mechanics (MM) force fields). The polarizable continuum model (PCM)34, 35 is a popular continuum solvent method. In this approach, the solvent effects are directly incorporated into the QM Hamiltonian. Although the TDDFT/PCM method36 has achieved some success in reproducing absorption and fluorescence spectra in solution, it has had difficulty describing the effect of microscopic interactions. Optical properties of solutes in a solvent environment are sensitive to the nature of the solvent-solvent and solute-solvent interactions, especially those that involve hydrogen-bonding interactions. The effective-fragment potential (EFP) method37, 38 is a sophisticated explicit solvent approach. The EFP method is based in quantum mechanics and provides a semi-classical description of the solvent-solvent and solute-solvent interactions including self-consistent induction, Coulombic interactions expressed as distributed multipoles, and exchange repulsion + charge transfer. In the initial implementation of EFP (EFP1), the latter term is fitted to either the Hartree-Fock (EFP1/HF) or DFT (EFP1/DFT) water dimer potential. The EFP1/DFT implementation is employed in the present work. EFP can be combined with DFT to incorporate solvent effects on QM ground-state properties39, 40 and with TDDFT to describe solvent effects on QM excited-state properties.41–44 A nonlinear-response TDDFT formalism, with the inclusion of EFP contributions for solvent effects, is presented here as recently implemented in the popular quantum-chemistry computer package GAMESS (General Atomic and Molecular Electronic Structure System).45, 46 The organization of this paper is as follows. A description of the density matrix formulation for nonlinear response TDDFT is followed by a brief overview of the EFP method and the manner in which the EFP method is interfaced with the nonlinear TDDFT method, in order to describe solvent effects on nonlinear properties. This is followed by a numerical analysis of the computational scheme that is developed and presented here. This is followed by a summary and conclusions.
J. Chem. Phys. 140, 18A523 (2014)
II. DENSITY-MATRIX BASED LINEARAND NONLINEAR-RESPONSE TDDFT
The present work employs the density matrix based approach to nonlinear-response TDDFT.47, 48 The method of Hirata and Head-Gordon47 is extended to the second-order response in the applied electric field. A density-matrix approach to second-order-response TDDFT was previously developed by using slightly different formalisms.48, 49 Although the linear-order response within this formalism recovers the well-known TDDFT matrix equations,47 it is presented here for consistency of the overall presentation. The system is initially in the ground state and its reduced one-electron density matrix ρ(r, r ) can be expanded in the Kohn-Sham spin-orbitals φp (r) as ρ(r, r ) =
Ppq φp (r)φq∗ (r ).
(1)
pq
The expansion coefficients Ppr (the discrete-index representation of the density matrix) satisfy the idempotence relation
Ppq Pqr = Ppr .
(2)
q
A simple derivation of this form of the idempotence relation for closed-shell systems is given in Appendix A. The electron density ρ(r) can easily be recovered from the density matrix ρ(r, r ): ρ(r) = ρ(r, r) =
Ppq φp (r)φq∗ (r).
(3)
pq
In addition to the idempotence relation, Ppr obeys the Heisenberg equation-of-motion (see p. 440 of Ref. 50): ∂Ppr . {Fpq Pqr − Ppq Fqr } = i ∂t q
(4)
For a non-hybrid density-functional, the Fock-like matrix F is Fpq =
−ZA 1 − ∇2 + |r − RA | 2 A ρ(r ) δE XC [ρ] φq (r)d r. + + d r |r − r| δρ(r) φp∗ (r)
(5)
The right-hand side of Eq. (5) contains the electronic kinetic, nuclear attraction, Coulomb repulsion, and exchangecorrelation energy components. For a hybrid densityfunctional, F is
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.120 On: Tue, 20 May 2014 20:29:47
18A523-3
F. Zahariev and M. S. Gordon
J. Chem. Phys. 140, 18A523 (2014)
⎞ 1 2 −ZA − ∇ + ⎟ ⎜ 2 |r − RA | ⎟ ⎜ A ⎟ φq (r)d r, = φp∗ (r) ⎜ N ⎜ ∗ XC φi (r )(1 − cH P (r, r ))φi (r ) δE [ρ] ⎟ ⎠ ⎝+ d r + (1 − cH ) |r − r| δρ(r) i=1 ⎛
Fpq
where P (r, r ) is the permutation operator that creates the exact Hatree-Fock (HF) exchange, and cH is a coefficient that mixes in HF exchange. The cH = 0 limit recovers the pure DFT expression of Eq. (5). The cH = 1 limit recovers the pure HF Fock matrix. Now, suppose there is an external perturbation, due to a time-dependent electric field, with matrix elements
1 gpq = λ hpq e−iωt + h∗qp eiωt , 2 where λ controls the strength of the perturbation and 3 Eξ μξ,pq hpq = φp∗ (r) (E · μ) φq (r)d r =
(7)
(8)
is a matrix element that couples the dipole-moment operator μˆ ξ = −ˆrξ (r1 = x, r2 = y, r3 = z), to the external electric field Eξ . For simplicity of notation, only one frequency ω is used but it is not difficult to generalize to the case of multiple frequencies. Both the density matrix and the Fock matrix can be expanded in powers of λ, thereby controlling the strength of the external perturbation (Eqs. (9) and (10)):
Fpq =
(0) Fpq
+
(1) λFpq
+λ
2
(2) Fpq
+λ
3
(3) Fpq
+ ....
(10)
q
is trivially satisfied due to (12)
(the indices i, j, . . . are used for occupied orbitals, a, b, . . . for unoccupied orbitals, and p, q, . . . for all orbitals), and Pij(0) = δij .
(13)
The right-hand side of Eq. (13) is the Kroneker delta. The ground state Kohn-Sham equation, Eq. (14), is a special case of the equation-of-motion, Eq. (4), such that (0) (0) (0) (0) Fpq (14) Pqr − Ppq Fqr = 0, q
where (0) Fpq = δpq εp .
q
(2) (0) (1) (1) (0) (2) (2) Ppq Pqr + Ppq Pqr + Ppq Pqr = Ppr ,
(17)
q
(3) (0) (2) (1) (1) (2) (0) (3) (3) Ppq Pqr + Ppq Pqr + Ppq Pqr + Ppq Pqr = Ppr , (18) and so on. The higher orders of the equations-of-motion are (0) (1) (1) (0) (0) (1) (1) (0) Fpq Fpq Pqr − Ppq Pqr − Ppq Fqr + Fqr q
(15)
q
=i
(1) ∂Ppr
∂t
(19)
,
(1) (1) (2) (0) (0) (2) (1) (1) Fpq Fpq Pqr − Ppq Pqr − Ppq Fqr + Fqr q
q
+
(9)
In the following, λ = 1 is assumed for simplicity. The expansions in Eqs. (9) and (10) result in a sequence of idempotence and equation-of-motion relations. In the zeroth order, the idempotence relation (0) (0) (0) Ppq Pqr = Ppr (11)
(0) Pia(0) = Pai(0) = Pab =0
In Eq. (15) εp are the Kohn-Sham orbital energies. The higher orders of the idempotence relation can also be derived from Eqs. (2) and (8): (1) (0) (0) (1) (1) Ppq Pqr − Ppq Pqr = Ppr , (16)
q
ξ =1
(0) (1) (2) (3) + λPpq + λ2 Ppq + λ3 Ppq + ..., Ppq = Ppq
(6)
(2) ∂Ppr (0) (2) (2) (0) , Fpq Pqr − Ppq Fqr = i ∂t q
(20)
(0) (3) (1) (2) (2) (1) (0) (3) Fpq Fpq Pqr − Ppq Pqr − Ppq Fqr + Fqr q
q
=i
(3) ∂Ppr
∂t
(21)
,
and so on. In order to relate the above density-matrix formalism to the molecular (hyper-) polarizabilities, consider the timedependent induced polarization pξ (t) [for simplicity, assume ∗ , i.e., the external perturbation is of the form gpq (t) fpq = fqp = hpq (ω)cos(ωt)], where the index ξ ∈ {x, y, z} specifies the Cartesian component of the vector: pξ (t) = μ(0) ξ + pξ (0) + pξ (ω) cos (ωt) + pξ (2ω) cos (2ωt) + . . . ,
(22)
where μ(0) ξ is the permanent ground-state dipole-moment, and the first three Fourier amplitudes of the time-dependent part are 3 1 βξ ζ η (0; −ω, ω) Eζ (−ω)Eη (ω) + . . . , (23) pξ (0) = 4 ζ,η=1
pξ (ω) =
3
αξ ζ (−ω; ω) Eζ (ω) + . . . ,
(24)
ζ =1
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.120 On: Tue, 20 May 2014 20:29:47
18A523-4
F. Zahariev and M. S. Gordon
pξ (2ω) =
J. Chem. Phys. 140, 18A523 (2014)
3 1 βξ ζ η (−2ω; ω, ω) Eζ (ω)Eη (ω) + . . . , 4 ζ,η=1
(25) α ξ ζ is the polarizability tensor and β ξ ζ η is the first hyperpolarizability tensor. The polarizability tensor α ξ ζ (−ω; ω) has two variables. The frequency ω to the right of the semicolon is the frequency of the incoming electric field Eζ (ω). The outgoing frequency −ω to the left of the semicolon is the frequency of the induced polarization pξ (ω). The sum of all of the variables is taken to be zero. Similarly, the hyperpolarizability tensor β ξ ζ η (−2ω; ω, ω) has two frequencies ω to the right of the semicolon that correspond to the two incoming fields Eζ (ω) and Eη (ω). There is one “outgoing” frequency −2ω to the left of the semicolon that corresponds to the induced polarization pξ (2ω). Once again the sum of all of the variables inside the parentheses is assumed to be zero. β ξ ζ η (−2ω; ω, ω) describes the SHG process and β ξ ζ η (0; −ω, ω) describes the OR process. Note that pξ (0),pξ (ω), pξ (2ω), etc. are the Fourier amplitudes of (1) (2) μξ,pq Ppq + Ppq + ... . (26) pξ (t) − μ(0) ξ =
further requires that dij = 0 and dab = 0, i.e., the Fourier component of the density matrix elements are non-zero only for the occupied-unoccupied and unoccupied-occupied pairs of indices. The first-order response of the Fock matrix is due to the original external perturbation, as well as to the dependence of the Fock matrix on the density (1) = gpq + Fpq
∂Fpq r,s
∂Prs
Prs(1) .
The response of the Fock matrix due to the density matrix response is δFpq =
∂Fpq r,s
+
∂Prs
Prs +
1 ∂ 2 Fpq Prs Pr s 2! r,s ∂Prs ∂Pr s r ,s
∂ 3 Fpq 1 Prs Pr s Pr s + . . . , 3! r,s ∂Prs ∂Pr s ∂Pr s r ,s r ,s
p,q
(33)
Based on Eqs. (23)–(25), the components of the (hyper-) polarizability tensors can be obtained by differentiation ∂pξ (ω) αξ ζ (−ω; ω) = , (27) ∂Eζ E=0 ∂ 2 pξ (2ω) βξ ζ η (−2ω; ω, ω) = 4 , (28) ∂Eζ ∂Eη E=0 ∂ 2 pξ (0) βξ ζ η (0; −ω, ω) = 4 . (29) ∂Eζ ∂Eη E=0 (1) Only Ppq from the perturbative expansion on the righthand side of Eq. (26) contributes to α ξ ζ (−ω; ω), as both are (2) contributes to β ξ ζ η (−2ω; ω, ω) and linear terms and only Ppq β ξ ζ η (0; −ω, ω) as these quantities are all quadratic terms, for example, (1) ∂ μξ,pq Ppq p,q αξ ζ (−ω; ω) = ∂Eζ E=0
=
(32)
μξ,pq
p,q
(1) ∂Ppq ∂Eζ
.
(30)
E=0
III. LINEAR-RESPONSE TDDFT
where Ppq is the full density matrix from the left-hand side of Eq. (9) and Fpq is the full Fock matrix from the left-hand side of Eq. (10). In order to obtain the Fock matrix response to a given order, the density matrix expansion from the right-hand side of Eq. (9) must be inserted into Eq. (33), and the result terminated at the desired order. Next, the explicit forms of the frequency-dependent density matrix and Fock matrix up to first order, given in Eqs. (31) and (32), respectively, are inserted into the firstorder equation-of-motion given in Eq. (19), and all of the terms in front of e−iωt on both sides of the resulting equation are collected. These two steps give ∂Fpq (0) (0) (0) dst Pqr Fpq dqr − dpq Fqr + hpq + ∂P st q st ∂Fqr (0) − Ppq dst (34) hqr + = ωdpr . ∂Pst st A similar procedure for eiωt will produce the complex conjugate of Eq. (34). Denoting Xai = dai and Yai = dia , the following linear-response matrix equation is obtained A B 1 0 X (ω) h −ω = . (35) B A 0 −1 Y (ω) h The solutions X and Y are vectors in terms of the index pairs “ai”, and the matrix elements of A and B are given by
The first-order (linear) response of the density-matrix is expressed as 1 ∗ iωt dpq e−iωt + dqp . (31) e 2 This form follows from standard time-dependent perturbation theory arguments. The first-order idempotence relation
Aia,j b = (εa − εi ) δij δab + (ia|j b) − cH (ab|j i) XC + (1 − cH )fia,j b
(1) = Ppq
(36)
and XC Bia,j b = (ai|bj ) − cH (aj |bi) + (1 − cH )fia,j b.
(37)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.120 On: Tue, 20 May 2014 20:29:47
18A523-5
F. Zahariev and M. S. Gordon
J. Chem. Phys. 140, 18A523 (2014)
The matrix elements of the exchange-correlation kernel fXC are 2 XC δ E [ρ] ∗ XC φ (r)φi (r)φb∗ (r )φj (r )d rd r . fai,bj = δρ(r)δρ(r ) a (38) The matrix elements of the two-electron integral are ∗ φp (r)φq (r)φs∗ (r )φr (r ) (pq|sr) = d rd r. (39) |r − r | The εp are KS orbital energies. The external perturbation 3 Eξ μξ,ai (see Eq. (8)) is also a vector in terms of the hai = ξ =1
index pairs “ai”. The linear response of the density matrix is Pai(1) = Xai + Yai .
(40)
The matrix elements of the solution to Eq. (35) are used on the right-hand side of Eq. (40). To obtain the polarizabil ∂P (1) ity tensor α ξ ζ (−ω; ω) using Eq. (30), one needs ∂Eaiξ E=0
(1) instead of Pai itself. In order to obtain a response equation (1) ∂Pai ∂ for ∂Eξ the derivative ∂Eξ is taken on both sides of E=0 E=0 Eq. (35). Referring to Eq. (8), one can see that after the differentiation is carried out, the dipole moment μξ becomes the perturbation on the right-hand side of Eq. (35). Since the matrices A and B do not depend on the electric field, the differentiation on the left-hand side of Eq. (35) acts only on the X and Y tensors that depend implicitly on the electric field. Hence, if a component of the dipole moment μξ is used as a perturbation in the linear-response equation, then A B 1 0 Xξ (ω) μξ =− . −ω Yξ (ω) B A 0 −1 μξ
(41) In Eq. (41), the dependence of the solution on the frequency ω and on the ξ -component of the dipole moment μ is explicitly noted. The polarizability tensor can be obtained from the solutions Xξ , ai (ω) and Yξ , ai (ω): αξ ζ (−ω; ω) = μξ,ai (Xζ,ai (ω) + Yζ,ai (ω)). (42)
The eigenvalue ωn corresponds to the nth excitation energy. The eigenvalues of Eq. (44) give the excitation spectrum and the eigenvectors form an orthonormal basis. The solutions of the inhomogenous matrix equation (44) can be expressed in terms of the eigenvectors and eigenvalues of the homogenous matrix equation by employing the spectral resolution method (see Appendix C). If the spectral resolution method is used, the averaged polarizability becomes 3
α¯ (−ω; ω) =
ξ =1
n
f0n , ωn − ω
(45)
IV. QUADRATIC-RESPONSE TDDFT
The second-order response of the density-matrix is 1 (2,−2) −i2ωt (2) (2,2) i2ωt (2,0) = + dpq e . (46) + dpq Ppq dpq e 2 The first superscript index of the frequency-independent (2,−2) (2,2) density matrix elements (Fourier coefficients) dpq , dpq (2,0) and dpq refers to the order of the response (2 means second order). The second superscript index refers to the frequency (2,−2) is the second order coin the exponent; for example, dpq i(−2)ωt (2,0) efficient multiplying e and dpq is the coefficient muli0ωt (2,−2) (2,2)∗ (2,0) (2,0)∗ tiplying e = 1. Also, dqp = dpq and dqp = dpq . (2,0) Since dpq contributes only to the zero-frequency component of the dipole moment and hence to β ξ ζ η (0; −ω, ω), it will not be considered here. The second-order idempotence relation determines the exact form of the occupied-occupied and unoccupied-unoccupied blocks of the second order density matrix directly in terms of the first-order density matrix components: (1) (1) Pia Paj (47) Pij(2) = − a
and (2) Pab =
Pai(1) Pib(1) .
(48)
i
The second-order response of the Fock-matrix is (2) = Fpq
∂Fpq r,s
∂Prs
Prs(2) +
where ∂ 2 Fpq = ∂Prs ∂Pr s
1 ∂ 2 Fpq P (1) P (1) , 2! r,s ∂Prs ∂Pr s rs r s
(49)
r ,s
E=0
The polarizability tensor α ξ ζ (−ω; ω) diverges at resonances. It is therefore assumed that when ω approaches an excitation frequency (resonance) ωn , the external perturbation fai can be adiabatically removed due to the resonant nature of the excitation. Then, the right-hand side of Eq. (41) is zero, and Eq. (41) becomes 1 0 A B Xn Xn = ωn . (44) Yn 0 −1 B A Yn
f0n is an oscillator strength for the 0 → n transition. Equation (45) reveals the divergent nature of the averaged polarizability near the excitation frequencies ωn . In practice, an eigenvalue matrix equation of half the dimension of Eq. (44) is solved (see Appendix B).
ai
Equation (42) is a direct consequence of Eqs. (30) and (40). Note that (1) ∂Ppq = Xξ,ai (ω) + Yξ,ai (ω) . (43) ∂Eξ
αξ ξ (−ω; ω) =
δ 3 EXC [ρ] φ ∗ (r)φq (r) δρ(r)δρ(r )δρ (r ) p
× φr∗ (r)φs (r)φr∗ (r)φs (r)d rd r d r .
(50)
Equation (50) follows from Eqs. (5) and (6). (2) Using the explicit form of Fpq from Eq. (49) and inserting it into Eq. (20) results in an explicit form of the secondorder equation-of-motion:
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.120 On: Tue, 20 May 2014 20:29:47
18A523-6
F. Zahariev and M. S. Gordon
J. Chem. Phys. 140, 18A523 (2014)
⎡
∂Fpq
∂Fqr
⎤
− + hpq + dst dqr − dpq hqr + dst ⎥ ∂Pst ∂Pst ⎥ ⎢ st st ⎥ ⎢ ⎛ ⎞ ⎛ ⎞ ⎥ ⎢ ⎢ ⎥ ⎢ ⎜ ⎟ ⎜ ⎟⎥ 2 2 ⎥ ⎢ ⎜ ∂Fpq (2) ⎟ ⎜ ⎟ ∂ F ∂F ∂ F 1 1 pq qr (2) qr q ⎢ (0) (0) ⎜ ⎥ Pqr dst + dst ds t ⎟ − Ppq dst + dst ds t ⎟ ⎥ ⎢+⎜ ⎜ ⎟ ⎜ ⎟ t 2! ∂Pst ∂Ps t ∂P 2! ∂P ∂P st st s ⎣ ⎝ s,t ∂Pst ⎠ ⎝ s,t ⎠⎦ (0) (2) ⎢ Fpq dqr
(2) (0) dpq Fqr
s,t s ,t
(51)
s,t s ,t
(2) = 2ωdpr .
(2) Recall that only the occupied-unoccupied dai and (2) unoccupied-occupied dia components of the solution are of interest. Using the notation
h(2) ai
=
(2) (2) Xai = dai ,
(52)
(2) Yai(2) = dia ,
(53)
hab +
b
∂Fab ck
∂Fab xck + yck ∂Pck ∂Pkc
(0) (0) and the diagonal form of Fpq and Ppq (see Eqs. (12), (13), and (15)). Equation (51) can be cast into a matrix equation form reminiscent of Eq. (35) (2) (2) A B h 1 0 X (ω) . =− − 2ω h(2) Y (2) (ω) B A 0 −1 (54) The components of h(2) are
xbi −
xaj hj i +
j
∂Fj i ck
∂ 2 Fai ∂ 2 Fai ∂ 2 Fai 1 xbj xb j + 2 xbj yb j + + 2! ∂Pbj ∂Pb j ∂Pbj ∂Pb j ∂Pbj ∂Pb j
∂Fj i xck + yck ∂Pck ∂Pkc ybj yb j ,
(55)
b,j b ,j
where x and y in Eq. (55) are the solutions to the nonhomogeneous linear-response Eq. (35). In Sec. III, the linear response solution was used to compute the polarizability based on Eq. (27) for the polarizability α ξ ζ (−ω; ω). In a similar manner, the second-order response solution gives rise to the SHG hyperpoliarizability β ξ ζ η (−2ω; ω, ω) based on Eq. (28) and to the OR hyperpolarizability β ξ ζ η (0; −ω, ω) based on Eq. (29). The most general form of the hyperpolarzability β ξ ζ η (ωξ ; ωζ , ωη ) can be obtained using an appropriate generalization of Eqs. (28) and (29). The dependence of the hyperpolarizability on the second-order response solution can be further transformed into a dependence on just the linear response solution by using the “(2n+1) theorem”:1 βξ ζ η (ωξ ; ωζ , ωη ) = A + B + C,
(56)
A=
perm.ξ,ζ,η
aij
Xξ,ai (ωξ ) −μζ,ij +
! +Yζ,ck (ωζ )) Yη,aj (ωη ) ,
C=
perm.ξ,ζ,η
iab
perm.ξ,ζ,η
H XC Xξ,ai (ωξ ) −μζ,ab + fab,ck (Xζ,ck (ωζ ) ck
(58)
XC gai,bj,ck (Xξ,ai (ωξ )+Yξ,ai (ωξ ))(Xζ,bj (ωζ )
aibj ck
+Yζ,bj (ωζ ))(Xη,ck (ωη ) + Yη,ck (ωη )) .
(59)
The “permutational sum” (perm. in Eqs. (57)–(59)) of a given quantity Qξ , ζ , η is
Qξ,ζ,η = Qξ,ζ,η + Qξ,η,ζ + Qζ,η,ξ + Qζ,ξ,η
perm.ξ,ζ,η
+Qη,ξ,ζ + Qη,ζ,ξ ,
H XC fij,ck (Xζ,ck (ωζ )
ck
! +Yζ,ck (ωζ )) Yη,bi (ωη ) ,
where any of the indices ξ , ζ , η can be x, y, or z and A, B, and C in Eq. (56) are defined as
B =−
(60)
and (57)
H XC fab,cd = (ia|j b) − cH (ab|j i) + (1 − cH )f XC ia,j b
(61)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.120 On: Tue, 20 May 2014 20:29:47
18A523-7
F. Zahariev and M. S. Gordon
J. Chem. Phys. 140, 18A523 (2014)
are the matrix elements of the combined Hartree-exchangecorrelation kernel and δ 3 E XC [ρ] XC gai,bj,ck = φ ∗ (r)φi (r)φb∗ (r ) δρ(r)δρ(r )δρ(r ) a × φj (r )φc∗ (r )φk (r )d rd r d r
lim (2ω − ωn ) βξ ζ η (−2ω; ω, ω) = σξ ζ (ωn )0|rη |n. (63)
2ω→ωn
Hence, the TPA cross-section itself is σξ ζ (ωn ) = A + B + C.
(64)
The three terms on the right-hand side of Eq. (64) are A=− Xξ,ai (ωξ ) − μζ,ij (1 − δζ n ) +
aij
H XC fij,ck (Xζ,ck (ωζ )
! + Yζ,ck (ωζ )) Yn,aj , (65)
ck
"
B=
perm.ξ,ζ,n
+
Xξ,ai (ωξ ) −ηζ,ab (1 − δζ n )
iab
H XC fabck (Xζ,ck (ωζ )
#
(66) C=−
perm.ξ,ζ,n
XC gai,bj,ck (Xξ,ai (ωξ ) + Yξ,ai (ωξ ))
aibj ck
×(Xζ,bj (ωζ ) + Yζ,bj (ωζ ))(Xn,ck
+ Yn,ck ) .
(67)
The value of the TPA cross section matrix element in Eq. (64) depends on the choice of coordinate system. However, orientationally averaged TPA transition moment quantities can be constructed from the matrix elements of the TPA. For a linear polarization of the light such an orientationally averaged TPA transition moment is 2Df + 4Dg ,
In the ground-state DFT/EFP1 method,39 the interaction between the DFT solute and the EFP water molecules39, 40 is represented by an effective one-body potential that has a Coulomb, polarization, and remainder component. The remainder term, which is fitted to a DFT water dimer potential, contains exchange repulsion, charge transfer, and short-range correlation contributions. The effective potential v EF P 1 (r) is added to the usual Kohn-Sham potential, but since only the EFP polarization term v P OL [ρ](r) depends on the electron density,39 only the EFP1 polarization term is updated during the SCF cycles. A key ingredient in the TDDFT linear response equation XC , is the exchange-correlation kernel f XC (r, r ) = δv δρ(r[ρ](r) ) which is a functional derivative of the exchange-correlation potential v XC [ρ] (r) with respect to the density. Since only the polarization component of the EFP1 effective potential v P OL [ρ] (r) has a non-zero functional derivative, the EFP1 correction to the linear-response TDDFT matrix equation41–44 is given by the following replacement: (73)
The last term in Eq. (73) is the EFP polarization kerP OL . The EFP polarization kernel nel f P OL (r, r ) = δv δρ(r[ρ](r) ) f P OL (r, r ), the derivative of the polarization potential with respect to the density, is expressed through the EFP polarizability tensor.41, 42 For applications of the combined TDDFT/EFP method, see Refs. 43 and 44. In order to extend the TDDFT/EFP1 scheme to nonlinear-response TDDFT, the replacement as given in Eq. (73) is performed in the expressions for the first hyperpolarizabilty and TPA cross-section. Since the EFP polarization potential v P OL [ρ](r) has only a linear dependence on the density, the EFP polarization kernel f P OL (r, r ) has no effective density dependence and hence g P OL (r, r , r ) P OL [ρ](r,r) = δf δρ(r = 0. Therefore, there is no need to correct ) g XC (r, r , r ) =
δf XC [ρ](r,r ) δρ(r )
in Eqs. (59) and (67).
(68) VI. COMPUTATIONAL DETAILS
where Df and Dg are defined as Df =
1 σξ ξ σζ ζ 30 ξ ζ
(69)
Dg =
1 σξ ζ σξ ζ , 30 ξ ζ
(70)
and
while the orientationally averaged TPA transition moment for a circular polarization is −2Df + 6Dg .
(72)
V. THE EFP METHOD
f XC (r, r ) → f XC (r, r ) + f P OL (r, r ).
+ Yζ,ck (ωζ )) Yn,bi ,
ck
(−Df + 3Dg ) . (Df + 2Dg )
(62)
are the matrix elements of the third functional derivative of the exchange-correlation energy density functional. The TPA absorption cross-section is proportional to the following residue:
perm.ξ,ζ,n
The polarization ration is
(71)
As an illustration, the difficult to model hyperpolarizability and two-photon absorption of bulk water will be computed by the new nonlinear-response TDDFT/EFP method as implemented in GAMESS. A system of 400 EFP molecules was equilibrated using a molecular dynamics (MD) simulation in the NVT ensemble. Upon equilibration, a smaller spherical sub-system of 150 EFP waters with a radius of about 9 A was extracted and one water molecule approximately in the middle of the sub-system was selected to be described by B3LYP/DH(d,p). A harmonic restraint potential with a force field of
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.120 On: Tue, 20 May 2014 20:29:47
18A523-8
F. Zahariev and M. S. Gordon
J. Chem. Phys. 140, 18A523 (2014)
The computed average TPA transition moment to the first excited state by two photons of equal energy is 9.9 a.u. with a solvent shift of 4.9 a.u. with respect to the gas phase value of 5.0 a.u. In Ref. 56, the same quantities were computed with a combined coupled cluster/molecular mechanics method to be, respectively, 7.2 a.u., 2.2 a.u., and 5.0 a.u. Reference 56 discusses the favorable qualitative comparison of the computed liquid phase TPA values with experiment57, 58 but there is no experimental basis to validate the computed TPA solvent shifts. VII. CONCLUSION
FIG. 1. The first excited state energy, second harmonic generation (SHG) hyper-polarizability (at ω = 0.043 a.u.) and the two-photon absorption transition moment and polarization ratio corresponding to the first excited state of water as a function of the simulation time (blue). The time averages are in green.
10 kcal/mol/A2 , effectively keeping the system in the NPT ensemble, was applied at the spherical boundary in order to prevent the evaporation of the water droplet.51, 52 The combined DFT/EFP system was equilibrated again for 55 ps and the next 45 ps production run was used to collect and average the optical properties of interest. All of the MD runs were done at 298 K using the velocity-Verlet algorithm with a time step of 0.5 fs. TDDFT (B3LYP)/EFP linear and nonlinear optical properties of water were computed on snapshots from the production run of the ground state DFT/EFP MD simulation. The d-aug-cc-pVDZ basis set was used for a proper description of the NLO properties. The snapshots from the MD trajectory were chosen with a separation of 250 fs to allow enough time for de-correlation between the snapshots to occur. Figure 1 shows both the running values (blue) and their time-averages (green) for several properties to illustrate the type and magnitude of statistical fluctuations during the MD simulation. All of the gas-phase optical values were computed with TDDFT using B3LYP/d-aug-cc-pVDZ. The average first excited state (11 B1 ) excitation energy is computed to be 7.72 eV with a solvent shift of 0.5 eV with respect to the gas phase first excited state energy of 7.22 eV. The corresponding experimental values are 8.2 eV for the liquid phase with a solvent shift of 0.8 eV relative to the 7.4 eV gas phase value.53 According to the experimental data, the sign of the SHG hyperpolarizability of liquid water changes relative to that of gas phase water. The average computed SHG hyperpolarizability of liquid water (11.8 a.u.) at the frequency ω = 0.043 a.u. compares reasonably well with the experimental value of 15.5 a.u. at the same frequency.54 The gas phase TDDFT/B3LYP/d-aug-cc-pVDZ SHG hyperpolarizability of water at ω = 0.0656 is predicted to be −5.9 a.u., again in reasonable agreement with the experimental value of −10.9 a.u. at the same frequency.55 Importantly, the calculations presented here predict the observed sign change.
A nonlinear response TDDFT computational scheme has been described that can successfully be combined with the EFP method to account for solvent effects. The combined method is able to correctly predict both gas phase values and aqueous solvent shifts for several important nonlinear properties. More details on the prediction of solvent effects (including non-aqueous solvents) within the presented TDDFT/EFP scheme will be described in subsequent publications. ACKNOWLEDGMENTS
This research is supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences, through the Ames Laboratory Chemical Physics program. The Ames Laboratory is operated for the U.S. Department of Energy by Iowa State University under Contract No. DE-AC02–07CH11358. APPENDIX A: THE IDEMPOTENCE RELATION OF THE DENSITY MATRIX
Expand the density-matrix in the basis of natural spinorbitals ρ(r, r ) =
N
φi (r)φi∗ (r ).
(A1)
i=1
Using this expansion it is easy to see that ρ r, r = ρ(r, r )ρ(r , r )d r
(A2)
holds. This is the idempotence relation of the density matrix in the coordinate representation. In order to get the form used in the text, insert the expansion of the density matrix in terms of the Kohn-Sham spinorbitals in the above equation. APPENDIX B: THE MODIFIED LINEAR-RESPONSE TDDFT EQUATION AND THE BIORTHOGONALITY CONDITION
The matrix equation can be simplified by a unitary transformation, which is equivalent to forming new block rows by adding and subtracting the original two block rows and thus forming two matrix equations of half-size, i.e., ( A + B) (X + Y ) = ω (X − Y )
(B1)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.120 On: Tue, 20 May 2014 20:29:47
18A523-9
F. Zahariev and M. S. Gordon
J. Chem. Phys. 140, 18A523 (2014)
and ( A − B) (X − Y ) = ω (X + Y ) .
(B2)
Multiplying both sides of the first matrix equation on the left by ( A − B) and using the second matrix equation results in a matrix eigenvalue equation ( A − B) ( A + B) (X + Y ) = ω2 (X + Y ) .
(B3)
Equation (B3) is half the size of the original matrix eigenvalue equation. This is still a non-Hermitian problem though. (X + Y )i and (X − Y )i are the right- and left-eigenvectors of the non-Hermitian matrix ( A − B) ( A + B)and the usual biorthonormality condition (X + Y )i · (X − Y )j = δij
(B4)
is imposed upon them, which condition could be also expressed as X i · X j − Y i · Y j = δij .
(B5)
The matrix equation can be modified further by multiplying both sides on the left by ( A − B)−1/2 and redefining the eigenvectors ( A − B)1/2 ( A + B) ( A − B)1/2 X˜ + Y˜ = ω2 X˜ + Y˜ , (B6) where X˜ + Y˜ = ( A − B)−1/2 (X + Y ). This is a Hermitian eigenvalue problem. The actual linear-response TDDFT implementation in GAMESS uses the Davidson algorithm based on both the Hermitian and non-Hermitian eiqenproblem formulations (see Ref. 59 for details). APPENDIX C: THE SPECTRAL RESOLUTION METHOD
The eigenvectors φi (r) of a Hermitian operator L, i.e., Lφi (r) = λi φi (r)
(C1)
form an orthogonal basis, so any function, say u(r), can be expanded in it, i.e., u(r) =
∞
αi φi (r).
(C2)
i=1
The action of L on u(r) becomes Lu(r) =
∞
αi λi φi (r).
(C3)
i=1
This representation is convenient for obtaining a function of L, i.e., f (L) u(r) =
∞
αi f (λi ) φi (r)
(C4)
i=1
and in the particular case of f (L) = (L − λ)−1 , it becomes −1
(L − λ)
u(r) =
∞ αi φi (r) i=1
λi − λ
.
(C5)
For additional information see Ref. 60. The solutions of the non-homogenous matrix equation in Eq. (41) can be expressed in terms of the eigenvectors and
eigenvalues of the homogenous matrix equation in Eq. (44) by employing the spectral resolution method: Yn Xn + (C6) [μξ · (Xn + Yn )] Xξ (ω) = ωn − ω ωn + ω n and Yξ (ω) =
n
[μξ · (Xn + Yn )]
Yn Xn + . (C7) ωn − ω ωn + ω
The scalar products on the right-hand side of Eqs. (46) and (47) are expressed in terms of the ai-matrix elements. The above expansions have a formal similarity with the sum-overstate (SOS) method, but the basis of the usual SOS expansion is different. 1 D.
Pugh, “Electric multipoles, polarizabilities and hyperpolarizabilities,” in Chemical Modeling: Applications and Theory, edited by A. Hinchliffe (RSC Publishing, 2000), Vol. 1. 2 M. Göppert-Mayer, Ann. Phys. 401, 273 (1931). 3 W. Kaiser and C. G. B. Garret, Phys. Rev. Lett. 7, 229 (1961). 4 I. D. Abella, Phys. Rev. Lett. 9, 453 (1962). 5 F. A. Moscatelli, Am. J. Phys. 54, 52 (1986). 6 J. D. Bhawalkar, G. S. He, and P. N. Prasad, Rep. Prog. Phys. 59, 1041 (1996); C. R. Shea, Y. Hefetz, R. Gillies, J. Wimberly, G. Dalickas, and T. Hasan, J. Biol. Chem. 265, 5977 (1990); P. Lenz, Photochem. Photobiol. 62, 333 (1995); J. D. Bhawalkar, N. D. Kumar, C. F. Zhao, and P. N. Prasad, J. Clin. Laser Med. Surg. 15, 201 (1997); W. G. Fisher, W. P. Partridge, C. Dees, and E. A. Wachter, Photochem. Photobiol. 66, 141 (1997). 7 D. Tunega and J. Noga, Theor. Chem. Acc. 100, 78 (1998). 8 A. B. Kumar, N. Vaval, and S. Pal, Chem. Phys. Lett. 295, 189 (1998). 9 O. Christiansen, A. Halkier, H. Koch, P. Jorgensen, and T. Helgaker, J. Chem. Phys. 108, 2801 (1998). 10 O. Chrisitiansen, J. Gauss, and J. F. Stanton, Chem. Phys. Lett. 292, 437 (1998). 11 P. B. Rozyczko and R. J. Bartlett, J. Chem. Phys. 109, 9201 (1998). 12 B. Fernandez, C. Hattig, H. Koch, and A. Rizzo, Chem. Phys. Lett. 288, 677 (1998). 13 P. Cronstrand, Y. Luo, and H. Agren, Chem. Phys. Lett. 352, 262 (2002). 14 P. Cronstrand, Y. Luo, and H. Agren, J. Chem. Phys. 117, 11102 (2002). 15 M. Drobizhev, Y. Stepanenko, Y. Dzenis, A. Karotki, A. Rebane, P. N. Taylor, and H. L. Anderson, J. Am. Chem. Soc. 126, 15352 (2004). 16 S. Delysse, P. Raimond, and J. M. Nunzi, Chem. Phys. 219, 341 (1997). 17 B. J. Orr and J. F. Ward, Mol. Phys. 20, 513 (1971). 18 C. Lambert, E. Schmalzlin, K. Meerholz, and C. Brauchle, Chem. – Eur. J. 4, 512 (1998). 19 U. S. Choi, T. W. Kim, S. W. Jung, and C. J. Kim, Bull. Korean Chem. Soc. 19, 299 (1998). 20 F. W. Vance, L. Karki, J. K. Reigle, J. T. Hupp, and M. A. Ratner, J. Phys. Chem. A 102, 8320 (1998). 21 Y. K. Lee, S. J. Jeon, and M. H. Cho, J. Am. Chem. Soc. 120, 10921 (1998). 22 P. Wang, C. G. Wang, P. W. Zhu, and C. Ye, Sci. China, Ser. B: Chem. 41, 156 (1998). 23 F. Grozema, R. Telesca, H. Jonkman, L. Siebbeles, and J. Snijders, J. Chem. Phys. 115, 10014 (2001). 24 F. Grozema, R. Telesca, J. Snijders, and L. Siebbeles, J. Chem. Phys. 118, 9441 (2003). 25 F. C. Grozema and P. T. van Duijnen, J. Phys. Chem. A 102, 7984 (1998). 26 A. Zangwill and P. Soven, Phys. Rev. A 21, 1561 (1980). 27 A. Dreuw, J. L. Weisman, and M. Head-Gordon, J. Chem. Phys. 119, 2943 (2003). 28 A. Dreuw and M. Head-Gordon, J. Am. Chem. Soc. 126, 4007 (2004). 29 D. J. Tozer, J. Chem. Phys. 119, 12697 (2003). 30 O. Gritsenko and E. J. Baerends, J. Chem. Phys. 121, 655 (2004). 31 M. E. Casida, C. Jamorsky, K. C. Casida, and D. R. Salahub, J. Chem. Phys. 113, 8918 (2000). 32 A. Dreuw and M. Head-Gordon, Chem. Rev. 105, 4009 (2005). 33 P. Elliott, F. Furche, and K. Burke, “Excited states from time-dependent density functional theory,” in Reviews in Computational Chemistry, edited by K. B. Lipkovitz and T. R. Cundari (Wiley, Hoboken, NJ, 2009), p. 91.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.120 On: Tue, 20 May 2014 20:29:47
18A523-10 34 J.
F. Zahariev and M. S. Gordon
Tomasi and M. Persico, Chem. Rev. 94, 2027 (1994). Tomasi, B. Mennucci, and R. Cammi, Chem. Rev. 105, 2999 (2005). 36 G. Scalmani, M. J. Frisch, B. Mennucci, J. Tomasi, R. Cammi, and V. Barone, J. Chem. Phys. 124, 094107 (2006). 37 M. S. Gordon, L. Slipchenko, L. Hui, and J. H. Jensen, Annu. Rep. Comput. Chem. 3, 177 (2007). 38 M. S. Gordon, M. A. Freitag, P. Bandyopadhyay, V. Kairys, J. H. Jensen, and W. J. Stevens, J. Phys. Chem. A 105, 293 (2001). 39 I. Adamovic, M. A. Freitag, and M. S. Gordon, J. Chem. Phys. 118, 6725 (2003); H. Li, H. M. Netzloff, and M. S. Gordon, ibid. 125, 194103 (2006); P. N. Day, J. H. Jensen, M. S. Gordon, S. P. Webb, W. J. Stevens, M. Krauss, D. Garmer, H. Basch, and D. Cohen, ibid. 105, 1968 (1996). 40 J. H. Jensen and M. S. Gordon, J. Chem. Phys. 108, 4772 (1998). 41 S. Yoo, F. Zahariev, S. Sok, and M. S. Gordon, J. Chem. Phys. 129, 144112 (2008). 42 N. Minezawa, N. De Silva, F. Zahariev, and M. S. Gordon, J. Chem. Phys. 134, 054111 (2011). 43 S. Sok, S. Y. Willow, F. Zahariev, and M. S. Gordon, J. Chem. Phys. A 115, 9801 (2011). 44 A. DeFusco, N. Minezawa, L. V. Slipchenko, and M. S. Gordon, J. Phys. Chem. Lett. 2, 2184 (2011). 45 M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. J. Su, T. L. Windus, M. Dupuis, J. A. Montgomery, J. Comput. Chem. 14, 1347 (1993). 35 J.
J. Chem. Phys. 140, 18A523 (2014) 46 M.
S. Gordon, M. W. Schmidt, in Theory and Applications of Computational Chemistry: The First Forty Years, edited by C. E. Dykstra, G. Frenking, K. S. Kim, and G. E. Scuseria (Elsevier, Amsterdam, 2005), Chap. 41, pp. 1167–1189. 47 S. Hirata and M. Head-Gordon, Chem. Phys. Lett. 314, 291 (1999). 48 F. Furche, J. Chem. Phys. 114, 5982 (2001). 49 F. Wang, C. Y. Yam, and G. Chen, J. Chem. Phys. 126, 244102 (2007). 50 R. McWeeny, Methods of Molecular Quantum Mechanics (Academic Press, 1992). 51 D. Beglov and B. Roux, J. Chem. Phys. 100, 9050 (1994). 52 C. H. Choi, S. Re, M. Feig, and Y. Sugita, Chem. Phys. Lett. 539–540, 218 (2012). 53 W. F. Chan, G. Cooper, and C. E. Brion, Chem. Phys. 178, 387 (1993). 54 B. F. Levine and C. G. Bethea, J. Chem. Phys. 65, 2429 (1976). 55 J. W. Ward and K. Miller, Phys. Rev. A 19, 826 (1979). 56 M. J. Paterson, J. Kongsted, O. Christiansen, K. V. Mikkelsen, and C. B. Nielsen, J. Chem. Phys. 125, 184501 (2006). 57 D. N. Nikogosyan, A. A. Oraevsky, and V. I. Rapusov, Chem. Phys. 77, 131 (1983). 58 C. L. Thomsen, D. Madsen, S. R. Keiding, J. Thogersen, and O. Christiansen, J. Chem. Phys. 110, 3453 (1999). 59 R. E. Stratmann and G. E. Scuseria, J. Chem. Phys. 109, 8218 (1998). 60 J. P. Keener, Principles of Applied Mathematics: Transformation and Approximation (Perseus Books, Cambridge, MA, 2000).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.120 On: Tue, 20 May 2014 20:29:47
