Representing an Infinite Solvent System with a Rectangular Finite

In the RIC method, water molecules outside a rectangular box are treated ... The total electrostatic energy of the system is given by The first term i...
0 downloads 0 Views 243KB Size
Download PDF
J. Phys. Chem. B 2002, 106, 2973-2982

2973

Representing an Infinite Solvent System with a Rectangular Finite System Using Image Charges Pei-Kun Yang,† Shwu-Huey Liaw,†,‡ and Carmay Lim§,|,* Institute of Molecular Medicine, College of Medicine, National Taiwan UniVersity, Taipei, Taiwan, R.O.C., Department of Life Science, National Yang Ming UniVersity, Taipei, Taiwan, R.O.C., Institute of Biomedical Sciences, Academia Sinica, Taipei 11529, Taiwan, R.O.C., and Department of Chemistry, National Tsing Hua UniVersity, Hsinchu 300, Taiwan, R.O.C. ReceiVed: July 26, 2001; In Final Form: NoVember 30, 2001

We have developed a novel solvent boundary potential for nonspherical charged systems based on a Rectangular Image Charge (RIC) method. This method is similar to the Kirkwood theory in that both methods solve for the potential in one dielectric cavity immersed in another dielectric material. The key difference is that the RIC method employs rectangular coordinates as opposed to spherical coordinates in the Kirkwood theory and hence can be applied to nonspherical systems. The RIC method yields electrostatic reaction fields and energies in accord with those obtained by numerical solution to a Poisson equation. It also yields accurate cation-water oxygen radial distribution functions as well as absolute free energies. Furthermore, the absolute free energy of an ion was verified to be independent of its location in the cubic box, indicating that the present method can extend the finite cubic system into an infinite solvent system.

Introduction Most chemical and biological processes occur in solution. Solvent effects play a critical role in reaction rates, equilibrium structures, and electronic spectra. Consequently, over the past decade various methods have been developed to treat solvent effects on the static and dynamical properties of molecules. These approaches generally fall under the following three groups. The most detailed approach is to employ molecular dynamics (MD) simulations with an explicit solvent using periodic boundary conditions (PBC).1 At the other end of the spectrum are implicit solvent models in which solvation effects are treated using continuum dielectric models2-5 or integral equations.6 A hybrid approach is to include a small number of solvent molecules explicitly, and approximate the influence of the remaining bulk solvent by an effective solvent boundary potential or reaction field.7-12 Application of PBC to nonperiodic, inhomogeneous, and nonequilibrium systems presents several difficulties. Although MD simulations using PBC can generally incorporate solvation effects on complex biomolecules, they are very computer intensive as many solvent molecules are needed to mimic bulk solutions. This is especially the case in simulations of highly charged species such as DNA or RNA in solution, where it is necessary to include a large enough number of water molecules to neutralize the repulsion among the negatively charged phosphate groups. Furthermore, it has been difficult to obtain stable trajectories without imposing added restraints or artificially modifying the charges on the phosphate.13 When the longrange Coulombic interactions are summed over an infinite periodic array using the Ewald method, the instantaneous fluctuations are replicated throughout the infinite system; hence * Author to whom correspondence should be addressed. † National Taiwan University. ‡ National Yang Ming University. § Academia Sinica. | National Tsing Hua University.

artificial counterions have to be added to neutralize the system.14 This presents difficulties in computing free energies of charged species. Implicit solvent models have a distinct advantage over MD simulations using PBC, especially for large biomolecular systems, in speeding up computation considerably by not taking solvent molecules into account explicitly. In particular, continuum dielectric methods, in which the solvent is approximated as a structureless continuous medium characterized by a dielectric constant, can yield relatively accurate electrostatic solvation free energies of small solutes.15-18 However, neglecting all atomic and structural details of the solvent molecules would not be appropriate in cases where solute-solvent hydrogen bonding interactions play a critical role in determining the thermodynamic and/or kinetic stability of the solute. Hybrid “microscopic” and “macroscopic” approaches provide a good balance between computational cost and accuracy. Several theoretical formulations of a solvent boundary potential that takes into account the influence of the bulk solvent surrounding a spherical region consisting of the solute hydrated with a finite number of explicit solvent molecules have been proposed, e.g., the mean force field approximation,8 the surface constrained all-atom solvent (SCAAS) model,19 reaction field with exclusion,20 image charge approximation,10,21 and the spherical solvent boundary potential.9,11 These solvent boundary potentials are based on the assumption that the simulation system is spherical on average. However, the instantaneous shape of the molecular surface is not perfectly spherical. Furthermore, many biomolecules such as DNA or RNA are not spherical and even proteins are often ellipsoidal rather than globular in shape, thus using a spherical solvent boundary would not be computationally effective for these nonspherical systems. Consequently, we have developed a novel solvent boundary potential for nonspherical systems based on a Rectangular Image Charge (RIC) method. This method is similar to the Kirkwood theory,22 which has been used to represent the reaction field in

10.1021/jp012900n CCC: $22.00 © 2002 American Chemical Society Published on Web 02/26/2002

2974 J. Phys. Chem. B, Vol. 106, No. 11, 2002

Yang et al.

lim(z f 0-) Ex ) lim(z f 0+) Ex

(1a)

lim(z f 0-) Ey ) lim(z f 0+) Ey

(1b)

lim(z f 0-) 1Ez ) lim(z f 0+) 2Ez

(1c)

