## Jane S Murray Peter Politzer

The quest for improved methods for elucidating and predicting the reactive behavior of molecules and other chemical species is a continuing theme of theoretical chemistry. This has led to the introduction of a variety of indices of reactivity; some are rather arbitrary, while others are more or less directly related to real physical properties. They have been designed and are used to provide some quantitative measure of the chemical activities of various sites and/or regions of the molecule.

In this chapter our focus is on one of these indices, the electrostatic potential V(r) that is created in the space around a molecule by its nuclei and electrons. V(r) can be computed rigorously, given the electronic density function p(r), by Eq. (3.1).

ZA is the charge on nucleus A, located at RA. V(r) is the potential created by the static charge distribution of the molecule; it can also be regarded as the exact interaction energy of this charge distribution with a proton situated at the point r. Of course such a proton would, in reality, produce some polarization of the molecule's electronic density, which is not taken into account if p(r) corresponds to the unperturbed ground state, as is normally the case. Despite this inconsistency, the pioneering work of Scrocco, Tomasi, and their collaborators (Bonaccorsi, Scrocco, and Tomasi 1970; Bonaccorsi, Scrocco, and Tomasi 1971) demonstrated the usefulness of the electrostatic potential as a guide to molecular interactive behavior, a role in which it is now well-established (Berthier et al. 1972; Berthod and Pullman 1975; Berthod and Pullman 1978; Bonaccorsi et al. 1972a; Bonaccorsi et al. 1972b; Bonaccorsi et al. 1975; Giessner-Prettre and Pullman 1975; Lavery, Corbin and Pullman 1982; Lavery, Pullman, and Pullman 1980; Lavery and Pullman 1981; Perahia and Pullman 1978; Politzer and Daiker 1981; Politzer, Laurence, and Jayasuriya 1985; Politzer and Murray 1990; Politzer and Murray 1991; Politzer and Truhlar 1981; Pullman and Berthod 1976; Pullman and Pullman 1980; Pullman and Pullman 1981b; Scrocco and Tomasi 1973, Scrocco and Tomasi 1978). Until a few years ago, its extensive applications have focused upon either (a) the extrema of V(r), its most negative and (more recently) positive values, or (b) the qualitative pattern of the electrostatic potential's positive and negative regions, plotted either on planes through the molecule or on its surface. More recently, a third approach has been pursued, which involves quantifying certain key features of the overall pattern of the electrostatic potential and relating them to macroscopic properties (Murray et al. 1994; Murray and Politzer 1994). [For a current exposition of a variety of applications of molecular electrostatic potentials, see Murray and Sen (Murray and Sen 1996).]

Unlike many of the other quantities used now and earlier as indices of reactivity, the electrostatic potential V(r) is a real physical property, one that can be determined experimentally by diffraction methods as well as computationally (Politzer and Truhlar 1981). It is through this potential that a molecule is first "seen" or "felt" by other chemical species. Equation (3.1) shows that the potential at any point r is the sum of a positive term, representing the contribution of the nuclei, and a negative one, reflecting the distribution of electrons. The sign of V(r) is determined by the term that dominates at the point in question. For a neutral atom, the electrostatic potential is positive everywhere (Politzer and Murray 1991; Sen and Politzer 1989): negative regions arise only in the space surrounding a molecule or an anion.

Our discussion in this chapter will focus on the use of the electrostatic potential as a means to understanding and predicting chemical interactions. First, we will examine some of its properties and important features. Next, we will discuss methodology. Finally we will review some recent applications of the electrostatic potential in areas such as hydrogen bonding, molecular recognition, and understanding and prediction of a variety of physio-chemical properties related to molecular interactions. Our intent has not been to provide a complete survey of the ways in which the potential has been used, many of which are described elsewhere (Politzer and Daiker 1981; Politzer, Laurence, and Jayasuriya 1985; Politzer and Murray 1990; Politzer and Murray 1991; Politzer and Truhlar 1981; Scrocco and Tomasi 1973), but rather to focus on some diverse examples.

### 3.1 Background

