J. Phys. Chem. B 2007, 111, 12263-12274
12263
Implicit Nonpolar Solvent Models Chunhu Tan, Yu-Hong Tan, and Ray Luo* Department of Molecular Biology and Biochemistry, UniVersity of California, IrVine, California 92697-3900 ReceiVed: May 3, 2007; In Final Form: July 24, 2007
We have systematically analyzed a new nonpolar solvent model that separates nonpolar solvation free energy into repulsive and attractive components. Our analysis shows that either molecular surfaces or volumes can be used to correlate with repulsive free energies of tested molecules in explicit solvent with correlation coefficients higher than 0.99. In addition, the attractive free energies in explicit solvent can also be reproduced with the new model with a correlation coefficient higher than 0.999. Given each component optimized, the new nonpolar solvent model is found to reproduce monomer nonpolar solvation free energies in explicit solvent very well. However, the overall accuracy of the nonpolar solvation free energies is lower than that of each component. In the more challenging dimer test cases, the agreement of the new model with explicit solvent is less impressive. Nevertheless, it is found that the new model works reasonably well for reproducing the relative nonpolar free energy landscapes near the global minimum of the dimer complexes.
Introduction Continuum or implicit solvent models are frequently used to calculate solvation properties of organic molecules and biomolecules for their high efficiency. A typical modern implicit solvent model consists of both electrostatic (polar) and nonelectrostatic (nonpolar) components. With a continuum dielectric model, the polar component can be calculated by solving the Poisson-Boltzmann equation.1-4 The nonpolar component of the solvation free energy is often estimated from the solventaccessible surface area (SAS, i.e., the surface traced out by the solvent probe sphere center as it rolls over the molecule) with a uniform surface tension coefficient,5-9
∆Gnp ) γ‚SAS + c
(1)
where the surface tension coefficient γ represents the contribution to the nonpolar solvation free energy (∆Gnp) per unit surface area. The constant offset is the solvation free energy for a point solute (SAS ) 0). This approach is often referred as the SA model. The use of a single γ is based on the observation that for linear alkanes, solvation free energy approximately increases linearly with SAS.10-12 This approach can be traced back to the scaled particle theory originally intended for large solute molecules13,14 but has widely been used for small solute molecules, as well. However, the correlation between solvation free energies computed by the SA solvent model and those computed by the explicit solvent model has been found to be poor for more general nonpolar organic molecules, as shown in Figure 1. The SA model, regardless of its poor accuracy, has been widely used for molecular mechanics simulations until the pioneer works by Levy and co-workers at the turn of this century.15-19 The new theoretical framework is originated from the following observation in liquid-state physical chemistry: ∆Gnp consists of two components of highly different interactionss ∆Grep, repulsive free energy; and ∆Gatt, attractive free energys and they must be modeled separately,15-23
∆Gnp ) ∆Grep + ∆Gatt
(2)
Figure 1. Correlation between implicit nonpolar solvation free energy (∆Gnp_imp) computed by eq 1 as implemented in Amber53 and explicit nonpolar solvation free energies (∆Gnp_exp) for 42 tested small molecules from the Amber force field database.52 See Theory and Computational Details on the tested small molecules and free energy simulations.
where ∆Grep is the solvation free energy from the solute-solvent repulsive interactions and the formation of solute cavity (the excluded volume effect). ∆Gatt is the free energy for establishing the solute-solvent attractive interactions, but may also include solvent-solvent reorganization component. To understand eq 2, one can envision the solvation process of a nonpolar solute with two steps. First, a solute-sized cavity is created in the solvent by turning on solute-solvent repulsive interactions, but with the absence of the attractive interactions between solute and solvent. Second, the attractive interactions between solute and solvent are turned on.24,25 Free energy contributions of the two steps to solvation can be studied independently.26 When decomposed in the way of eq 2, ∆Grep was found to have a very good correlation with SAS.15,27
∆Grep ) γ‚SAS + c
(3)
Use of molecular volumes, for example, solvent-accessible volume (SAV, i.e., the molecular volume enclosed by SAS), has also been proposed to correlate with ∆Grep,28-33
∆Grep ) p‚SAV + c
10.1021/jp073399n CCC: $37.00 © 2007 American Chemical Society Published on Web 10/05/2007
(4)
12264 J. Phys. Chem. B, Vol. 111, No. 42, 2007
Tan et al.
where p is a solvent pressure parameter. Use of both molecular surfaces and volumes has also be tested.33 Free energy simulations of spherical cavities have shown, however, that for small cavities, ∆Grep correlates with the cavity volume, whereas for large cavities, ∆Grep correlates with the cavity surface.14,34,35 The crossover occurs around the cavity radius of 10 Å, as shown in recent studies.30,36 On the other hand, with the decomposition of eq 2, efficient computation of ∆Gatt is a new challenge. Fortunately, according to Chandler and co-workers,24,37 ∆Gatt can be approximated by the van der Waals attractive interaction potential energy between solute (u) and solvent (v):
∆Gatt ≈
or is often simplified as
uLJ(r) )
(5)
Ns
∑ ∫Faw(raw)Vatt(raw) draw a)1
(6)
Here, the sum is over all solute atoms (Ns), and the integration is over the solvent-occupied volume. Faw(raw) is a solvent distribution function around solute atom “a” at a given solutesolvent distance, raw. Vatt(raw) is the attractive van der Waals potential in a decomposition scheme; for example, the WCA scheme.20,38 In implicit solvents, eq 6 has to be further approximated because Faw(raw) cannot be known a priori without equilibrium simulations in explicit solvent. As a first approximation, a uniform distribution (i.e., constant density) can be used. In this study, we intend to address three questions related to the applications of the new nonpolar solvation model (eq 2) to biomolecules. First, how well does the repulsive free energy (often termed cavity free energy) correlate with molecular surface or volume for molecules of biological interest? If a highquality linear correlation can be found between ∆Grep and the molecular surface or volume, the computation of ∆Grep is equivalent to that of molecular surface or volume, which has been thoroughly studied.33,39-51 Second, given eq 6, how do we calculate the attractive interaction energy between solute and solvent in implicit solvent? And finally, how good is the new nonpolar solvation model, that is, how good is eq 2 in reproducing total nonpolar solvation free energies? To answer these questions, we have chosen two test sets to analyze the new approach. The first set consists of 42 molecules derived from the Amber force field database that cover all templates/ residues for proteins and nucleic acids.52 The second set consists of two groups of C-G base pairs, derived from hydrogenbonded and stacked conformations in B-DNA. For ease of analysis, the assessment of accuracy in this study is based on simulated nonpolar solvation free energies in explicit solvent, and the TIP3P explicit solvent is used.
Vrep(r) )
1. Decomposition Schemes. In most molecular mechanics force fields, such as Amber,53 CHARMM,54 and OPLS,55 the classical Lennard-Jones 6-12 potential is often used to model nonpolar interactions in molecular mechanics:
[(σr ) - (σr ) ]
uLJ(r) ) 4
12
A r12
(9)
B r6
(10)
Vatt(r) ) -
Another well-known decomposition scheme is the WeeksChandler-Andersen (WCA) scheme:20,38
Vrep(r) )
{
Vatt(r) )
uLJ(r) + r ermin r > rmin 0
(11)
{
(12)
- r e rmin uLJ(r) r > rmin
where rmin ) 21/6σ. In the present study, a new decomposition scheme, termed the σ scheme, is proposed to improve the overall accuracy in modeling nonpolar solvation free energies with continuum solvents.
Vrep(r) ) Vatt(r) )
{ {
uLJ(r) r eσ 0 r>σ
(13)
0 r eσ uLJ(r) r > σ
(14)
The σ decomposition scheme is similar to the BarkerHenderson theory, in which the potential is split on the basis of the sign of energy, whereas in the WCA scheme, the potential is split on the basis of the sign of force.56,57 Note that the σ decomposition scheme defined by eqs 13 and 14 is continuous but not smooth (i.e., without continuous first derivative) at r ) σ. This may cause the divergence theorem (to be applied below) to be inapplicable. To overcome the limitation, eqs 13 and 14 can be redefined in the following way: For a given ξ, a polynomial is defined for each of Vrep and Vatt so that each term smoothly goes to zero from uLJ in the range [σ - ξ, σ + ξ] (the polynomial and its first derivative are continuous at both ends). For example, the following polynomial can be used for Vrep/att within [σ - ξ, σ + ξ]:
Vrep/att ) C0 + C1/r3 + C2/r6 + C3/r12
Theory and Computational Details
6
(7)
(8)
with A ) 4σ12 and B ) 4σ.6 Here, and σ are the well-depth and radius parameters of the Lennard-Jones potential, respectively. These parameters are the van der Waals coefficients as specified in a force field. The Lennard-Jones potential is typically decomposed into the repulsive and attractive regions according to a decomposition scheme. The most straightforward decomposition scheme is the 6/12 scheme, in which18,21,22
This has subsequently been confirmed in simulations by Levy and co-workers.15 In general, the solute-solvent van der Waals interaction energy can be analytically expressed as the following volume integral:
Uuv att )
B A r12 r6
(15)
For Vrep, the coefficients are
C0 )
1 [B(p6 + q6 + p3q3) - 3A] T
(16)
C1 )
2 [2A - B(p6 + q6)](p3 + q3) T
(17)
Implicit Nonpolar Solvent Models
C2 )
J. Phys. Chem. B, Vol. 111, No. 42, 2007 12265
q3 [B(q9 + p6q3 + p3q6 + 3p9) - 6Ap3] T (18)
C3 ) -
q9 [A(q3 - 2p3) + Bp9] T
(19)
For Vatt, the coefficients are
C0 )
1 [3A - B(p6 + q6 + p3q3)] T
(20)
C1 )
2 [B(p6 + q6) - 2A](p3 + q3) T
(21)
C2 )
p3 [6Aq3 - B(p9 + p6q3 + p3q6 + 3q9)] T (22)
C3 )
p9 [A(p3 - 2q3) + Bq9] T
(23)
Here, A and B are the van der Waals parameters, p ) σ - ξ, q ) σ + ξ, and
T ) [(p9 - q9) - 3p3q3(p3 - q3)](p3 + q3)
(24)
Note that the sum of Vrep and Vatt so defined is equal to uLJ, as in the other two decomposition schemes. The performance of all three decomposition schemes will be presented in Results and Discussion. 2. Repulsive Free Energy and Forces. In the new nonpolar solvent model, computation of repulsive free energy is straightforward because a linear relation with either molecular surface or volume can be found. To further quantify the correlations between repulsive free energy and molecular surface/volume, we have performed extensive thermodynamic integration (TI) simulations for 42 training molecules from the Amber force field database.52,53 All repulsive free energies are calculated by the difference between total nonpolar free energies from TI and attractive free energies as approximated by eq 5. Once repulsive free energy can be modeled as linearly proportional to either molecular surface or volume, repulsive forces can be computed given derivatives of molecular surface/ volume as documented in the literature.33,44,46,48,49,58-63 3. Attractive Free Energy. As reviewed in the Introduction, the free energy contribution of adding the attractive interaction upon solvation may be approximated by the attractive interaction energy between the solute and the solvent,15,24,30,37 which can be calculated by eq 6. To simplify the calculation, Faw(raw) can be expressed in terms of the bulk water density, Fw, and a correlation function, g(raw), as
Faw(raw) ) Fwg(raw)
(25)
Figure 2. Geometrical relationships used in the derivation of nonpolar attractive free energy and force. Dab is the vector pointing from the center of atom a to the center of atom b. Ps is the position of area element dσs with respect to the center of atom b. ras is the vector pointing from the center of a to the area element dσs, and nˆ s is the outward normal vector associated with dσs.
However, it is necessary to avoid having to incorporate more complicated and expensive integral equation theories into the models.64,65 Floris and Tomasi21 have shown that under the continuum approximation, the volume integral in eq 6 can be transformed into an integral over the solute-accessible surface with a simple application of the divergence theorem, Ns
Uuv att ) Fw
{
Here Cs is the solute-occupied volume. This corresponds to the so-called “homogeneous or continuum approximation”: the average solvent number density outside the solute-occupied volume is constant. Clearly, this is a potentially severe approximation and is questionable for a condensed-phase system.
b
Θatt(ras)ras‚nˆ s dσs
(27)
where indexes a and b run over all solute atoms. The surface integration is over the SAS of atom b (Sb). Θatt(ras) is a function defined on the SAS of atom b, and nˆ s is the outward normal vector associated with the SAS element, dσs (see Figure 2). In the case of Lennard-Jones potential with the 6/12 decomposition scheme, as shown by eqs 9 and 10, the corresponding Θatt in eq 27 can be shown to be21,23,66
Θatt ) -
Baw
(28)
3|ras|6
For the WCA decomposition scheme (eq 11 and 12),
Θatt )
Θatt )
[
Aaw Baw 1 (ras3 - rmin3) + 3 9 3ras 3rmin rmin3 Aaw 9|ras|
12
-
Baw
]
if ras e rmin
if ras > rmin
3|ras|6
(29) (30)
Similarly, for the σ decomposition scheme (eq 13 and 14),
Θatt ) Θatt )
(26)
∑ ∑ ∫S
a)1 b)1
The first-order approximation of g(raw) is a step function:21,22
0 for raw ∈Cs g(raw) ) 1 for raw ∉Cs
Ns
[
]
Aaw Baw 1 - 3 3 3|ras| 3σ9 σ Aaw 9|ras|
12
-
Baw 3|ras|6
if ras eσ
if ras > σ
(31)
(32)
Here, Aaw and Baw are the van der Waal coefficients in eq 8 for interactions between solute atoms (a) and water oxygen atoms (w). Formulation of attractive forces based on eq 27 for numerical calculations and analysis of the numerical procedure can be found in the Appendix.
12266 J. Phys. Chem. B, Vol. 111, No. 42, 2007
Tan et al.
4. Simulation Details. Explicit SolVent Simulations. All molecules in the present study were solvated by a TIP3P water box with a buffer of 11 Å in XLEAP of Amber 9.53 The system temperature was coupled to 300 K, and the system pressure was coupled to 1 bar by the Berendsen’s coupling algorithms.67 SHAKE was turned on for bonds containing hydrogen atoms68 so that a time step of 2fs was used in all molecular dynamics simulations. Particle Mesh Ewald with a real space cutoff of 9 Å was used for nonbonded interactions.69 van der Waals interactions outside the cutoff distance were computed on the basis of the continuum approximation.53 For easy comparison with implicit solvent calculations that were performed with rigid conformations, all solute atoms were restrained to their initial positions with a harmonic potential with a force constant of 50 kcal/mol-Å2. Before free energy simulations were started as discussed below, all systems were fully relaxed with the PMEMD program in Amber 953 for at least 4 ns or until there was no systematic drift in the running average potential energies. Total nonpolar solvation free energies were calculated by performing the thermodynamic integration method both in explicit solvent and in vacuum,70
∆G )
∫0
1
〈 〉
dH(λ) dλ dλ λ
(33)
where H is a parametrized Hamiltonian, and λ is the coupling coefficient. In practice, the above integral is estimated by a numerical method in which the integrand is determined at a set of discrete λ’s (or “windows”). The SANDER program in Amber 853 was used to compute the integrands. Because only the nonpolar free energies were of interest, all atoms of a studied molecule were uncharged in the initial state of the system. The van der Waals interaction between the solute and the solvent was decoupled by increasing λ from 0 to 1 in steps of ∆λ. A total of 60 λ values (from 0 to 0.5 in steps of 0.025 and from 0.5 to 1 in steps of 0.0125) were used. To handle the integrand divergence problem for “dummy” atoms at λ ) 1 (where atoms “disappear” in the perturbed state), KLAMBDA is set to 6 in SANDER to use the nonlinear mixing rule of potential function,71
V(λ) ) (1 - λ)kV0 + [1 - (1 - λ)k]V1
(34)
where V0 is the potential with the original Hamiltonian, and V1 is the potential with the perturbed Hamiltonian. In the final state, all atoms become “dummies”, which have neither electrostatic nor van der Waals terms. For each intermediate state (or each “window”), 40 ps of equilibration was run, followed by another 40 ps of data collection. To guarantee convergence, both equilibration and production simulation times were doubled until the free energy difference was below the cutoff of 0.25 kcal/ mol. Figure 3 shows an example convergence behavior in the TI calculations: the nonpolar solvation free energy profiles for both hydrogen-bonded and stacked C-G base pairs. Attractive solvation free energies were approximated by the averages of attractive interaction energy between the solute and solvent as shown by eq 5, which were calculated with the saved MD equilibration trajectories. The repulsive free energies were then obtained by the difference between the total nonpolar free energies and the attractive free energies. The standard errors were used as statistical uncertainties for simulated free energies.72 For each window in TI, the standard error of 〈dH/dλ〉 is δ〈dH(λ)/dλ〉λ ) (σ/xNe), where σ is the standard deviation of 〈dH/dλ〉, Ne is the “independent” sampling
Figure 3. Convergence of free energy simulations in explicit solvent with respect to the length of each TI window. The free energy profiles are for two sets of C-G base pairs, including (A) hydrogen-bonded and (B) stacked. See Figure 4 for conformations used.
number, which can be calculated by the correlation time, τ, of 〈dH/dλ〉 and the data collection time Tdc: Ne ) Tdc/(2τ).72 Usually, for the molecules studied here, τ < 0.3ps. Implicit SolVent Simulations. As described above, repulsive free energy is modeled as a term linearly proportional to SAS/ SAV. Solvent-excluded surface (SES, i.e., the surface traced by the inward-facing surface of the solvent probe sphere as it rolls over the molecule) and solvent-excluded volume (SEV, the molecular volume enclosed by SES) are also investigated. In this study, SAS, SAV, and SEV were calculated with PBSA73 in Amber 9.53 In the PB procedure, the SAS area is computed numerically by representing each atomic surface by spherically distributed dots.74 These dots are first checked against covalently bonded atoms to see whether any of the dots are buried. The exposed dots from the first step are then checked against a nonbonded pair list with a cutoff distance of 9 Å to see whether any of them are buried. The exposed dots of each atom after the second step then represent the solvent-accessible portion of the atom surface and are used to compute the SAS area of the atom. The molecular SAS area is simply a summation of the atomic SAS areas. In the PB procedure, SEV was numerically represented by a cubic grid enclosed by SES. The SEV grid is generated in three steps. First, all grid points within SAV are marked. This limits further checking to be within SAV only. Second, all grid points within the van der Waals volume (i.e., the union of all atomic volumes of the molecule) are marked as within SEV. These grid points constitute the contact region of the SEV. In the meantime, the grid points outside the van der Waals volume but within the contact region of SAV are marked as outside of the SEV. Thus, the only grid points remaining to be checked are the reentry region of SAV. These grid points are checked in the third step: the grid points within the reentry region of SAV are checked to see whether they are within any solvent probe located on the solvent-accessible arcs. If a grid point passes the checking, it belongs to the reentry region of SEV. The solvent-accessible arcs are numerically represented with a resolution of half-grid spacing. Here a grid spacing of 0.25 Å was used in this study. The SEV grid would be used for dielectric assignment in a typical PB calculation.73 SES was calculated with the MSMS program.45 Atomic van der Waals radii (Rmin) in the ff99 Amber force field were used. This will be further discussed below.
Implicit Nonpolar Solvent Models
Figure 4. Two sets of tested C-G base pairs: (A) hydrogen-bonded conformations and (B) stacked conformations. For each set, a total of 11 different conformations (of different pairwise distances) are used. These conformations are generated by keeping G (shown as licorice) fixed while moving C (shown with line) to the right in A and upward in B.
The attractive free energy and forces were computed numerically as a surface integration over the SAS, as described in Appendix. Thus, the accuracy of both attractive free energy and forces depends on the numerical resolution of the SAS. This is also discussed in the Appendix. Tested Molecules. Two sets of testing systems were considered. The first set consists of 42 building blocks for proteins and nucleic acids in the Amber force field database.52,53 The second set consists of two groups of C-G base pairs, initially in hydrogen-bonded and stacked conformations. These initial conformations were generated by NUCGEN of Amber 9.53 Subsequent modified conformations were generated by keeping G fixed in space while moving C to the right in hydrogenbonded conformations and upward in stacked conformations. All base pair conformations are shown in Figure 4. Results and Discussion 1. Optimization of Repulsive Free Energy. For each decomposition scheme, we have systematically explored all existing approaches to achieve an optimal agreement between explicit and implicit solvents. As reviewed in the Introduction, both surface and volume estimators have been proposed to correlate with the repulsive free energies. Here, four commonly used molecular surfaces and volumes,s SES, SAS, SEV, and SAVsare analyzed. In addition, two sets of van der Waals radii, σ or Rmin, can be used in molecular surface or volume calculation. To find the optimal correlation, we have systematically searched the water probe radius (from 0.2 to 1.6 Å with a resolution of 0.02 Å) for each of the eight combinations of estimators and radius sets with respect to the repulsive free energies of the 42 small molecules in the monomer test set. It was found that the Rmin radius set yields slightly better correlations with repulsive free energies for all tested estimators and was used. Note also that the correlations only weakly depend on the water probe radii in the monomer test set; that is, the correlation coefficients change only at the fourth digit after the decimal point. However, the water probe does play an important role in the following dimer test set (see Section 4, below). For example, in the σ decomposition scheme, the best correlation between repulsive free energies and SAV is 0.9979 for the 42 small molecules with a water probe radius of 0.94 Å, and the correlation is 0.9977 with the water probe radius of 1.30 Å
J. Phys. Chem. B, Vol. 111, No. 42, 2007 12267
Figure 5. Correlations of repulsive free energies for the 42 tested small molecules between explicit (∆Grep_exp) and implicit (∆Grep_imp) solvents with the four tested estimators (the 6/12 decomposition scheme). Corresponding parameters can be found in Table 1. Standard errors are plotted as error bars. Note that the error bars are too small to see. (A) SES: CC, 0.987; rmsd, 1.07 kcal/mol; rms relative deviation, 0.054. (B) SEV: CC, 0.933; rmsd, 2.41 kcal/mol; rms relative deviation, 0.123. (C) SAS: CC, 0.986; rmsd, 1.13 kcal/mol; rms relative deviation, 0.058. (D) SAV: CC, 0.987; rmsd, 1.06 kcal/mol; rms relative deviation, 0.054.
Figure 6. Same as Figure 5, but for the WCA decomposition scheme. Corresponding parameters can be found in Table 2. (A) SES: CC, 0.996; rmsd, 0.39 kcal/mol; rms relative deviation, 0.032. (B) SEV: CC, 0.985; rmsd, 0.75 kcal/mol; rms relative deviation, 0.080. (C) SAS: CC, 0.996; rmsd, 0.40 kcal/mol; rms relative deviation, 0.033. (D) SAV: CC, 0.996; rmsd, 0.38 kcal/mol; rms relative deviation, 0.029.
optimized in the dimer test set. Thus, only correlation data with the water probe radii optimized from the dimer test set are shown here. Figures 5-7 show the correlations between the four tested estimators and repulsive free energies for all three decomposition schemes, respectively. Interestingly, for both the WCA and σ decomposition schemes, all four estimators work very well (SEV works less well than the other three) in correlation with the explicit repulsive free energies in the monomer test set. However, the correlations are all lower for the 6/12 decomposition scheme. All corresponding surface tensions, pressures, and constant offsets can be found in Tables 1-3. Fitting quality measures can be found in Figures 5-7. 2. Optimization of Attractive Free Energy. Similarly, there are two choices of radius sets (σ or Rmin) for SAS with which the attractive free energies were computed. It was found that both radius sets yield attractive free energies in very high correlation with those in explicit solvent for the 42 small molecules in the monomer test set: the correlation coefficients are higher than 0.999 with both the WCA and σ decomposition schemes, and the correlation coefficient is higher than 0.99 with the 6/12 decomposition scheme. To be consistent with the repulsive component, the Rmin radius set is used. In addition, two more parameters in the attractive component are still needed
12268 J. Phys. Chem. B, Vol. 111, No. 42, 2007
Tan et al.
Figure 7. Same as Figure 5, but for the σ decomposition scheme. Corresponding parameters can be found in Table 3. (A) SES: CC, 0.997; rmsd, 0.30 kcal/mol; rms relative deviation, 0.026. (B) SEV: CC, 0.987; rmsd, 0.69 kcal/mol; rms relative deviation, 0.082. (C) SAS: CC, 0.997; rmsd, 0.30 kcal/mol; rms relative deviation, 0.026. (D) SAV: CC, 0.998; rmsd, 0.27 kcal/mol; rms relative deviation, 0.022.
TABLE 1: Surface Tensions (γ), Pressures (p), Constant Offsets (c), and Optimized Water Probe Radii (probe) for Repulsive Free Energies Using the Four Tested Estimators (the 6/12 Decomposition Scheme) estimator (probe, Å)
γ or p (kcal/mol-Å2 or 3)
c (kcal/mol)
SES (1.70) SEV (0.20) SAS (0.50) SAV (1.80)
0.1730 0.1421 0.1374 0.0480
-0.9332 4.6277 -3.1479 -3.2655
TABLE 2: Same as Table 1, but for the WCA Decomposition Scheme estimator (probe, Å)
γ or p (kcal/mol-Å2 or 3)
c (kcal/mol)
SES (1.30) SEV (1.60) SAS (0.31) SAV (1.37)
0.1116 0.0990 0.0941 0.0387
-0.2794 2.8340 -1.0268 -0.7834
TABLE 3: Same as Table 1, but for the σ Decomposition Scheme estimator (probe, Å)
γ or p (kcal/mol-Å2 or 3)
c (kcal/mol)
SES (1.30) SEV (1.80) SAS (0.28) SAV (1.30)
0.1053 0.0928 0.0894 0.0378
-0.3146 2.6363 -0.8824 -0.5692
to achieve an optimal agreement with explicit solvent: the water probe radius and the effective water number density (Few). For the σ decomposition scheme, the water probe radius was found to be 0.557 Å to remove the intercept in the linear regression between explicit and implicit attractive free energies. Few was found to be 1.128 68 Fw, where Fw ) 0.033 33/Å3 is the bulk water number density at 300 K and 1 bar to achieve a slope of 1.0 in the above linear regression. For the WCA scheme, the water probe was found to be 0.406 to remove the intercept and achieve a slope of 1.0 at the same time. However, for the 6/12 decomposition scheme, the intercept cannot be removed, no matter which water probe radius was used. Thus, the probe radius is set to 1.06 Å, with which the highest correlation coefficient can be achieved. Finally, the attractive free energies in explicit solvent can be reproduced by the implicit solvent with very high correlations for all three decomposition schemes, as shown in Figure 8 and Table 4. 3. Overall Performance: Monomer Nonpolar Solvation Free Energies. Given each of the repulsive and attractive components optimized with respect to explicit solvent, we go ahead to investigate how well nonpolar solvation free energies
Figure 8. Correlation of attractive free energies for the 42 tested small molecules between explicit (∆Gatt_exp) and implicit (∆Gatt_imp) solvents. Corresponding parameters can be found in Table 4. Standard errors are plotted as error bars. Note that the error bars are too small to see. (A) 6/12 decomposition scheme: CC, 0.9989; rmsd, 0.36 kcal/mol; rms relative deviation, 0.018. (B) WCA decomposition scheme: CC, 0.9997; rmsd, 0.13 kcal/mol; rms relative deviation, 0.008. (C) σ decomposition scheme: CC, 0.9995; rmsd, 0.16 kcal/mol; rms relative deviation, 0.010.
TABLE 4: Optimized Parameters of Attractive Free Energies for the Three Tested Decomposition Schemesa scheme 6/12 WCA σ
probe (Å)
Few
intercept
1.060 0.406 0.557
1.472 1.000 1.129
1.071 0.000 0.000
a Effective water density, Few, in the unit of bulk water density, Fw (0.033 33 Å-3), fitting intercept in kcal/mol.
for the 42 small molecules in the monomer test set can be reproduced with the three decomposition schemes. Figures 9-11 show the correlations of total nonpolar solvation free energies between explicit and implicit solvents with the three decomposition schemes, respectively. All parameters can be found in Tables 1-3. With the WCA and σ decomposition schemes, SES, SAS, and SAV perform very well, but SEV generates significant discrepancy between explicit and implicit solvents in the monomer test set. With the 6/12 decomposition scheme, none of the estimators works as well as those in the WCA and σ decomposition schemes, although each component (repulsive and attractive) is reasonably well reproduced after optimization, as shown in previous sections. It should be pointed out that the σ decomposition scheme works the best in the monomer test set. The reason for its best performance turns out to be quite simple. Note that the total nonpolar free energies are very close to zero for the tested molecules, whereas the magnitudes of both repulsive and attractive components are large but with opposite signs (Figures 5-8). Thus, the accuracy of total nonpolar free energies is in part determined by how well the errors in the respective calculations cancel one another. With the σ decomposition scheme, the average absolute value of each component is decreased by nearly 60% from that with the 6/12 scheme
Implicit Nonpolar Solvent Models
Figure 9. Correlations of total nonpolar solvation free energies for the 42 tested small molecules between explicit (∆Gnp_exp) and implicit (∆Gnp_imp) solvents with the four tested estimators (the 6/12 decomposition scheme). Standard errors are plotted as error bars. Note that the error bars for ∆Gnp_imp are too small to see. (A) SES: CC, 0.889; rmsd, 1.36 kcal/mol; UAVG, 1.09 kcal/mol. (B) SEV: CC, 0.728; rmsd, 2.54 kcal/mol; UAVG, 1.79 kcal/mol. (C) SAS: CC, 0.886; rmsd, 1.42 kcal/ mol; UAVG, 1.14 kcal/mol. (D) SAV: CC, 0.883; rmsd, 1.35 kcal/mol; UAVG, 1.07 kcal/mol.
Figure 10. Same as Figure 9, but for the WCA decomposition scheme. (A) SES: CC, 0.980; rmsd, 0.38 kcal/mol; UAVG, 0.31 kcal/mol. (B) SEV: CC, 0.881; rmsd, 0.74 kcal/mol; UAVG, 0.54 kcal/mol. (C) SAS: CC, 0.981; rmsd, 0.38 kcal/mol; UAVG, 0.31 kcal/mol. (D) SAV: CC, 0.982; rmsd, 0.35 kcal/mol; UAVG, 0.27 kcal/mol.
Figure 11. Same as Figure 9, but for the σ decomposition scheme. (A) SES: CC, 0.981; rmsd, 0.33 kcal/mol; UAVG, 0.26 kcal/mol. (B) SEV: CC, 0.891; rmsd, 0.67 kcal/mol; UAVG, 0.50 kcal/mol. (C) SAS: CC, 0.984; rmsd, 0.31 kcal/mol; UAVG, 0.25 kcal/mol. (D) SAV: CC, 0.986; rmsd, 0.28 kcal/mol; UAVG, 0.21 kcal/mol.
and by ∼10% from that with the WCA scheme. Thus, the σ scheme behaves the best in this regard. It is interesting to compare the parameters optimized here with those used by Baker and co-workers,33 since both their study and this study are intended for the Amber force field in
J. Phys. Chem. B, Vol. 111, No. 42, 2007 12269 the TIP3P solvent. Nevertheless, the comparison can only be qualitative, since different radius sets, different training sets, and different benchmarks (free energies versus mean forces) are used. Additionally, in the formulation of Baker and co-workers, the repulsive component was estimated with both SAS and SAV but without the constant offset.33 For the repulsive component, a 1.14 Å water probe was used for the WCA scheme by Baker and co-workers.33 In this study, the water probes are uniformly larger with SAV: 1.30 Å for the σ scheme and 1.37 Å for the WCA scheme (Tables 1-3). Interestingly, the water probes are uniformly smaller with SAS: 0.28 Å for the σ scheme and 0.31 Å for the WCA scheme (Tables 1-3); however, the smaller water probe radii are determined only by the dimer test sets in this study. The monomer test sets would perform equally excellently with the water probes found by Baker and co-workers.33 For the attractive component, 0.80 Å was used for the WCA scheme by Baker and co-workers, but here, 0.557 Å is used for the σ scheme and 0.406 Å is used for the WCA scheme (Table 4). Nevertheless, the value reported here is within the 99% confidence interval of the fitted water probe radius reported by Baker and co-workers,33 after considering the difference in radius sets used (Rmin versus σ). 4. Overall Performance: Dimer Nonpolar Potential of Mean Force Profiles. It turns out that a more challenging test set for the new nonpolar solvation scheme is the dimer nonpolar potential of mean force (PMF) profiles. As described in Theory and Computational Details, two groups of C-G base pair conformations were used in the validation: hydrogen-bonded and stacked, as found in B-DNA. For each conformation in Figure 4, the nonpolar solvation free energy for “annihilating” base C was first calculated with TI. The sum of nonpolar solvation free energy and the van der Waals interaction energy between bases C and G for each conformation was then used to construct the relative PMF profiles between bases C and G. Here, the reference state was set as the complex conformation in each PMF profile. Even if three estimators with both the WCA and σ decomposition schemes work very well in the monomer test set, none of them can reproduce the base-pair association energetics consistently, no matter which water probe radius is used for the repulsive component. (Fortunately, the performance of the nonpolar solvation scheme in the monomer test set is highly insensitive to the water probe radii used. Thus, we can optimize its performance primarily on the basis of the dimer test set.) Specifically, the new nonpolar solvation scheme mostly fails to reproduce the nonpolar solvation free energy differences (∆∆Gnp) between the free energy changes of fully separated conformations and those of complex conformations, no matter which water probe radius is used. Tables 5 and 6 show two versions of ∆∆Gnp based upon base pair association calculated by implicit solvent: in Table 5, the dimer repulsive free energy is computed with one offset; in Table 6, it is computed with two offsets. Here, the implicit ∆∆Gnp are calculated with the σ decomposition scheme. (Results from the WCA decomposition scheme are very similar (data not shown). Not surprisingly, the results from the 6/12 decomposition scheme are much worse than those of the WCA and σ decomposition schemes.) Except SAS in Table 6 with two offsets, all ∆∆Gnp in implicit solvent are too positive when compared with those in explicit solvent. Bearing this limitation in mind, we have decided to optimize the water probe radius so that the rms deviations between the explicit and implicit nonpolar PMF profiles for both the hydrogen-bonded and stacked dimers are minimized. In addition,
12270 J. Phys. Chem. B, Vol. 111, No. 42, 2007
Tan et al.
TABLE 5: ∆∆Gnp upon Binding Calculated by Explicit Solvent and Implicit Solvent with the Four Tested Estimators (SES/SEV/SAS/SAV) for Repulsive Free Energiesa explicit SES (1.80 Å) SEV (0.20 Å) SAS (1.80 Å) SAV (1.80 Å)
H-bonded (kcal/mol)
stacked (kcal/mol)
-2.84 0.92 -0.60 -0.29 0.63
-2.39 0.75 2.13 -0.62 1.29
a For each estimator, the water probe radius is scanned from 0.20 to 1.80 Å to find the best agreement of ∆∆Gnp between implicit and explicit solvent. The optimized water probe radii are shown in parentheses. The corresponding surface tension (γ) or pressure (p) and the constant offset (c) are as follows: SES: γ ) 0.106 kcal/mol/Å2, c ) -0.346 kcal/mol. SEV: p ) 0.086 kcal/mol/Å3, c ) -0.174 kcal/ mol. SAS: γ ) 0.061 kcal/mol/Å2, c ) -6.407 kcal/mol. SAV: p ) 0.029 kcal/mol/Å3, c ) -1.764 kcal/mol.
Figure 13. Same as Figure 12, but for the WCA decomposition scheme. The minimized rmsd’s are 0.94 kcal/mol for SES, 6.22 kcal/ mol for SEV, 1.32 kcal/mol for SAS, and 1.39 kcal/mol for SAV.
TABLE 6: Same as Table 5, but Using Two Offsets for the Dimer Repulsive Free Energy Calculationa explicit SES (1.80 Å) SEV (0.20 Å) SAS (0.80 Å) SAV (1.80 Å)
H-bonded (kcal/mol)
stacked (kcal/mol)
-2.84 0.58 2.50 -2.51 -1.14
-2.39 0.40 5.23 -2.29 -0.47
a The corresponding surface tension (γ) or pressure (p) and the constant offset are the same as for Table 5 (because the same optimal water probe radii can be used), except for SAS: γ ) 0.078 kcal /mol/ Å2, c ) -2.929 kcal/mol.
only eight data points (covering distances from the closest dimer separation to the second minimum) were used. Finally, only relative free energies were considered, and the global minimum was used as the reference to compute the relative free energies. Figures 12-14 show the comparisons between explicit and implicit solvents for the three decomposition schemes. Consistent with the observation in the monomer test set, the 6/12 decomposition scheme is still the worst of the three after
Figure 12. Relative nonpolar potential of mean force profiles for both hydrogen-bonded and stacked C-G base pairs in explicit and implicit solvents with the four tested estimators (the 6/12 decomposition scheme). Standard errors are plotted as error bars. Note that the error bars for implicit solvents are too small to see. (A) SES and SEV on the hydrogen-bonded conformations. (B) SAS and SAV on the hydrogen-bonded conformations. (C) SES and SEV on the stacked conformations. (D) SAS and SAV on stacked conformations. The water probe radius is optimized by minimizing the rmsd between explicit and implicit PMF for both sets of conformations (see text for more information). The minimized rmsd’s are 4.18 kcal/mol for SES, 13.91 kcal/mol for SEV, 3.85 kcal/mol for SAS, and 1.49 kcal/mol for SAV. The optimized water probe radii, surface tensions, pressures, and constant offsets can be found in Table 1. See Figure 4 for tested conformations.
Figure 14. Same as Figure 12, but for the σ decomposition scheme. The minimized rmsd’s are 1.22 kcal/mol for SES, 6.11 kcal/mol for SEV, 1.35 kcal/mol for SAS, and 1.34 kcal/mol for SAV.
optimization. Similarly, in both the WCA and σ decomposition schemes, SEV does not perform well for the dimer test set, no matter which water probe radius is used (within the range of 0.20-1.80 Å). However, the three other estimators in both the WCA and σ decomposition schemes can reproduce the dimer nonpolar PMF with a certain degree of success, as described below. SES performs reasonably well for the hydrogen-bonded dimer: the new scheme can reproduce both the height and location of the barrier (all relative free energies with respect to the global minimum) (see Figures 13A, 14A). However, it fails to reproduce the barrier height for the stacked dimer, even if it can capture the barrier location (see Figures 13C, 14C). Unfortunately, no second minimum presents in any implicit nonpolar PMF profile. Its behavior is much like PB in our previous studies of the polar PMF profiles of the same base pairs: only the major barriers and the overall trends can be reproduced, but detailed features due to explicit solvent structure are absent.52 SAS moves the barrier locations a bit closer to the global minima, although both barrier heights are reproduced with acceptable accuracy (Figures 13B, 14B and 13D, 14D). Similar to SES, the overall PMF trends can be reproduced even if no second minimum presents. SAV can reproduce the barrier locations better than SAS, but fail to reproduce the barrier heights (see Figures 13B, 14B and 13D, 14D). Similarly to both SES and SAS, no second minimum presents, even if the overall trend of the PMF profiles can be reproduced.
Implicit Nonpolar Solvent Models Concluding Remarks We have systematically analyzed a new approach to model nonpolar solvation for biomolecules. The new method proposes to model total nonpolar solvation free energy by two separate components: repulsive and attractive. Our analysis shows that all tested molecular surfaces (SES/SAS) and volumes (SEV/ SAV) work very well in correlation with explicit repulsive free energies for the tested molecules in the monomer set (with correlation coefficient 0.99 or higher). In addition, the explicit attractive free energies (mostly van der Waals dispersion free energy) can also be reproduced by the new approach with a correlation coefficient higher than 0.999. Surface- or Volume-Dependence for Repulsive Component? A surprising conclusion from our current analysis is that both molecular surfaces and volumes can be used as estimators of repulsive (cavity) solvation free energies with very similar high accuracies, even if the tested monomer molecules are all within the previously reported switching region (spherical radii around 10 Å) from volume-dependence to surface-dependence.30,36 A probable reason for the apparent discrepancy between our conclusion and the literature is that the simulations are on “realistic” molecules in this study, whereas the simulations were performed on ideal spherical cavities in previous studies. Excellent performance in the Monomer Test Set and Decomposition Schemes. Given each component optimized with respect to explicit solvent, the new nonpolar solvent model is tested in reproducing on monomer nonpolar solvation free energies. Except SEV as the estimator for repulsive free energy, SES, SAS, and SAV all work very well for the tested molecules. However, the overall accuracy of the total nonpolar solvation free energies is lower than that of each component, due to apparent imbalances in the cancellation of error between the two components. In addition, it seems that the decomposition scheme of solvent/solute van der Waals potential also plays a role in the overall accuracy of the new approach. It is found that the σ decomposition scheme works better than the WCA and 6/12 decomposition schemes. Difficulties in the Dimer Test Set. A more challenging test for the new nonpolar solvent model is the dimer nonpolar potential of mean force. We have used two sets of C-G base pair conformations as test cases. Unfortunately, none of the estimators can reproduce the base-pair association energetics consistently, even if three estimators work very well in the monomer test cases. Specifically, the new nonpolar solvent model mostly fails in reproducing nonpolar solvation free energy differences (∆∆Gnp) between fully separated conformations and complex conformations: ∆∆Gnp with most estimators are far more positive than those in explicit solvent. Nevertheless, it is found that three estimators work reasonably well for reproducing the free energy landscapes near the global minimum of the tested dimer complexes when both the WCA and σ decomposition schemes are used. Specifically, SES performs reasonably well for the hydrogen-bonded dimer: the new model can reproduce both the height and location of the barrier. However, it fails to reproduce the barrier height for the stacked dimer, even if it can capture the barrier location. Interestingly no second minimum presents in any implicit nonpolar PMF profile. Its behavior is much like the PoissonBoltzmann solvent in our previous studies of the polar PMF profiles of the same base pairs: only the major barriers and the overall trends can be reproduced, but detailed features due to explicit solvent structure are absent. With SAS, the barrier locations are shifted closer to the global minima, although both
J. Phys. Chem. B, Vol. 111, No. 42, 2007 12271 barrier heights are reproduced with acceptable accuracy. Similar to SES, the overall PMF trends can be reproduced, even if no second minimum presents. With SAV, the barrier locations are reproduced better than SAS, but the barrier heights are not reproduced well. Similar to both SES and SAS, no second minimum presents, even if the overall trend of the PMF profiles can be reproduced. Future Directions. Even with the limited success in the dimer test cases, the new nonpolar solvent model (eq 2) still performs much better than the traditional SA model (eq 1) in reproducing ∆Gnp. This is apparent when Figure 1 is compared with Figures 9-11. In the new nonpolar solvent model, even the worst case, SEV with the 6/12 decomposition scheme, performs better than the traditional SA model. A remaining question is whether the new nonpolar solvent model, with an improved decomposition scheme and with SES, SAS, or SAV as the estimator for the repulsive component, works equally well on larger and more complex biomacromolecules. Second, if volume and surface work equally well as estimators for the repulsive solvation free energies for small solutes discussed here, is it necessary to switch from volume to surface when the solute size is increased to the sizes of typical biomolecules? Research in this direction will be addressed in our next report. Acknowledgment. We are grateful for critical reading of the manuscript by Drs. R. Levy, D. Case, M. K. Gilson, and N. Baker. This work is supported in part by NIH GM069620. Appendix Numerical Calculation of Attractive Forces. To calculate attractive forces, we first define
Γ(ras) ) Θ(ras)ras‚nˆ s
(1)
Thus, the partial derivatives of Uuv att with respect to the coordinates of solute atom J (RJ, R ) x, y, z) are
∂Uuv att
Ns
) Fw
∂RJ
Ns
∂
∑ ∑ ∫S a)1 b)1 ∂R
b
Γ(ras) dσs
(2)
J
We can calculate the partial derivative of the integral by definition,
∂ ∂RJ
∫S
∫S
∆
Γ(ras) dσs ) lim
b
∆RJ
∆RJf0
b
Γ(ras) dσs (3)
where
∆
∫S
b
∫S (R + ∆R ) Γ|R +∆R dσs ∫S (R ) Γ|R dσs ) ∫S (R ) [Γ|R +∆R - Γ|RJ] dσs + ∫S′(R )∆R Γ|R + ∆R dσs + O(∆RJ2)
Γ(ras) dσs ) b
b
J
J
J
J
J
b
b
J
J
J
J
J
J
J
J
(4)
In deriving eq 4, we have used the following relation:
Sb(RJ + ∆RJ) = Sb(RJ) + S′b(RJ)∆RJ Thus, we have
(5)
12272 J. Phys. Chem. B, Vol. 111, No. 42, 2007
Tan et al.
Γ|RJ+∆RJ - Γ|RJ
∂D ˆ aJ (DaJ)R ∂DaJ ∂|DaJ| D ˆ + ) D ˆ aJ + |DaJ| ) ∂RJ ∂RJ ∂RJ |DaJ| aJ Rˆ |DaJ| - (DaJ)RD ˆ aJ ) Rˆ (14) |DaJ| 2 |DaJ|
∂ ∂RJ
∫S
b
Γ(ras) dσs ) lim
∆RJf0
{∫ }
∫S′(R )∆R Γ|R +∆R dσs b
J
J
J
J
∆RJ
∆RJ
Sb(RJ)
)
dσs +
∂Sb ∂Γ dσs + Γ h |RJ (6) b ∂R ∂RJ J
∫S
Here, we have computed the second integral with the intermediate value theorem for integration. The average value of Γ (on ∆Sb ) S′b(RJ)∆RJ), Γ h , can then be taken out of the integral. Finally, (∂Sb/∂RJ) is the derivative of SAS and can be found in the literature.58 The first term in eq 6 can be computed numerically from a discrete summation over the uniform surface area elements, δsi on atom “b”:
∫S
b
∂Γ ∂RJ
dσs )
∑
δsi ∈ Sb
∂Γ ∂RJ
δsi
3. Otherwise,
∂Dab )0 ∂RJ
∂Uuv att ∂RJ
∂ras ∂nˆ s ∂ (ras‚nˆ s) ) ‚nˆ s + ras‚ ∂RJ ∂RJ ∂RJ
(11)
Note that the outward normal vector, nˆ s, of every element, dσs, is a constant vector with respect to the atomic center so that we have
∂ras ∂ (r ‚nˆ ) ) ‚nˆ ∂RJ as s ∂RJ s
i
∂Γ ∂RJ
b
δsi + Γ h
|
∂Sb
RJ
]
(17)
∂RJ
(18)
(19)
Equations 17-19 are for the 6/12 decomposition scheme. For the WCA decomposition scheme, eq 19 is replaced by eqs 20 and 21.
(rˆ as)R
(
Aaw
)
Baw
- rmin3 ras‚nˆ s |ras| 3rmin rmin3 Baw (nˆ s)R Aaw + (ras3 - rmin3) if ras e rmin (20) 3 9 3|ras| 3rmin rmin3
M)
M)
4
[
(
Baw
3|ras|
6
9
-
-
Aaw 9|ras|
12
) ( ) (nˆ s)R 4Aaw
3|ras|13
(12)
2Baw |ras|7
]
-
(rˆ as)Rras‚nˆ s if ras > rmin (21)
For the σ decomposition scheme, eq 19 is replaced by eqs 22 and 23.
As shown in Figure 2,
∂ras ∂Dab ∂Ps ∂Dab ) + ) ∂RJ ∂RJ ∂RJ ∂RJ
[
∑∑ ∑ a)1 b)1 δs ∈S
{
Note that
for any vector ras. To compute the second part of ∂Γ/∂RJ (eq 8), we can utilize the following identity:
Ns
M, if a * b, J ) a ∂Γ ) - M, if a * b, J ) b ∂RJ 0, otherwise Baw(nˆ s)R 2Baw(rˆ as)R ras‚nˆ s M) 3|ras|6 |ras|7
(9)
(10)
Ns
) Fw
where
as
∂ras ∂|ras| ) rˆ as‚ ∂RJ ∂RJ
(16)
where (DaJ)R is the R component of vector DaJ, and Rˆ is the vector directed along axis R. With the above preparation, we have
(8)
From eq 28 in Theory and Computational Details and using the 6/12 decomposition as an example, we obtain
∂Θatt 2Baw ∂|ras| ) ∂RJ |r |7 ∂RJ
∂DJb ∂|DJb| ∂D ˆ Jb (DJb)R D ˆ + ) D ˆ Jb + |DJb| ) ∂RJ ∂RJ ∂RJ |DJb| Jb - Rˆ |DJb| + (DJb)RD ˆ Jb ) -Rˆ (15) |DJb| 2 |DJb|
(7)
The derivative ∂Γ/∂RJ can be computed with the strategy of Cossi et al.:59
∂(ras‚nˆ s) ∂Γ ∂Θ ) ras‚nˆ s + Θ ∂RJ ∂RJ ∂RJ
2. a * b, J ) a
(13)
Dab is the vector pointing from the center of atom a to the center of atom b. Ps is the position of area element dσs with respect to the center of atom b. Here, Ps is a constant vector with respect to the atomic center, just like nˆ s. Now we need to distinguish three different cases: 1. a * b, J ) b
M)
(
Aaw
1
M)
|ras| 3σ 3
(
Baw
3|ras|6
-
9
-
)[
]
Baw (rˆ as)R (nˆ s)R ras‚nˆ s 3 |ras| 3 σ
Aaw
)
if ras e σ (22)
(nˆ s)R 9|ras|12 4Aaw 2Baw (rˆ as)Rras‚nˆ s if ras > σ (23) 7 |ras| 3|ras|13
(
)
Implicit Nonpolar Solvent Models
Figure 15. Relative errors in attractive forces (Γ term) and attractive free energy for propane versus the numerical resolution of the solventaccessible surface (represented by N surface dots) in the proposed numerical procedure. The relative errors were computed with respect to the brute-force numerical solution in MATLAB for propane.
Figure 16. Relative errors in attractive forces (S term) for propane versus the numerical resolutions of solvent-accessible surface arcs represented by arc dots separated by arcres in the proposed numerical procedure. The relative errors were computed with respect to the bruteforce numerical solution in MATLAB.
2. Numerical Accuracy in Attractive Free Energy and Forces. The attractive free energy and forces calculated with the proposed numerical procedure is compared with the bruteforce numerical integration procedure in MATLAB using propane as a test case. Since the numerical force procedure consists of two terms (a Γ term that is related to ∂Γ and the S term that is related to ∂S (see eq 17)), the performance of the two terms will be discussed separately below. The accuracy of the Γ term is directly related to the numerical surface resolution, represented as the number of surface area elements per atom, N, in this study. We have studied the relative errors between the proposed numerical procedure and the MATLAB procedure with six different N values, 400 × 2n, (n ) 0, 1, 2, 3, 4, 5), and the results are shown in Figure 15 for propane. From this comparison, it can be found that N g 2500 is necessary for a relative error