where Ex, Ey, and Ez are, respectively, the x, y, and z components of the electrostatic field E. The image charge Q′ is given by24

 1 - 2 Q′ ) Q 1 + 2

(2)

Thus, the potential φ(r) can be evaluated by Q and Q′ in a homogeneous dielectric 1:

φ(r) )

Q Q′ + 4π1r1 4π1r2

(3)

Extension of the Method of Images to a System with Two Infinite Boundaries. Consider a point charge Q(x0) in a dielectric a bounded by planes a1 and a2, which separate the medium from another dielectric w (Figure 2a). The electric potential φ(x) in a can be evaluated by placing image charges Ql at positions xl, as shown in Figure 2b:

φ(x) )

∞

Ql

l)-∞

4πa|x - xl|

∑

where Figure 1. (a) A point charge Q in a dielectric 1, a distant d away from a plane interface, which separates the medium from another semiinfinite dielectric 2. (b) The potential φ in 1 is solved using the method of images by replacing the dielectric material 2 by 1, and placing an image charge Q′ at a distance d away from the plane interface.

MD simulations of spherical systems.10-12 Both methods solve for the potential in one dielectric cavity immersed in another dielectric material. The key difference is that the Kirkwood theory solves this problem using spherical coordinates, while the RIC method employs rectangular coordinates and hence can be applied to nonspherical systems. In the RIC method, water molecules outside a rectangular box are treated implicitly by including their electrostatic and van der Waals (vdW) effects on atoms in the simulation region. The formulas for treating the electrostatic interactions of infinite bulk solvent are given in the next section. To validate the RIC method, electrostatic fields and energies were computed and compared to respective values obtained by solving a Poisson equation using finite difference methods. Furthermore, MD and free energy simulations were carried out for Na+ and K+. The structural and thermodynamic properties of these ions were compared to respective experimental data and simulation results obtained using a spherical solvent boundary potential8,11 and the minimum image convention.23 Theory Consider a point charge Q in a semi-infinite dielectric 1, a distance d from a planar interface that separates this medium from another semi-infinite dielectric 2 (Figure 1a). The surface is taken as the plane z ) 0. In the method of images, the electric potential at a distance r1 from the charge can be evaluated by placing an image charge Q′ in 2 at a distance d to the infinite plane. The image charge Q′ is chosen to satisfy the following boundary conditions at z ) 0 (Figure 1b) :

Ql ) Q

(

)

a - w  a + w

(4)

|l|

(5)

and

xl ) (-1)l x0 + la

(6)

In eq 6, a ) a2 - a1, the distance between the two boundaries, and the center of the system corresponds to a1 + a2 ) 0. Extension of the Method of Images to a 3-Dimensional System. Consider a point charge Q(x0,y0,z0) in a dielectric a bounded by a rectangular box of lengths a, b, and c, which separates the medium from another dielectric w. The origin is set equal to the center of the box. Using the image charge method, the electric potential φ(x,y,z) or φ(R) in a can be evaluated by

φ(R) )

∞

Ql,m,n

l, m, n ) -∞

4πa|R - Rl,m,n|

∑

where

Ql,m,n ) Q

(

)

a - w  a + w

(7)

|l|+|m|+|n|

(8)

and

Rl,m,n ) {xl, ym, zn} ) {(-1)lx0 + la, (-1)my0 + mb, (-1)nz0 + nc} (9) Extension of the Method of Images to a Multi-Charged System. Consider N point charges, Q1(x1, y1, z1), Q2(x2, y2, z2) ... QN(xN, yN, zN), in a dielectric a bounded by a rectangular box of lengths a, b, and c, which separates the medium from

Representing an Infinite Solvent System

J. Phys. Chem. B, Vol. 106, No. 11, 2002 2975

Figure 2. (a) A point charge Q(x0) in a dielectric a; the boundaries a1 and a2 separate the medium from another dielectric w. (b) The electric potential φ(x) in a can be evaluated by the method of images: the image charges Qi+1 and Q-(i+1) (i ) 0, 1, 2, ...) are induced by Q-i and Qi to satisfy the boundary conditions at the a2 and a1 planes, respectively.

another dielectric w. Using the image charge method, the electric potential at a distance Ri ) {xi, yi, zi} in a can be evaluated by

Qj

N

φ(Ri) )

+∑ ∑ j)1 4π |R - R | j)1

j*i

a

i

j

N

∞

N

term is the interactions between each charge and all image charges. The force at atom i is given by

Qj,l,m,n

∑ l,m,n)-∞ l2+m2 +n2 * 0

Fi(R))

4πa|Ri - Rj,l,m,n|

∑ j)1 j*i

QiQj(Ri - Rj)

∑ ∑ j)1 l,m,n)-∞

where

(

|l|+|m|+|n|

(11)

N

Rj,l,m,n ) {xj,l, yj,m, zj,n} ) {(-1) xj + la, (-1) yj + mb, (-1) zj + nc} (12) n

The total electrostatic energy of the system is given by N

U)

N

QiQj

+ ∑ ∑ i)1 j)i+1 4π |R - R | a

i

1

∑ j)1 j*i

m

j N

∞

∑ ∑ 2i,j)1 l,m,n)-∞ l2+m2+n2*0

QiQj,l,m,n 4πa|Ri - Rj,l,m,n|

(13)

The first term in eq 13 is the direct electrostatic interactions between the real charges in the simulation region. The second

4πa|Ri - Rj,l,m,n|3

(14)

while the electric field at atom i is given by