The expression for V(r), given as Eq. (3.1), follows from the definition of electrical potential, which will be reviewed here. Any distribution of electrical charge creates a potential V(r) in the surrounding space. For an assembly of point charges Q. located at positions r;, this electrical potential is simply a sum of Coulombic potentials, as given in Eq. (3.2).

Qt may be either positive or negative in sign. If the charge distribution is continuous, then integration replaces the summation in Eq. (3.2), giving Eq. (3.3).

D(r) is the total charge density; its sign can vary as a function of r.

If we consider a molecule as having a static but continuous distribution of electronic charge around a rigid nuclear framework, then its electrical or "electrostatic" potential will have a term similar to Eq. (3.2), with Q. being the positive charges of the nuclei, ZA, and a term similar to Eq. (3.3), with D(r) being replaced by the electronic density function p(r). Since p(r) is customarily defined as a positive function [unlike Z)(r)], the second term on the right side of Eq. (3.1) comes in with a negative sign. The net result is Eq. (3.1). Our purpose in reviewing this background is to show explicitly that Eq. (3.1) follows from basic electrostatics.

It should be noted that the electrostatic potential V(r) is related rigorously to the total charge density D(r) through Poisson's equation, Eq. (3.4).

V2 V(r) = -4T7D(r) = 4irp(r) - 4ir£zA8(r-RA) (3.4)

V(r) accordingly plays a key role in the very fundamental and rapidly developing area of density functional theory (Parr 1983; Parr and Yang 1989). An important aspect of this has been the development of relationships, both exact and approximate, between the energies of atoms and molecules and the electrostatic potentials at their nuclei (Levy, Clement, and Tal 1981; Politzer 1980; Politzer 1981; Politzer 1987). More recently, we have shown that the position of the minimum potential along a bond provides us with a realistic set of co-valent radii (Wiener et al. 1996) and that its magnitude is related to the bond dissociation energy (Wiener et al. 1997). The significance of these findings is discussed elsewhere (Politzer and Murray 1996).

The sign of the electrostatic potential V(r) in any particular region around a molecule, which depends upon whether the effects of the nuclei or electrons are dominant, is a key to assessing its reactivity there. Regions in which V(r) is negative in sign are those to which electrophiles are initially attracted, in particular to the most negative potentials, Vmin. These Vmin are typically between 1 and 2 A from the nuclear framework and associated with one of the following three molecular features: (a) heteroatoms with lone pairs, such as O, N, F, S, P, CI, Se, As, and Br; (b) pi regions, such as are found in aromatic and double- and triple-bonded systems; and (c) strained bonds (Politzer and Murray 1990).

Sites susceptible to nucleophilic attack can also be identified and ranked by means of positive electrostatic potential regions, but it is necessary to analyze the latter at distances at least 1 to 2 A away from the nuclei, e.g., in planes removed from the molecular framework (Murray, Lane, and Politzer 1990; Politzer, Abrahmsen, and Sjoberg 1984; Politzer et al. 1984) or on molecular surfaces (Murray et al. 1991b; Murray and Politzer 1991; Murray and Politzer 1992; Pullman, Perahia, and Cauchy 1979; Sjoberg and Politzer 1990). This is because the electrostatic potentials of atoms and molecules have local maxima only at the nuclei (Pathak and Gadre 1990). To identify sites for nucleophilic attack, it is necessary, accordingly, to look for the most positive values in planes or on surfaces that are at some distance from the nuclei. (These are, of course, not true local minima.)

As mentioned earlier, the electrostatic potential around a free neutral atom is positive everywhere (Politzer and Murray 1991; Sen and Politzer 1989), due to the very highly concentrated positive charge of the nucleus in contrast to the dispersed negative charges of the electrons. It is when atoms interact to form molecules that regions of negative potential may and usually do develop as a consequence of the subtle electronic rearrangements that accompany the process.

Figures 3.1 and 3.2 show calculated electrostatic potentials for guanine (1), in the plane of the molecule and on the molecular surface, respectively. Looking at Figures 3.1 and 3.2, it can be seen that there are negative potentials associated with N3, N? and the carbonyl oxygen, with the latter two overlapping to form a strong and extensive negative region on one side of the molecule. Both Figs. 3.1 and 3.2 allow us to rank N? as the site most suscepti-

