Quantum Calculation of Molecular Energies and Energy Gradients in Solution by a

Quantum Calculation of Molecular Energies and Energy Gradients in Solution by a

(Parte 1 de 3)

Quantum Calculation of Molecular Energies and Energy Gradients in Solution by a Conductor Solvent Model

Vincenzo Barone and Maurizio Cossi* Dipartimento di Chimica, UniVersita Federico I, Via Mezzocannone 4, I-80134 Napoli, Italy

ReceiVed: May 21, 1997; In Final Form: October 17, 1997

A new implementation of the conductor-like screening solvation model (COSMO) in the GAUSSIAN94 package is presented. It allows Hartree-Fock (HF), density functional (DF) and post-HF energy, and HF and DF gradient calculations: the cavities are modeled on the molecular shape, using recently optimized parameters, and both electrostatic and nonelectrostatic contributions to energies and gradients are considered. The calculated solvation energies for 19 neutral molecules in water are found in very good agreement with experimentaldata; the solvent-inducedgeometryrelaxationis studiedfor some closed and open shell molecules, at HF and DF levels. The computational times are very satisfying: the self-consistent energy evaluation needs a time 15-30% longer than the corresponding procedure in vacuo, whereas the calculation of energy gradients is only 25% longer than in vacuo for medium size molecules.

1. Introduction

Many chemical reactions of biological and technological relevance take place in a condensed medium, in particular in liquid solutions, and often an accurate theoretical treatment of such processes cannot leave a realistic description of the environmental effects aside. The recent advances in quantum calculationtechniques are making ab initio calculationspossible on systems of ever-increasingsize: thus there is a growing need of ab initio procedures providing reliable treatments of solutesolvent interactions.

Among the several approaches proposed to describe the solvent effect at the ab initio level, continuum models are quite popular,1-4 due to their flexibility and efficiency. In such models the solute molecule, possibly supplemented by some solvent moleculesbelongingto the first solvationshell, is placed in a cavity surrounded by a polarizable continuum, whose reaction field modifies the energy and the properties of the solute. In the most advanced ab initio models (e.g., PCM,5,6 SCIPCM,7 SCRF,8,9 COSMO,10,1 and GCOSMO12) the cavities are of molecular shape, and the reaction field is described in terms of apparentpolarizationcharges or reactionfield factors included in the solute Hamiltonian, so that it is possible to perform iterative procedures leading to the self-consistence between the solute wave function and the solvent polarization. Quantum calculations considering the solvent reaction field are not limited to monodeterminant wave functions: for example, Mikkelsen et al.13,14 and Karlstrom15 have proposed multiconfigurational procedures where the solvent is described by a continuum dielectric.

In addition, some models also describe nonelectrostatic solute-solvent interactions (which are usually referred to as cavitation, dispersion, and repulsion energies). Moreover, the geometry relaxation induced by the solvent on the solute molecules cannot often be neglected; thus an efficient solvation model must provide energy gradients and allow geometry optimizations in solution. In other words, it is desirable that both direct (i.e., polarization) and indirect (i.e., relaxation) solvent effects are treated with the same accuracy.

In this work we report the novel implementation in the package GAUSSIAN9416 of an ab initio solvation model based on the conductor-likesolvation model (COSMO), first proposed by Klamt and Schuurmann for classical calculations10 and then extended to quantum mechanical systems.1,12 The present implementation allows Hartree-Fock (HF), density functional (DF), and post-HFenergycalculationsand HF and DF geometry optimizations in solution; moreover, all the manipulations of the molecular wave function provided by GAUSSIAN94 for isolated systems are available for solvated molecules also.

The COSMO approach describes the solvent reaction field by means of apparent polarization charges distributed on the cavity surface, which are determined by imposing that the total electrostaticpotential cancels out on the surface. This boundary condition, suited for cavities in conductor media, can describe (i) the interaction between molecules and metals (e.g., in the simulationof electrodicprocesses)and (i) the solvation in polar liquids.