Ei(R) )

and

l

QiQj,l,m,n(Ri - Rj,l,m,n)

l2+m2+n2*0

)

 a - w  a + w

∞

N

(10)

Qj,l,m,n ) Qj

+

4πa|Ri - Rj|3

Qj(Ri - Rj)

+

4πa|Ri - Rj|3 N

∞

∑ ∑ j)1 l,m,n)-∞ l2+m2+n2*0

Qj,l,m,n(Ri - Rj,l,m,n) 4πa|Ri - Rj,l,m,n|3

(15)

Figure 3 illustrates the image charges induced by a system containing 10 water molecules in a 23.6 Å simulation cube for -1 e {l, m, n} e 1. In computing the image charges using eq 11, a was set equal to 1 and w to 80. As water molecules approach the implicit solvent boundary (see Figure 4), the second term in eq 13 will go to infinity, hence vdW interactions between water oxygen atoms in the simulation and implicit solvent regions were included to prevent water molecules from approaching too close to the implicit solvent boundary.

2976 J. Phys. Chem. B, Vol. 106, No. 11, 2002

Yang et al.

Figure 3. A two-dimensional illustration of a 23.6 Å box containing 10 water molecules and their image charges.

vdW Boundary Force. Since the distance to the first peak of the water oxygen-water oxygen radial distribution function (rdf) is ∼2.8 Å25 (see below) the implicit solvent region starts at 2.8 Å from the edge of the box. Atoms can move in the simulation region in Figure 4; i.e., in the box of length (L + 5.6) Å, but are prohibited from entering the implicit solvent region. The force on an explicit water molecule is assumed to stem from the average water structure in the implicit solvent region. Therefore, the water oxygen at position r in the simulation region (Figure 4) was subjected to a boundary force of the form8

F(r) )

∑ Fg(rOO) FvdW(rOO) dV

(16)

In eq 16, the summation is over all implicit water molecules in the hatched region in Figure 4, F is the water density, g(rOO) is the water oxygen-water oxygen rdf, and FvdW is the vdW force between the water oxygen atoms. Details of the derivation and implications of eq 16 can be found in ref 8. Calculations The reaction field components (eq 15) and energies (eq 13) were computed with -51 e {l, m, n} e 51 for atoms at various positions in a cube and compared to corresponding quantities obtained by solving a Poisson equation using finite difference methods15,26,27 to verify the accuracy of the RIC method. In addition, the ith image shell reaction energy (eq 13 with {l, m, n} ) (i) was calculated for various particle positions in a cube to elucidate the convergence of eq 13. This led to an approximation for the total reaction energy and force in terms of the 1st image shell reaction energy and force (eqs 13 and 14 with {l, m, n} ) (1), as described in the Results section (see eqs 20 and 21 below). Next, a pure water conventional MD simulation was carried out using PBC with Ewald summation to compute the water oxygen-water oxygen rdf, g(rOO), for evaluating the vdW force (eq 16). It was found that the vdW forces from eq 16 could be simply approximated by eq 23 (see Results). The latter approximation in conjunction with eqs 20 and 21 below for the electrostatic energy and force, respectively, was incorporated into the CHARMM28 program.

Figure 4. The RIC system. Water molecules in the simulation region cannot enter the implicit solvent region (hatched). The distance between the edge of the box (light solid) and the implicit solvent boundary (bold solid) is 2.8 Å, the distance to the first peak position of the water oxygen-water oxygen rdf. The boundary force on a water oxygen in the simulation region was calculated by summing its vdW interactions with water oxygen atoms in the implicit solvent region (see the Calculations section).

Simulations using the RIC method were carried out for 33 and 267 water molecules in a 10 and 20 Å cube, respectively, to verify the water density along the X, Y, and Z axes. In addition, RIC simulations were performed for Na+ or K+ solvated by 33, 113, and 267 water molecules in a 10, 15, and 20 Å cube, respectively. To verify the ion-solvent structures and energetics resulting from simulations using the RIC method, reference simulations were carried out for Na+ and K+ in a 20 and 40 Å cube containing 267 and 2137 water molecules, respectively, using PBC as well as in a 6, 12, and 20 Å-radius sphere containing 33, 267, and 1118 water molecules, respectively, using spherical boundary conditions.7,28 The latter employs a spherical solvent boundary potential similar to the vdW boundary force described above in that the force on an explicit water molecule is assumed to stem from the average water structure in the implicit solvent region. Finally, free energy simulations using the RIC method were performed for Na+ or K+ solvated by 33, 113, and 267 water molecules in a 10, 15, and 20 Å cube, respectively, and the resulting solvation free energies were compared to those obtained in previous work using spherical boundary conditions.11 Reaction Field from the Rectangular Image Charge Method. A finite (instead of an infinite) number of image charges was used in eqs 13 and 15 to compute the reaction energy and field, respectively. This approximation was tested for a cube of length 1.0 Å containing two ions of unit charge 1e, Q1(x1, y1, z1), and Q2(x2, y2, z2). Various geometries in the cube were generated by placing Q1 and Q2 at positions ranging from (-0.4 Å, -0.4 Å, -0.4 Å) to (0.4 Å, 0.4 Å, 0.4 Å) in increments of 0.2 Å. For a given geometry of the ions in the cube, the electrostatic energy at R1 ) (x1, y1, z1) due to Q2 at (x2, y2, z2) was calculated using eq 13 with -51 e {l, m, n} e 51 and -1 e {l, m, n} e 1 to yield U(R1) and U1(R1), respectively. Reaction Field from a Finite-Difference Solution to the Poisson Equation. The program DELPHI26,27,29 was used to