ble to electrophilic interactions, with N3 and the carbonyl oxygen being fairly similar but less negative than N?. This is consistent with the experimental observation that N7 is the favored site for the protonation and alkylation of guanine (Fiskin and Beer 1965; Lawley 1957), with some alkylation also occurring at the oxygen (Friedman, Mahapatra, and Stevenson 1963), which is more accessible than N3 (Figs. 3.1 and 3.2). Focusing next on positive regions of potential, it is clear that the surface shown in Fig. 3.2 is superior to the contour map in Fig. 3.1 in revealing the relative magnitudes of the positive potentials associated with the hydrogens. Specifically, Fig. 3.2 shows that the hydrogens of the amine group and the one bonded to Nx are more positive than those on the five-membered ring; the former would accordingly be predicted to be more favored for hydrogen bonding. This is indeed found to be the case; guanine hydrogen bonds to negative sites on cytosine through an amine hydrogen and that onNr

Uses of the electrostatic potential that will be emphasized in this chapter will be both qualitative and quantitative. It is important to recognize however that these applications are based on the following exact interpretations of V(r):

1. Given a point charge ± Q located at the point r, then ± QV(r) is equal to the electrostatic interaction energy between the unpolarized molecule and the point charge.

2. In a perturbation theory treatment of the total (not just electrostatic) interaction between the molecule and the point charge, ±QV(r) is the first-order term in the expression for the total interaction energy (which would include polarization and other effects).

3. The negative gradient of ±QV(r) equals the electrostatic force that is exerted by the molecule's unperturbed charge distribution upon the point charge ±Q.

As just mentioned, ± QV{r) is an energy quantity. Even though V(r) itself is a potential, not an energy, it is customary to express V(r) in units of energy (e.g., kcal/mole). This is actually QV(r) with Q equal to +1.

Since the electrostatic potential is closely related to the electronic density, it may be useful to discuss how the information that can be obtained from V(r) differs from that provided by the p(r). Both are real physical properties, related by Eqs. (3.1) and (3.4). An important difference between V(r) and p(r) is that the electrostatic potential explicitly reflects the net effect of all of the nuclei and electrons at each point in space, whereas the electron density directly represents only the concentration of electrons at each point. A molecule's interactions with another chemical system is affected by its total charge distribution, both positive and negative, and thus can be better understood in terms of its electrostatic potential than its electronic density alone. Examples illustrating this point have been discussed elsewhere (Politzer and Daiker 1981; Politzer and Murray 1991).

### 3.2 Methodology

Although Eq. (3.1) is an exact formula for the electrostatic potential due to a set of nuclei (ZA) and an electronic density p(r), the latter [p(r)] is generally obtained computationally from an ab initio or semi-empirical wave function or, more recently, from density functional procedures and is therefore necessarily approximate. It follows that the resulting V(r) is also approximate.

Within this framework, given any particular p(r), V(r) can be evaluated rigorously (all of the integrals in Eq. (3.1) being calculated exactly) or approximately (Politzer and Daiker 1981). With the ready availability of software packages like the Gaussian series, the former is certainly the more widely used procedure. We will discuss some aspects of this first. We will then briefly examine some approximate methods for obtaining V(r) from p(r); detailed analyses of these procedures are found elsewhere (Politzer and Daiker 1981; Politzer, Laurence, and Jayasuriya 1985; Tomasi 1982).

3.2.1 Dependence of Rigorously Evaluated V(r) Upon Computational Level

How does a rigorously calculated electrostatic potential depend upon the computational level at which was obtained p(r)? Most ab initio calculations of V(r) for reasonably sized molecules are based on self-consistent field (SCF) or near Hartree-Fock wavefunctions and therefore do not reflect electron correlation in the computation of p(r). It is true that the availability of supercomputers and high-powered work stations has made post-Hartree-Fock calculations of V(r) (which include electron correlation) a realistic possibility even for molecules with 5 to 10 first-row atoms; however, there is reason to believe that such computational levels are usually not necessary and not warranted. The M0ller-Plesset theorem states that properties computed from Hartree-Fock wave functions using one-electron operators, as is V(r), are correct through first order (M0ller and Plesset 1934): any errors are no more than second-order effects.

