8050
J. Phys. Chem. 1993,97, 8050-8053
Frozen Density Functional Approach for ab Initio Calculations of Solvated Molecules Tomasz Adam Wesolowskit and Arieh Warshel' Department of Chemistry, University of Southern California. bs Angeles, California 90089- 1062 Received: March 25, 1993; In Final Form: May 13, 1993
A new density functional method for treatment of chemical processes in solution is presented. The method is based on freezing the electron density of the solvent molecules, while solving the eigenvalue problem for the solute Hamiltonian, which includes the effective potential of the solvent molecules. The method is developed and examined in the simple case of one solvent and one solute molecule. The results are encouraging and the deviation between the unfrozen and frozen treatments can be attributed to the kinetic energy functional used. The method can be implemented in ab initio calculations of solvation free energies, following a recent pseudopotential approach [Vaidehi et al., 19921. I. Introduction Quantum mechanical calculations of chemical processes in solutions and proteins present a major challenge due to the enormoussize of the reacting systems.' Thus the common strategy in such calculations involves simplified models that describe quantum mechanically only a small part of the system, e.g. the solute, while treating the surrounding (e.g. the solvent) classically.2" Obtaining a more quantum mechanical description of the solvent is obviously far from simple and requires at present significant simplifications. Our attempt to progress in this direction resulted recently in the development of an approach that presented the solvent molecules by a pseudopotential, which was incorporated into the solute Hamiltonian at the ab initio level.' This approach, however, has not provided a rigorous prescription for the derivation of the relevant pseudopotentials, thus making it difficult to assess the approximations involved. A promising alternative, at last from the formal point of view, can be provided by the density functional theory (DFT). For example, one may try to use the recently developed method of Yang,* which divides the system into subsystems and solves the eigenvalue problem for each of them. However, at present this approach is too expensive to be used for calculations of solvation free energies or related problems. A more effective treatment of solvation problems might be developed by following the implicit assumption of our pseudopotential treatment and keeping the electron density of the solvent constant. Such an approximation reflects the fact that the density of the solvent electrons is approximately known and the solution of the eigenvalueproblem for the entire system is not essential. In this case, one can restrict the explicit DFT treatment to the solute molecule while representing the contribution of the solvent molecules as an effective potential. Such an approach can provide a useful way of improving pseudopotential treatments since the DFT formulation provides a clear picture for the meaning of different interaction terms. Furthermore, while parameters in pseudopotential methods depend on the type of the interacting atoms, a DFT-based approach might lead to more general interaction terms which do not require reparametrization for different systems. In view of the above discussion we consider in this work a new DFT approach which is based on the use of a frozen electron density for the surrounding of the reacting system. This approach is examined in a preliminary way by evaluating the interaction between a Li+ ion and a single water molecule (the generalization to many water molecules is straightforward). The calculations
provide encouraging indications that the frozen-DFT method can be used effectivelyin ab initio calculations of chemical processes in solutions.
II. Theoretical Method (a) Background. Our task is to start from a formal DFT treatment of the entire solute-solvent system and then to freeze the solvent density (possibly with induction corrections), while reproducing the energetics of the unfrozen system. Trying to approximate the properties of a complete system with frozen densities of subsystems presents a major challenge. However, because we are dealing here with noncovalent interactions, we expect such an approximation to give reasonable results. Before we start the development of the frozen density approach it is useful to consider the regular DFT approach for N electrons under the effect of an external potential V(r). For simplicity, we restrict all subsequent considerations to closed-shell systems. According to the Hohenberg-Kohn theorem9 there exist functionals of the total electron density ( p ) representing the total energy (E) and kinetic energy (T):
T = Tbl The ground-state energy corresponds to the minimum of E[p] with respect to all electron densities, obeying the constraint j p ( r ) dr
=N
(3)
The Euler equation for this problem islo (4) where p is a Lagrange multiplier associated with constraint (3) and 6F[p]/6p(r) denotes the variational derivative of functional F defined by
Kohn and Sham11 proposed a method for solving eq 4 using molecular orbitals (ql, ...,cpn). In this method, refereed to here as DFT/KS, the total electron density is represented by
* To whom correspondence should be addressed. t Permanent address: Department of Biophysics, ul. Zwirki i Wigury 93, 02-089 Warszawa, Poland. 0022-3654/93/2097-8050$04.00/00 1993 American Chemical Society
(6)
The Journal of Physical Chemistry, Vol. 97, No.30, 1993 8051
Frozen Density Functional Approach
and the effective potential can be expressed as
are defined as The functionals Ts[p]and Exc[p] N 2
Ts[PI = -2
vj1/2v21vj ) j= 1
Exhl = F b l
-
JJm P(r)P(r') dr dr'
TJp]
and the Euler equation (eq 4) now becomes
The corresponding Euler equation takes the form:
6Ts[Pl
P = Veff(d + Here Verr(r)is an effective potential defined as
Veff= ~ ( r + )
J&
d t + Vexc(r) The solution of eq 19 can be obtained analogous to the original KohnSham method by solving
where Vexcis the exchange-correlation potential given by
+
[-1/2V2 V#efr(r)]qi= tivi
Within the framework of the KohnSham method the molecular orbitals, vi, are obtained as a solution of the following equation:
+
[-1/2V2 Veff(r)]q = tipi
i = 1, 2,
...,N/2
= p 2 ( r ) + pl(r)
(13)
If pZ(r) is taken as a given function, which obeys the constraint J p 2 ( r ) dr = N - N'
where the new effective potential WE&), is given by
(12)
Equation 6 provides a link between the solution of eq 12 and the effective potential V&). Usually eqs 6 and 10-12 are solved in an iterative way. The next section presents a simplified scheme that freezes a part of the total electron density during the self-consistent calculations. (b) The Frozen Density Treatment. Let us divide our system into two subsystems S1 and Sz with N' and N - N' electrons respectively, so that the total density can be expressed as p(r)
i = 1,2, ...,N'/2 (20)
(14)
thenpl is the densitythat minimizes the energy of the total system
Equations 15,20, and 21 form a new scheme for self-consistent calculations, analogous to the DFT/KS method, which is referred to here as FDFT/KS, where "F" denotes the fact that part of the electronic density (e.g. solvent) is "frozen". Instead of solving the N-electron problem it is now possible to solve an "-electron problem (N'< N),with the modified effective potential of eq 21. It is important to point out that our approach for treating the kinetic energy is very similar to the approach introduced in a recent work of Cortona. (c) The Kinetic Energy Functional. In the FDFT/KS scheme the number of explicitly treated electronsin the system is reduced fromNto N'. Thus, insteadof calculatingN X Nmatrixelements of the form
V ( k l )= (VkIV.&)lQ) (22) one has to calculate N' X N' matrix elements of the form
fl(W) = ('pkIV#eff(r)IQ)= (VkIVefi(r)IQ) +
HPl. Following the idea of the KohnSham method, we propose to express P I as "12
Now we can define the effective potential, V&), in the same way as in eqs 9 and 10, since this potential depends on the total electron density p(r) and does not depend on the way this density is represented. However, thekineticenergy can not be represented in the same way as in eq 7 as long as the wave function that corresponds to theentire system (solvent plus solute) is not known. Thus for p2 we use thedensity dependent kinetic energy functional, T[p2],and for PI a wave function dependent functional, Ts[pl], while considering the remaining nonadditive term Tsnadd.That is, we express the total kinetic energy as T,[PI= T , [ P+ ~ p21 = T , [ P ~+I T,[P,I+ T ? ~ ~ [p21 P~,
(16) With this definition we can write the variational derivative 6Ts[Pl/6Pl(r)as
(23)
(q{-b)
Obviously, the computational efficiency of the FDFT/KS method is superior to that of the DFT/KS method, only if the matrixelements (q ~ G T d " d d [ p l , p 2 ] / 6 p l ( r ) ( ( harenot ) much harder to calculate than the original (qklVe!dr)le). A very simple analytical form of Tudd[pl,pz] can be obtained by applying the original Thomas-Fermi theory, in which kinetic energy depends on local values of the electron density TTF[p]= cTFJPs'3 dr
(24)
where CTF= 2.871, In this case the nonadditive kinetic energy functional is given by TTFnadd[pl, P21 = CTFJ{(P1+ P2)'I3
-
- P25/'} dr (25)
The corresponding variational derivative is given by
Another analytical form of T,,&[Pl,pz] can be obtained when the kinetic energy functional contains terms that also depend on
Wesolowski and Warshel
8052 The Journal of Physical Chemistry, Vol. 97, No. 30, 1993
3
"I$ d
dft
+
fdft-tf
fdft-add
0
ldft-tfw
-
c
dn fdft(H2O) fdft(Li+)
in
-165"'
'
16
'
'
7-1
'
'
'
26
'
31
'
'
36
41
Figure 1. Dependance of the total energy of the Li+-H20 complex on the distance between the Li+ ion and the 0 atom. The figure compares energiescalculatedwith the regular DFT/KS method toenergiesobtained by the FDFT/KS method using frozen electron density for the Li+ ion. The FDFT/KS calculationsrepresent the nonadditivekinetic energy term by the following functionals: additive (0), Thomas-Fermi (TF), ThomasFermi with gradient correction (TFW), and scaled Thomas-Fermi
(TF-1.4).
the electron density gradient.12
+
&s%f P
dr
The nonadditive kinetic energy functional is now given by
where
while the corresponding variational derivative is given by 6TTFwnaddbl,
6P I
P21 = 6TTFnaddb1, P21 6Pl
'
1
1.6
2.1
2.6
3.1
3.6
4.1
R(O-Li+)[A]
R(O-Li+)[A]
TTFw[p] = TTF[p]
-50
+
6TWmddb1,
P21
6P 1
where
III. Results In order to examine the performance of the FDFT/KS method in studiesof intermolecularinteractionswe consider the interaction between an H20 molecule and a Li+ ion located at the bisector of the H-O-H angle. The calculations were performed with Whitten's basis set.14 The kinetic and exchange-correlation functionals, together with their variational derivatives, were obtained by means of numerical integration over Cartesian grids, centered on atomic nuclei. The spacing of the grid decreased discontinuously, from 0.5 to 0.05 A, with increasing proximity to the corresponding atomic nucleus. The exchange-correlation functional proposed by Vosko, Wilk, and Nursairls was used in all calculations. Figure 1 compares DFT/KS energies with those obtained by FDFT/KS with different analytical forms of the functionalsused for Tnadd[p1, ~ 2 1 . All FDFT/KS calculations correspond to the use of frozen density for the Li+ ion. As can be easily seen from the figure the neglect of Tmdd[Pl, p2] leads to totally wrong results for short distances. Nonadditive treatment using a ThomasFermi functionalgives reasonable results, while a more elaborated functional that contains both Thomas-Fermi terms and Weizsacker corrections does not lead to significant improvement of results. For short distances, energies are even worse. The
Figure 2. Dependance of the total energy of the Li+-H20 complex on the Li+-O distance for calculationsthat use the regular DFT/KS method and the FDFT/KS method with a scaled Thomas-Fermi kinetic energy functional. The FDFT/KScalculationsinvolve the frozen electron density for (a) Li+ and (b) H20.
agreement between the DFT/KS and FDFT/KS energies can be improved significantly using scaled Thomas-Fermi functionals for Tnadd[P1, p ~ ] .The best results are obtained with a scaling factor of 1.4, where the difference between the DFT/KS and FDFT/KS interaction energies is only 0.002 kcal/mol at the minima of the corresponding interaction curves. The difference between FDFT/KS and DFT/KS results can be attributed to either introduction of the frozen density constraint or to the accuracy of the nonadditive kinetic energy functional (Tna[pl, ~21).In this respect it is important to point out that in the DFT treatment (in contrast to pseudopotential treatment) one does not have to worry about orthogonality of the wave functions of the solute-solvent system, provided the "correct" nonlocal functional is used. In the case studied in Figure 1 one can attribute all differences between the FDFT/KS and DFT/ KS methods to the functional used. That is, in this case we represented the electronic configuration of Li+ ion by a single ls-type orbital. Consequently,the electron density of the Li+ ion is constant in both the DFT/KS and FDFT/KS calculations. The situation is, however, quite different in calculationsthat freeze the density of the water molecule, where both inductive effects and quality of the functional might be important. Fortunately, the same functional optimized in frozen Li+ calculations can be used in frozen H2O calculations, allowingone to explore the effect of the frozen density constraint. This point is illustrated in Figure 2, which compares the FDFT/KS interaction energy of the Li+.-H20 system for two sets of calculations: one which uses frozen electron density for Li+ and another which uses frozen density for the water molecule. Both calculations use a nonadditive kinetic energy term based on a scaled Thomas-Fermi functional, with the same 1.4 scaling factor (obtained in frozen Li+calculations). As seen from Figure 2, frozen H20 calculations have higher interaction energiesthan both frozen Li+ calculations and DFT/KS calculations. In this case the difference between the FDFT/KS calculations that use frozen electron density for the water molecule and DFT/KS calculations can be attributed to changes in the electron density of the water molecule upon approaching a Li+ ion (these changes are, of course, neglected in the FDFT/KS calculations). As seen in Figure 3 the above difference is proportional to R-4 when the distance between water oxygen and the Li+ exceeds 2 A. One can easily improve the accuracy of the FDFT/KS calculationsfor this range by including electronic polarization in the model. It is important to point out that the scaled Thomas-Fermi functional is far from being the perfect functional. While interaction energies are quite accurate, the correspondingelectron densities are not accurate (see also ref 16). Several strategies could be used in the future to reduce this problem, but considering the impressive results obtained for the energetics, we believe that the present treatment can provide a starting point for calculations of solvation energies.
Frozen Density Functional Approach
The Journal of Physical Chemistry, Vol. 97, No. 30, 1993 8053
-- I
-12
u 5
I O
15
2.0
2.5
h(R) Figure3. Difference (diff) between the DFT/KS and FDFT/KS energies with electron density for HzO,at Li+-0 distances of more than 2 A. The correlation between In ( 4 i f f ) and In R has a linear cocffcient close to 4, indicating that the difference is due to the neglect of inductive effects.
IV. Discussion The present paper considers a new quantum mechanical approach for studies of chemical processes in solution. This approach is based on the use of a density functional treatment for the solute-solvent system, while freezing the electron density of the solvent molecules. The first and most crucial step in developingthe FDFT/KS approach is implemented and examined in the present work by considering a solute and one solvent molecule. It is found that FDFT/KS treatment gives reasonable results, even with simplekinetic energy functionals. It also appears that the FDFT/KS method provides a uniform formalism that can be applied to a wide range of distances between solute and solvent molecules. This formalism can also be easily extended to intermediate ranges by including electron polarization. The present treatment is not perfect and can probably be improved by using more sophisticated functionals. However, the use of density functional formulation provides a much clearer physical picture than our previous pseudopotential treatment. It is also very significantto note that we obtain a reasonable intermolecular potential surface, without use of adjustable parameters. The present approach might seem similar to the “divide and conquer” approach of Yang.8 However, the philosophy of the two approaches is quite different. Yang’s approach divides the system to subsystems and solves the eigenvalue problem for each of them. The FDFT/KS approach, on the other hand, follows the philosophy of our early hybrid quantum/classical method17 and freezesjhe electron density of the solvent molecules. This allows one to restrict explicit DFT treatment to the solute molecule while representing the effect of the solvent molecule as an effective potential (which is, however, rigorously defined in terms of the assumed electron density). While this idea (as any other simple
idea) might seem obvious,we are not aware of similar DFTstudies. In fact, our approach could be termed “freeze and conquer”. Some readers of this paper might not appreciate the computational advantagesof the FDFT method since this workconsiders only a simple solvent molecule. The true advantage of our approach is, of course, associated with the treatment of many solvent molecules. For example, a study of a system with 1000 solvent molecules only requires solution of the solute’s eigenvalue problem while that of Yang requires solution of an additional 1000 eigenvalue problems. The use of frozen density for solvent molecules is obviously a major approximation,even when inductive effects are taken into account. However, this approximation might be justified in studies of solvation problems since the interactions between the solute and the solvent molecules are much weaker than corresponding interactions between bonded fragments. Thus, one might obtain reasonable results even if the performance of the method at very short distances is not perfect. Furthermore, deviations at short distances can be improved by developingkinetic energy functionals that give a better descriptionof intermolecularinteractions. Such functionals may be obtained by fitting the KS kinetic energy term to a function of the corresponding overlapping densities of the interacting molecules. The implementation of the FDFT/KS method in actual calculationsof processes in solution will be reported in subsequent papers. These studies will evaluate solvation free energiesby the same free energy perturbation approach used in our pseudopotential study.7
Acknowledgment. This work was supported by NIH Grant GM24492 and by ONR Grant N00014-91-5-1318.
References and Notes (1) Warshel, A. ComputerModelingof ChemicalReactions In Enzymes AndSolutions; John Willey & Sons: New York, 1991. (2) Warshel, A. Curr. @in. S t m t . Biol. 1992, 2, 230. (3) Kolman, P. A.; Hayes, D. M. J. Am. Chem. Soc. 1981,103,2955. (4) Tapia, 0.;Anders, J.; Aullo, J. M.; Branden, C. I. J . Chem. Phys. 1985. 83. 4673. (S) ;an Duijnen, P.;Thole,Th.;Brcer, B.;Nieuwpoort,R.lnt. J. Quantum Chem. 1980, 17,651. (6) Bash, P. A.; Field, M. J.; Davenport, R. C.; Petsko, G. A.; Ringe, D.; Karplus, M. Biochemistry 1991,30, 5826. 17) Vaidehi.. N.:. Wesolowski.. T.:. Warshel. A. J. Chem. Phvs. 1992.97. , . 4264. ‘ (8) Yang, W. Phys. Reo. Lett. 1991, 66, 1438. (9) Hohenberg, P.; Kohn, W. Phys. R w . 1964, 136,9864. (IO) Parr, R. G.; Yang, W. Density-Functional Theory of Atom and Molecules; Oxford University Reas: New York, 1989; p 52. (11) Kohn, W.; Sham, L. J. Phys. Reu. 1965,110, A1133. 1121 Weizsaecker. C. F. 2.Phvs. 1935.96.431. (l3j Cortona, P. (1991) Phys.wReu.B 1991,448454. (14) Whitten, J. L. J. Chem. Phys. 1965, 44, 359. (1s) Vosko, S. H.; Wilk, L.; Nursair, M. Can. J . Phys. 1980.58, 1200. (16) Reference 10; p 135. (17) Warshel, A.; Levitt, M. J. Mol. Biol. 1976,103, 227. ~
