9116
J. Phys. Chem. 1996, 100, 9116-9126
Molecular Response Method for Solvated Molecules in Nonequilibrium Solvation Kurt V. Mikkelsen* Department of Chemistry, H. C. Ørsted Institute, UniVersity of Copenhagen, DK-2100 Copenhagen Ø
Kristian O. Sylvester-Hvid Department of Chemistry, Aarhus UniVersity, DK-8000 A° rhus C, Denmark ReceiVed: NoVember 30, 1995; In Final Form: February 28, 1996X
We present an extension of the molecular response method of solvated compounds. The extension treats in a proper fashion the effects of nonequilibrium solvation on the molecular properties of the solvated compound. The molecular properties evaluated within the response formalism are properties that are measured using high-frequency electromagnetic fields. We consider the changes in the molecular properties due to nonequilibrium solvation. The nonequilibrium solvation arises from changes in the molecular charge distribution induced by the high-frequency perturbation. The optical polarization vector is at all times assumed to be in equilibrium with the charge distribution of the solute, contrary to the inertial polarization, which is not necessarily in equilibrium with the charge distribution of the solute. We calculate shifts of excitation energies and frequency-dependent polarizabilities as a function of the optical and static dielectric constants.
1. Introduction The theoretical efforts for calculating molecular properties in the gas phase have successfully reached the level of reproducing, predicting, interpreting, and correcting many gas phase experiments.1,2 The obvious need for developing methods for investigating the effects of the condensed phase on molecular properties stems from the fact that the majority of experimental measurements occur in condensed phases: solutions, crystals, polymers, or glasses. The interactions between the molecular compound and the condensed phase can be of both long- and short-range nature. Solvent effects on molecular properties have been addressed by a large number of research groups. An improved understanding of how solvents affect the molecular properties and chemical reactions has been reached through these collective efforts starting from considerations of solvation energies,3-6 developing into an increasingly important research area within theoretical chemistry.7-32 For a more detailed discussion we refer to the recent review by Tomasi and Persico on molecular interactions in solution.33. Many-electron structure methods for describing solventinduced effects on molecular properties range from supermolecular models to continuum models, semicontinuum models being a compromise. Supermolecular models take account of the short-range interactions between the solute and the solvation shells, totally neglecting long-range interactions. On the other hand, continuum models disregard short-range interactions but take account of the long-range interactions. Semicontinuum models attempt to take the best from each of the two models by including a solvation shell around the solute and surrounding the resulting molecular complex by a dielectric medium. Models accounting for nonequilibrium solvation effects on timedependent molecular properties within a continuum model have not reached a satisfactory level. In fact, few models for solving the time-dependent quantum mechanical equations for a molecular system coupled to a solvent exist.34-43 Therefore, this article presents a novel method for treating nonequilibrium solvation states (the solvent is described as a dielectric medium), in conjunction with a recently established response method.41 X
Abstract published in AdVance ACS Abstracts, May 1, 1996.
S0022-3654(95)03540-4 CCC: $12.00
This involves solving the time-dependent response equations for a solvated molecule in a nonequilibrium solvation state due to an external, high-frequency perturbation. A continuum model for investigating the effects of the solvent on molecular properties and chemical reactions is one where a solute, being a molecule or a molecular complex, is surrounded by a continuous and polarizable medium. To account for the occurrence of nonequilibrium solvent states, optical and inertial polarization vectors are introduced and a free energy functional is derived for theses states.44-49 The medium is characterized by two dielectric constants op (the electronic response of the solvent molecules) and st (the full response of the solvent). The description of this model is, therefore, one of equilibrated quantized fields for the solute electrostatic charge distribution and for the solvent response to this charge distribution through an optical polarization vector. These two fields are, however, dependent on an additional arbitrary (classical) polarization vector field, which is slow to incorporate the changes in the electronic structure of the composite solute-solvent system. This field introduces the contribution of the inertial (slow or nuclear) modes of the solvent molecules to the electrostatic free energy of the total system. The optical polarization vector is evaluated as one being in equilibrium with the charge distribution of the solute which leads to the semiclassical nonequilibrium or equilibrium approaches.17,34,41,44,50-57 The development of electronic structure theory of the solutesolvent/dielectric model has occurred within both semiempirical and ab-initio frameworks, and the implementations have used either dipole fields, extended multipole expansion fields, or fully numeric (point charge) electrostatic fields to represent the molecular charge distribution. The shape of the cavity enclosing the solute is either spherical17,34,41,53,56-61 or nonspherical54,55,62,63 Recently, uncorrelated and correlated electronic structure methods, taking accurate account of nonequilibrium solvation, have been developed.64,65 However, they do not enable response calculations of time-dependent molecular properties. The present article concerns the incorporation of the optical and inertial polarization vectors into the response methodology where the molecules are surrounded by a dielectric medium. The formalism is an extension of previous work34-40 and © 1996 American Chemical Society
Solvated Molecules in Nonequilibrium Solvation
J. Phys. Chem., Vol. 100, No. 21, 1996 9117
generalizes the response formalism presented in ref 41 to handle nonequilibrium solvation states which occur during the application of external high-frequency laser fields to the solvated molecule. The energy functional of the nonequilibrium solvation state arising due to changes in the solute induced by the external high-frequency electromagnetic field is given by
2. Theoretical Section
Esol ) ∑Rlm(Ff,op)〈Tlm(Ff)〉 + lm
remains fixed during the period of interaction with the highfrequency electromagnetic field. The Results and Discussion section demonstrates the significance of introducing the nonequilibrium solvation states in the calculation of solute properties such as excitation energies and frequency-dependent polarizabilities.
Rlm(Fi,st,op)[2〈Tlm(Ff)〉 - 〈Tlm(Fi)〉] ∑ lm
The Hamiltonian of a solvated molecular system perturbed by an external field is
(1)
where Esol is the energy of the state that exists in nonequilibrium solvation. The external perturbation has induced changes in the solute giving rise to a new charge distribution of the solute denoted as Ff. The charge distribution of the solute, prior to the presence of the external high-frequency perturbation, is denoted Fi. For a spherical cavity the polarization fields due to the induced polarization charges in the dielectric medium are
Rlm(Ff,op) ) gl(op)〈Tlm(Ff)〉
(2)
Rlm(Fi,st,op) ) gl(st,op)〈Tlm(Fi)〉
(3)
where the gl factors are given in terms of the dielectric constants (st,op), the cavity radius (Rcav), and the multipole expansion parameter (l):
H ) H0 + Wsol + V(t)
(11)
where the first term, H0, is the Hamiltonian of an isolated molecule, and the second term is the interaction operator describing the interactions between the solute and the dielectric medium. The third term is the interaction operator, V(t), which accounts for the interaction between the molecule and the external field. The time-dependent perturbation operator can be written in terms of the Fourier integral:
V(t) ) ∫-∞dω Vω exp[(-iω + η)t] ∞
(12)
where η is a positive infinitesimal number which ensures that
V(t f -∞) ) 0
(13)
and where Vω is the Fourier transform of V(t). Due to η, the perturbation is applied adiabatically in the limit t f -∞. We optimize our reference state, |0〉, such that
(l + 1)( - 1) 1 gl() ) R-(2l+1) cav 2 l + (l + 1)
(4)
(H0 + Wsol)|0 〉 ) E0|0〉
gl(st,op) ) gl(st) - gl(op)
(5)
is valid. In the situation where the external time-dependent electromagnetic field is absent we recover the generalized Brillouin condition for the optimized wave functions |0〉:
The charge moments 〈Tlm(F)〉 are given as n e - 〈Tlm (F)〉 〈Tlm(F)〉 ) Tlm
(6)
n ) ∑Zitlm(Ri) Tlm
(7)
e ) ∑tlm Tlm pqEpq
(8)
i
〈0|[λ, H0 + Wsol]|0〉 ) 0
lm tlm pq ) 〈φp|t (r)|φq〉
(9)
(15)
where λ denotes the orbital and configurational variation parameters. The response of the molecular system to a time-dependent perturbation will, as in refs 34 and 41, be governed by Ehrenfest’s equation of motion:
pq
where the nuclear charge and position vector are given as Zi and Ri, respectively. Further, the electronic spherical moment integrals are given as
(14)
〈 〉
∂B d - i 〈[B,H]〉 〈B〉 ) dt ∂t
(16)
where H is given in eq 11 and B is an arbitrary operator. The expectation values are determined using the time-dependent state |0t〉, such that
〈...〉 ) 〈t0|...|0t〉
(17)
σ
|0t〉 ) exp(iκ(t)) exp(iS(t))|0〉
(18)
and aqσ are the creation and annihilation operators where for an electron in spin orbital φj. For further details see refs 34, 41, and 56. The theoretical section of this paper describes the mathematical derivation of the time-dependent equations taking account of the different relaxation times of the polarization vectors of the dielectric medium. The optical polarization vector changes instantaneously with respect to changes of the solute’s charge distribution, as opposed to the inertial polarization vector, which
where exp(iκ(t)) performs unitary transformations in the orbital space and where κ(t) is given as
and the electronic spin-free operators are given as † Epq ) ∑apσ aqσ † apσ
where
(10)
κ(t) ) ∑(κk(t)q†k + κ*k (t)qk)
(19)
q†k ) Epq; p > q
(20)
k
with
9118 J. Phys. Chem., Vol. 100, No. 21, 1996
Mikkelsen and Sylvester-Hvid
Similarly, the unitary transformations in configuration space are given by exp(iS(t)) with
S(t) ) ∑(Sn(t)R†n + S*n(t)Rn)
(21)
n
illustrates that expectation values of effective solvent Hamiltonians are not the solvent contribution to the total energy.41,44,64 2.2. Ehrenfest’s Equation of Motion. The operators to be considered in Ehrenfest’s equation of motion are the timetransformed operators T†,t defined as
()
with the state transfer operators Rn defined as
R†n
) |n〉〈0|
(22)
and where |n〉 belongs to the orthogonal complement space to |0〉. 2.1. Effective Solvent Hamiltonian. The effective solvent Hamiltonian must give the same expressions for the timeindependent properties as those obtained by differentiation of the energy expression. This is inherent to the utilization of Frenkel’s variation principle, which for approximate wave functions gives the same results for time-independent properties as energy differentiation or finite field calculations. Utilizing the same strategy as in ref 41, where the analogy with the gradient of the energy functional is used to define the effective solvent Hamiltonian Wsol, we obtain the following nonequilibrium expression: n 2 n 2 ) + ∑g(st,op)(Tlm ) + Tg(st,op) Wsol ) ∑gl(op)(Tlm lm
lm
(23)
The operator Tg(st,op) is the sum of two very similar contributions, Tg1(op) and Tg2(st,op), which take the forms
Tg1(op) ) ∑(-2∑gl(op)〈Tlm(Ff)〉tlm pqEpq)
(24)
Tg2(st,op) ) ∑(-2∑gl(st,op)〈Tlm(Fi)〉tlm pqEpq)
(25)
pq
pq
lm
lm
One may recognize the operator Tg1(op) as that referred to as Tg in the equilibium studies of ref 41, where the distinction between Fi and Ff is nonexistent. Considering the expectation value of the solvent Hamiltonian for a state corresponding to the charge distribution Ff leads to n 2 n 2 ) + ∑g(st,op)(Tlm ) + 〈Wsol〉 ) ∑gl(op)(Tlm lm
lm
〈Tg(op)〉 + 〈Tg(st,op)〉
n 2 e ) - 2〈Tlm(Ff)〉〈Tlm (Ff)〉} + ) ∑gl(op){(Tlm lm
∑ lm
n 2 gl(st,op)(Tlm )
- 2∑
e gl(st,op) 〈Tlm(Fi)〉〈Tlm (Ff)〉
lm
e (Ff)〉2) + ) ∑gl(op)(〈Tlm(Ff)〉 + 〈Tlm 2
lm
gl(st,op)(2〈Tlm(Ff)〉〈Tlm(Fi)〉 + ∑ lm
e (Ff)〉2) + ) ∑gl(op)(〈Tlm(Ff)〉2 + 〈Tlm
qtk ) exp(iκ(t))qk exp(-iκ(t))
(28)
† q†,t k ) exp(iκ(t))qk exp(-iκ(t))
(29)
Rtn ) exp(iκ(t)) exp(iS(t))Rn exp(-iS(t)) exp(-iκ(t))
(30)
† R†,t n ) exp(iκ(t)) exp(iS(t))Rn exp(-iS(t)) exp(-iκ(t))
(31)
We obtain the response functions by solving Ehrenfest’s equation for each order in the perturbation strength. When solving this equation, we need to evaluate the expectation value of the commutator [T†,t, H], which implies that
-i〈[T†,t,H]〉 ) -i〈[T†,t,H0]〉 - i〈[T†,t,Vpert(t)]〉 i〈[T†,t,Wsol]〉 (32) of which the first two terms have been investigated and implemented; see ref 66. We need to consider the terms arising from -i〈[T†,t,Wsol]〉, which can be written as
-i〈[T†,t,Wsol]〉 ) -i〈[T†,t,Tg1(op)]〉 - i〈[T†,t,Tg2(st,op)]〉 (33) Thus, the commutator involving Wsol gives rise to two terms containing the operators defined in eqs 24 and 25. The structure of the first term has already been discussed and implemented; see ref 41. Inserting the definitions of the operators Tg1(op) and Tg2(st,op), we obtain
gl(op)tlm ∑ pqEpq × lm,pq
lm t n 〈 0|Ep′q′|0t〉 - Tlm )}]|0t〉 (∑ tp′q′ p′q′
i〈t0|[T†,t,{-2
∑gl(st,op)(2〈Tlm(Ff)〉〈Tlm(Fi)〉 -
(27)
where the individual operators are given as
-i〈[T†,t,Wsol]〉 ) -i〈t0|[T†,t,{2
n 2 n ) - 2〈Tlm(Fi)〉Tlm ) (Tlm
lm
qt Rt T†,t ) †,t q R†,t
t gl(st,op)〈Tlm(Fi)〉tlm ∑ pqEpq}]|0 〉 lm,pq
(34)
lm
〈Tlm(Fi)〉〈Tlm(Fi)〉 + 〈Tlm(Fi)〉2)
e e (Ff)〉2 + ∑gl(st,op)〈Tlm (Fi)〉2 ) Esol + ∑gl(op)〈Tlm lm
lm
which may be written in terms of three contributions, two of (a) (b) , are known from the previous equilibrium and Gsol which, Gsol studies from ref 41. Thus, we have
(26)
Esol is the nonequilibrium solvation energy given in eq 1. This
(a) (b) †,t (c) (T†,t) + Gsol (T ) + Gsol (T†,t) -i〈[T†,t,Wsol]〉 ) Gsol
(35)
Solvated Molecules in Nonequilibrium Solvation
J. Phys. Chem., Vol. 100, No. 21, 1996 9119 Inserting the definitions of the orbital rotation operator κ(t) yields the following in the case of qtk:
with the following definitions: (a) Gsol (T†,t) ) -i〈t0|[T†,t,{-2
∑ gl(op)tlmpqTlmn Epq}]|0t〉
(36)
lm,pq
(c,1) t Gsol (qk) ) 2
∑ gl(op,st)〈Tlm(Fi)〉{(〈0|[qk,Epq]|0R〉 +
lm,pq
lm 〈0L|[qk,Epq]|0〉)tlm pq + Qpq〈0|[qk,Epq]|0〉} (46)
(b) †,t Gsol (T ) )
-i〈t0|[T†,t,{2
lm t t t gl(op)tlm ∑ pqtp′q′〈 0|Ep′q′|0 〉Epq}]|0 〉 lm,pq,p′q′
(c) Gsol (T†,t) )
-i〈t0|[T†,t,{-2
∑
(37)
and an analogous expression results in the case of q†,t k . we obtain For the operators Rtn and R†,t n
∑ gl(op,st)〈Tlm(Fi)〉 × {(Sn(t)〈0|Epq|0〉 +
(c,1) t (Rn) ) 2 Gsol t gl(st,op)〈Tlm(Fi)〉tlm pqEpq}]|0 〉 (38)
lm,pq
lm 〈n|Epq|0R〉)tlm pq + Qpq〈n|Epq|0〉} (47)
lm,pq
2.3. Contribution of Nonequilibrium Solvation Terms to Linear Response Equations. We introduce the following states
|0R〉 ) -∑SnR†n|0〉 ) -∑Sn|n〉 n
and (c,1) †,t (Rn ) ) 2 Gsol
(39)
gl(op,st)〈Tlm(Fi)〉{(S′n(t)〈0|Epq|0〉 ∑ lm,pq lm 〈0L|Epq|n〉)tlm pq - Qpq〈0|Epq|n〉} (48)
n
and
〈0L| ) ∑〈0|(S′nR†n) ) ∑S′n〈n| n
(40)
n
We also introduce the one-index transformed solvent integrals as lm Qlm ∑ pqEpq ) ∑[κ(t),tpqEpq] pq pq
(41)
Although calculated in previous work,41 we repeat in appendix (a) (b) A the results for the terms Gsol and Gsol . (a) (b) The addition of the linear contributions from Gsol , Gsol , and (c) Gsol is done most efficiently by defining the following operators:
Tg ) -2
{gl(op)〈Tlm〉 + gl(op,st)〈Tlm(Fi)〉}tlm ∑ pqEpq lm,pq
Tyo ) -2
{gl(op)〈Tlm〉 + gl(op,st)〈Tlm(Fi)〉}Qlm ∑ pqEpq lm,pq
(50)
where
gl(op)〈Qlm〉tlm ∑ pqEpq lm,pq
(51)
e e gl(op){〈0|Tlm |0R〉 + 〈0L|Tlm |0〉}tlm ∑ pqEpq lm,pq
(52)
Txo ) 2 lm lm Qlm pq ) ∑[κpstsq - tps κsq]
(42)
s
(c) for the qtk component of the vector T†,t Consideration of Gsol leads to
(c) (qk) ) -i〈0|exp(-iS(t)) exp(-iκ(t))[exp(iκ(t))qk × Gsol
exp(-iκ(t)),{-2
gl(st,op)〈Tlm(Fi)〉tlm ∑ pqEpq}] × lm,pq
exp(iκ(t)) exp(iS(t))|0〉 (43) Linear response calculations require terms linear in either κ(t) or S(t). Expanding the exponentials yields the following firstorder contribution: (c,1) (qk) ) 2 Gsol
∑ gl(st,op)〈Tlm(Fi)〉tlmpq ×
Txc ) 2
Next, using these effective operators, one may collect the (a,1) (b,1) (c,1) , Gsol , and Gsol terms for the individual elements of the Gsol vector T†,t such that we arrive at
-i〈[qk,Wsol]〉 ) -(〈0L|[qk,Tg]|0〉 + 〈0|[qk,Tg]|0R〉 + 〈0|[qk,Txc|0〉) - 〈0|[qk,Tyo+Txo|0〉 (53) -i〈[q†k ,Wsol]〉 ) -(〈0L|[q†k ,Tg]|0〉 + 〈0|[q†k ,Tg]|0R + 〈0|[q†k ,Txc|0〉) - 〈0|[q†k ,Tyo + Txo|0〉 (54) -i〈[Rn,Wsol]〉 ) -(〈n|Tg|0R〉 + 〈0|Tg|0〉Sn(t) + 〈n|Txc|0〉) 〈n|Tyo + Txo|0〉 (55)
lm,pq
{〈0|[S(t),[qk,Epq]]|0〉 + 〈0|[qk,[κ(t),Epq]]|0〉} (44) For the operator q†k we obtain a similar expression by exchanging qk with q†k . In the case of the configurational operator Rtn, we obtain (c,1) t (Rn) ) 2 Gsol
(49)
∑ gl(st,op)〈Tlm(Fi)〉tlmpq ×
-i〈[R†n,Wsol]〉 ) -(〈0L|Tg|n〉 + 〈0|Tg|0〉S′n(t) - 〈n|Txc|0〉) + 〈n|Tyo + Txo|0〉 (56) The implementation of these equations allows us to investigate the effects of nonequilibrium solvation on the molecular properties of the solute. We use the same implementation strategy as in ref 41.
lm,pq
{