It has been shown that the electrostatic potentials of formamide calculated at near-Hartree-Fock (HF/6-31G*) and post-Hartree-Fock (MP2/6-31G*) levels are qualitatively similar (Politzer and Murray 1991). Both computational approaches predict the oxygen to be the preferred site for electrophilic attack (Seminario, Murray, and Politzer 1991). It is further noteworthy that SCF results obtained with minimal basis sets (e.g., HF/STO-3G and HF/STO-5G) are also in good agreement with those calculated at the higher computational levels.

We feel justified in restating our earlier conclusion that varying the ab initio computational level does not greatly affect the overall pattern of the electrostatic potential for a given molecule (Politzer and Daiker 1981; Politzer and Murray 1991). There are, however, certainly differences in detail, the exact locations and magnitudes of the potential minima may change to some degree, and not always predictably (Luque, Illas, and Orozco 1990; Politzer and Daiker 1981; Politzer and Murray 1991; Seminario, Murray, and Politzer 1991). The key point is that a generally reliable picture of the electrostatic potential can be obtained with an SCF wavefunction, even if only of minimum basis set quality (Boyd and Wang 1989; Daudel et al. 1978; Gatti, MacDougall, and Bader 1988; Luque, Illas, and Orozco 1990; Politzer and Daiker 1981; Politzer and Murray 1991; Seminario, Murray, and Politzer 1991). However, we have found that the inclusion of polarization functions for molecules with second-row atoms is recommended, even at the minimum basis set level. There are in dications that the semi-empirical MNDO and AM 1 methods also yield qualitatively reliable electrostatic potentials (Ferenczy, Reynolds, and Richards 1990; Luque, Illas, and Orozco 1990; Luque and Orozco 1990).

Within the past ten years, density functional procedures (Dahl and Avery 1984; Labanowski and Andzelm 1991; Parr and Yang 1989; Seminario and Politzer 1995) have emerged as an extremely promising alternative to the more traditional ab initio and semi-empirical procedures for computing molecular properties. Density functional theory is based on the Hohenberg-Kohn theorem (Hohenberg and Kohn 1964), according to which all of the electronic properties of a chemical system, including the energy, are determined by the electronic density. Important features of this approach are that it takes account of electron correlation but nevertheless requires considerably less computer time and space than do comparable ab initio techniques. The effectiveness of density functional procedures for computing molecular electrostatic potentials and other molecular properties is still being explored (Labanowski and Andzelm 1991; Laidig 1994; Murray et al. 1992; Seminario and Politzer 1995; Sola et al. 1996). Electrostatic potentials obtained by a local density functional method were shown to be similar to those from SCF calculations (Murray et al. 1992). More recently, analyses of density distributions obtained by density functional techniques (Laidig 1994; Sola et al. 1996) suggest that they provide generally reliable distributions of charge. These results are encouraging in regard to the use of density functional methods for obtaining electrostatic properties.

3.2.2 Approximate Evaluation of V(r)