Representing an Infinite Solvent System

J. Phys. Chem. B, Vol. 106, No. 11, 2002 2977

TABLE 1: Reaction Field Components (in kcal/mol/e/Å) and Energies (in kcal/mol) from Finite-Difference Solution to the Poisson Equation and the RIC Method DELPHI atom

X

Y

Z

charge

1 2 3

-3.0 3.0 3.0

0.0 0.0 3.0

0.0 0.0 0.0

1.0 2.0 -2.0

1 2 3

-3.0 3.0 3.0

0.0 0.0 3.0

0.0 0.0 0.0

1 2 3

-3.0 3.0 3.0

0.0 0.0 3.0

1 2 3

-3.0 3.0 3.0

0.0 0.0 3.0

Ex

image charges

Ey

Ez

U

Ex

Ey

Ez

U

1.00 0.33 0.65

1.10 1.59 2.02

0.00 0.00 0.00

-21.3

0.95 0.34 0.63

1.07 1.52 1.92

0.00 0.00 0.00

-20.7

1.0 0.0 0.0

0.94 0.58 0.61

0.00 0.00 0.04

0.00 0.00 0.00

-15.8

0.89 0.56 0.59

0.00 0.00 0.04

0.00 0.00 0.00

-15.6

0.0 0.0 0.0

0.0 2.0 0.0

-1.16 -1.89 -1.64

0.00 0.00 0.30

0.00 0.00 0.00

-63.1

-1.13 -1.78 -1.56

0.00 0.00 0.28

0.00 0.00 0.00

-62.2

0.0 0.0 0.0

0.0 0.0 -2.0

1.22 1.64 1.68

1.10 1.59 1.68

0.00 0.00 0.00

-67.3

1.18 1.56 1.60

1.07 1.52 1.60

0.00 0.00 0.00

-66.3

calculate the reaction field and electrostatic energy for atoms in a 20 × 20 × 20 Å3 cube. The cubic cavity was constructed by placing dummy atoms of radius 1 Å and zero charge at positions ranging from (-9.0 Å, -9.0 Å, -9.0 Å) to (9.0 Å, 9.0 Å, 9.0 Å) in increments of 1.0 Å. Charges were then assigned to atoms at specific positions in the cubic box, which was treated as a low dielectric medium with an internal dielectric constant of unity surrounded by a high dielectric solvent. The ionic strength was set to zero. The Poisson equation was solved by finite difference methods for the electrostatic potentials using a 65 × 65 × 65 grid with an initial spacing of 1.40 Å, and refined with a spacing of 0.35 Å.15 The difference between the electrostatic field components in solution (out ) 80 for water) and in the gas phase (out ) 1) was calculated as well as the total electrostatic energy. The resulting electrostatic energies and fields were compared to the respective values derived from eqs 13 and 15 with -51 e {l, m, n} e 51. vdW Boundary Force. To compute the force on a water oxygen at position r in the simulation region (eq 16), the water oxygen-water oxygen rdf, g(rOO), was calculated from a room temperature 300 ps simulation of a 28 Å cube containing 729 TIP3P30 water molecules using PBC with Ewald summation. The resulting water oxygen-water oxygen rdf compares reasonably well with that obtained from neutral diffraction data,25 and reproduces the distribution of nearest neighbors centered at 2.8 Å. The vdW force between the water oxygen atoms, FvdW(rOO), was computed using TIP3P30 oxygen Lennard-Jones parameters. F(r) (eq 16) was computed using a grid in the implicit solvent region with various grid spacings and cutoffs to ensure that the forces are converged with respect to these two parameters. It was found to be approximated by eq 23 (see the Results section), which was hence employed in the MD simulations. For the reference simulations, the spherical solvent boundary potential was approximated by a quartic function, as implemented in the CHARMM program.28 Simulation Protocol. The simulations were carried out in an NVE ensemble using the CHARMM program28 modified to include the RIC method (see the Theory section above). The TIP3P30 model of water was employed and the bonds with water hydrogen atoms were fixed with the SHAKE algorithm.31 Each ion (Na+ or K+) was fixed at the center of a cubic box of length L containing water molecules at an experimental density of 0.0334 molecules/Å3. The Lennard-Jones parameters of the ions were taken from previous work.11 For the RIC simulations no nonbond cutoff was applied, whereas for the reference simulations using spherical boundary conditions or PBC a nonbond cutoff equal to the radius of the sphere or half the cube length was employed. All atoms were propagated according to New-

ton’s equations using the leapfrog Verlet integrator and a timestep of 2 fs at a mean temperature of 298 K. Each system was first minimized for 1000 steps, then equilibrated for 20 ps and subjected to 200 ps of production dynamics, during which coordinates were stored every 20 fs for analysis. Free Energy Perturbation with Image Charges. Solvation free energies of the ions were obtained from free energy simulations using the same simulation protocol as described in the previous section. The free energy for charging an ion was computed using free energy perturbation methods with λ ) 0 corresponding to an uncharged dummy atom (i.e., zero charge and zero Lennard-Jones parameters) and λ ) 1 to the fully charged ion (q ) 1). The Helmholtz free energy ∆A is given by32,33

∆A )

∑i ∆Ai(λi f λi+1)

(17)

where the free energy ∆Ai between two states, λi and λi+1, is