For the latter applications, the conductor-like model is physically less founded than dielectric models; nevertheless, the conductor approach is attractive, since its boundary condition is computationally simpler, especially in the expression of the energy gradients. Some authors have pointed out that the conductor model well reproduces the solute energies and properties obtained with the dielectric approach, using the dielectric constants characteristic of polar solvents.12,17 Our results confirm that, using cavity parameters optimized for the well-known polarizable continuum model (PCM), the COSMO proceduregives hydrationenergies in very good agreementwith the experimental results.

The present implementationalso providesthe nonelectrostatic contributions to the solute free energy and the first derivatives of these contributions with respect to the nuclear coordinates: some examples of geometry optimizations are reported, using both the electrostatic and the nonelectrostatic terms. To our knowledge, these are the first examples of complete geometry optimizations for a conductor model. In fact the derivatives of the nonelectrostatic solute-solvent energies require an accurate description of the changes induced on the cavity shape by the nuclear motion, which is provided by a recently developed procedure.18 On the other hand, the electrostatic energy gradients also are improved by considering these geometrical contributions properly.

S1089-5639(97)01699-X C: $15.0 © 1998 American Chemical Society Published on Web 02/20/1998

2. Theoretical Background

2.1. Definition of the Cavity. The solute molecules are embedded in cavities formed by interlocking spheres centered on the solute atoms or atomic groups. The surface is smoothed by addingsome otherspheres,not centeredon atoms,to simulate the so-called solvent-excluding surface, following the GEPOL procedure.19-21 Then, the cavity surfaceis partitionedinto small domains, called tesserae, obtained by projecting on the surface the faces of polyhedra inscribed in each sphere; the tesserae completely buried into other spheres are discarded, and those partially cut are replaced by suitable polygons.18

The present implementation exploits the recently developed

PolyGen procedure,2,23 allowing a wide choice of polyhedra, so that the surface can be covered with finer and finer tessellations, if needed. A long experience with PCM calculations indicates that a tessellation using 60 tesserae per sphere (obtained by projecting the faces of inscribed pentakisdodecahedra) is a reasonable compromise between accuracy and efficiency; we shall show that this tessellation can be used in standard COSMO calculations, too. Anyway, it is a good rule to verify the invariance of the results with respect to this parameter, at least in test cases. 2.2. Molecular Free Energy. Under the influence of the solvent, the molecular Hamiltonian is perturbed:

where Hö 0is the Hamiltonian of the isolated solute; the operator Vö , describing the electrostatic solute-solvent interactions, linearly depends on the solute wave function j¾〉: in this case, it has been shown24 that the SCF procedure leads to the variational minimization of the free energy of the solute, G:

The interaction operator Vö is written in terms of apparent polarization charges: in each tessera i a charge qi appears, according to the conductor-like boundary condition:

where V and Vq are the electrostatic potential due to the solute and to the polarization charges, respectively, and rb is a point on the surface. The vector of the conductor-like polarization charges, Q, can be determined by the equation where the vector V contains the electrostatic potential due to the solute on the tesserae. The elements of the matrix A are

where Si is the area of tessera i; the expression for Aii has been found by Klamt and Schuurmann in the case of a single sphere

partitioned into a variable number of identical tesserae10 and also used by Truong and Stefanovich in their own GCOSMO implementation.12

If the COSMO model is used to simulate a solvent with dielectric constant , the polarization charges have to be scaled so that the total polarization charge will obey the Gauss law: a very effective way to do this is multiplying each charge by the factor ( -1)/ , so that the actual charges are

In analogy with the boundary element method formulation of PCM,24 it is convenient to separate the potential due to the solute nuclei, VN, and that due to the electrons, Ve, defining two sets of charges:

In a finite basis matrix formulation, the molecular Hamiltonian in the presence of the solvent can be written where is the Hamiltonian for the isolated molecule, and j, y, X, and

UNN represent the interactions between qe and the solute nuclei, qN and the solute electrons, qe and the electrons, and qN and the nuclei,respectively. The explicitexpressionsfor these terms have been presented elsewhere for the PCM model;6 they can be used in this framework without changes. The corresponding Fock matrix is