With the continuing surge of development in computing capabilities, approximate methods for evaluating V(r) are now generally used only for very large molecular systems, such as those studied in nucleic-acid, protein, and other biomolecular research. Historically, the most widely used approximate procedures for computing the electrostatic potential have been those based upon multipole expansions (Etchebest, Lavery, and Pullman 1982; Politzer and Daiker 1981; Rabinowitz, Namboodiri, and Weinstein 1986; Tomasi 1982; Williams 1988; Williams 1991). Such representations can approach the rigorously computed V(r) to varying degrees, depending on the number of terms (i.e., quadrupole, octa-pole, etc.) that are included. Terminating the expansion after the monopole terms [which corresponds to using a set of point charges to obtain V(r)] is the simplest possibility, the results of which obviously depend on the number, locations, and magnitudes of the point charges (Politzer and Daiker 1981). Overall, this approach has had only limited success. (On the other hand, there continues to be considerable interest in using the molecular electrostatic potential as a basis for obtaining physically meaningful atomic charges (Besler, Merz, and Kollman 1990; Breneman and Wiberg 1990; Chirlian and Francl 1987; Francl et al. 1996; Williams 1991; Williams and Yan 1988; Woods et al. 1990).) Expansions through the quadrupole terms have been shown to yield V(r) comparable with that obtained rigorously from the same p(r) (Murray et al. 1990; Rabinowitz, Namboodiri, and Weinstein 1986). This success has stemmed from the recognition that an electronic density function written in terms of a Gaussian basis set can be expressed as a finite multicenter expansion (Rabinowitz, Namboodiri, and Weinstein 1986), with the centers not limited to the nuclei.

Another methodology for computing the electrostatic potential that has been of interest for a number of years involves representing V(r) of large systems as a combination or superposition of contributions from their constituent units or fragments (Bonaccorsi et al. 1980; Nagy, Angyán, and Náray-Szabó 1987; Náray-Szabó 1979; Politzer and Daiker 1981;

Pullman and Pullman 1981a; Scrocco andTomasi 1973;Tomasi 1981). Breneman, quite recently, has developed a method in which V(r) is computed from densities obtained through "transferable atom equivalents" (Breneman 1996); the resulting electrostatic potentials are of ab initio quality.

3.3 Some Applications

### 3.3.1 Analysis of Noncovalent Interactions

Noncovalent interactions, both inter- and intramolecular, are of considerable importance in determining the physical properties of molecules. Such interactions can be classified as hydrogen-bonding or non-hydrogen-bonding. In this section we will explore some recent uses of the electrostatic potential in the analysis of both types.

3.3.1.1 Family-Independent Relationships Between Computed Electrostatic Potentials on Molecular Surfaces and Solute Hydrogen Bond Acidity/Basicity

In view of the well-established importance of the electrostatic component in hydrogen bonding (Benzel and Dykstra 1983; Buckingham and Fowler 1985; Kollman 1977; Legon and Millen 1987; Lin and Dykstra 1986; Umeyama and Morokuma 1977), it is not surprising that the molecular electrostatic potential V(r) has been found to be an effective means for analyzing and correlating hydrogen-bonding interactions (Espinosa et al. 1996; Hagelin et al. 1995; Kollman et al. 1975; Leroy, Louterman-Leloup and Ruelle 1976c; Murray and Politzer 1991; Murray and Politzer 1992; Murray, Ranganathan, and Politzer 1991; Politzer and Daiker 1981; Politzer and Murray 1991). For example, it has been used successfully to predict the sites and directionality of hydrogen bonds in a variety of systems, including many hydrogen-bonded dimers (Kollman et al. 1975; Leroy, Louterman-Leloup, and Ruelle 1976a; Leroy, Louterman-Leloup, and Ruelle 1976b; Leroy, Louterman-Leloup, and Ruelle 1976c). Specifically, the positions of the most negative potentials, Vmin, associated with the hydrogen-bond accepting heteroatoms of isolated gas-phase molecules were shown to be effective for predicting the sites and therefore directionality of hydrogen bonds to that particular heteroatom (Kollman et al. 1975; Leroy, Louterman-Leloup, and Ruelle 1976a; Leroy, Louterman-Leloup, and Ruelle 1976b; Leroy, Louterman-Leloup, and Ruelle 1976c). In addition, a good correlation was found between calculated hydrogen bond energies and the value of V(v) at a fixed distance from the hydrogen-bond accepting molecule in a series of complexes between HF and various acceptors (Kollman et al. 1975).