∆Ai (λi f λi+1) ) -kBT ln〈e-[∆H(λifλi+1)]/kBT〉i

(18)

In eq 18, kB is Boltzmann’s constant, T is the temperature, 〈...〉i indicates averaging over the ensemble of configurations representative of state i, and the Hamiltonian varies from Hi to Hi+1 as follows:

∆H(λi f λi+1) ) {H(Iλi+1, Wλi) - H(Iλi, Wλi)} + {H(Iλi+1, I′λi+1) - H(Iλi, I′λi)} + {H(Iλi+1, W′λi) H(Iλi, W′λi)} + {H(Wλi, I′λi+1) - H(Wλi, I′λi)} (19) The first term in parentheses involves the interactions between the ion I and explict water molecules W, as in conventional free energy simulations. The second and third terms in parentheses involve the interactions between the ion I and its image I′ as well as water images W′, while the last term corresponds to the interactions between water W and the image of the ion I′. The solvation free energy of an ion was computed by perturbing the ion to a dummy atom of zero charge and radius using 10 windows with λi varying from 0.05 to 0.95 in increments of 0.1. At each window, the system was equilibrated for 5 ps simulation, followed by 20 ps of production dynamics. Perturbation energies were collected every 2 fs for computing the relative free energies. The reverse perturbations were performed from a configuration that is independent from the forward perturbation runs. The free energies for the forward and reverse runs were then averaged.

2978 J. Phys. Chem. B, Vol. 106, No. 11, 2002

Yang et al.

Figure 5. The reaction energy of two ions of unit charge at various positions n representing {x1, y1, z1, x2, y2, z2} in a cube of length 1 Å. (a) The reaction energy from the first shell of image charges; i.e., eq 13 with l ) m ) n ) (1. (b) The reaction energy from the (2nd + 3rd), (4th + 5th), and (6th + 7th) image shells separately and from the sum of the 2nd to 51st image shells.

TABLE 2: The Mean and Variance of the Reaction Energy (in kcal/mol) Computed Using Eq 13 with Various Image Charge Layers, {l, m, n} |i|a

mean

variance

1 2+3 4+5 6+7 8+9 10 + 11 2-51

-323.6 34.2 9.1 4.0 2.2 1.4 54.2

111.0 2.9 0.35 0.11 0.04 0.02 3.4

a The reaction energy was computed with eq 13 with l ) m ) n ) (|i|.

Results Comparison between the Rectangular Image Charge and Poisson Methods. Table 1 gives the reaction field components and energy of charged atoms at specific positions in a 20 Å cube (see the Calculations section). The reaction field components and energies computed using eqs 15 and 13, respectively, with -51 e {l, m, n} e 51 differ from the corresponding values obtained by solving a Poisson equation by less than 6.0% (average 4.3%) and 2.9% (average 1.8%). The slight discrepancy between the two sets of results may be due to the nonperfect rectangular cavity in the Delphi calculations. Overall, the results in Table 1 show that the RIC method can yield reasonably accurate reaction field components and energies. Reaction Energy from a Finite Number of Image Charges. Figure 5 shows the ith image shell reaction energy (eq 13 with {l, m, n} ) (i) of two ions of unit charge at various positions in a cube of length 1 Å. The reaction energy from the first shell of image charges, U1(R) (eq 13 with -1 e {l, m, n} e 1) is large and strongly dependent on the particle positions (Figure 5a and Table 2). The reaction energy corresponding to each successive shell of image charges dramatically decreases and becomes less dependent on the particle positions (Figure 5b and Table 2). When the reaction energies from the 2nd to 51st image shells in eq 13 were summed, the resulting energy was found to be only slightly dependent on the particle positions (Figure 5b and Table 2). This implies that the reaction energy can be

Figure 6. (a) Comparison of the reaction energy U computed using -51 e {l, m, n} e 51 in eq 13 (y axis) with -1 e {l, m, n} e 1 (x axis). (b) Comparison of the force F(R) computed using -51 e {l, m, n} e 51 in eq 14 (y axis) with -1 e {l, m, n} e 1 (x axis).

approximated by the first image shell energy and a positionindependent term. Since the reaction energy U(R) (eq 13) converges quickly with the number of image charge layers surrounding the primary cell, it was computed using eq 13 with a finite (as opposed to infinite) number of image charges, viz., -51 e {l, m, n} e 51, and compared to U1(R) (see the Calculations section). Figure 6a shows that the infinite series in eq 13 can be approximated by summing image charges Qj,l,m,n for {l, m, n} ) -1 to +1 since U(R) can be fitted by a linear function of U1(R):

U(R)/(kcal/mol) ) U1(R)/(kcal/mol) +

54.19(Q/e)2 L/Å + 5.6

(20)

Eq 20 implies that

F(R) ≈ F1(R)

(21)

The validity of eq 21 is verified in Figure 6b, which shows that a plot of F(R) (eq 14 with -51 e {l, m, n} e 51) vs F1(R) (eq 14 with -1 e {l, m, n} e 1) gives a straight line of unit slope and zero intercept. Thus, the electrostatic forces were approximated by F1(R) using eq 14 and -1 e {l, m, n} e 1 in the MD simulations. vdW Boundary Force. Figure 7 shows the boundary force along the x direction, Fx, exerted on a water oxygen in the simulation region (Figure 4) for a 20 Å cube in the xy plane (see eq 16 and the Calculations section). The origin is set equal to the center of the cubic box. When the |x| coordinate of the water oxygen is e10 Å, |Fx| increases slightly, but when it is outside the cubic box |Fx| is a parabolic function of x and reaches a maximum at the implicit solvent boundary; i.e., at 12.8 Å (Figure 7a). When the |y| coordinate of the water oxygen is e10 Å, Fx is independent of y, but when it is outside the cubic box, |Fx| decreases linearly (Figure 7b). In most simulations,

