Thomas A. Darden

and the other charge q2 is at position r (x, y, z), a distance r away [r (x2 y2 z2)1/2]. Then the force F (a vector) on charge q2 due to charge q1 is given by

whereas the force on q1 due to q2 is F.

In Eq. (1), rˆ r/r is the unit vector in the direction from charge q1 to q2, and k is the Coulomb force constant, which depends on the units being used. In many texts k is written as 1/4πε0, where ε0 is the permittivity of a vacuum. In Gaussian or cgs units the constant k is 1. If q1 and q2 are in units of elementary charge, so that for example the charge of a sodium ion is 1, distance is measured in angstroms, and energy is measured in kilocalories per mole (kcal/mol), k is approximately 332. For convenience k will be taken to be equal to 1 in the rest of this chapter.

The electrostatic energy of the pair of charges is the work that is required to move the second charge q2 against the electrostatic force, from infinity (where the pair experiences no interaction) along a path to the point r. This work is the integral of the dot product F dr (we are now dealing with vectors) of the force F with the infinitesmal displacement dr along the path from infinity to r, and is given by

Note that we have not defined which path we take from infinity to r, but it turns out not to matter. This is because F is curl-free away from the origin.

The force should be minus the gradient of the energy, F U, where the gradient, written U, is a vector whose direction is that of the maximal increase in U and whose length is equal to its rate of increase. In Cartesian coordinates the gradient can be calculated by taking partial derivatives with respect to components: U (∂U/∂x, ∂U/∂y, ∂U/∂z). Applying this to U given by Eq. (2), we see that indeed F U, where F is given by Coulomb’s law in Eq. (1).

The electric field E at charge q2 due to charge q1 is defined by

and can be thought of as the force on a unit test charge at the position of charge q2 due to charge q1 (i.e., set q2 1). Note that the force on q2 is now given by F q2E. The electric field due to q1 can also be defined by Eq. (3) at arbitrary points r ≠ 0, i.e., there need not be a charge present at the point r.

The electrostatic potential φ(r) at the position r due to the charge q1 is the energy of the pair when the second charge is a unit test charge:

Note that φ(r) is the work required to move a unit test particle from infinity to r when q1 is at the origin. Following the above argument, this work is obtained by taking the integral of the field E dotted with displacement, and thus, as above, the field E is the negative of the gradient of the potential φ, or E φ. As with the field, the potential φ(r) can be defined at a point r even when there is no charge at the point.

having total charge Q, enclosed in a volume V with closed surface ∂V, the surface integral over ∂V of the field due to the collection of charges is given by 4πQ.

If instead of discrete charges the charge is described by a smooth charge density ρ(r), the total charge Q contained within a volume V is given by the integral of ρ. Noting that the electric field is the negative of the gradient of φ, we can write Gauss’ law as

If we are considering electrostatics in dielectric media instead of in vacuo, which is necessary for a continuum treatment of solvation, Gauss’ law and Poisson’s equation must be generalized. In dielectric media, Gauss’ law is expressed in terms of the electric displacement D in place of the electric field E. In a linear isotropic dielectric medium, D εE, where ε is the permittivity of the medium. Gauss’ law in such a medium states that the surface integral of D over any closed surface is given by 4π times the total charge contained within.

As an example of the application of Gauss’ law in a dielectric medium, let us consider a simple continuum model of an ion in water. The ion is modeled as a point charge q at the center of a sphere of radius a that is immersed in a dielectric continuum. The interior of the sphere has dielectric constant 1, whereas the continuum has a larger dielectric constant ε. In the dielectric continuum, by spherical symmetry and Gauss’ law, the electric field at any point r distant r a from the center of the ion is given by E(r) (q/εr2) rˆ. Thus the electrostatic potential there is given by φ(r) q/εr. Inside the sphere the electric field is given for r ≠ 0 by E(r) (q/r2)rˆ. The electrostatic potential at r is given by its value on the surface of the ionic sphere plus the work to move the charge inside against the field, or φ(r) q/εa q/r q/a. If we remove the work required to move the test charge in from infinity to r against the unscreened vacuum electrostatic field due to the ion, which is given by q/r, we have the work performed by the dielectric medium in moving the test charge, or φdiel(r) (1 1/ε)q/a. This result can be used to calculate the work of charging the ion in the dielectric medium, arriving at the Born solvation free energy,

Although the continuum model of the ion could be analyzed by Gauss’ law together with spherical symmetry, in order to treat more general continuum models of electrostatics such as solvated proteins we need to consider media that have a position-specific permittivity ε(r). For these a more general variant of Poisson’s equation holds:

which is derived in the same way as Eq. (10). In this equation ρfree refers to the free charge density (for example, the charge density of the protein in a continuum treatment), as opposed to the charge density induced in the dielectric continuum or at dielectric boundaries.

Early in the twentieth century physicists established that molecules are composed of positively charged nuclei and negatively charged electrons. Given their tiny size and nonclassical behavior, exemplified by the Heisenberg uncertainty principle, it is remarkable (at least to me) that Eq. (1) can be considered exact as a description of the electrostatic forces acting between the atomic nuclei and electrons making up molecules and molecular systems. For those readers who are skeptical, and perhaps you should be skeptical of such a claim, I recommend the very readable introduction to Jackson’s electrodynamics book [1].

Thus if electrons as well as nuclei could be treated classically, molecular simulation would be much simpler. Since the nonelectrostatic forces operating on the system (strong nuclear force, weak force, etc.) are negligibly weak over molecular scales, the evaluation of the potential energy of a conformation would be a straightforward application of Eq. (6), where the charged particles would now be atomic nuclei or electrons, whose charges are known precisely. There would be no need for empirical parameters, and our only problem would be to generate a sufficient number of conformations of nuclei and electrons to average over. Unfortunately, electrons cannot be described classically in terms of configurations like the nuclear configurations; electron density is the relevant physical observable. However, it can be said that all molecular interactions can be derived from Coulomb’s law and the principles of quantum mechanics, together expressed in the Schro¨dinger equation.

The only problem with the foregoing approach to molecular interactions is that the accurate solution of Schro¨dinger’s equation is possible only for very small systems, due to the limitations in current algorithms and computer power. For systems of biological interest, molecular interactions must be approximated by the use of empirical force fields made up of parametrized terms, most of which bear no recognizable relation to Coulomb’s law. Nonetheless the force fields in use today all include terms describing electrostatic interactions. This is due at least in part to the following facts.

1.At long range, complex interactions having their origin in quantum mechanics are negligible, and molecular interactions can be described by classical electrostatics of nuclei and electron density.

2.At long range, the charge density need not be known precisely, and the interactions are well approximated by using a simplified representation, such as partial charges at the nuclei.

3.Although electrostatic interactions between pairs of molecules may be weak in many cases, a consequence of their long-range nature is that in large systems the energy due to electrostatic interactions must be calculated between all the pairs of the system and thus will dominate all other interactions. It is essential to include them.

As motivation for the rest of the chapter, a few further observations can be made. First, the calculation of full electrostatics is expensive. A typical molecular mechanics potential function is of the form

where E is the total molecular mechanics energy, Eb and Eθ are harmonic terms describing bond and angle vibrations, Eφ is a torsion term to describe the energetics of rotations about bonds, and Evdw and Ees are non-bond terms to describe interactions between atom pairs that are not part of a common bond, valence angle, or torsion angle. The van der Waals interactions in Evdw are described by dispersion and exchange repulsion terms originating in quantum mechanics, whereas Ees are the Coulombic interactions given by Eq. (6).

In large systems the computer time required to calculate the potential energy of a particular molecular conformation is dominated by the cost of calculating the non-bond interactions. This is due to the fact that the number of non-bond pairs is so much larger than the number of terms involved in the bond, angle, and torsion interactions. For example, in a system of 104 atoms, which is typical in current biomolecular simulations, there are about 104 bond terms and roughly the same number of angle and torsion terms (actually there are much fewer torsions if the system is mainly water, which is often the case). In contrast, there are N(N 1)/2 or about 5 107 non-bond pairs to be calculated. If these were calculated in a straightforward way the simulation would be very expensive.

The non-bond dispersion and repulsion terms decay rapidly, and thus the calculation of Evdw can be approximated by restricting the sum over all non-bond pairs to that over

˚

neighboring pairs (i.e., interactions are typically cut off at 8–10 A), reducing the number of interaction pairs by more than a factor of 10 in our example. However, electrostatic interactions do not decay rapidly, and thus the above cutoff approximation is quite inaccurate for Ees. Until recently, in order to speed the calculation, cutoffs were applied to electrostatic interactions as well as to the other more rapidly decaying non-bond interactions. Although this approximation may not always cause severe artifacts in the simulation, it is to be avoided whenever possible, especially if the simulation system contains mobile charged entities. A simple example of the problems that can occur was given by Auffinger and Beveridge [2]. They simulated a mix of sodium and chloride ions in water, under periodic boundary conditions. When cutoffs were applied to the electrostatic interactions they observed a strong tendency for ions of the same charge to be separated by distances close to the cutoff distance even when large cutoffs were used.