In these terms it is possibleto perform an iterativecalculation, formally equal to the usual SCF procedure for isolated molecules, leading to self-consistentpolarization charges and solute wave function. Notice that both Hartree-Fock and density functional, as well as hybrid HF-DF, Hamiltonians can be modified according to eq 10. Moreover, the present implementation allows us to perform perturbative many-body (MPn), configuration interaction, and coupled cluster calculations using the polarization charges determined at the HF level or iterating the calculation to get a complete self-consistence between the polarization charges and the post-Hartree-Fock wave function, following the same procedure illustrated for the PCM method, e.g., in ref 25.

The total polarization charges appearing on the cavity surface are subject to Gauss’ law:

where Zn and Ne are the atomic number of nucleus n and the tesserae

Aij ) 1 jrbi - rbjj (6) qiN ) qGaussN )- nuclei Zn (14)

∑i qie ) qGausse )

1996 J. Phys. Chem. A, Vol. 102, No. 1, 1998 Barone and Cossi number of solute electrons, respectively. In practice, the conditions 14 and 15 are not fulfilled, both for numerical and for physical reasons: the numerical error arises from the approximated description of the polarization charge density in terms of discrete point charges; the physical error is due to the little amount of solute charge escaping from the cavity because of the exponential decay of the electronic tails. The numerical error is present for both qN and qe, approximately at the same extent, whereas the escaped charge error, usually heavier than the former one, affects qe only. Some methods have been proposed to correct such effects, both in dielectric26 and in conductor models;27 however, it is known that the conductor approach is less sensitive to these sources of error than the dielectricone, and actuallywe foundthat the calculatedCOSMO charges for neutral solutes sensibly satisfy eqs 14 and 15 (see below). Then we deferred the treatment of numerical and electronic errors to a further work, devoted to ionic solutes.

According to eq 2 the molecular free energy in solution can be written

where the subscript “es” recalls that we are considering electrostaticinteractions only: for brevity, Vi indicates the value of the electrostatic potential on the tessera i.

Recalling eqs 4 and 7, the solute potential at tessera i can be related to the set of polarization charges

and ¢Ges can be put in an equivalent form,

that will be useful in the formulation of energy gradients.

Taking into account the nonelectrostatic interactions, the solute free energy becomes

The free energy associated with the formation of the cavity in the continuummedium, Gcav, is calculatedwith the expression derived by Pierotti from the hard sphere theory,28 adapted to the case of nonsphericalcavities.25 The dispersionand repulsion terms, Gdis, Grep, are calculated following Floris and Tomasi’s procedure,29,30 with the parameters proposed by Caillet and

Claverie;31 the implementation details are described in ref 18. Clearly the nonelectrostatic terms are exactly the same in the PCM and the COSMO models, provided equal cavities are employed. 2.3. Free Energy Derivatives. The free energy derivatives can be written

where the superscript R indicates the partial derivative with respect to the nuclear coordinate R.

The expressions for the derivatives of the nonelectrostatic terms have been given elsewhere;18,32 as said in the Introduction, they heavily depend on the derivatives of the tesserae geometrical elements (area, position of the vertices), which are neglected, totally or in part, in the preceding COSMO implementations.10,1,3 The derivative of the electrostatic contribution is

[〈¾jHö 0j¾〉]R is the usual energy derivative calculated by standard programs; from eq 19 we have

Recalling eq 18, the terms depending on the derivatives of the polarization charges cancel out, leaving

This expression for ¢GesR is computationally very efficient, since it doesn’t require the calculation of charge derivatives:

notice that in dielectric models an analogous expression does not exist.

The free energy derivatives (eq 21) can be used in standard optimizationprocedures(e.g., the Berny algorithmimplemented in GAUSSIAN94, used throughout this work), in order to calculate the solvent-induced molecular relaxation.

3. Computational Details

The calculationswere performedat the HF, MP2, and density functional levels; for the latter we resorted to the hybrid B3LYP functional, which combines HF and Becke34 exchange terms with the Lee-Yang-Parr correlation functional;35 the 6-31G- (d) and 6-311+G(d,p) basis sets were used.

The solvent was water at 25 °C, with dielectric constant ) 78.4.

(Parte 1 de 3)