Representing an Infinite Solvent System

J. Phys. Chem. B, Vol. 106, No. 11, 2002 2979

Figure 7. The x-component of the solvent boundary force, Fx, in the xy plane for a simulation cube of length 20 Å. (a) Fx as function of x for various values of y; (b) Fx as a function of y for x ) 12.8 Å.

the center of the water oxygen will be confined to the region |y| < 11 Å (see below and Figure 8), so the force can be approximated by the curve for y ) 0-10 Å in Figure 7a. Since the potential energy at the edge of the cube is about -0.25 kcal/ mol, a Fx equal to 0.05 kcal/mol/Å can be assumed when 5 < |x| < 10 Å. Hence, the force along the x direction can be described as follows:

Fx/(kcal mol-1Å-1) ) 0.0 S x/Å e 5

(22a)

Fx/(kcal mol-1Å-1) ) 0.05 S 5 < x/Å e 10 (22b) Fx/(kcal mol-1Å-1) ) (x/Å - 12.8)2 - 2.82 + 0.05 S 10 < x/Å < 12.8 (22c) In the general case of a cubic box of length L and L/2 > 5 Å, the force should approach zero for all explicit water molecules more than 5 Å from the simulation boundary. For water molecules near the boundary, the behavior of a force component should be similar to that shown in Figure 7, and is independent of the cube size. Because x, y, and z are interchangeable, Fy and Fz can adopt the same form as Fx. Thus the solvent boundary force can generally be given by:

Fx,y,z/(kcal mol-1Å-1) ) 0.0 S x/Å e (L/Å)/2 - 5 (23a) Fx,y,z/(kcal mol-1Å-1) ) 0.05 S (L/Å)/2 - 5 < x/Å e (L/Å)/2 (23b) Fx,y,z/(kcal mol-1Å-1) ) [x/Å - (L/Å)/2 - 2.8]2 - 2.82 + 0.05 S (L/Å)/2 < x/Å < (L/Å)/2 + 2.8 (23c) The simulation results described below were obtained using eq 23 to describe the solvent boundary force on the explicit water molecules.

Figure 8. The water density along the X axis (dotted line), Y axis (dashed line), and Z axis (solid line) for (a) a 10 Å cubic box containing 33 water molecules, and (b) a 20 Å cubic box containing 267 water molecules.

Water Simulations with Image Charges. During the 200 ps simulation of the small 10 Å cubic box containing 33 water molecules, the total, kinetic, and potential energy remained roughly constant with an average (variance) of -237.6 (0.06), 56.8 (4.6), and -294.4 (4.6) kcal/mol, respectively, while the temperature fluctuated around a mean of 297.8 (24.2) K after 20 ps. With the larger 20 Å cubic box containing 267 water molecules, the fluctuations are reduced, as evidenced by a mean temperature of 297.5 (8.9) K, and average values of the total, kinetic, and potential energies of -2025.7 (0.2), 471.8 (14.1), and -2497.5 (14.1) kcal/mol, respectively. Figures 8a and 8b show the water density along the X, Y, and Z axes for the 10 and 20 Å cubic boxes, respectively. The water density is close to its bulk value, but near the boundary it decays from one at about L/2 - 0.75 Å to zero around L/2 + 0.75 Å. The latter means an absence of water molecules within 2 Å of the implicit solvent boundary, hence artificial divergences in the forces and energies are avoided in the simulations with image charges. Ion Simulations with Image Charges. Ion-SolVent Structure and Energetics. Figures 9 and 10 show the ion-water oxygen rdf, running coordination number, and electrostatic potential as a function of distance from K+ and Na+, respectively. They show that the computed quantities are accurate up to half the cube length as expected. The first peak positions of the K+-water oxygen rdf (2.70 ( 0.03 Å) and Na+-water oxygen rdf (2.31 ( 0.03 Å) agree with the respective experimental numbers34 (2.80 ( 0.08 Å for K+, and 2.36 ( 0.06 Å for Na+). The first-shell hydration numbers for K+ (6-8) and Na+ (6) are also in accord with the respective values derived from X-ray and neutron diffraction studies.34 Comparison between RIC and Reference Simulations. Figures 11 and 12 compare the ion-water oxygen rdfs computed using the RIC method for L ) 10, 15, and 20 Å with the respective rdfs obtained using the minimum image convention23

2980 J. Phys. Chem. B, Vol. 106, No. 11, 2002

Yang et al.

Figure 11. K+-water oxygen rdfs as a function of distance from the center of the ion. (a) The dashed, solid and dotted curves correspond to the rdfs computed using the minimum image convention with L ) 40 Å, the RIC method for L ) 20 Å, and spherical boundary conditions with a radius of 20 Å, respectively. (b) The black solid, black dashed, black dotted, grey solid, and grey dotted curves correspond to the rdfs computed using the RIC method with L ) 20, 15, and 10 Å, and spherical boundary conditions with a radius of 12 and 6 Å, respectively. Figure 9. K+-water oxygen rdf (a), running coordination number (b) and electric potential (c) as a function of distance from the center of the ion calculated by the RIC method for L ) 10 Å (dotted line), 15 Å (dashed line), and 20 Å (solid line).