Ions of the same charge repel each other unless they are further apart than the cutoff distance, in which case they do not interact at all. Many workers had assumed that the dielectric screening of water would diminish artifacts due to the use of cutoffs. However, normal dielectric screening does not occur when cutoffs are applied. Bader and Chandler [3] observed a strong attractive well for the potential of mean force for the Fe2 –Fe3 ion pair when cutoffs were applied, whereas the Ewald summation (see below) led to normal dielectric screening of the repulsion. Thus cutoffs are not a viable option for accurate simulations, and therefore efficient calculation of full electrostatics is the focus of intense effort by a number of groups. Describing the currently popular methods is one goal of this chapter.

A second observation is that owing to the slow decay of Coulombic interactions with distance r, it is not straightforward to calculate the actual value of the Coulombic energy or potential in certain circumstances. This observation is by now well known for lattice sums in periodic boundary conditions, which we treat later in this chapter. However, it is also true for simple situations in nonperiodic boundary conditions. A simple example is the problem of calculating the electrostatic potential φ of a neutral atom, e.g., neon, at the center of an isotropic bath of water molecules. For simplicity, assume that the neon atom is at the origin of our coordinate system and the water molecules surrounding it are modeled with partial charges at the oxygen and hydrogen positions. We thus have a collec-

tion of charges, and we wish to calculate the electrostatic potential due to them at the origin. The electrostatic potential due to one of the charges qi at distance ri from the origin is simply qi/ri. If the water bath is infinite, we cannot immediately sum these individual contributions; instead, we must devise a limiting process that involves a sequence of finite calculations. For example, for 0 r ∞, let S1(r) denote the sum of qi/ri over all water atoms with ri r, and let S2(r) denote the sum of qi/ri over all water atoms belonging to water molecules whose oxygen is closer than r to the origin. Then consider the limits of S1(r) and S2(r) as r → ∞. These two converge rapidly, but to different limits! In fact, the limit of S1(r) is positive, whereas that of S2(r) is negative [4]. Thus although the electrostatic potential at the neon atom is clearly a physically reasonable quantity, it is not immediately clear how to calculate it. Hummer et al. [5] and Ashbaugh and Wood [6] have argued that the atom-based summation S1(r) leads to the correct result, whereas Aqvist and Hansson [7] have argued that a molecule-based summation such as S2(r) is correct. Let us consider this dilemma carefully.

The difficulty in calculating the potential at the center of the water bath is due to the slow decay of the Coulombic interactions, leading to the conditional convergence of infinite Coulombic sums. To explain the notion of conditional versus absolute convergence, consider the simple infinite series 1 1/2 1/3 1/4 1/5 . An infinite series converges absolutely if the sum of the absolute values of its terms converges. In this case that would mean that the sum 1 1/2 1/3 1/4 1/5 converges, which is not true. Thus this series is not absolutely convergent. The original alternating series, however, does converge, i.e., the sequence of its successive partial sums converges to a real number. However, if the series is rearranged as 1 1/3 1/2 1/5 1/7 1/4 , that is, two positive terms followed by a negative term, the partial sums now converge to a different limit [8]. The series is said to be conditionally convergent, meaning that the result is conditional on the method of summation. To summarize so far, it is possible for an infinite series to converge but only conditionally, meaning that the answer you obtain depends on the algorithm for calculating it. This strange behavior is not possible if a series converges absolutely.

If, in our example of the water bath, the individual terms qi/ri were replaced by their absolute values |qi |/ri, the resulting sum would be infinite. Thus the Coulombic sum is not absolutely convergent, and it is no surprise that S2(r) and S1(r) converge to different limits. If Coulombic interactions decayed faster than the inverse third power of distance, as do Lennard-Jones non-bond terms, their sums would converge absolutely and calculating them would be straightforward. However, since Coulombic interactions are longrange, care is needed in their calculation in order to arrive at a physically correct result.

Because of the spherical symmetry of the system, we can use Gauss’ law to deduce the correct value of the potential of a neon atom at the center of an infinite water bath. The electrostatic potential of the neon atom is the work required to bring a unit test charge from infinity, against the field E produced by the water molecules, to the origin. For 0 r ∞, let Q(r) be the sum of qi over all water atoms with ri r. By Gauss’ law the surface integral of E over the sphere of radius r is 4πQ(r) (note that we are performing the calculation in vacuo instead of in a dielectric medium, i.e., the water molecules will provide the intrinsic dielectric screening), and by spherical symmetry the field at any point on the surface of the sphere is Q(r)/r2. Using this result and simple calculus we can calculate the potential φ at the origin. The result agrees with the limit of S1(r), and the dilemma is resolved. The details and application to ionic charging free energies are found in Ref. 9.