1578
J. Phys. Chem. 1996, 100, 1578-1599
A Comprehensive Analytical Treatment of Continuum Electrostatics Michael Schaefer* and Martin Karplus* Department of Chemistry, HarVard UniVersity, 12 Oxford Street, Cambridge, Massachusetts 02138, and Laboratoire de Chimie Biophysique, Institut le Bel, UniVersite´ Louis Pasteur, 4, rue Blaise Pascal, 67000 Strasbourg, France ReceiVed: July 28, 1995; In Final Form: October 12, 1995X
An analytical treatment of the continuum model for electrostatic solvation is presented and shown to be in good agreement with numerical finite-difference continuum calculations for the same system. Because the analytic method, called ACE, is relatively fast and yields a pairwise potential from which gradients and second derivatives can be obtained, its use with molecular mechanics/molecular dynamics programs will be possible. For calculating the self-energies of molecules in aqueous solution, a new analytic method is introduced. It is based on a continuum electrostatics model in which charges are defined as Gaussian charge distributions, rather than point charges, and in which the solute interior is described by a continuous density function. The self-energy potential is a pairwise energy function which only depends on charges, interatomic distances, and the effective, solvent-inaccessible volumes of the atoms. The potential is derived by integrating the energy density of the electric field and termed the integrated field (IF) method. A combination of the analytic self-energy potential with the generalized Born equation for charge-charge interactions,1-3 termed analytic continuum electrostatics (ACE) potential, makes possible an analytic treatment of the electrostatic contribution to solvation. The focusing of field lines is included. The IF method is applied to the calculation of self-energies of all charges on a series of small organic compounds and the protein bovine pancreatic trypsin inhibitor (BPTI). The agreement with self-energies as calculated using a finite-difference program is excellent for small molecules and satisfactory for BPTI. Application of the combined ACE potential to the calculation of solvation energies of small compounds shows good agreement with results of finite-difference calculations. The potential on a plane through the protein superoxide dismutase using ACE and a finitedifference program is in good agreement. For the calculation of the solvation energies of small molecules (20 atoms), ACE is more than 3 orders of magnitude faster than the program UHBD91 which was used for all numerical continuum electrostatics calculations; for BPTI, the speedup is more than 2 orders of magnitude.
1. Introduction The importance of solvation effects on polar and ionic molecules has been recognized for many years.4-10 The solvation free energies of such molecules, which tend to be dominated by electrostatics, are of the order of 10 and 100 kcal/ mol, respectively.11,12 Particularly for macromolecules of biological interest, such as proteins and nucleic acids, it is essential to include electrostatic solvation in investigations of their stability and their functions, e.g., enzymatic processes, conformational alterations, molecular association, and substrate binding, all of which can involve important changes in the charge distribution and in solvation. Water influences the electrostatic free energy of a solute in two ways. It solvates each individual charge of the molecule, which determines the self-energies of the charges, and it screens the interaction between charge pairs. The two effects are related because, in a qualitative sense, the better a charge is solvated the more it is likely to be screened from interacting with other charges. Based on the assumption of a linear response of the polarizable solute-solvent system, the electrostatic energy of a solute with N charges q1, ..., qN can be decomposed into N ∝ qi2 and N(N - 1)/2 interaction energy self-energy terms Eself i int terms Eij ∝ qiqj. Both the self-energies of single charges and the interaction energies of charge pairs depend on the geometry of the solute. In molecular dynamics (MD) calculations in vacuum, a distance-dependent dielectric function was frequently employed to account, at least qualitatively, for the screening of X
Abstract published in AdVance ACS Abstracts, December 15, 1995.
0022-3654/96/20100-1578$12.00/0
charge-charge interactions while the self-energies of charges were not considered.13-18 It is the solvent effect on the interaction energies which leads to the fact that the solvation free energy of a molecule cannot be calculated as a sum of independent charge contributions.19-22 Since the self-energy of a single charge depends only on its location and on the solute geometry, it is independent of the magnitudes and locations of other charges. There are three types of approaches for describing electrostatic solvation in a more or less quantitative manner: (i) empirical methods related to the solvent-accessible surface of a solute, (ii) microscopic methods, and (iii) methods based on the continuum electrostatics model. Methods from all three categories have been incorporated into molecular mechanics programs.14,23-25 The empirical methods generally assume that the solvation energy is a sum of atom or group contributions. Each group contribution is approximated by a linear function of the solventaccessible surface of the group25-30 or the volume occupied by solute atoms within an appropriately defined first solvation shell.22,31-34 Although such an approach is theoretically flawed because it omits the solvent screening of interactions between charges, it has been applied with some success to the calculation of solvation energies of small organic molecules with as few as five different atom types and solvation parameters for amino acids.22,27 Despite the conceptual simplicity of these methods, their incorporation into molecular mechanics programs requires time-consuming calculation of gradients and second derivatives of the accessible surface area.25,30,35 This increases the com© 1996 American Chemical Society
Analytical Continuum Electrostatics Potential putation time to be nearly equal to that of an explicit solvent calculation in many cases. Interestingly, the application of an occupied volume method, which makes use of a Gaussian atom pair potential to approximate the fraction of solute volume within the first solvation shell around each charge, was shown to increase the time required for molecular dynamics simulation of the bovine pancreatic trypsin inhibitor (BPTI) by only 20%,22 relative to a simulation of the same protein in Vacuo. In the microscopic methods, the system of interest is generally embedded in a finite volume filled with simple spheres36 or more appropriate models for water molecules.37-40 The calculation of solvation energies on the basis of these models requires the performance of free energy simulations, either by a Monte Carlo (MC) sampling procedure40-43 or by molecular dynamics.44-47 Several studies have demonstrated good agreement for the solvation free energies calculated from microscopic methods with experimental data and with the results of continuum electrostatics methods.21,48 Both the MC and MD sampling methods are computationally several orders of magnitude more time-consuming than continuum methods. Clearly, cases where specific water molecules play an essential role (e.g., in enzymatic reactions) require explicit solvent calculations. The use of the microscopic methods limits drastically the time frame that is accessible by simulation or the number of configurations that can be investigated, compared to an analogous study of a system in Vacuo. For example, Stouten et al.22 report that MD simulation of BPTI in a water box with about 2350 water molecules increased the total number of atoms by a factor of 13 relative to the number of protein atoms and the required CPU time by a factor of 30 relative to a vacuum simulation of BPTI, despite use of a relatively short (12 Å) nonbonded cutoff distance. Approaches are being developed (e.g., fast multipole methods, parallel programs) to reduce the required computer time.49-52 Another method that accounts for solvent effects at the microscopic level is based on the solution of integral equations. The most widely used technique is the reference interaction site model (RISM).53 The “extended” XRISM method, which is applicable to polar and ionic molecular systems,21,54-56 has been applied to the calculation of solvation energies of small organic molecules. Although only approximate results are obtained,57 it is a very useful approach because it is much faster than explicit solvent calculations and because the mean field two-site potential commonly derived within RISM accounts for both electrostatic and dispersion interactions between solute atoms and the solvent. In the continuum electrostatics models, the solute interior and the solvent are defined as regions with different dielectric constants. For systems with simple solute geometry (e.g., a sphere as a schematic model of a protein), the Poisson equation or the Poisson-Boltzmann equation can be solved analytically, as in the original work of Kirkwood.4,6,58 A simplified approximation to the Kirkwood solution is the method of image charges.59,60 The Kirkwood solution has played a very important role in theoretical studies involving charge-charge interactions in proteins,61-64 despite the obvious loss of detail in the description of the molecular geometry. In the TanfordKirkwood theory of protein titration62 and in modified versions of that theory63,65,66 the self-energy contribution to the intrinsic pKa of the titratable groups in proteins is ignored. In fact, there is no obvious way to apply the analytical solutions of the Poisson-Boltzmann equation to the calculation of self-energies of charges within molecular structures that are treated as composed of atoms. This inadequacy is, in part, a consequence of the point-charge concept which leads to a divergence of the self-energy at the boundary between the solute and the solvent and, in part, a consequence of the dependence of self-energies
J. Phys. Chem., Vol. 100, No. 5, 1996 1579 on the somewhat arbitrary procedure of mapping a given molecular structure onto a spherical model system. The problem of mapping can be circumvented by introducing a spherical simulation region that contains the solute and a suitable number of explicit solvent molecules.67 This mixed microscopicmacroscopic approach has the problem of controlling the number of solvent molecules and pressure in the simulation region. However, due to the possibility of restricting the number of explicit solvent molecules that need to be included in the simulation, it is computationally less demanding than fully microscopic methods. With the advent of modern computers, numerical techniques could be employed to solve the Poisson-Boltzmann equation for molecular boundaries of arbitrary geometry, e.g., by use of the finite-difference algorithm.10,68-72 In the protein dipoleLangevin dipole (PDLD) method of Warshel and co-workers,14,73 the dielectric response of water is modeled by a grid of Langevin dipoles. Since the use of a discrete grid of Langevin dipoles is physically equivalent to the use of a finite-difference grid within numerical continuum electrostatics programs, the PDLD method falls under the category of macroscopic, rather than microscopic, electrostatic methods. The numerical methods have the advantage of treating the continuum electrostatics model rigorously, given sufficient resources in computer time and memory. The effect of the solvent dielectric on selfenergies (single charges) as well as on interaction energies (charge pairs) is accounted for by these methods. The numerical continuum methods have been applied successfully to a variety of questions of biological interest, including substrate binding,74 solvation energies,12,75 pK shifts and titration curves,76,77 and redox potentials.78 The results suggest that the continuum model provides a useful description at the molecular level for cases where the electrostatic contribution dominates. An important disadvantage of the numerical continuum electrostatics methods is that it is conceptually difficult to incorporate them into molecular mechanics programs, mainly because of the problem of assigning forces related to the dielectric boundary (“dielectric boundary pressure”) to individual atoms.24,79,80 Other problems are the convergence of the results as a function of the resolution of the solute-solvent boundary and charge representation81,82 and computational speed. Despite considerable progress in the implementation of numerical continuum electrostatics algorithms over the past years,83,84 a single calculation of energy and forces of an average-size protein typically requires of the order of seconds of CPU time on a Convex C-2 computer, which is 1-2 orders of magnitude slower than the effective solvation term developed by Sander and coworkers,22 for example. However, some combined numerical continuum electrostatics-MD calculations have been made, showing that it is possible to speed up the method significantly if a Poisson-Boltzmann solution is calculated only every 200 or 250 dynamics steps, as proposed by Sharp24 and by Niedermeier and Schulten.79 Because of the difficulties of incorporating numerical continuum electrostatics methods into molecular mechanics programs, the former are frequently used in a separate, postprocessing step for the study of structures that were generated by use of the latter;85,86 i.e., different physical models and, consequently, force fields are applied to the generation and to the evaluation of a given structure. To avoid this inconsistency and to include solvent effects on electrostatics in conformational searches in a time-efficient manner, the development of a simple and accurate analytical continuum electrostatic potential is required and has been the subject of several studies in recent years.3,20,60,87,88
1580 J. Phys. Chem., Vol. 100, No. 5, 1996 In this paper, we present an analytical treatment of continuum electrostatics (ACE) which is in good agreement with numerical continuum electrostatics methods. The new method can be used in a straightforward manner in molecular mechanics programs because it provides a direct way of calculating gradients and second derivatives of the electrostatic energy. Atomic selfenergies and charge-pair interaction energies are calculated in separate steps. For the calculation of self-energies in step I, we introduce a pairwise potential that is based on a continuum model in which charges are defined as Gaussian distributions and in which the solute volume is described by a continuous density function. Since the self-energy potential is derived by integrating the energy density of the electric field, it is termed the integrated field (IF) method. For step II, in which the charge-charge interaction energies in solution are calculated, we use the generalized Born (GB) equation1,2 as published in Still et al.3 According to the GB equation, the interaction energy between two charges depends on both the intercharge distance and the effective solvation radii of the charges as a measure of their solvent exposure. Thus, the GB equation accounts for solvent screening. The effective solvation radii of all charges can be determined from the charge self-energies calculated in step I. Several analytical continuum electrostatics methods for calculating self-energies have been proposed that are based on the concept of integrating the energy density of the electric field.3,20,87,88 The basic idea dates back to the work of Born,89 who derived an exact solution for the electrostatic solvation of ions, assuming spherical shape and a homogeneous surface charge density. The integration methods all approximate the electric field of a charge by the Coulomb field in an infinite, homogeneous dielectric. This approximation is equivalent to omission of the reaction field and is termed “Coulomb field approximation” in the following. It leads to an overestimation of the self-energy of surface charges (see section 2.5). In ref 87, a method for including the reaction field in the integration of the energy density is outlined which makes use of an iterative algorithm. Since this leads to significant complications, we neglect such effects in the present treatment. The numerical integration method for the Coulomb field approximation described in the supplementary material to ref 3 uses the solvent region as the integration volume, while the other self-energy methods are derived from an equivalent expression which includes an integral of the energy density over the solute volume. In the latter cases, this leads to analytical, pairwise self-energy potentials;20,87,88 every pairwise self-energy contribution involves one charged atom as the source of an electric field and another atom as part of the solute volume (integration volume). The main advantage of the self-energy potential introduced in this paper over the numerical integration method presented in the supplementary material to Still et al.3 is that it allows for the calculation of solvation forces and second derivatives of the self-energy. A related analytic approximation is the pairwise self-energy potential,88 which is proportional to 1/r4 at all distances r between the charge and the atom involved. This is identical to the long-range limit of the potential given by Schaefer and Froemmel.20 However, use of the functional form 1/r4 at all distances is likely to lead to an overestimate of the short-range contribution to the self-energy of changes. Together with the simplification of using a unique integration volume for all atom types in proteins, this may have been the reason for finding only a moderate correlation between calculated atomic self-energies and the results of finite-difference calculations. The self-energy potential in Schaefer and Froemmel20 does not diverge at zero charge atom distance and was derived for effective volumes that depend on atom type and
Schaefer and Karplus that fill spheres centered at the atomic positions. Despite these improvements, the assumption of spherical atomic integration volumes in conjunction with the volume overlap of bonded atoms leads to a discontinuous description of the solute volume and unphysical fluctuations of the self-energy potential in the interior of a molecule.90 This problem is overcome in the present implementation by the introduction of a Gaussian continuous function describing the atomic volumes. To test the present approach, the self-energy potential and its combination with the GB equation, termed ACE, are applied to a series of small organic molecules and to the proteins bovine pancreatic trypsin inhibitor (BPTI) and superoxide dismutase. To verify that the new analytical potential is a good approximation to the continuum electrostatics model, the results of the calculations are compared directly to data obtained from the finite-difference program UHBD (University of Houston Brownian Dynamics)91 when it is applied to the identical system with the same parameters. The comparison of the ACE potential with a numerical continuum electrostatics program makes possible a detailed and rigorous analysis of the accuracy of the new electrostatic method. Furthermore, this approach avoids the uncertainties that are present in a comparison with experimental data, such as the effect of the inclusion of nonpolar and hydrophobic contributions and the consideration of conformational equilibria in a comparison between calculated and experimental solvation energies. 2. Theory 2.1. Integral Formula for the Electrostatic Energy. In the continuum electrostatics model, a solute is defined as a region with low dielectric constant which is embedded in an infinite region with the dielectric constant of the bulk solvent; in the present treatment we are concerned with water, though other solvent could be used as well. The dielectric constants of the solute and the solvent are denoted by i and s, respectively, and we assume that there are no saturation effects in the dielectric response of the solute and solvent media.92 For simplicity, this study is restricted to zero ionic strength of the solution; i.e., the potential of the solute charges must satisfy the Poisson equation rather than the Poisson-Boltzmann equation, which applies to the general case of nonzero ionic strength. A generalization of the method for finite ionic strength is straightforward and will be the subject of a future publication. Finally, we use a fixed set of solute charges; in applications of the theory, we assume the solute charges are consistent with a particular choice for the solute dielectric constant, e.g., i ) 1 in conjunction with the CHARMM parameter sets param1993-95 and param22.96,97 Extension of the method to systems with titratable sites can be done in the same way as for numerical solutions of the Poisson-Boltzmann equation.76,77 Following the approach used to derive the Born solvation energy of ions,89 the electrostatic energy, Eel, of a solute with N charges q1 ... qN can be calculated by integrating the energy density u(x b) ) D B 2(x b)/(8π(x b)) of the electric field over all 98 space; i.e.,
Eel ) ∫R3u(x b) d3x )
1 ∫ DB2(xb) d3x + 8πs solvent 1 ∫ DB2(xb) d3x (1) 8πi solute
Here, (x b) is the position-dependent dielectric constant, and the dielectric displacement D B is a superposition of the contributions from the individual charges; i.e., D B ) ∑iD B i, where i ) 1, ..., N.
Analytical Continuum Electrostatics Potential
J. Phys. Chem., Vol. 100, No. 5, 1996 1581
By formally assigning a zero charge to uncharged solute atoms, we can assume in the following that the charge index i refers to the atom number (charged or uncharged) and that the number of solute atoms is equal to N, the number of charges. The symbol R3 denotes the infinite three-dimensional Cartesian space. By adding and subtracting the integral of D B 2/(8πs) over the solute volume, we can write
Eel )
1 ∫ DB2(xb) d3x + 8πτ ∫soluteDB2(xb) d3x; 8πs R3 τ)
1 1 - (2) i s
The dielectric displacement D B in both integrals of eq 2 must satisfy the boundary conditions on the electric field at the interface between solute and solvent; that is, the tangential components of the electric field, B E ) D B /, and the normal component of the dielectric displacement D B do not change when passing through the solute-solvent boundary. In the first integral term of eq 2, the dielectric displacement can be approximated by the Coulomb field. This has been shown to introduce a relative error of at most a few percent in the electrostatic energy.20 Qualitative considerations can also be used to justify this assumption. If the solute is small, all charges are highly exposed to solvent, so that there is a very small departure of the field from the Coulomb field. If, on the other hand, the solute is large, the dominant contribution to the electrostatic energy comes from the second integral in eq 2, so that the error introduced by approximating the first integral term is small on a relative scale. For point charges, the Coulomb field (dielectric displacement) b-b xi|. When integrating the square D B is given by D B ) ∑iqi/|x B j yield the of D B over all space, the off-diagonal terms 2D B i‚D B i‚D B j d3x ) qiqj/(srij), Coulomb interaction energy, (1/4πs)∫R3D whereas the diagonal terms, corresponding to the self-energy integrals of D B i2, diverge. To avoid the divergences, it is appropriate to introduce finite sizes for all charges, e.g., to assume that the charge qi is distributed homogeneously on the surface of a sphere with radius Ri, as first proposed by Born.89 Based on this assumption, the self-energy integral of charge i in a homogeneous dielectric with dielectric constant s is given B i2 d3x ) qi2/(2sRi). by (1/8πs)∫R3D With the first integral term in eq 2 given by the Coulomb and Born energies in a medium with the solvent dielectric constant, calculation of the electrostatic energy requires evaluation of the second integral over a finite volume, namely, the from solute volume. Distinguishing the self-energy terms Eself i interaction energy terms Eijint, one has
Eel ) ∑Eiself + ∑Eijint i
qi2 τ + ∫ D B 2(x b) d3x 2sRi 8π solute i
(4)
qiqj τ + ∫ D B (x b)‚D B j(x b) d3x srij 4π solute i
(5)
Eiself ) Eijint )
(3)
iC) )O -CH2-NH-CH3 )O/-OH/-NH2 -OH -NH3+ -SH -S-X*< -X*-C*-C*-
main chain amide nitrogen main chain CR; Cβ of Val, Ile, Thr; Cγ of Leu main chain carbonyl carbon; C of COOH, CONH2; Arg Cζ main chain carbonyl oxygen Cβ etc.; Gly CR Arg N methyl average O or N of COOH, CONH2; Arg Nη1,Nη2 aliphatic or phenolic hydroxyl Lys Nζ; N-terminus Cys Sγ Met Sδ; disulfide S ring Cγ; Tyr Cζ; Trp Cδ2, C2; Pro N, CR 5-ring C or N 6-ring C, except Trp Trp 6-ring C
13.4 15.4 8.7 15.6 25.5 25.5 39.0 21.1 19.8 34.4 39.2 25.7 8.7 18.8 20.0 22.1
Asterisk indicates ring atom; table taken from ref 20. b Volumes deduced from ref 103.
that are too small to accommodate a solvent probe sphere. In this work, we adopt the solvent-inaccessible volume definition because, with respect to electrostatics in aqueous solution, the main aspect of the solute volume is the absence of water; i.e., it creates a cavity with low polarizability. For a given, static solute structure, the property of low polarizability also applies to gaps between solute atoms that are too small to accommodate a solvent molecule. The solvent-inaccessible volume definition is the definition that is commonly implemented in numerical finite-difference programs to define the regions of solute and solvent dielectric constants;91,102 in particular, this is the definition employed in the UHBD program which we use for comparison. To evaluate the integral in eq 4, the solute integration volume can be formally subdivided into volumes associated with the individual atoms, A1 ... AN. We introduce a “molecular density b) describing the solute volume according to function” PS(x
PS(x b) )
{
x is inside solute volume 1, if b 0, otherwise
(11)
and rewrite the self-energy of charge i as a function that involves integration of the energy density τDi2/(8π) over all space
Eiself
qi2 τ ) + ∫ D B 2(x b) PS(x b) d3x 2sRi 8π R3 i
(12)
b) acts as a weight function on the energy density. where PS(x b), k ) 1, We now introduce an atomic density function Pk(x ..., N, describing the volume Ak associated with atom k. With respect to the electrostatic energy, the only condition applying to the atomic density functions is that their sum must be equal b) at any point b x, to the molecular density function PS(x
∑k Pk(xb) ) PS(xb)
(13)
This condition can be satisfied by many schemes for distributing the solvent-inaccessible volume among solute atoms. One possibility would be to construct Voronoy polyhedra, which have been used in many studies of the volume and packing of atoms in proteins.103-105 The advantage of using Voronoy polyhedra as atom volumes is that they do not overlap and that there are no voids between atoms. With respect to solute atoms at the surface, Voronoy polyhedra can be constructed by positioning solvent probe spheres on the solvent-accessible surface, as suggested by Finney.105 We emphasize that eq 13 b) of an individual does not imply that the density function Pk(x atom must be given by a step function, i.e., a function which is
equal to 1 inside Ak and zero otherwise. Other functional forms b) are allowed, especially functions of b x that are for Pk(x continuously differentiable, as long as the condition given in eq 13 is satisfied. Finney’s study of the volume of atoms in proteins,105 which used Voronoy polyhedra for the atoms, showed that for a given atom type the standard deviation of the volume distribution is only 10-15% of the mean value. This suggests that the solventinaccessible volume contribution of atoms in proteins is, to a first approximation, independent of the geometry and that effective volumes can be assigned to all atoms as a function of atom type. The volume contribution of an atom to the solventinaccessible volume is termed “effective volume” in the following. In the context of integrating the energy density of the electric field, the effective atom volume V ˜ k gives the size of the integration volume associated with atom k. From the effective atom volume, one can derive the radius R ˜ k of the ˜ k/4π)1/3; this quantity was atomic sphere according to R ˜ k ) (3V introduced in ref 20 and is termed the “effective atom radius”. In general, the effective radius of an atom is different from its van der Waals radius. It is introduced here solely as a measure of the effective volume and does not imply the use of spheres for the integration volumes Ak (see section 2.6). In Table 1, the effective volumes of the united atom types in the standard amino acids are listed; hydrogen volumes are not given because the publication on atom volumes and packing103 from which the table is derived was restricted to united atom types. The use of united atom types for effective atom volumes, i.e., the occurrence of (hydrogen) atoms with zero volume, does not lead to problems with the description of the solute volume and the theory presented in the following; however, the assignment of effective volumes to hydrogens will obviously lead to a more detailed description of the solute volume. To be able to assign effective volumes to atoms in solutes other than proteins, especially small organic compounds, we developed a utility program which uses the atomic coordinates and van der Waals radii and the radius of a solvent probe sphere as input. The program calculates solvent-inaccessible volume elements on a cubic grid and then assigns each volume element to the nearest atom, such that atom volumes can have very complicated shape. Since effective atom volumes are on the order of 10 Å3, a grid constant of 0.2 Å for the volume grid, leading to volume elements of 0.008 Å3, is sufficiently small to ensure that the error in the calculated atom volumes is less than 0.1%. Because of the assignment of volume elements to the nearest atom center, the algorithm implemented in the utility program is equivalent to the use of Voronoy polyhedra for calculating effective atom volumes. Furthermore, this method
Analytical Continuum Electrostatics Potential
J. Phys. Chem., Vol. 100, No. 5, 1996 1583
of determining interior volume elements is equivalent to the method that is commonly used in finite-difference programs, in particular the program UHBD which we use for comparison, to determine interior links between grid points.91,102 2.4. Pairwise Self-Energy Potential. The introduction of atomic density functions, eq 13, makes it possible to write the self-energy Eself i , eq 12, in the form of a pairwise atom potential:
Eiself Eikself )
)
qi2 2sRi
+ ∑Eikself
(14)
k*i
τ ∫ DB 2(xb) Pk(xb) d3x 8π R3 i
(15)
In eq 14, the self-energy of atom i is given by its Born selfenergy in the solvent dielectric s plus a sum of N - 1 integral terms Eself ik involving all other solute atoms k * i. The integral terms are proportional to a factor τ which is characteristic of a transfer from the dielectric s to i; i.e., Eself ik is the energy that is required to replace the solvent dielectric by the solute dielectric within the volume of atom k, provided that only the charge of atom i is present to generate an electric field. In the following, the charge of atom i will be referred to as “charge i”. In the usual case of a solvent with higher dielectric constant are positive. Since the than the solute, both τ and Eself ik dielectric displacement D B generally decreases with increasing distance from a charge, Eself ik represents an effective, repulsive interaction between charge i and atom k of a solute molecule. The term “effective interaction” is used to emphasize that atom k is involved solely because it contributes to the solute volume, preventing charge i from interacting with the solvent dielectric within the volume Ak; the absence of an energetically favorable interaction of charge i with the solvent leads effectively to a repulsive interaction with solute atom k. There are two limiting cases to the self-energy Eself i : the first corresponds to a single ion in solution and the second to a charged atom embedded in an infinite solute with uniform dielectric constant i. For an ion in solution, only the first term ) qi2/ remains in eq 14 which yields the correct result of Eself i (2sRi). In the second case, if the charge of atom i is distributed homogeneously on the surface of a sphere with van der Waals radius Ri, and the sum of the integration volumes ∑k*iPk represents all space except the van der Waals sphere of atom i, the sum ∑k*iEself ik in eq 14 must yield the Born transfer energy τqi2/(2Ri). When added to the first term qi2/(2sRi), this results in the correct Born self-energy qi2/(2iRi) within a dielectric i for the case of an infinite solute. In section 2.7, this limiting case will be used as a normalizing condition on the self-energy potential Eself ik (see eq 26). An interesting point about the case of a single ion in solution is that if it were possible to determine the individual volumes Ak exactly and to calculate the associated integrals exactly, it would not be necessary to exclude the contribution with k ) i from the sum in eq 14. For a single ion, the atomic integration volume Ai is given by the van der Waals sphere, while the field of a surface charge on the same sphere vanishes inside. This means that the self-contribution Eiiself is zero, so that the exclusion k * i is unnecessary. In what follows (see sections 2.6 and 2.7) we introduce an approximation for the effective potential Eself ik that is nonvanishing at zero distance rik because it is derived for a Gaussian describing the volume Ak of atom k and a Gaussian charge distribution of atom i. When employed for an ion in solution, such a potential Eself ik that is nonzero for
k ) i would lead to a departure of the self-energy from the Born energy. Consequently, we explicitly exclude the case k ) i from the sum in eq 14. 2.5. Coulomb Field Approximation. To evaluate the volume integral in eq 15, the dielectric displacement D B i(x b) at a point b x due to charge i at b xi is approximated by the Coulomb b) ) qi/|x b-b xi| for a point charge. This makes it field; i.e., D B i(x possible to integrate the energy density of the electric field without prior knowledge of the electrostatic potential; to introduce the latter would require an iterative procedure.87 The same approximation, termed Coulomb field approximation, was already used to simplify the first integral term in eq 2. Here, we address the question of its applicability to the self-energy integral in eq 15. Since the reaction field is defined as the non-Coulombic field contribution in an inhomogeneous dielectric, the Coulomb field approximation is equivalent to the assumption that the reaction field contribution to the energy density of the electric field can be neglected. In the usual case where the polarizability of the solvent is higher than the solute polarizability (s > i), the reaction field leads to the focusing of field lines toward the solvent and a decrease of the magnitude of the dielectric displacement in the solute interior. Therefore, omission of the reaction field when integrating D B i over Ak in the interior of a large solute can be expected to lead to an overestimation of the energy Eself ik . The Coulomb field approximation may be justified by the short-range nature of the electric field which decreases in a homogeneous dielectric with r-2 from any given charge, such that the change of the energy density of the electric field due to the reaction field can be expected to be small within the first solvation shell around a charge. Note that the first solvation shell around an atom with, e.g., radius Ri ) 2.0 Å accounts for about 58% of the self-energy and solvation energy, assuming 2.8 Å for the width of a solvation shell. To determine an estimate for the error in the self-energy, eq 4, that is due to the Coulomb field approximation, we consider the case of a charge with radius Ri at a distance d from an infinite, planar dielectric boundary; the dielectric constant on the charge side of the boundary and on the opposite side is set equal to the solute dielectric constant i ) 1 and the solvent dielectric constant s ) 80, respectively. For this geometry, the potential and electric field can be determined analytically from an image charge solution,98 while the field according to the Coulomb field b) ) qi/|x b-b xi|. By integrating the energy approximation is D B i(x density of the electric field over all space on the charge side of the boundary (the “solute volume” in this example), excluding points that are within a distance r < Ri from the charge, it is possible to calculate the relative increase of the self-energy integral due to the Coulomb field approximation; if the charge is in contact with the boundary, d ) Ri, the self-energy integral increases by 59%, while the relative increase for a charge that is separated by one atom layer from the boundary, d ) 3Ri, is only 9.4%.106 The case of a charged atom in contact with an infinite, planar boundary can be seen as an extreme example of a surface charge, for which the magnitude of the reaction field relative to the Coulomb field is substantially larger than for a charge in a small solute or for a charge in a macromolecule. The calculated error of 59% is, therefore, an estimate for the upper bound to the relative error in the integral term of eq 4 that is due to the Coulomb field approximation. In applications to biological molecules, in particular small solutes and solutes with a high degree of solvent exposure (including unfolded proteins), the error is expected to be much smaller. 2.6. Atom Charge and Volume Description. In general, the exact subdivision of the solute volume into atomic integra-
1584 J. Phys. Chem., Vol. 100, No. 5, 1996
Schaefer and Karplus When calculating the self-energy of q by integrating the energy density of the electric field with the volume density acting as a weight function (see eq 12), the error ∆E in the self-energy caused by the density fluctuation can be estimated to be
∆E ≈
Figure 1. Wire graph of the protein BPTI. The intersection of atom spheres with the paper plane shown as a gray-shaded area. Bonds between atoms below the paper plane are drawn with broken lines; effective atom radii R ˜ k are derived from effective volumes as given in Table 1. The use of spheres to describe atom volume gives rise to overlaps between bonded atoms and gaps between nonbonded atoms in the protein interior. The cross section is the same as in Figures 2 and 3.
tion volumes and subsequent integration of the energy density of the electric field over these volumes is computationally difficult. It is, therefore, necessary to make a simplifying assumption for the geometry of the individual atomic integration volumes. Based on the Coulomb field approximation, an analytical solution for the integral in eq 15 was presented ref 20, assuming the effective atom volumes to be spherical and that the atom charges were uniform surface charges on these spheres. When the resulting self-energy for proteins of known structure was compared with the results of numerical continuum electrostatics methods, it was found that the calculated selfenergy of a test charge as a function of position inside a protein fluctuates significantly (ref 90; see also Figure 9 in ref 20). The main reason for the fluctuations was that the density function b) was discontinuous and strongly fluctuating in the interior ∑kPk(x of a molecule, in contradiction to the expected constant (unit) density, eq 11. This arose from the assumption of spherical atom volumes in conjunction with the large difference between bonded and nonbonded nearest-neighbor distances. In the case of a fluctuating density, the sum of atom contributions to in eq 14 is no longer equivalent to the integral of the field Eself i energy over the solute volume in eq 4. To obtain an error estimate, we consider the case of two overlapping atom spheres, each with uniform density, with an overlap volume of V ) 4 Å3 (see Figure 1). Let the overlap volume be centered at a point which is at a distance of R ) 1.5 Å from a charged atom with charge q ) 0.5e, where e is the elementary charge. Inside the overlap volume, the volume density is equal to 2. We assume the overlap volume to be compensated by a gap (zero volume density) of equal size V ) 4 Å3, centered at a distance of D ) 2 Å from the overlap volume and at a distance of D + R ) 3.5 Å from the charge. A distance of D ) 2 Å is a typical distance between adjacent maxima and minima for density fluctuations inside a protein (see Figure 3a).
(
)
τV 2 τVq2 1 1 ≈ (D1 - D22) ) 8π 8π R4 (R + D)4 2.7 kcal/mol (16)
where D1 and D2 denote the average dielectric displacements at the overlap volume and at the gap volume, respectively; we assume dielectric constants of i ) 1 and s ) 80 such that τ ) 0.9875. The error is given relative to an integral of the energy density where a constant volume density of 1 is used for both the overlap and the gap volume. Thus, fluctuations of the order of (1 in the molecular density function ∑kPk can lead to an error of the order of 3 kcal/mol in the calculated self-energy of a charge. Since the typical distance between a maximum and a minimum of the volume density is D ) 2 Å, the error gives rise to a force of ∆E/D ) 1.5 kcal/(mol‚Å), which is of the same magnitude as the force between two partial charges q1 ) q2 ) 0.25e at a distance of r12 ) 3.9 Å in vacuum; i.e., fluctuations of the order of (1 in the density function describing the solute interior lead to significant unphysical forces acting on charges in the interior of a macromolecule. To overcome the problem of discontinuities and density fluctuations in the function ∑kPk, we introduce a rigorous continuum model in which both the atom charge distributions b) and the density functions Pk(x b) describing atom volume Fi(x are given by Gaussians; i.e., we write
b) ) qi π-3/2ai3 exp(-ai2(x b-b x i)2); ai ) Fi(x b) ) Pk(x
(π/2)1/2 Ri
4 exp(-(x b-b x i)2/(RR ˜ k)2) 3π1/2R3
(17) (18)
As before, Ri is the van der Waals radius of charged atom i. The relation between the parameter ai and the van der Waals radius Ri is derived from the condition that the self-energy of the charge distribution Fi in a homogeneous dielectric must ) (1/2)∫R3Fiφi d3x ) be equal to the Born energy, Eself i 2 ˜ k is the effective radius of atom k introduced in qi /(2Ri). R section 2.3. The definition of the atomic density Pk is such that the volume associated with it is equal to the effective b) d3x ) 4πR ˜ k3/3 ) V ˜ k. The effective volume volume, ∫R3Pk(x V ˜ k is independent of R, which is a smoothing parameter introduced to be able to vary the width of the volume density function Pk. Other choices for the descriptions of atomic charge and volume distributions are possible, e.g., a homogeneous surface charge distribution on the surface of a sphere instead of the Gaussian charge distribution3,89 or a cubic switching function describing the atomic volume density as a function of distance from the atom position. However, the requirements that must be met by any atom charge and volume description are (i) that the size of the charge distribution is nonzero and leads to the Born self-energy when integrating the energy density of the electric field over all space and (ii) that the sum of all solute b) is a good approximation to the function atom densities ∑kPk(x b) defining the solute volume according to eq 11. The PS(x requirement of nonzero size for charges is necessary to avoid the occurrence of infinite self-energies. The condition on the b) is necessary to ensure that the selfvolume density ∑kPk(x energy calculated as a sum of atom contributions, eq 14, is a
Analytical Continuum Electrostatics Potential
Figure 2. Molecular density function ∑kPk(x b) in a two-dimensional slice through the protein BPTI using Gaussian atom volumes, eq 18, with effective volumes taken from Table 1. Smoothing factor R equal to 1.0 (lower plot) and 1.8 (upper plot); same cross section as in Figures 1 and 3, molecule rotated by 90° around an axis perpendicular to the plane of the cross section.
good approximation to the self-energy given as an integral of the energy density of the electric field over the solute volume, eq 4. The scale parameter R in the definition of the volume distribution Pk of atoms controls the density fluctuations inside a molecule. A large R leads to a smooth density function ∑kPk in the interior of the solute at the expense of a less detailed description of the solute-solvent boundary. Figure 2 compares the molecular volume density in a two-dimensional slice through BPTI for two values of the smoothing parameter R. The fluctuations of the density ∑kPk at R ) 1.0 are of the order of (0.4 relative to the mean value of 1. The average density fluctuation for the model of overlapping atom spheres of Schaefer and Froemmel20 is also 0.4 and was found to lead to significant errors in the self-energy of a charge in the interior of a macromolecule (see eq 16). Therefore, the density function ∑kPk with R ) 1.0 is not suited for describing the solute volume. In contrast, the density function with R ) 1.8 leads to a much smoother density; the variation now is (0.1, relative to the mean value of 1. At the same time, the contour or “shape” of the protein as a whole is not greatly altered. Thus, use of a smoothing parameter in the neighborhood of R ) 1.8 leads to a density function that is a good approximation to the molecular density function PS in eq 11. Figure 3 gives a comparison between the solvent-inaccessible volume of a molecule and the molecular density function ∑kPk for R ) 1.0 and R ) 1.8. In Figure 3a,b, the same solventinaccessible volume in a two-dimensional slice through BPTI is shown as a gray-shaded area. In the calculation of the solventinaccessible volume, we approximated the solvent-accessible surface by the Shrake and Rupley dot surface,107 probing 800 surface dots per solute atom. The molecular density ∑kPk was calculated using Gaussian distribution of the atom volumes, eq 18. Effective volumes of all atoms were taken from Table 1. The density function with R ) 1.8 shows very little variation in the interior of BPTI when compared to the density with R ) 1.0. As one may expect from the definition of atomic volume density by a Gaussian, the contour of constant density at 1/e ≈ 0.368 follows closely the surface of the inaccessible volume, except that single atoms protruding from the molecule are no longer resolved by the density function with R ) 1.8. Since the volume density ∑kPk varies continuously between the solute interior and the solvent, as opposed to the mathematical definition of the solute-solvent boundary given in eq 11, one cannot directly deduce an error from the differences between
J. Phys. Chem., Vol. 100, No. 5, 1996 1585 the contour at 1/e and the contour of the solvent-inaccessible volume. The objective of resolving the molecular surface at the atomic level is clearly in conflict with the objective of reducing fluctuations in ∑kPk by increasing R. In section 4.1, a quantitative assessment of the error is given by determining the root-mean-square difference between self-energies calculated by integration, involving the continuous volume density ∑kPk at different values of R, and self-energies calculated using a finite-difference program in which the solute-solvent boundary is defined by the step function PS, eq 11. Of course, the best boundary to define the electrostatic properties of a physical system is not certain; e.g., in recent continuum calculations, Davis and McCammon81 have used a smoothing function at the boundary for the transition between the solute and solvent dielectric constants. This leads to better grid independence of the calculated potential, even though the smoothing algorithm itself is grid dependent because it interpolates between the solute and the solvent dielectric constants only for grid links that are passing through the solute boundary. As a result, the function describing the dielectric response approaches a step function when the grid constant is reduced (e.g., in a focusing calculation). 2.7. Self-Energy Approximation. In a homogeneous dielectric , the potential φi of the Gaussian charge distribution Fi, eq 17, is given by
φi(x b) ) qi
erf(ai|x b-b x i|)
(19)
|x b-b x i|
where “erf” denotes the error function. Based on the Coulomb field approximation, the dielectric displacement of a Gaussian charge i is
b) ) -qi∇ B D B i(x
erf(ai|x b-b x i|)
(20)
|x b-b x i|
Insertion of eq 20 and of eq 18 for the Gaussian volume distribution of atom k into eq 15 leads to the following expression for the energy Eself ik which atom k contributes to the self-energy of charge i:
Eikself
)
(
τqi2
∫ R
3/2 3 R3
6π
∇ B
) (
erf(ai|x b-b x i|) |x b-b x i|
2
exp -
)
(x b-b x k)2 (RR ˜ k)
2
d3x (21)
While there is no analytical solution to the integral in eq 21, one can derive an accurate approximation. The functional form chosen for the approximation is based on the assumption that the effective interaction Eself ik between charge i and atom k is finite and monotonically decreasing as a function of the charge atom distance rik. These assumptions fail at the point charge ˜ k) is above a critical limit, Ri f 0, and if the ratio x ) Ri/(RR value (xc ≈ 1) where one parameter involved in the approximation becomes undefined (see discussion below and Figure 5). The range of validity of the analytical approximation is, therefore, restricted to charge and effective atom radii satisfying is a nonzero e Ri e RR ˜ k, where Rmin the conditions Rmin i i minimum charge radius. To ensure the applicability of the analytical approximation, we formally assign charge radii vdW is the van der ,Rmin according to Ri ) Max(RvdW i i ), where Ri Waals radius of atom i, and we introduce an atom-pair˜ k) such that the dependent scale parameter Rik ) Max(R,Ri/R ˜ k is satisfied in all cases. For example, the condition Ri e RikR smallest effective volume in Table 1 is that of the carbonyl carbon (symbol “>C)”) and corresponds to an effective radius of R ˜ k ) 1.3 Å. If one uses the parameter set param19 of CHARMM to assign van der Waals radii to charges and a
1586 J. Phys. Chem., Vol. 100, No. 5, 1996
Schaefer and Karplus
Figure 3. Comparison of the solvent-inaccessible volume (gray-shaded area) with the molecular volume density function ∑kPk, eq 18, in a twodimensional slice through BPTI; effective atom radii R ˜ k derived from effective volumes as given in Table 1. (a) Smoothing parameter R ) 1.0, (b) R ) 1.8. Lines of constant volume density at ∑kPk ) 1/e (dotted), 0.5 (continuous), 1.0 (dashed), and 1.5 (- ‚ -). Solvent-inaccessible volume calculated by initially assigning every point within RvdW + Rprobe from any protein atom to the interior and then reassigning points within Rprobe from any dot on the Shrake and Rupley dot surface (800 dots per atom) to the outside; probe sphere radius Rprobe ) 1.4 Å.
smoothing parameter R ) 1.8, the only charge type with a radius ) larger than RR ˜ k ) 2.34 Å is the CR atom type with RvdW i 2.365 Å, leading to a scale parameter Rik ) 1.82; all other atom ˜ k. pairs in this example satisfy the condition Ri e RR By numerically integrating eq 21, we find that for atom pairs ˜ k the self-energy contribution Eself (ik) with Ri e RR ik decreases ˜ k) and then asymptotilike a Gaussian at short range (rik < RR ˜ k/ cally approaches the long-range limit which is equal to τqi2V decreases monotonically, leading to (8πrik4). The energy Eself ik a repulsive interaction between charge i and atom k at all distances rik. On the basis of this behavior, we make the following Ansatz for Eself ik :
Eikself )
(
)
τqi2V ˜k τqi2 rik3 exp(-rik2/σik2) + ωik 8π r 4 + µ 4 ik ik
4
(22)
The parameters ωik and σik determine the height and width of the Gaussian that approximates Eself ik in the short-range domain. The first term in eq 22 becomes negligible for large rik, and the second (long-range) term vanishes at rik ) 0 due to introduction of the parameter µik; it asymptotically approaches the limit ˜ k/(8πrik4) for large distances. τqi2V Other functions with the correct long-range limit, in particular functions of the form 1/P(rik) with P(x) a polynomial of fourth order, were investigated but failed to reproduce the self-energy potential eq 21 as well as eq 22. For example, the simple function f(rik) ) 1/(a4 + rik4) is not applicable because it is smaller than the (long-range) function 1/rik4 at all distances; this leads to a large departure from the self-energy potential in the intermediate distance range rik ≈ σik (see Figure 4 and related text). A study to establish the applicability of Pade´ approximants, i.e., fractions of two polynomials in rik with the polynomial in the nominator and denominator of order O(N) and O(N + 4), respectively, will be reported separately. To determine the parameters in eq 22, we use the fact that it is possible to evaluate the integral expressions for Eself ik and for 2 /∂r analytically at r ) 0 (x bi ) b xk the second derivative ∂2Eself ik ik ik in eq 21).108 Since the second term of the approximation to Eself ik , eq 22, as well as its second derivative with respect to rik vanishes at rik ) 0, this makes it possible to determine the
Figure 4. Contribution Eikself of volume V ˜ k to the self-energy of charge qi as a function of distance rik. Continuous curves, numerical integration of eq 21; long dashes, analytical approximation, eq 22; short dashes, ˜ k/(8πrik4) of the analytical approximation. long-range limit τqi2V Parameters: top curve, radii Ri ) R ˜ k ) 1.0 Å, smoothing parameter R ) 2.0, charge qi ) 31/2; middle curve, Ri ) R ˜ k ) 1.0 Å, R ) 1.0, qi ) 1; bottom curve, Ri ) 0.1 Å, R ˜ k ) 1.0 Å, R ) 1.0, qi ) 0.116; i ) 1 and s ) 80 such that τ ) 0.9875. The charges were chosen to obtain values of Eikself in the same energy range.
parameters ωik and σik according to
1 4 1 ) (Q - arctan Qik) ωik 3πR 3 ik RikR ˜k
(23)
ik
σik2 )
3(Qik - arctan Qik)
R 2R ˜ 2; (3 + fik)Qik - 4 arctan Qik ik k Qik ) fik )
qik2 (2qik2 + 1)1/2
;
( )
˜k 1 2 π RikR 2 ; q ) ik 2 Ri qik2 + 1 2qik2 + 1
Rik ) Max(R,Ri /R ˜ k); Ri ) Max(RivdW,Rimin)
2
(24) (25)
Both ωik and σik depend only on the charge radius Ri, the effective atom radius R ˜ k, and on the smoothing parameter R.
Analytical Continuum Electrostatics Potential
J. Phys. Chem., Vol. 100, No. 5, 1996 1587
To determine the long-range parameter µik, one can make use of the normalizing condition that in an infinite, homogeneous ˜ k, medium of atoms with volume V ˜ k and density of atoms 1/V the self-energy of charge i, eq 14, must be equal to the Born ) energy with the dielectric constant of the solute, i.e., Eself i qi2/(2iRi). By equating this energy with the right-hand side of eq 14 and subtracting qi2/(2sRi) on both sides, the condition becomes
τqi2 2Ri
∞
) ∑Eikself )
[
k*i
1 V ˜k
∫R Eikself(rik) d3rik 3
(
2 2 V ˜k τqi2 rik3 ∞ 2 exp(-rik /σik ) ) 4π∫0 rik + V ˜k ωik 8π r 4 + µ 4 ik ik
)] 4
drik (26)
where we changed the sum over the infinite number of selfenergy contributions Eself ik into the equivalent integral over all space. The factor τqi2 cancels out such that µik, too, depends ˜ k and the smoothing parameter only on the atomic radii Ri and R R. We solved the integral in eq 26 using the program Mathematica109 and obtained the following result for the parameter µik:
µik )
77π(21/2)
Ri
512(1 - 2π3/2σik3Ri /(ωikV ˜ k))
(27)
To test the accuracy of the analytical approximation, Figure 4 compares the energy Eself ik , determined from eq 22, as a function of rik with Eself ik obtained by numerically integratig eq ˜ k) of 0.1, 0.5, 21. The three cases shown are for ratios Ri/(RR ˜ k of the and 1.0, which all satisfy the validity condition Ri e RR analytical approximation. The approximation is exact at rik ) 0, and the absolute error is vanishingly small in the long-range domain, as expected. The absolute error is largest in the region of transition between the short- and long-range domains, which is characterized by rik ≈ σik, with σik given in eq 24. Relative to the maximum of Eself ik at rik ) 0, we found the error of the analytical approximation at any distance rik to be less than 5%. For comparison, Figure 4 also shows the curves for the long˜ k/(8πrik4) which was employed by Gilson and range limit τqi2V 88 Honig to calculate Eself ik at all distances rik. Since all solute atoms k * i contribute to the self-energy of a given charge i, including atoms that are bonded to atom i, the range of distances is relevant is rik > 1 Å. From the curves at which Eself ik ˜ k ) 1.0 Å in Figure 4 (top curves, R ) calculated for Ri ) R 2.0, σik ) 2.48 Å; middle curves, R ) 1.0, σik ) 1.85 Å), it follows that use of the long-range limit as an approximation to Eself ik at all distances leads to an overestimation of self-energy contributions and self-energy forces at short range, in particular for bonded atom pairs. In an intermediate distance range which is characterized by σik < rik < 2σik, the long-range function is considerably smaller than the effective self-energy potential Eself ik derived in this work. To analyze the dependence of the parameters ωik, σik, and ˜ k, we introduce the dimensionless µik on the radii Ri and R ˜ k), and yµ ) µik/(RR ˜ k) quantities yω ) ωik/(R3Ri), yσ ) σik/(RR ˜ k) of the charge which all depend only on the ratio x ) Ri/(RR radius and the scaled effective atom radius. In Figure 5, yω, yσ, and yµ are plotted as functions of x in the range 0 < x < 1.6. All three parameters yω, yσ, and yµ and, consequently, ωik, σik2, and µik are positive and well-defined in the range 0 < x < ˜ k, corresponding to 1.123. At the point charge limit Ri , RR the limit x f 0, the quantities yω, yσ, and yµ are finite, with values (3/2)π1/2 ≈ 2.66, 1, and 77(2π)1/2/256 ≈ 0.754, respectively. Since 1/ωik ) 1/(yωR3Ri), the first (short range) term of
Figure 5. Quantities yω ) ωik/(R3Ri) (continuous curve), yσ ) σik/ (RR ˜ k) (long dashes), and yµ ) µik/(RR ˜ k) (short dashes) as functions of the ratio x ) Ri/(RR ˜ k) of the charge radius and the scaled effective atom radius.
the self-energy potential, eq 22, diverges for zero charge radius Ri, leading to an undefined energy Eself ik . The quantities µik ) ˜ k and σik ) yσRR ˜ k diverge at critical values xc ≈ 1.123 yµRR and x′c ≈ 1.532. Because a divergent parameter µik leads to a second (long-range) term in eq 22 that vanishes at all distances rik, it is necessary to require that x < xc ) 1.123 to ensure the applicability of the analytical approximation. To guarantee that all parameters of the approximation vary slowly as a function of x, we set the upper bound for the ratio of the radii equal to x ) 1.0. In combination with the requirement that Eself ik is e R e RR ˜ for all finite, this leads to the condition Rmin i k i is charge atom pairs. The minimum charge radius Rmin i addressed in section 2.9. Inserting eq 14 for the self-energy into eq 9, the atomic solvation energy becomes
∆Eiself
)
Eiself
-
qi2
)-
2iRi
τqi2 2Ri
+ ∑ Eikself
(28)
k*i
where Eself ik is given in eq 22. The atomic solvation energy (transfer from a low dielectric to a high dielectric, i < s) is always negative because for every charge i, one has ∑k*i 2 Eself ik e τqi /(2Ri). In the last expression, the equal sign applies only to the limiting case of an infinite solute, where the sum ∑k*i Eself ik is equivalent to the integral of the energy density τD B 2i /(8π) over all space (see the normalizing condition on self Eself ik , eq 26). For a finite solute, the sum ∑k*i Eik corresponds to the integral of the energy density over the solute volume, i.e.,
∑ Eikself ) k*i
τ 8π
∫soluteDBi2 d3x