Figure 12. As in Figure 11 for Na+.

Figure 10. As in Figure 10 for Na+.

for L ) 40 Å as well as spherical boundary conditions for a 6, 12, and 20 Å sphere. The ion-water oxygen rdf derived from

a 20 Å cube (267 water molecules) simulation using the RIC method appears indistinguishable from the respective rdf obtained from the much larger 40 Å cube (2137 water molecules) simulation using the minimum image convention or a 20 Å sphere (1118 water molecules) with spherical boundary conditions (see Figures 11a and 12a). Note that the first and second peak of the ion-water oxygen rdf derived from simulations using spherical boundary conditions for a 6 Å sphere (33 water molecules, gray dotted curve) and a 12 Å sphere (267 water molecules, gray dashed curve) differ from those obtained from simulations using a larger number of water molecules (Figures 11b and 12b). In contrast, the first and second peak of the ion-water oxygen rdf computed using the RIC method for the same number of water molecules (33 water molecules, black dotted curve; 267 water molecules, black solid curve in Figures

Representing an Infinite Solvent System 11b and 12b) are similar to those obtained from simulations with a larger number of water molecules. This shows that the image charge potential can improve the rdf near the implicit solvent boundary. Table 3 compares the average electrostatic potential energy of Na+ and K+ computed using eq 20 with corresponding values from periodic boundary simulations with the minimum image convention.23 The Na+ and K+ potential energies computed from RIC simulations of a 20 Å cube containing 267 water molecules are similar to those obtained from PBC simulations of a 40 Å cube with 2137 water molecules. Thus, the RIC method can yield accurate cation-water oxygen radial distribution functions as well as potential energies, but with fewer explicit water molecules compared to PBC simulations. Ion Solvation Free Energy. Table 4 compares the Na+ and + K solvation free energies computed from this work with those obtained using spherical boundary conditions11 and with experimental solvation free energies. The absolute and relative solvation free energies obtained using the RIC method with a cubic cell length of 10, 15 and 20 Å are in close agreement with the corresponding experimental values.35 In particular, they agree with the respective values derived using spherical boundary conditions for the same set of Lennard-Jones ion-water parameters.11 Since the image charges can extend the finite cubic system into an infinite solvent system, the solvation free energy of an ion should be independent of its position in the cubic box. To verify this, the solvation free energy of K+ was evaluated at various grid points in a 10, 15, and 20 Å3 cubic box. The resulting values listed in Table 5 are similar and agree with the value calculated at the origin to within 2.3%.

J. Phys. Chem. B, Vol. 106, No. 11, 2002 2981 TABLE 3: The Average Potential Energy (kcal/mol) of K+ and Na+ from Simulations Using the RIC Method and PBC PE (kcal/mol) method

no. of waters

K+

Na+

RIC RIC PBC PBC

33 267 267 2137

-129.4 -150.3 -125.8 -143.9

-171.0 -192.0 -165.2 -185.8

TABLE 4: Comparison between Experimental and Computed Solvation Free Energies (kcal/mol) method

no. of waters

-∆G(K+)

-∆G(Na+)

∆G(K+) ∆G(Na+)

experimenta RICb,c RICb,d RICb,e SBCf SBCf SBCf

33 113 267 25 50 100

