J. Phys. Chem. B 2007, 111, 8965-8973
8965
Degenerate Four-Wave Mixing in Solution by Cubic Response Theory and the Polarizable Continuum Model Lara Ferrighi, Luca Frediani,* and Kenneth Ruud Centre for Theoretical and Computational Chemistry, Department of Chemistry, UniVersity of Tromsø, N-9037 Tromsø, Norway ReceiVed: March 16, 2007; In Final Form: April 19, 2007
We have derived and implemented the solvent contribution to the cubic response function for the polarizable continuum model in its integral equation formulation. The present formulation is valid both at the HartreeFock and at the Kohn-Sham density functional levels of theory, because both bear the same formal description of the solvent contribution through an additional additive term in the Hamiltonian operator. This new implementation is used to compute the solvent effect for the degenerate four-wave mixing process on three different classes of heteroaromatic chromophores: (I) triazine, benzoxazole, benzimidazole, and benzothiazole; (II) three benzothiazole derivatives; and (III) three tri-s-triazine derivatives. The results are first discussed in terms of the solvent effects on the calculated property and then compared to experimental data for the substrates of class (I) and (II).
I. Introduction Degenerate four-wave mixing (DFWM) is the generation of a wave at frequency ω due to the interaction of three waves at the same frequency. It was first observed by Maker et al.1 and later by Hellwarth.2 Its effect can be viewed as an intensitydependent refractive index3 and can be exploited in building ultrafast “all-optical” switches,4,5 where the refractive index variation is due to an oscillating field, as opposed to other phenomena such as the Kerr effect,3 where the modulation of the refractive index is triggered by a static field. At the molecular level, the DFWM process corresponds to the determination of the frequency-dependent second hyperpolarizability γ(-ω; ω,-ω,ω). Within the framework of quantum mechanics, γ(-ω; ω,-ω,ω) can be obtained by solving the time-dependent Schro¨dinger equation to the fourth order. A viable approach for obtaining the solution is given by response theory,6 which allows for the theoretical modeling of virtually all nonlinear optical phenomena through the same formalism: the investigated process is expressed in terms of one or more sum-over-states (SOS) expressions involving the different electronic (and possibly also vibrational) states of the system and the operators corresponding to the applied fields. Each SOS expression is then converted into its response function counterpart. The important feature of the response formalism is that it does not require any knowledge about the excited-state spectrum of the molecule. Indeed, only the actual state (generally the ground state) wave function or density is in general required. In order to employ response theory for real-world applications, it is important to be able to perform calculations on systems of the same size as the substrates used in experimental observations (hundreds or even thousands of atoms in the molecular skeleton are nowadays common) and to provide a faithful description of the environment where the chromophore is probed. One approach for describing the surrounding environment is given by continuum models,7 whose basic ansatz is to describe the environment effect on the substrate through an effective contribution to the Hamiltonian operator derived from the * Corresponding author. E-mail:
[email protected].
electrostatic interaction between the chromophore and the surrounding environment. Among continuum models, the polarizable continuum model8,9 (PCM) in the integral equation formalism10 (IEF) has been widely employed for the calculation of molecular properties. Its main advantage compared to other continuum methods is the possibility to employ a moleculeshaped cavity instead of spheres or ellipses: this is an important issue when large systems with a possibly nonuniform shape are considered because, in these cases, simpler shapes would not yield a faithful description of the electrostatic interactions. The IEF-PCM formulation of linear11 and quadratic response12 theory has previously been presented and applied.13-17 For large molecular systems, an efficient parallel formulation has also been developed.18 In this work, the formalism for the IEF-PCM contribution to the cubic response function is derived at the selfconsistent field (SCF) level of theory. The approach is thus applicable both to Hartree-Fock (HF) and density functional theory (DFT), because the formulation of the solvent terms in the Hamiltonian does not differ for these two models. Therefore, once the cubic response formalism is available at the DFT level,19 the inclusion of solvent effects is identical to the HF case. We apply the present formalism to the theoretical determination of the DFWM process. The measured quantity is the third-order optical susceptibility χ(3)(-ω; ω,-ω,ω), which for a DFWM process corresponds to the microscopic hyperpolarizability γ(-ω; ω,-ω,ω). The development of methods to deal with fourth-order timedependent molecular properties, including environment effects, has already been proposed by other authors. Cammi et al.20,21 have implemented the coupled-perturbed Hartree-Fock (CPHF) formalism for a molecular solute that was subsequently applied by Mennucci et al.22 to obtain the first and second hyperpolarizabilities of urea in water; Champagne et al.23 employed the same CPHF approach on a number of polyacetylene chains including also local field effects. The first implementation of cubic response theory with solvent effects described as a dielectric continuum employing a spherical cavity model was presented by Sylvester-Hvid et al.,24 who derived the formalism both at the HF and at the multiconfigurational self-consistent
10.1021/jp0721191 CCC: $37.00 © 2007 American Chemical Society Published on Web 07/12/2007
8966 J. Phys. Chem. B, Vol. 111, No. 30, 2007
Ferrighi et al.
field (MCSCF) levels. Other authors have recently derived the response equations for a continuum solvation model at the coupled cluster (CC) level25 up to the fourth order (cubic response), although, to the best of our knowledge, applications have only been presented using quadratic response theory,26 whereas the second hyperpolarizability was obtained using finite difference calculations.27 Another possibility for including environment effects in the calculation is provided by quantum mechanics/molecular mechanics methods (QM/MM).28 This has recently been done by Kongsted et al. at the CC level for the wave function,29 linear,30 and nonlinear response theory.26,27 Although such an approach is in principle superior to the use of a continuum solvation model, it is also significantly more expensive, because the CC property calculations have to be performed on a (potentially large) set of configurations obtained from molecular simulations. This prevents at present its use on bigger molecules such as those studied in the present work. Jensen et al.31-33 have developed a similar approach by making use of DFT instead. They have developed the underlying equations up to and including quadratic response properties and have calculated fourth-order properties by means of finite differences.33 The remainder of the paper is organized as follows: In section II, we will present the equations for the IEF-PCM contribution to the CR function, in section III, we will give the details of our calculations, in section IV, we will discuss the results obtained, and in section V, we will give some concluding remarks.
Yˆ is the interaction of the nuclei with the apparent charges induced by the electrons, Xˆ (0) is the interaction of the electrons with the electronic apparent charges, and the last term, UNN, is the nuclear solvation energy. In the second quantization formalism, the terms can be written as36
II. Theory
In the case of a time-dependent perturbation, the expression for the free energy has a different form to take into account the dynamics of the solvent.37,38 The electronic term Xˆ is then divided into a so-called fast component Xˆ d(0), always in equilibrium with the solute, and a slow component Xˆ i(Fi) related to the initial density Fi. The free energy expression then becomes
A. Equilibrium and Nonequilibrium Solvation. The interaction energy between a charge distribution F(r), placed in a cavity C in a dielectric medium, and the polarization of the dielectric itself is given by34
Eint )
∫ΓV(s)σ(s) ds
(1)
F(r) dr ΓC|r - s|
(2)
with
V(s) )
∫
and the polarization being represented by the apparent charge density σ(s) defined on the surface ΓC of the cavity. Because σ(s) is dependent on the total potential, V(r) + Vσ(r), it is not possible to determine the interaction analytically. This problem is solved in the IEF-PCM10 approach by discretizing the cavity surface ΓC into a set of tesserae and representing the polarization σ(s) by a set of charges q. At the representative point sτ of each tessera τ, the charge qτ ) σ(sτ)aτ is defined, where aτ is the area of the tessera. The interaction energy can now be approximated as a discrete sum
Eint ) V‚q )
Vτ[D()]-1 ∑τ Vτqτ ) ∑ τσ Vσ ) τσ V‚(D()-1V) (3)
where the integral operator appearing in eq 2 has become a matrix D(), which only depends on the geometry of the cavity and the appropriate dielectric constant(s) of the medium. The solute model Hamiltonian can be written as35,36
H ˆ )H ˆ 0 + Jˆ + Yˆ + Xˆ (0) + UNN
(4)
where H ˆ 0 is the unperturbed Hamiltonian, Jˆ represents the interaction of the electrons with the nuclear apparent charges,
Jˆ ) V ˆ e‚qN )
Vτ,pqqNτ Eˆ pq ) ∑JpqEˆ pq ∑τ Vˆ eτ qNτ ) ∑ τ,pq pq
(5)
Yˆ ) VN‚qˆ e )
e VNτ qτ,pq Eˆ pq ) ∑YpqEˆ pq ∑τ VNτ qˆ eτ ) ∑ τ,pq pq
(6)
Xˆ (0) ) V ˆ e‚〈0|qˆ e|0〉 )
∑ Vτ,pqqτ,rsEˆ pq〈0|Eˆ rs|0〉 )
τ,pqrs
)
∑XpqrsEˆ pq〈0|Eˆ rs|0〉
(7)
pqrs
where Vˆ eτ is the operator that yields the electrostatic potential due to the electronic density at the representative point of tessera τ and qˆ eτ ) [D()-1V]τ is the corresponding charge operator. The nuclear term can be written as UNN ) VN‚qN ) ∑τVNτ qNτ . Because Jˆ and Yˆ are formally equivalent and UNN does not contribute to the wave function, the free-energy functional of the solute can be written as
〈|
|〉
1 ˆ 0 + Jˆ + Xˆ (0) 0 G ) 〈0|G ˆ |0〉 ) 0 H 2
(8)
〈|
|〉
1 G ) 〈0|G ˆ |0〉 ) 0 H ˆ 0 + Jˆ + Xˆ i(Fi) + Xˆ d(0) 0 2
(9)
B. The Cubic Response Equations. Using the notation of ref 6, we will consider the time evolution operator eκ(t),
κ(t) )
∑µ κµq†µ + κ′µqµ
(10)
where q†µ is the orbital excitation operator and κµ is the corresponding time-dependent amplitudes. The orbital excitation q†µ and de-excitation qµ ) (q†µ)† operators are defined through the singlet (+) and triplet (-) elementary excitation operators as † † † † ajR ( ajβ aiβ q†µ ) Eˆ ij(() ) aiR
i>j
(11)
and the CR-SCF terms can then be derived within the PCM formalism following the scheme already presented for linear11 and quadratic12 response. Following the formalism of Olsen and Jørgensen,6 we define the free energy derivative of order n as [n+1] [n+1] [n+1] ) Ej,l + Vj,l ‚[qN + qi] + Gj,l 1,l2,...,ln 1,l2,...,ln 1,l2,...,ln n
[n-k+1] d[k] Vj,l ∑ ,l ,...,l ‚ql ,l ,...l k)0 k+1 k+2
n
1 2
k
(12)
[n+1] [n+1] where Ej,l is the gas-phase energy term, Vj,l is the 1,l2,...,ln 1,l2,...,ln vector collecting the electrostatic potentials at the tesserae, qN is the nuclear apparent surface charge (ASC) vector, qi is the
Solvent Contribution to the Cubic Response Function
J. Phys. Chem. B, Vol. 111, No. 30, 2007 8967
inertial ASC vector, and qld[k] is the dynamic ASC vector. 1,l2,...lk The subdivision of the electronic ASCs into a static and a dynamic component is a consequence of the nonequilibrium formalism that has been used; the corresponding equilibrium expression is obtained by assuming that the inertial charges are zero. For details about nonequilibrium solvation within the PCM, see refs 11,38. The PCM contribution to the cubic response function can then be derived expanding eq 12 for n ) 3 as N i d [4] [4] [3] d[1] G[4] + jlmp ) Ejlmp + Vjlmp‚(q + q + q ) + Vjlm‚qp d[2] [1] d[3] V[2] jl ‚qlmp + Vj ‚qlmp (13) d[3] We will here only address V[4] jlmp and qlmp explicitly, because the expressions for the remaining terms have already been given in previous works.11,12 The expression for the perturbation expansion of the potentials [n+1] and charges ql[n] appearing in eq 12 are for the Vj,l 1,l2,...,ln 1,l2,...,ln SCF state for n ) 3
1 † V[4] jlmp ) - 〈0|[Oj ,[Ol, [Om, [Op, V]]]]|0〉 6
(14)
1 d qd[3] lmp ) - 〈0|[Ol, [Om, [Op, q ]]]|0〉 6
(15)
Figure 1. Molecules of groups I and II. In the first row, triazine (1.a), benzoxazole (1.b), benzimidazole (1.c), and benzothiazole (1.d). In the second row, 2-phenyl-benzothiazole (2.a) and 2-pyrid-2-yl-benzothiazole (2.b). In the third row, 2-thiophen-2-yl-benzothiazole (2.c) in the two employed conformations (cis and trans labels refer to sulfur atoms).
The operators Oj are defined as the general transformations of the elementary orbital excitations T ) (q†, q)T of eq 11 by a unitary matrix X12
O ) TX
(16)
The exponent of the time evolution operator (see eq 10) can then be rewritten as
κ(t) )
∑j OojRj
(17)
where R ) X†β, and β is a column vector collecting the coefficients of the orbital operator. For a detailed explanation of the formalism, the reader is referred to refs 11 and 12. By making use of eq 17 and contracting eq 14 with three excitation vectors, we obtain e [4] VjlmpNlNmNp
1 ) - 〈0|[qj, V(3κ, 2κ, 1κ)]|0〉 6
(18)
The equation for the transformed charges qd[3] lmp contracted with three vectors and transformed into the elementary basis is given as e d[3] qlmp NlNmNp
1 ) - 〈0|qd(3κ, 2κ, 1κ)|0〉 6
(19)
Identical expressions are obtained for the de-excitation operator qj instead of q†j . III. Computational Details To test the implementation of the formalism presented in the previous section, we have carried out a series of calculations of the static and frequency-dependent second-order hyperpolarizability γ(-ω; ω,-ω,ω) for a set of heteroaromatic molecules. The molecules investigated in this work can be classified in three groups, which we will refer to as I, II and III. The molecular structures of the molecules of group I and II are shown in Figure 1, whereas the molecules of group III are shown in Figure 2.
Figure 2. Molecules of group III. Tri-s-triazine derivatives: triethenyltri-s-triazine (3.a), triaminoethenyl-tri-s-triazine (3.b), and triazido-tris-triazine (3.c).
Group I. Group I contains small molecules with one or two aromatic rings: triazine (1.a), benzoxazole (1.b), benzimidazole (1.c), and benzothiazole (1.d). Group II. Group II contains three benzothiazole derivatives linked with three different aromatic rings in position 2: 2-phenyl-benzothiazole (2.a), 2-pyrid-2-yl-benzothiazole (2.b), and 2-thiophen-2-yl-benzothiazole (2.c). Group III. Group III contains three tri-s-triazine derivatives, namely, triethenyl-tri-s-triazine (3.a), triaminoethenyl-tri-striazine (3.b), and triazido-tri-s-triazine (3.c). Experimental data are available for groups I and II for the DFWM process at a wavelength of 602 nm in tetrahydrofuran (THF) solution,39,40 whereas for the molecules of group III, only theoretical HF values in gas phase are available in the literature.41 For molecules in groups I and II, the DFWM process corresponding to γ(-ω; ω,-ω,ω) has been calculated using the experimental frequency39 (ω ) 0.076 a.u., corresponding to λ ) 602 nm). The results obtained are discussed in section IV.A. To the best of our knowledge, there are no experimental data available for the molecules of group III, and we have, therefore, performed a more detailed study of the dispersion of the second hyperpolarizability employing the frequencies ω )
8968 J. Phys. Chem. B, Vol. 111, No. 30, 2007 0.0, ω ) 0.043 a.u. (corresponding to λ ) 1064 nm), and ω ) 0.076 a.u. in a set of commonly used solvents. The results are discussed in section IV.B. The DFWM process can be affected by resonances close to ω and 2ω. For molecules in groups I and II, experimental evidence suggests that there is no one-photon or two-photon absorption for DFWM at the wavelength of 602 nm,39 whereas no such experimental information is available for the molecules in group III. Because the γ(-ω; ω,-ω,ω) obtained using response theory can only be used off resonance unless special care is taken,42-44 we have calculated a few roots of the excitation spectrum to avoid unreliable calculations of the frequency-dependent polarizabilities, thus ensuring that the frequencies employed do not fall too close to a resonance. When near-resonance effects occur, we do not report the results of the response calculations. The molecules in groups I and II have been optimized in the gas phase and THF solution at the HF and DFT/B3LYP levels of theory using the cc-pVDZ basis set, while the third-order polarizabilities were calculated at the HF, B3LYP, and CAMB3LYP levels of theory using the aug-cc-pVDZ basis set. The group III molecules have only been optimized in the gas phase at the HF and B3LYP levels using the 6-31G* basis set; the hyperpolarizabilities have been calculated in gas phase and solvent using HF and DFT (B3LYP and CAM-B3LYP functionals) with the cc-pVDZ basis set. Group I and III molecules are kept planar, whereas group II molecules have been optimized without constraints to investigate the rotation between the two rings. For (2.a), we have obtained a planar geometry both at the HF and B3LYP levels. For (2.b), we have found that the configuration displayed in Figure 1 is planar both for HF and for B3LYP. The opposite configuration (cis with respect to the nitrogen atoms) is found to display a dihedral angle (N-CC-N) of 25° for HF and 5° for B3LYP. However, these configurations lie significantly higher in energy (both for HF and B3LYP the energy difference is more than 5 kcal/mol), thus leading to a negligible population of this conformation. For (2.c), we have obtained two planar configurations: their energies are similar (the cis configurations with respect to the sulfur atoms are less stable then their trans counterparts by around 1.5 kcal/mol) and the transition state is low enough (around 5 kcal/mol with respect to the trans conformer) to assume an equilibrium distribution. We have therefore employed both conformers of (2.c) in the property calculations. For the solvent calculations, a cavity defined by a set of interlocking spheres has been employed. The spheres are centered on the heavy atoms and the following set of radii has been employed: R(C,CH) ) 2.28 Å, R(N,NH) ) 2.16 Å, R(O) ) 2.04 Å, and R(S) ) 2.22 Å for groups I and II, and R(N,C,CH,CH2) ) 2.28 Å and R(Nazido) ) 2.04 Å for group III. The choice of radii is based on our previous work12,16 and on the van der Waals radii.45 To keep a uniform cavity within each ring/fused ring, we have here chosen to employ the same value for all identical nuclei belonging to the same ring/fused-ring, independent of the nature of the atoms to which they are bonded. This choice, tested on (1.a), (1.b), and (3.a) with respect to variations of the radii, was confirmed to have a negligible impact on the hyperpolarizability. The dependency of polarizabilities and hyperpolarizabilities on the choice of the radii has also been studied by Cammi et al.,20 who have obtained a variation of γ smaller than 5% within the spanned interval of radii. The static and optical dielectric constants of the different solvents used are reported in Table 1.
Ferrighi et al. TABLE 1: Static and Optical Dielectric Constants for the Solvents Employed in the Present Work C6H12 CHCl3 THF H2O
stat
opt
2.02 4.90 7.58 78.39
2.02 2.09 1.97 1.78
The hyperpolarizability results are presented as the orientationally averaged scalar γ, which can be compared to the experimental measurements. The orientationally averaged quantity is obtained from the second hyperpolarizability tensor γijkl as
〈γ〉 )
1 15
(γiijj + γijji + γijij) ∑ i,j
i, j ) x, y, z (20)
Units of 10-35 esu are used for the hyperpolarizability. The conversion factor is 1 au ) 5.0367 × 1040 esu. All calculations have been performed using a locally modified version of the Dalton program.18,46 IV. Results and Discussion A. Groups I and II. In this section, we will report and discuss the results obtained for the molecules in groups I and II. For all molecules in these groups, three sets of calculations have been carried out: DFWM gas-phase calculations with gas-phase optimized geometries and DFWM solvent calculations both with gas phase and with solvent-optimized geometries. All results are reported in Tables 2 and 3 for the two groups, respectively. For group I, we have employed three different functionals for the DFT calculations: the commonly used B3LYP functional, and the recent Coulomb-Attenuated Method B3LYP functional (CAM-B3LYP) with two different parametrizations. Using the same notation as in ref 47, the two sets of parameters that have been employed for CAM-B3LYP are R ) 0.19, β ) 0.46, and µ ) 0.33 and R ) 0.19, β ) 0.81, and µ ) 0.30. For the sake of brevity, we will refer to these parametrizations as CAMB3LYP and CAM-B3LYP(2), respectively, in the rest of this paper. The first parametrization is the one suggested in ref 47, whereas the second one fulfills the requirement for the exact long-range Coulomb interaction (R + β ) 1) and has been successfully employed in the calculation of the second harmonic generation process.16 In Table 2, we report the calculated values for the molecules of group I. The direct solvent effect, defined by comparing gas phase and solvent results at the gas-phase geometry, appears to be significant in all cases, being around 35% for triazine and between 40 and 50% for the other molecules. In contrast, the so-called “indirect” solvent effect, defined comparing values calculated in solvent with the gas phase and solvent reoptimized geometries, respectively, is always less than 1%, indicating that geometry variations of these substrates upon solvation are very small, as expected for rigid ring systems, and do not significantly contribute to the overall solvent effect. The introduction of electron correlation through DFT leads to higher second hyperpolarizabilities with respect to the HF values, as previously observed in the gas phase in a number of studies.19,48,49 This increase is further amplified by the solvent, and we have previously noted this synergy between correlation and solvation effects also in the case of the first hyperpolarizability.12,16 For instance, for B3LYP or CAM-B3LYP, the enhancement of γ in the solvent compared to the gas-phase data is up to 10% larger than the enhancement observed with Hartree-Fock.
Solvent Contribution to the Cubic Response Function
J. Phys. Chem. B, Vol. 111, No. 30, 2007 8969
TABLE 2: Values of γ(-ω; ω,-ω, ω) for DFWM of Group I Molecules in Gas Phase and THF Solution at the Wavelength of 602 nma (1.a)
HF CAM-B3LYP(2) CAM-B3LYP B3LYP
(1.b)
(1.c)
(1.d)
gas phase
THF (sg)
THF (gpg)
gas phase
THF (sg)
THF (gpg)
gas phase
THF (sg)
THF (gpg)
gas phase
THF (sg)
THF (gpg)
0.31 0.46 0.50 0.61
0.42 0.62 0.68 0.84
0.42 0.62 0.67 0.83
1.05 1.22 1.36 1.72
1.48 1.75 1.96 2.48
1.47 1.74 1.95 2.47
1.34 1.50 1.68 2.17
1.91 2.16 2.43 3.11
1.89 2.15 2.42 3.09
1.34 1.61 1.80 2.29
1.95 2.38 2.69 3.47
1.94 2.37 2.68 3.45
a Geometries have been optimized in gas phase and THF with the cc-pVDZ basis. For the solvent calculations, both the results with gas-phase geometry (gpg) and the results with solvent optimized geometry (sg) are reported. The second hyperpolarizabilities have been calculated with the aug-cc-pVDZ basis set. The nonequilibrium formalism for the calculations in solvent has been employed. Values are reported in 10-35 esu.
TABLE 3: Values of γ(-ω; ω,-ω, ω) for DFWM of Group II Molecules in Gas Phase and THF Solution at the Wavelength of 602 nma (2.a)
HF CAM-B3LYP(2) CAM-B3LYP B3LYP
(2.b)
(2.c) cis
(2.c) trans
gas phase
THF (sg)
THF (gpg)
gas phase
THF (sg)
THF (gpg)
gas phase
THF (gpg)
gas phase
THF (gpg)
5.46 9.14 12.93
8.42 15.01 22.94 23.34
8.36 14.91 22.66 21.87
5.47 10.67 19.52
8.22 17.08 41.93
8.22 17.12 42.65
5.79 10.68 14.06 23.26
9.08 18.06 25.28 45.97
5.65 10.78 14.10 25.28
8.84 17.82 25.93 49.13
a The geometries have been optimized in the gas phase and THF with the cc-pVDZ basis. For (2.a) and (2.b) the THF results are reported with both the gas-phase geometry (gpg) and the solvent optimized geometry (sg). For (2.c) both the cis and the trans conformations (with respect to the sulfur atoms) have been used, employing the gas-phase geometry only. The second hyperpolarizabilities have been calculated with the aug-ccpVDZ basis set. The nonequilibrium formalism for the calculation in solvent has been employed. Values are reported in 10-35 esu.
TABLE 4: Excitation Energies in eV in the Gas Phase and THF for the First Three Roots of Symmetry A′ and A′′ for the Benzothiazole Derivatives (Group II)a HF
CAM-B3LYP(2)
A′ root
gas
A′′ THF
gas
A′ THF
gas
CAM-B3LYP
A′′ THF
A′
B3LYP A′′
A′
A′′
gas
THF
gas
THF
gas
THF
gas
THF
gas
THF
4.30 4.60 4.93
4.20 4.58 4.90
5.02 5.61 5.81
5.07 5.67 5.85
3.98 4.09 4.56
3.89 4.07 4.53
4.68 5.19 5.35
4.74 5.26 5.43
4.25 4.48 5.13
4.18 4.46 5.14
4.70 5.02 5.43
4.73 5.08 5.52
3.84 3.96 4.61
3.83 3.91 4.68
4.28 4.63 5.06
4.30 4.70 5.14
1 2 3
4.81 5.42 5.89
4.74 5.41 5.56
6.07 6.17 6.23
6.10 6.19 6.27
4.46 4.81 5.07
4.37 4.79 5.04
5.17 5.73 6.10
(2.a) 5.22 5.80 6.09
1 2 3
4.86 5.39 5.77
4.78 5.39 5.74
5.93 6.19 6.21
5.97 6.22 6.25
4.43 4.72 5.28
4.36 4.71 5.28
4.89 5.19 5.62
(2.b) 4.92 5.24 5.70
1 2 3
4.59 5.42 5.73
4.49 5.42 5.70
5.94 6.08 6.23
5.98 6.07 6.29
4.16 4.76 5.27
4.05 4.75 5.21
5.16 5.42 5.77
(2.c) cis 5.24 4.01 5.46 4.53 5.82 5.10
3.90 4.52 5.05
5.00 5.19 5.45
5.09 5.23 5.60
3.74 3.98 4.71
3.63 3.97 4.65
4.63 4.65 5.08
4.67 4.74 5.14
1 2 3
4.61 5.42 5.73
4.52 5.43 5.70
6.01 6.10 6.26
6.05 6.12 6.33
4.14 4.76 5.23
4.04 4.75 5.18
5.19 5.60 5.98
(2.c) trans 5.27 3.90 5.66 4.52 5.96 5.01
3.90 4.52 5.01
5.12 5.50 5.74
5.12 5.50 5.74
3.74 3.99 4.65
3.63 3.98 4.61
4.67 4.94 5.22
4.76 5.00 5.30
a For the calculations in THF, solvent optimized geometries were used for (2.a) and (2.b), whereas gas-phase geometries have been used for the two conformers of (2.c). Experimental DFWM frequency is ω ) 2.07 eV (2ω ) 4.14 eV).
The benzothiazole derivatives of group II were studied experimentally in ref 39, because benzothiazole had the highest value of γ among the group I molecules, as also confirmed by our calculations. We have taken into account the compounds obtained by bridging a benzene, a pyridine, or a thiophene ring to the five-membered ring of benzothiazole, as shown in Figure 1. It has been found experimentally that the second hyperpolarizability is strongly enhanced by the presence of the extended π system.39 The calculated values of γ(-ω; ω,-ω, ω) for group II in the gas phase and in solvent are reported in Table 3. For molecule (2.a) and (2.b), we have considered both the direct and the indirect solvent effects, as for group I, whereas for (2.c), the indirect solvent effect has not been reported because the
results are in agreement with the observations for the other molecules of this group, and we limit the discussion for this molecule to the conformational analysis. To discuss the results obtained, it is necessary to recall that DFWM is a multiphoton process that can be influenced by resonances at ωi = ω or ωi = 2ω, where ωi is an excitation energy from the ground state. The first excitation energies for the molecules of group II are reported in Table 4. Both symmetries are coupled to the ground state because x and y belong to the A′ symmetry class and z belongs to A′′. Nevertheless, the transitions to the A′ states, which are coupled to the ground states through operators lying in the molecular plane, have transition moments one or more orders of magnitude larger than transitions to the A′′ states.
8970 J. Phys. Chem. B, Vol. 111, No. 30, 2007
Ferrighi et al.
TABLE 5: Calculated Values of γ(0; 0,0,0) in Different Environments at the HF, CAM-B3LYP, and B3LYP Level of Theory for Triethenyl-tri-s-triazine (3.a), Triaminoethenyl-tri-s-triazine (3.b), and Triazido-tri-s-triazine (3.c)a (3.a) Gas phase C6H12 CHCl3 THF H2O
(3.b)
(3.c)
HF
CAM-B3LYP
B3LYP
HF
CAM-B3LYP
B3LYP
HF
CAM-B3LYP
B3LYP
1.66 2.58 3.55 3.91 4.57
2.71 4.32 6.18 6.85 8.22
3.56 5.93 8.65 9.63 11.70
3.64 5.97 8.63 9.59 11.60
4.63 7.80 11.91 13.49 16.92
5.52 9.72 15.04 17.09 21.62
1.00 1.54 2.12 2.32 2.73
1.25 2.12 3.12 3.49 4.25
1.69 2.90 4.29 4.80 5.87
a Geometries are optimized in the gas phase with the 6-31G* basis set. Second hyperpolarizabilities are calculated with the cc-pVDZ basis set. Values are reported in 10-35 esu.
The energy of the incoming radiation is 2.07 eV, corresponding to a wavelength of 602 nm. There are, therefore, no problems connected to the resonance at ωi = ω. The resonance at ωi = 2ω (4.14 eV) is also not a problem at the HF level, but it affects the DFT calculations significantly. Because the DFT calculations give lower excitation energies than HF, the DFT results for (2.a) and (2.b) yield excitation energies that are very close to the resonance condition. In particular, we have verified that when one or more excitation energies are below 2ω, unphysical hyperpolarizability values can be obtained. As a consequence, the set of results reported in this paper is not complete, because in several cases the divergence prevented a meaningful solution of the response equations, and even for the reported cases, where a solution has been achieved, the observed trends (correlation effect, solvation effect, etc.) are strongly influenced by the proximity of the resonance and a detailed discussion of these results is, therefore, not meaningful. We will, therefore, limit ourselves to a couple of general considerations. First of all, we note that the HF values are very similar for the three molecules and also that a very similar solvent effect is obtained with an increase of about 50-60% of γ(-ω; ω,-ω, ω), which is slightly larger than that observed for the molecules of group I. The DFT values are much larger than the HF values and the solvent effect is also enhanced. Although this had already been observed for the group I molecules, there is, for the group II molecules, a significant lack of consistency that is certainly due to the closeto-resonance condition. For instance, γ(-ω; ω,-ω, ω) of (2.b) with CAM-B3LYP is exceedingly large compared to (2.a) and (2.c). For the two different conformers of (2.c), we have obtained no significant differences in γ(-ω; ω,-ω, ω) at any level of theory. B. Group III. The octupolar triazines of group III differ from each other by the lateral substituent in the tri-s-triazine core, having, respectively, -CHdCH2 (3.a), -CHdCH-NH2 (3.b), and -N3 (3.c). It has been found that the triazine-based compounds give higher polarizabilities compared to the benzene counterparts,50,51 and they have in addition higher thermal stability. For the molecules of group III, we have performed both gas phase and solvent calculations, studying an extended set of solvents because there are no experimental data to compare our results with. For all three molecules, we have calculated the static second-order hyperpolarizability γ(0; 0,0,0) in gas phase, C6H12, CHCl3, THF, and H2O at the HF and DFT levels, using the B3LYP and CAM-B3LYP functionals. The results are shown in Table 5. For all calculations of the group III molecules, only the gas-phase optimized geometries have been employed, because the molecules of group III are rigid aromatic systems and we can, thus, safely assume negligible indirect solvation effects, as already observed for the molecules of groups I and II. The direct solvent effect increases significantly with increasing dielectric constant of the solvent. For triethenyl-tri-s-triazine at the HF level, the enhancement of the second hyperpolariz-
TABLE 6: Calculated Values of γ(-ω; ω,-ω, ω) at the HF, CAM-B3LYP, and B3LYP Levels of Theory at Different Frequencies and in Different Environments for Triethenyl-tri-s-triazine (3.a)a HF
CAM-B3LYP
B3LYP
gas gas gas frequency phase THF H2O phase THF H2O phase THF H2O 0.000 0.043 0.076
1.66 2.01 3.44
2.55 2.42 2.71 3.12 2.94 3.56 5.43 5.08
4.23 3.97 3.56 5.84 5.46 5.38
5.79 5.40 9.04 8.41
a Geometries are optimized in the gas phase with the 6-31G* basis set. Second hyperpolarizabilities are calculated with the cc-pVDZ basis set. Values are reported in 10-35 esu.
TABLE 7: Calculated Values of γ(-ω; ω,-ω, ω) Using HF, CAM-B3LYP, and B3LYP at Different Frequencies and Different Environments for Triaminoethenyl-tri-s-triazine (3.b)a HF gas freq. phase THF
CAM-B3LYP gas H2O phase THF
B3LYP
gas H2O phase THF
H 2O
0.000 3.64 6.05 5.76 4.63 7.20 6.42 5.52 8.96 8.00 0.043 4.62 7.90 7.52 7.06 12.33 11.10 10.83 19.93 18.16 0.076 9.11 17.51 16.70 a Geometries are optimized in the gas phase with the 6-31G* basis set. Second hyperpolarizabilities are calculated with the cc-pVDZ basis set. Values are reported in 10-35 esu.
ability is between 55% for C6H12 and 175% for water, which, respectively, are the solvents with the lowest and highest stat. At the DFT level, both the B3LYP and the CAM-B3LYP functionals give higher gas-phase hyperpolarizabilities γ(0; 0,0,0) and also higher enhancements upon solvation than observed for HF. The CAM-B3LYP results are, as expected, between the HF and B3LYP values: the CAM-B3LYP results give an enhancement between 59% for cyclohexane and 200% for water, whereas B3LYP yields an amplification of γ(0; 0,0,0) between 66 and 228%. Triaminoethenyl-tri-s-triazine (3.b), where one hydrogen of each side chain of (3.a) has been replaced by an amino group, displays a higher hyperpolarizability than (3.a). Moreover, the observed solvent effect relative to the gas-phase results is larger, leading to an increase of γ from two to four times with respect to the gas phase. Similar trends are found for the triazido-tris-triazine (3.c), even though the calculated γ(0;0,0,0) are generally smaller compared to (3.a) and (3.b), in agreement with the CPHF values of ref 41. For (3.a) and (3.b), we have also carried out a frequencydependent study for the DFWM process in gas phase, THF and H2O employing HF, CAM-B3LYP, and B3LYP. The results are shown in Table 6 for (3.a) and in Table 7 for (3.b). The frequencies ω ) 0.0, 0.043, and 0.076 a.u. have been employed. The static values (ω ) 0.0) reported in these tables are calculated using the nonequilibrium formalism to discuss the dispersion
Solvent Contribution to the Cubic Response Function
J. Phys. Chem. B, Vol. 111, No. 30, 2007 8971
consistently. As for the two other groups of molecules, we have reported only the nonresonant values of γ. The dispersion of the hyperpolarizability of (3.a) at the HF level is overall similar in all selected environments (gas phase, THF, water) leading to a 20% increase of γ for ω ) 0.043 a.u. and to a 105-110% increase for ω ) 0.076 a.u. as compared to the zero-frequency limit. At the DFT level, the near-resonance condition prevented us from obtaining values at the highest frequency (ω ) 0.076). By comparing the results at ω ) 0.043 a.u. with the zero-frequency limit, two observations can be made: (a) correlation increases the dispersion, leading to a gasphase enhancement of 30% for CAM-B3LYP and 50% for B3LYP; (b) solvation leads to a slightly higher dispersion than that observed in the gas phase, with an increase of the hyperpolarizability of about 38% for CAM-B3LYP, to be compared with 30% in the gas phase and 56% instead of 50% for B3LYP, thus showing again a “synergy” effect between solvation and correlation. Similar considerations can be drawn for (3.b). We will, therefore, focus on the main differences with respect to (3.a). First of all, we note that the dispersion is larger: the gas-phase values for ω ) 0.043 a.u. are for instance 1.25, 1.5, and 2.0 times bigger than the corresponding zero-frequency limits for HF, CAM-B3LYP, and B3LYP, respectively. Moreover, solvation leads to an enhancement of the dispersion also for HF, although the solvation effect on the dispersion is still further accentuated by the correlated methods. As for (3.a), the DFT results are available only for ω ) 0.043 a.u. The study of dispersion on (3.c) has not been considered because calculations of γ(0; 0,0,0) showed that among the three selected molecules, (3.c) displays the smallest values in addition to being less sensitive to correlation effects. We do not, therefore, expect significant variations to the observations made for (3.a) and (3.b) as far as the dispersion of the hyperpolarizability is concerned. C. Comparison to Experimental Data. The molecules of groups I and II have been selected for the present study because experimental measurements of second hyperpolarizabilities were available.39 Nevertheless, the comparison to experimental data requires some extra care. The results reported in ref 39 have been obtained by converting the macroscopic third-order susceptibility χ[3] into the microscopic averaged hyperpolarizability γ by making use of the equation [3]
χ
) f (Nsγs + Ncγc) 4
(21)
where f is the Lorentz local field factor (f ) (n2 + 2)/3, n being the refractive index), N is the number density, γ is the microscopic second hyperpolarizability, and the subscripts s and c refer, respectively, to the solvent and the chromophore. In ref 39, γc is obtained by fitting eq 21 to measurements performed at different concentrations. The γ values obtained in this way are in principle directly comparable to gas-phase calculations because the solvation effect is accounted for through the factor f. This prevents us from comparing our calculated values, which include solvation effects in a more accurate way, to the results reported by Zhao et al.39 To be able to compare our calculations to the experimental results, it is first necessary to formulate the problem in a rigorous thermodynamical context.52 Equation 21 can be rewritten as
χ′ )
χ[3]V ) γ˜ sns + γ˜ cnc NA
(22)
where NA is Avogadro’s number, V is the system volume, ns and nc are the moles of the solvent and the chromophore, respectively, and γ˜ x ) f4γx (x ) s,c). χ′ is an extensive property of the system, as shown by the linear dependence on nc and ns, and the hyperpolarizabilities γ˜ c and γ˜ s can, therefore, be defined as
γ˜ s )
( ) ∂χ′ ∂ns
, γ˜ c )
P,T,nc
( ) ∂χ′ ∂nc
(23)
P,T,ns
In principle, both γ˜ c and γ˜ s depend on all thermodynamic parameters of the system (pressure, temperature, and concentrations). To derive molecule-specific properties, one has to refer to a standard state that, for a solution, is commonly taken as the infinite dilution limit. By dividing eq 22 by the total number of moles ns + nc, we obtain
χ′′ )
χV ) γ˜ sxs + γ˜ cxc NA(ns + nc)
(24)
and it follows that γ˜ s ) χ′′ (xc ) 0). Equation 24 can be differentiated to obtain
dχ′′ ) γs dxs + γc dxc + xs dγs + xc dγc
(25)
By thermodynamical arguments,52 one can show that xs dγs + xc dγc ) 0, and by recalling that for a binary system dxc ) -dxs we may write
dχ′′ ) γc - γs dxc
(26)
The above equation is valid for all concentration ranges. Taking the infinite dilution limit (xc f 0) and converting the equation to a molarity scale, we obtain
dχ′′ dχ′′ dMc dχ′′ ) lim ) xcf0 dxc xcf0 dMc dxc dMc lim
|
Mc)0
Fs Ws
(27)
where Fs is the solvent density, Ws is the solvent molecular weight, and their ratio corresponds to limxcf0 dMc/dxc. dχ′′/ dMc|Mc)0 is the slope of the fitting obtained by the experimental data extrapolated to Mc ) 0. By recalling that
γ˜ s ) χ′′ (xc ) 0) )
Ws [3] χ (Mc ) 0) N AF s
(28)
the final expression for γ˜ c reads
γ˜ c )
[ |
1 dχ[3] NA dMc
Mc)0
-
]
Ws [3] χ (Mc ) 0) Fs
(29)
According to Cammi et al.,53 once γ˜ c has been obtained from experimental measurements through this procedure, it can be directly compared to the calculated values obtained by a calculation, including the reaction field and cavity field effects. In section II, we have described how the reaction field effect is introduced at the CR level, and the results have been discussed in section IV. The cavity field effect is due to the modification of the incoming fields as a result of the presence of the cavity. We have implemented these corrections for the cavity field effect following the work of Cappelli et al.,54 to which the reader is referred for a detailed description. Here we limit ourselves to noting that the inclusion of the cavity field effect is obtained by making use of a modified dipole operator µ˜ for the incoming
8972 J. Phys. Chem. B, Vol. 111, No. 30, 2007
Ferrighi et al.
TABLE 8: Second Hyperpolarizability without (γ) and with (γ˜ ) the Local Field Correction for Molecules of Group Ia HF
CAM-B3LYP(2)
CAM-B3LYP
B3LYP
HF
CAM-B3LYP(2)
(1.a) γ γ˜ exp. γ γ˜ exp.
0.42 0.72
1.91 3.32
0.62 1.05
2.16 3.69
CAM-B3LYP
B3LYP
1.96 3.31
2.48 4.19
2.69 4.49
3.47 5.77
(1.b) 0.68 1.14
0.84 1.41
1.48 2.57
1.75 2.96
1.27
4.9 ( 0.99
(1.c)
(1.d)
6.2 ( 0.93
2.43 4.16
3.11 5.32
1.95 3.38
2.38 4.00
8.1 ( 1.05
a
For each molecule, the experimental result is also reported. Experimental data are taken from ref 40 for triazine and ref 39 for the other molecules. All experimental data are multiplied by the Lorentz factor f4 (f ) (n2 + 2)/3). Values are reported in 10-35 esu.
TABLE 9: Second Hyperpolarizability without (γ) and with (γ˜ ) the Local Field Correction for Molecules of Group IIa HF
CAM-B3LYP(2)
CAM-B3LYP
B3LYP
(2.a) γ γ˜ exp.
8.42 11.97
15.01 20.41 13.3 ( 1.73
γ γ˜ exp.
8.22 11.59
17.08 22.95 14.6 ( 1.46
22.95 30.62
23.34 30.22
(2.b)
V. Summary
41.93 54.49
(2.c) cis γ γ˜
9.08 13.18
18.06 25.32
γ γ˜ exp.
8.84 12.92
17.82 25.01 16.7 ( 3.4
25.28 35.17
45.97 62.21
25.03 34.84
49.13 67.07
(2.c) trans
a For each molecule, the experimental result is also reported. Experimental data are taken from ref 39. All experimental data are multiplied by the Lorentz factor f4 (f ) (n2 + 2)/3). Values are reported in 10-35 esu.
fields.55 The appropriate response quantity corresponding to DFWM is thus obtained as
γ˜ ijkl ) 〈〈µi; µ˜ j, µ˜ k, µ˜ l〉〉ω,-ω,ω
agreement with experimental findings, whereas DFT leads to a significant overestimation of γ˜ . In view of our findings for the group I molecules, we are inclined to consider the agreement between experiment and the HF results to be fortuitous and that the experimental data are actually somewhat too low. Unfortunately, the problems with resonances in the DFT calculations make it difficult for us to explore this further.
(30)
where the operators corresponding to the three incoming waves bear the cavity field correction. In Tables 8 and 9, we report our calculated results with (γ˜ ) and without (γ) the cavity field effect compared to the experimental measurements multiplied by the factor f4. The refractive index n for THF is 1.407, thus f4 ) 3.10. We notice that all gas-phase values are significantly underestimated, and the inclusion of the solvent brings our results into better agreement with experiment. The triazine results obtained using CAM-B3LYP and B3LYP are in very good agreement with experiment, whereas the HF and CAM-B3LYP(2) underestimate the experimental data. For the other molecules of group I, a general underestimation is observed. Nevertheless, the theoretical values at the DFT level are all within three times the experimental error bars for (1.b) and (1.c), whereas for (1.d) B3LYP is just outside this confidence interval, and the other methods all underestimate the experimental observations. The γ˜ values for the group II molecules are collected in Table 9 together with available experimental data. As for the discussion of the theoretical results alone, we limit ourselves here to a few general considerations due to the proximity of the experimental frequencies to resonances in the DFT results. We note, in particular, that the HF results are in surprisingly good
In this paper, we have presented the theory and implementation of the solvent contribution to the cubic response function of a molecular solute. The equations are derived at the SCF level and are, therefore, valid both for HF and DFT. The solvent is described as a macroscopic dielectric surrounding the molecule through the IEF-PCM. The present implementation has been tested on three different classes of molecules. The first two sets have been chosen because of the availability of experimental results, whereas the third group has been chosen to show the applicability of the method to bigger substrates. The comparison to experimental data has been carried out by including also cavity field effect contributions through a formally correct implementation within the PCM, instead of making use of local field factors that are only strictly valid for spherical cavities. The cavity-field-corrected hyperpolarizability calculated in this way must then be compared with the corresponding experimentally derived infinite dilution susceptibility without the inclusion of the Lorentz field factor. The overall solvent effect has been found to be quite large and dependent on the substrate, but also on the chosen method: in general, the solvent effects are more pronounced at the DFT level than at the HF level. The comparisons of our results with experiments have shown that for the smallest substrates of the group I molecules, HF underestimates the properties, whereas DFT shows a rather good agreement. For group II, the calculated hyperpolarizabilities are significantly larger than what has been observed experimentally, and in this case HF, therefore, yields results in better agreement with experimental data. A possible explanation for this is the resonance conditions observed for the DFT calculations. The proximity of the energy of the incoming light to resonances in the molecules is also documented experimentally, although it does not seem to affect the experimental measurements.39 The observed behavior seems to show that the calculation of γ in the molecules of group II cannot be considered satisfactory, neither from a purely theoretical perspective nor with respect to a comparison with experimental data. Response theory in the present formulation is not able to overcome these problems. The inclusion of a finite lifetime approximation for the electronic excited states also in the calculation of cubic response functions, as previously implemented for linear11 and quadratic12 response
Solvent Contribution to the Cubic Response Function functions, appears highly desirable for the study of processes involving the second hyperpolarizability. The different functionals employed have been used to investigate the effect of a better description of the long-range Coulomb interactions. The second hyperpolarizability is clearly reduced along the series B3LYP/CAM-B3LYP/CAM-B3LYP(2), but a definitive assessment of the quality of the results compared to experimental findings cannot be given because all functionals yield results within the generous experimental error bars. Acknowledgment. This work has received support from the Norwegian Research Council through a Center of Excellence Grant (Grant No. 179568/V30), a Strategic University Program in Quantum Chemistry (Grant No. 154011/420), a YFF grant to KR (Grant No. 162746/V00), and a Ph.D. grant to L.F. through the Nanomat program (Grant No. 158538/431), as well as through a grant of computer time from the Norwegian Supercomputing Program. References and Notes (1) Maker, P. D.; Terhune, R. W.; Savage, C. M. Phys. ReV. Lett. 1964, 12, 507. (2) Hellwarth, R. W. Prog. Quantum Electron. 1977, 5, 1. (3) Butcher, P. N.; Cotter, D. The Elements of Nonlinear Optics; Cambridge University Press: Cambridge, U.K., 1990. (4) Jin, Z. H.; Li, Z. Y.; Kasatani, K.; Okamoto, H. Chin. Phys. Lett. 2006, 23, 2786. (5) Dawes, A. M. C.; Illing, L.; Clark, S. M.; Gauthier, D. J. Science 2005, 308, 672. (6) Olsen, J.; Jørgensen, P. In Modern Electronic Structure Theory; Yarkony, D. R., Ed.; World Scientific: New York, 1995; p 857. (7) Tomasi, J.; Persico, M. Chem. ReV. 1994, 94, 2027. (8) Miertus, S.; Scrocco, E.; Tomasi, J. Chem. Phys. 1981, 55, 117. (9) Tomasi, J.; Mennucci, B.; Cammi, R. Chem. ReV. 2005, 105, 2999. (10) Cances, E.; Mennucci, B. J. Math. Chem. 1998, 23, 309. (11) Cammi, R.; Frediani, L.; Mennucci, B.; Ruud, K. J. Chem. Phys. 2003, 119, 5818. (12) Frediani, L.; Ågren, H.; Ferrighi, L.; Ruud, K. J. Chem. Phys. 2005, 123, 144117. (13) Pecul, M.; Lamparska, E.; Cappelli, C.; Frediani, L.; Ruud, K. J. Phys. Chem. A 2006, 110, 2807. (14) Ruud, K.; Mennucci, B.; Cammi, R.; Frediani, L. J. Comp. Meth. Sci. Eng. 2004, 4, 381. (15) Ruud, K.; Frediani, L.; Cammi, R.; Mennucci, B. Int. J. Mol. Sci. 2003, 4, 119. (16) Ferrighi, L.; Frediani, L.; Cappelli, C.; Salek, P.; Ågren, H.; Helgaker, T.; Ruud, K. Chem. Phys. Lett. 2006, 425, 267. (17) Frediani, L.; Rinkevicius, Z.; Ågren, H. J. Chem. Phys. 2005, 122, 244104. (18) Ferrighi, L.; Frediani, L.; Fossgaard, E.; Ruud, K. J. Chem. Phys. 2006, 125, 154112. (19) Jansik, B.; Salek, P.; Jonsson, D.; Vahtras, O.; Ågren, H. J. Chem. Phys. 2005, 122, 054107. (20) Cammi, R.; Cossi, M.; Tomasi, J. J. Chem. Phys. 1996, 104, 4611. (21) Cammi, R.; Cossi, M.; Mennucci, B.; Tomasi, J. J. Chem. Phys. 1996, 105, 10556. (22) Mennucci, B.; Cammi, R.; Cossi, M.; Tomasi, J. J. Mol. Struct.: THEOCHEM 1998, 426, 191.
J. Phys. Chem. B, Vol. 111, No. 30, 2007 8973 (23) Champagne, B.; Mennucci, B.; Cossi, M.; Cammi, R.; Tomasi, J. Chem. Phys. 1998, 238, 153. (24) Sylvester-Hvid, K. O.; Mikkelsen, K. V.; Jonsson, D.; Norman, P.; Ågren, H. J. Phys. Chem. A 1999, 103, 8375. (25) Christiansen, O.; Mikkelsen, K. V. J. Chem. Phys. 1999, 110, 8348. (26) Kongsted, J.; Osted, A.; Mikkelsen, K. V.; Christiansen, O. J. Chem. Phys. 2003, 119, 10519. (27) Kongsted, J.; Osted, A.; Mikkelsen, K. V.; Christiansen, O. J. Chem. Phys. 2004, 120, 3787. (28) Warshel, A.; Levitt, M. J. Mol. Biol. 1976, 103, 227. (29) Kongsted, J.; Osted, A.; Mikkelsen, K. V.; Christiansen, O. Mol. Phys. 2002, 100, 1813. (30) Kongsted, J.; Osted, A.; Mikkelsen, K. V.; Christiansen, O. J. Chem. Phys. 2003, 118, 1620. (31) Jensen, L.; van Duijnen, P. T.; Snijders, J. G. J. Chem. Phys. 2003, 118, 514. (32) Jensen, L.; van Duijnen, P. T.; Snijders, J. G. J. Chem. Phys. 2003, 119, 3800. (33) Jensen, L.; van Duijnen, P. T.; Snijders, J. G. J. Chem. Phys. 2003, 119, 12998. (34) Jackson, J. D. Classical Electrodynamics; Wiley Science: New York, 1975. (35) Cammi, R.; Tomasi, J. J. Comput. Chem. 1995, 16, 1449. (36) Cammi, R.; Frediani, L.; Mennucci, B.; Tomasi, J.; Ruud, K.; Mikkelsen, K. V. J. Chem. Phys. 2002, 117, 13. (37) Mennucci, B.; Cammi, R.; Tomasi, J. J. Chem. Phys. 1998, 109, 2798. (38) Cammi, R.; Tomasi, J. Int. J. Quantum Chem. 1995, 56, 465. (39) Zhao, M. T.; Samoc, M.; Prasad, P. N.; Reinhardt, B. A.; Unroe, M. R.; Prazak, M.; Evers, R. C.; Kane, J. J.; Jariwala, C.; Sinsky, M. Chem. Mater. 1990, 2, 670. (40) Keshari, V.; Wijekoon, W. M. K. P.; Prasad, P. N.; Karna, S. P. J. Phys. Chem. 1995, 99, 9045. (41) Zheng, W. X.; Wong, N. B.; Li, W. K.; Tian, A. M. J. Chem. Theory Comput. 2006, 2, 808. (42) Norman, P.; Ruud, K.; Helgaker, T. J. Chem. Phys. 2004, 120, 5027. (43) Krykunov, M.; Kundrat, M. D.; Autschbach, J. J. Chem. Phys. 2006, 125, 194110. (44) Norman, P.; Bishop, D. M.; Jensen, H. J. A.; Oddershede, J. J. Chem. Phys. 2005, 123, 194103. (45) Bondi, A. J. Phys. Chem. 1964, 68, 441. (46) Dalton, an ab initio electronic structure program, Release 2.0; http:// www.kjemi.uio.no/software/dalton/dalton.html, 2005. (47) Peach, M. J. G.; Helgaker, T.; Salek, P.; Keal, T. W.; Lutnaes, O. B.; Tozer, D. J.; Handy, N. C. Phys. Chem. Chem. Phys. 2005, 8, 558. (48) Rizzo, A.; Cappelli, C.; Jansik, B.; Jonsson, D.; Salek, P.; Coriani, S.; Wilson, D. J. D.; Helgaker, T.; Ågren, H. J. Chem. Phys. 2005, 122, 234314. (49) Rizzo, A.; Cappelli, C.; Jansik, B.; Jonsson, D.; Salek, P.; Coriani, S.; Ågren, H. J. Chem. Phys. 2004, 121, 8814. (50) Ray, P. C.; Das, P. K. Chem. Phys. Lett. 1995, 244, 153. (51) Zhu, W. H.; Wu, G. S. J. Phys. Chem. A 2002, 106, 7216. (52) Denbigh, K. Principles of chemical equilibrium: With applications in chemistry and chemical engineering; Cambridge University Press: Cambridge, U.K., 1964. (53) Cammi, R.; Mennucci, B.; Tomasi, J. J. Phys. Chem. A 2000, 104, 4690. (54) Cappelli, C.; Corni, S.; Mennucci, B.; Cammi, R.; Tomasi, J. J. Phys. Chem. A 2002, 106, 12331. (55) Cammi, R.; Mennucci, B.; Tomasi, J. J. Phys. Chem. A 1998, 102, 870.