In order to further explore the relationship between the magnitude of the V irl in the vicinity of a hydrogen-bond accepting heteroatom and its tendency to form a hydrogen bond, we studied the relationship between Vmin and the solvent hydrogen-bond-accepting parameter (3 (Murray, Ranganathan, and Politzer 1991). (3 is one of the "solvatochromic parameters" introduced by Kamlet et al. (Kamlet et al. 1983; Kamlet, Abboud, and Taft 1981; Kamlet et al. 1979; Kamlet, Solomonovici, and Taft 1979; Kamlet and Taft 1976) in the course of an extended effort to separate, identify, and quantify various types of solvent effects upon experimentally measurable solution properties (e.g., rate constants, equilibrium constants, and IR, NMR, ESR, and UV/vis absorption maxima and intensities). It is inter preted as providing a measure of a solvent's ability to accept a proton in a solute-to-solvent hydrogen bond (Kamlet et al. 1983). We showed that there exist good correlations between P and Vmm, computed at the HF/STO-5G*//HF/STO-3G* level, for four series of oxygen-or nitrogen-containing molecules (Murray, Ranganathan, and Politzer (1991): azines, primary amines, alkyl esters, and molecules containing double-bonded oxygens.

The success of this simple approach led us to consider correlating some value of the electrostatic potential associated with a hydrogen-bond-donating molecule with the solva-tochromic parameter a. The latter is viewed as indicative of a solvent's ability to donate a proton in a solute-solvent hydrogen bond (Kamlet et al. 1983). Because there are no true maxima associate with any regions of a molecule's V(r) away from its nuclei, we chose to compute the electrostatic potential on molecular surfaces defined by a contour of the electronic density (e.g., the 0.002 au or 0.001 au contour). Good correlations were found between a and the most positive value of the electrostatic potential on the surface, max, for a group of -OH and a group of alkyl-hydrogen-bond donors (Murray and Politzer 1992). These calculations were also carried out at the HF/STO-5G*//HF/STO-3G* level.

In the preceding studies, we had focused upon the hydrogen-bond basicity and acidity of solvents. Our next step was to investigate whether our calculated Vmin and Vs max would also correlate with Abraham et al.'s more recently developed scales of solute hydrogen-bond basicity and acidity (Abraham et al. 1988; Abraham et al. 1989a; Abraham et al. 1990; Abraham et al. 1989b), designed (35? and a^, respectively, ai? and had been obtained from equilibrium constants for the formation of 1:1 complexes between a solute molecule and a given reference base or acid, respectively, in CC14. We found good correlations between Fs max and a" and Vmin and , for groups of molecules (Murray and Politzer 1992). The relationships are very similar to those found for the solvent parameters a and p (Murray, Ranganathan, and Politzer, 1991). These findings confirmed that the calculated electrostatic potential, which refers to the molecule in the gas phase, can be quantitatively related to its tendency to form hydrogen bonds in solution, whether as a part of the solvent interacting with a solute or as a solute molecule forming a 1:1 complex with a reference system.

The correlations that have been described are all family dependent, a different one applying to each different group of structurally related molecules. Our next objective, accordingly, was to ascertain whether they could perhaps be made more general by improving the quality of the wave functions used to calculate the electrostatic potentials and, in the case of the hydrogen-bond-acceptor molecules, by using surface electrostatic potential minima (Vs min) instead of the three-dimensional spatial minima (Vlnin). Optimized structures and surface electrostatic potentials were computed for eighteen hydrogen-bond-donating and thirty-three hydrogen-bond-accepting molecules, at the HF/6-31G* level (Hagelin et al. 1995).

The eighteen hydrogen-bond-donors are listed in Table 3.1 along with their experimentally derived and statistically corrected1 a^1 values and our calculated Fs max. At the HF/6-31G* level, an excellent general correlation was found between ai? (corrected) and max; the latter is invariably associated with the hydrogen(s) to be donated in the hydrogen bond (Hagelin et al. 1995). This relationship, given as Eq. (3.5), is shown in Fig. 3.3. The correlation coefficient is 0.991 and the standard deviation is 0 .04.

For the group of thirty-three hydrogen-bond-accepting molecules in Table 3.2, we chose to seek correlations directly with the equilibrium constant, KHB, for 1:1 complexation of

Molecule |
Ui2 |
(kcal/mole) |

## Post a comment