80.6 79.2 ( 0.4 80.9 ( 0.2 81.4 ( 0.3 79.1 80.4 81.5

98.2 99.8 ( 0.2 100.6 ( 0.3 100.9 ( 0.1 100.7 100.8 102.6

17.6 20.6 19.7 19.5 21.6 20.4 21.1

a Experimental values from Burgess, 1978.35 b Parameters from Beglov and Roux, 1994.11 c Computed using image charges in a cubic box of length 10 Å. d Computed using image charges in a cubic box of length 15 Å. eComputed using image charges in a cubic box of length 20 Å. f Net solvation free energy, including charging and cavity free energies, from Beglov and Roux, 1994,11 using spherical boundary conditions.

TABLE 5: Solvation Free Energy of K+ (in kcal/mol) Calculated at Various Positions in a Cubic Box of Length L Using the RIC Method -∆G(K+)

X, Y, Z L ) 10 Å

Concluding Discussion The RIC method was developed as a cost-effective and accurate means of simulating nonspherical, highly charged systems. It has potential use not only in DNA, RNA, or nonglobular protein simulations, but also in folding/unfolding simulations where the denatured state of the macromolecule is often nonspherical. Even with as few as 33 water molecules, the RIC method can yield accurate ion-water oxygen rdfs up to the second hydration shell (Figures 11 and 12) as well as absolute free energies (Table 4). Unlike PBC the RIC method does not induce artificial periodicity nor require charge neutrality when Ewald summation is used to treat Coulombic interactions. In contrast to Friedman’s image charge approximation which assumes that 2 . 1, such an assumption is not needed here (see the Theory section) and thus the RIC method is also suitable for simulations in nonaqueous solution with known dielectric constant. However, MD simulations employing the present RIC method encounter some limitations that are also inherent in simulations employing spherical boundary conditions. As with other methods that treat implicit solvent molecules as a dielectric continuum, there is uncertainty in the dielectric boundary. In this work the distance from the edge of box to the implicit solvent boundary was assumed to be equal to 2.8 Å (Figure 4). As rationalized in the Introduction section, this distance was chosen as it is the distance to the first peak of the water oxygen-water oxygen rdf. In analogy with spherical boundary conditions that apply strictly to a perfect sphere such as the Kirkwood theory22 the RIC method assumes a perfect planar boundary. However, errors stemming from a nonperfect planar boundary appear to be negligible from the results in Table 1 (see the Results section). The current formulation of the RIC method neglects the density and volume fluctuations of bulk water at constant pressure. This

0, 0, 0 0, 0, 2 0, 2, 2 2, 2, 2 L ) 15 Å 0, 0, 0 0, 0, 4 0, 4, 4 4, 4, 4 L ) 20 Å 0, 0, 0 0, 0, 6 0, 6, 6 6, 6, 6

79.2 ( 0.4 78.7 ( 0.3 (0.6%) 78.0 ( 0.1 (1.5%) 77.6 ( 0.1 (2.0%) 80.9 ( 0.2 80.7 ( 0.1 (0.2%) 79.6 ( 0.2 (1.6%) 79.2 ( 0.3 (2.1%) 81.4 ( 0.3 81.2 ( 0.2 (0.2%) 80.7 ( 0.2 (0.9%) 79.5 ( 0.2 (2.3%)

can be taken into account by including the cavity free energy contribution, as described in previous works.11 Acknowledgment. We are grateful to Dr. C. Satheesan Babu for stimulating discussions and many helpful suggestions. We thank Professor Martin Karplus for the CHARMM program. This work was supported by Academia Sinica, the National Science Council (NSC87-2312-B002-030), and the National Center for High-Performance Computing, Taiwan. References and Notes (1) de Leauw, S. W.; Perram, J. W.; Smith, E. R. Proc. Royal Soc. London 1980, 373, 27-56. (2) Still, W. C.; Tempczyk, A.; Hawley, R. C.; Hendrickson, T. J. Am. Chem. Soc. 1990, 112, 6127-6129. (3) Honig, B.; Sharp, K.; Yang, A. J. Phys. Chem. 1993, 97, 11011109. (4) Cramer, C. J.; Truhlar, D. Chem. ReV. 1999, 99, 2161-2200. (5) Babu, C. S.; Lim, C. J. Phys. Chem. B 1999, 103, 7958-7968. (6) Hansen, J. P.; McDonald, I. R. Theory of simple liquids; Academic Press: New York, 1986.

2982 J. Phys. Chem. B, Vol. 106, No. 11, 2002 (7) Berkowitz, M.; McCammon, J. A. Chem. Phys. Lett. 1982, 90, 215217. (8) Brooks, C. L., III; Karplus, M. J. Chem. Phys. 1983, 79, 63126325. (9) Alper, H.; Levy, R. M. J. Chem. Phys. 1993, 99, 9847-9852. (10) Wallqvist, A. Mol. Simul. 1993, 10, 13-17. (11) Beglov, D.; Roux, B. J. Chem. Phys. 1994, 100, 9050-9063. (12) Wang, L.; Hermans, J. J. Phys. Chem. 1995, 99, 12001-12007. (13) Mcconnell, K. J.; Nirmala, R.; Young, M. A.; Ravishanker, G.; Beveridge, D. L. J. Am. Chem. Soc. 1994, 116, 4461-4462. (14) Ravishanker, G.; Auffinger, P.; Langley, D. R.; Jayaram, B.; Young, M. A.; Beveridge, D. L. ReV. Comput. Chem. 1997, 317-372. (15) Lim, C.; Bashford, D.; Karplus, M. J. Phys. Chem. 1991, 95, 56105620. (16) Jayaram, B. J. Phys. Chem. 1994, 98, 5773-5777. (17) Sitkoff, D.; Sharp, K. A.; Honig, B. J. Phys. Chem. 1994, 98, 19781988. (18) Reddy, M. R.; Erion, M. D.; Agarwal, A.; Viswanadhan, V. N.; McDonald, D. Q.; Still, W. C. J. Comput. Chem. 1998, 19, 769-780. (19) King, G.; Warshel, A. J. Chem. Phys. 1989, 91, 3647-3661. (20) Rullman, J. A.; van Duijnen, P. T. Mol. Phys. 1987, 61, 293-311. (21) Friedman, H. L. Mol. Phys. 1975, 29, 1533-1543. (22) Kirkwood, J. G. J. Chem. Phys. 1934, 2, 351-361.

Yang et al. (23) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1990. (24) Jackson, J. D. Classical Electrodynamics; Wiley: New York, 1975. (25) Narten, A. H.; Thiessen, W. E.; Blum, L. Science 1982, 217, 10331034. (26) Gilson, M. K.; Honig, B. H. J. Comput. Chem. 1988, 9, 327-335. (27) Sharp, K. A.; Honig, B. Annu. ReV. Biophys. Biophys. Chem. 1990, 19, 310-332. (28) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187-217. (29) Nicholls, A.; Honig, B. J. Comput. Chem. 1991, 12, 435-445. (30) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926-935. (31) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (32) Beveridge, D. L.; DiCapua, F. M. Free energy Via molecular simulation: A primer ; ESCOM: Leiden, 1989. (33) McQuarrie, D. A. Statistical Mechanics; Harper and Row: New York, 1976. (34) Marcus, Y. Chem. ReV. 1988, 88, 1475-1498. (35) Burgess, M. A. Metal ions in solution; Ellis Horwood: Chichester, England, 1978.