Continuum Polarizable Force Field within the ... - ACS Publications

J. Phys. Chem. B 2008, 112, 7675–7688

7675

Continuum Polarizable Force Field within the Poisson-Boltzmann Framework Yu-Hong Tan, Chunhu Tan, Junmei Wang, and Ray Luo* Department of Molecular Biology and Biochemistry, UniVersity of California, IrVine, California 92697-3900 ReceiVed: NoVember 21, 2007; ReVised Manuscript ReceiVed: March 28, 2008

We have developed and tested a complete set of nonbonded parameters for a continuum polarizable force field. Our analysis shows that the new continuum polarizable model is consistent with B3LYP/cc-pVTZ in modeling electronic response upon variation of dielectric environment. Comparison with experiment also shows that the new continuum polarizable model is reasonable, with accuracy similar to that of B3LYP/ccpVTZ in reproduction of dipole moments of selected organic molecules in the gas phase. We have further tested the validity to interchange the Amber van der Waals parameters between the explicit and continuum polarizable force fields with a series of dimers. It can be found that the continuum polarizable model agrees well with MP2/cc-pVTZ, with deviations in dimer binding energies less than 0.9 kcal/mol in the aqueous dielectric environment. Finally, we have optimized atomic cavity radii with respect to experimental solvation free energies of 177 training molecules. To validate the optimized cavity radii, we have tested these parameters against 176 test molecules. It is found that the optimized Poisson-Boltzmann atomic cavity radii transfer well from the training set to the test set, with an overall root-mean-square deviation of 1.30 kcal/mol, an unsigned average error of 1.07 kcal/mol, and a correlation coefficient of 92% for all 353 molecules in both the training and test sets. Given the development documented here, the next natural step is the construction of a full protein/nucleic acid force field within the new continuum polarization framework. Introduction Most commonly used force fields are additive or nonpolarizable in computational simulations of structures, dynamics, and functions of biomolecules. They incorporate a relatively simple potential energy function:

V(r) )

∑

bonds

∑

torsions

kb(b - b0)2 +

∑

kθ(θ - θ0)2 +

angles

kφ[cos(nφ + δ) + 1] +

∑

nonbond pairs

[

Aij Bij qiqj + 12 - 6 rij r r ij

ij

]

(1) The first three summations are over bonds (1-2 interactions), angles (1-3 interactions), and proper/improper torsions (1-4 interactions). The final summation (over pairs of atoms i and j) excludes 1-2 and 1-3 interactions and often uses different parameters for atoms separated by three covalent bonds (i.e., the 1-4 interactions) from those for atoms separated by more than three covalent bonds. It describes electrostatics interactions that are represented by a Coulombic potential, and dispersion and exchange repulsion interactions that are represented by a Lennard-Jones 6-12 potential. Many such force fields have been developed for biomolecular simulations, such as Amber,1–6 CHARMM,7,8 GROMOS,9 and OPLS.10–12 In force field applications, the covalent bonding terms are not that critical in biomolecular simulations, as long as bond and angle parameters are taken from small molecule structures and vibrational frequencies, and torsion terms are fitted to highlevel quantum mechanical potential energy surfaces. The van der Waals parameters are not very different among various force fields either, mostly derived from Jorgensen’s OPLS parameter * Corresponding author. E-mail: [email protected].

set.10 The most important differences among various force fields are the treatments of electrostatics. To develop parameters for electrostatics, all force fields face the problem of how to effectively wrap polarization into the total potential energy as physically and cleanly as possible. Amber based its parameter development on a set of HF/6-31G* calculations to effectively increase molecular dipole moments,3 to an amount similar to those observed in the TIP3P and SPC water models.13,14 In CHARMM charges were derived primarily from fits to solute-water dimer energetics. In addition to fitting the dimer interaction energies, charges for model compounds were adjusted to obtain dipole moments somewhat larger than experimental or gas-phase ab initio values. OPLS placed a strong emphasis on deriving nonbonded interactions by comparison to liquid-phase thermodynamics.11,12 Although the OPLS parameters were principally derived with reference to condensedphase simulations, comparisons to gas-phase peptide energetics also show good results.15 A cornerstone for accurate biomolecular simulations is the molecular mechanical force field. Accuracy of the force field parameters ultimately determines that of the simulations. Arguably, polarizable force fields of various forms offer a more accurate and balanced treatment of solute-solute and solute-solvent interactions for biomolecular simulations, at least for the environment-sensitive electrostatic components. Indeed, it is widely recognized that the use of an electrostatic model based on fixed charges as in nonpolarizable force fields has a significant shortcoming: the model is unable to respond directly to the molecular environment. Nevertheless, it should be pointed out that the need for explicit inclusion of polarization effects in force fields is still a matter of current discussion in the literature.16,17 Three basic methods for including polarization have been studied: fluctuating charge, Drude oscillator, and induced dipole models. Fluctuating charge models use the principle of elec-

10.1021/jp7110988 CCC: $40.75  2008 American Chemical Society Published on Web 05/29/2008

7676 J. Phys. Chem. B, Vol. 112, No. 25, 2008 tronegativity equalization to produce a set of point charges that optimize the total electrostatic energy. In theory, intermolecular charge transfer is then handled by requiring conservation of charge for the whole system. In applications, charge conservation is often enforced for each individual molecule.18–21 Drude oscillator methods, also referred to as shell models, use a harmonic restraint potential to tether a mobile point charge to an interaction site. In a general Drude oscillator model, an atom carries a charge fixed at the nucleus and a second restrained charge of variable position. The charge magnitudes and harmonic force constants are fit to atomic and molecular polarizability data and experimental energies. As with the fluctuating charge model, the Drude oscillator model is most often used with an extended Lagrangian treatment of the variable charges during a molecular dynamics simulation.22 Probably the best-studied method for handling polarization is the use of induced multipole moments. While higher-order multipole polarization and hyperpolarization can be included in force fields,23 only induced dipoles are usually considered. To date, most polarizable force fields are still under development with a few published applications limited to those from the development groups. The Amber ff02 potential represents an initial polarizable member in the Amber force fields with most of its components retained from nonpolarizable force fields.24,25 The Amoeba force field is a fully polarizable force field with its own polarizable water potential.26,27 The CHARMM fluctuating charge force field is another fully polarizable force field.28,29 This force field has been tested with seven representative globular proteins in explicit polarizable solvent represented by the TIP4P-FQ potential.28,29 The Drude polarizable CHARMM force field incorporates induced polarization into the potential function on the basis of the classical Drude oscillator model.30–34 The early versions of PROSA force fields were first based on a fluctuating charge formalism.35 This model was later extended via the addition of induced dipoles on some sites.36 The SIBFA force field has been applied to a series of protein structure and molecular recognition problems.37–39 The energy function is based on a rather direct decomposition of ab initio SCF computations, including multipole electrostatics, repulsion, and dispersion-like terms.37–39 Notably, separate terms are used for polarization and charge transfer.37–39 Very recently, the NEMO potential for a capped glycine residue has been reported.40 This study devoted considerable attention to the intramolecular potential.40 The SDFF force field is unique among current polarizable protein force field efforts in its extensive use of cross-terms and emphasis on valence potential functions.41 With ever-improving accuracy in implicit solvents,42–45 an interesting question in the development of polarizable force fields is how to incorporate electronic polarization in implicit solvents. Maple et al. directly incorporated an explicit polarizable force field in their generalized Born implicit solvent.46 Schneiders et al. combined a Poisson-Boltzmann solvent and the explicit polarizable Amoeba force field successfully.47 Alternatively, it is also reasonable to treat electronic polarization in a continuum manner since implicit solvents have treated solvent polarization in a continuum manner. However, adoption of such a strategy does not argue for its capability to handle atomic-detailed polarization effects within a molecular environment. Instead, it is intended to offer an efficient and selfconsistent approach in treating polar interactions in biomolecular simulations more satisfactorily than existing additive force fields with implicit solvents. Continuum treatment of electronic polarization was often used in many calculations related to solvation free energies.48,49 Typical implicit solvents, such as

Tan et al. numerical Poisson-Boltzmann approaches50,51 and generalized Born approaches,52,53 were developed to deal with a nonvacuum solute interior. However, a molecular mechanics force field that is designed to be consistent with such a treatment of electronic polarization just began to emerge.54 In a previous study we explored the feasibility of treating electronic polarization in a continuum manner.54 We found that an electrostatic model with a single set of parameters can be used in different environments and conformations with such a continuum treatment of electrostatic polarization.54 In this work we want to investigate the highly important nonbonded components for such a force field. Methods Continuum Polarizable Force Field. In our proposed continuum polarizable model, a continuum dipole moment density within the solute interior is used to represent electronic polarization. The dipole moment density, termed polarization (P) in classical electrostatics, is related to the electric field (E) as

P)

ε-1 E 4π

(2)

where ε is the solute dielectric constant; a value of 4 was found to be optimal for tested molecules.54 Note that this is different from the commonly used high protein dielectric constants >20 or the widely accepted low protein dielectric constant of 2. This discrepancy can be better understood after a brief review of protein dielectric constants as discussed by Warshel and co-workers.55 Warshel and coworkers have pointed out that the protein dielectric constant is 1 when all electrostatic details, such as permanent dipoles and induced dipoles, are represented and sampled explicitly. This is the case for explicit polarizable force fields. Incidentally, this also argues for why existing additive force fields that are designed with a protein dielectric constant of 1 are already fundamentally flawed. They have also pointed out that the protein dielectric constant is 2 or higher when all except induced dipoles are represented and sampled. This is the case for continuum polarizable force fields. In other words, the solute dielectric in the continuum polarizable force fields should be similar to its dielectric at the high-frequency limit, ε∞, as in the Debye equation for the frequency-dependent complex dielectric constant: ε ) ε′ + iε′′, where ε′ ) ε∞ + [(εs - ε∞)/(1 + ω2τ2)] and ε′′ ) [(εs - ε∞)ωτ]/(1 + ω2τ2). Here εs is the static dielectric value at zero frequency, ω is angular frequency, and τ is the relaxation time for molecular orientation. It should be pointed out that ε∞ values are different for different molecules, ranging from about 2.0 for cyclohexane to 4.0 for N-methylacetamide (NMA) and 5.2 for liquid water at room temperature. Thus, this value is very different from the measured or simulated protein macroscopic dielectric constant, which can be 20 and higher. The much higher value is to account for rotational polarization, the exchange of protein and water in the solventexposed region of protein, and also the average over ionization states of acids and bases. This in part explains why a value as high as 20 must be used in many previous applications of implicit solvents in protein stability and pKa predictions where no explicit representation or sampling of any polarization details are present.55–57 The current model adopts an empirically optimized “average” solute dielectric constant of 4 in all molecular mechanics calculations. This value was shown to give the best agreement in included dipole moments between the continuum electronic polarizable model and the B3LYP/ccpVTZ theory for tested model compounds in tested continuum solvents.54

Continuum Polarizable Force Field

J. Phys. Chem. B, Vol. 112, No. 25, 2008 7677

Once polarization is known from eq 2, the polarization energy (Vpol) can then be calculated from induced charges (Find) due to polarization and potential generated by solute point charges (φperm):

Vpol ) -

1 2

∫ Findφperm dυ

Find ) -

1 ∇·P 4π

(3) (4)

Note that the use of eqs 2, 3, and 4 requires the electrostatic field, E. In the proposed continuum polarizable force field, E () -∇φ) can be obtained by solving the Poisson equation with a given set of solute point charges (Fperm, permanent charges):54

∇ · (ε∇φ) ) -4πF

perm

(5)

In the continuum polarizable force field, atomic charges of each molecule are fitted with respect to the quantum mechanical electrostatic potential (ESP) simultaneously in vacuo, a lowdielectric environment (ε ) 4), and a high-dielectric aqueous environment (ε ) 78.39) with a modified RESP methodology that was developed to accommodate the proposed continuum electronic polarization treatment.54 The modified RESP was based on a least-squares fit between the molecular mechanical (MM) electrostatic potential and the quantum mechanical (QM) electrostatic potential as in the original RESP.58,59 The procedure can be summarized in the following steps: (1) Atomic charges Q are initialized by fitting EQM, i.e., the electrostatic potential from QM, without considering electronic polarization. (2) Induced electronic dipole moment density (P) is calculated by eq 2 and the Poisson equation (∇ · εE ) 4πFperm) with given Fperm, i.e., atomic charges Q. (3) Eind, the electrostatic potential from the induced electronic dipole moment density (P), is computed. (4) New atomic charges Q are refitted to the difference between the electrostatic potential from QM and that from the induced electronic dipole moment density (EQM - Eind). (5) Go to step 2 until convergence is achieved. As discussed in our previous study, when the environment is not in vacuo, there are two contributions to polarization: the first is electronic and the second is solvent. What complicates the picture is that there is, in general, a coupling between the two polarization effects. Thus it is not straightforward to separate the final polarization effect into two isolated components: solvent and electronic polarization effects. The MM electrostatic field from solving the Poisson equation includes all the effects. However, the QM electrostatic field does not include the contribution of solvent polarization due to the way that ESP is computed in the Gaussian program.60 Therefore, the contribution of solvent polarization has to be subtracted from the included electrostatic field to achieve a consistency between QM and MM as discussed in our previous study.54 In this study, the quantum mechanical electrostatic potential was obtained in B3LYP/cc-pVTZ with the Gaussian 03 program.60 The PCM solvent61,62 was used to model polarization responses to different dielectric environments. The charge fitting calculations were performed with a revised Amber 9 package.63 The electrostatic field was calculated with a numerical Poisson-Boltzmann (PB) solver64,65 (igb)10) in the SANDER module. Atomic radii were set to be consistent with those in the PCM solvent61,62 in Gaussian 03.60 The probe radius was set to be 0.6 Å.44 Solvent excluded surface was used to assign the dielectric constant in the solute interior. The dielectric constant for the solute region was set to be 4 as discussed above,

Figure 1. Thermodynamic cycle for solvation free energy calculation.

and that for the solvent region was set to be 1 for vacuum, 4 for a low-dielectric environment, and 78.39 for a high-dielectric aqueous environment to be consistent with the PCM solvent in Gaussian 03.60 The ionic strength was set to 0 mM. A grid spacing of 0.2 Å was used in the finite-difference PB solver. To ensure accurate boundary potential on the finite-difference grid, the longest dimension of the finite-difference grid was set to be 4 times that of solute. The relative convergence criterion was 0.0001 for the finite-difference solver. No cutoff for electrostatic interactions was used. Parameter Development Strategy. Since covalent components vary little upon treatment of noncovalent components, it is safe to retain these parameters from previous generations of force fields, at least at this stage of development. However, the two noncovalent components, namely, electrostatics and van der Waals interactions, need careful refinement. Accurate treatment of noncovalent interactions were pioneered by OPLS developers in their efforts to pursue accurate reproduction of the thermodynamic properties of pure organic liquids.11,66 Indeed, the OPLS noncovalent parameters have become the basis for all other modern force fields. The accuracy of noncovalent parameters is also crucial for the development of any polarizable force fields. In this study we assumed that there is an overall consistency between implicit and explicit polarizable models, at least for intramolecular electrostatic interactions, because both models are designed to be consistent with the same benchmark ab initio training data, the B3LYP/cc-pVTZ electrostatics.54 Due to the consistency in electrostatics, intramolecular noncovalent parameters can be made interchangeable, at least approximately, between implicit and explicit polarizable models in Amber. Any future update of van der Waals parameters for the explicit polarizable model can be plugged into the continuum polarizable force field without further changes. This strategy apparently simplifies future refinements of all force fields based on the same benchmark ab initio quantum mechanical theory. The only issue left in the development of noncovalent terms in the continuum polarizable force field is reproduction of solvation free energies of model compounds in a given solvent. In this study, we adopted the widely used thermodynamic cycle shown in Figure 1 to compute solvation free energies of small organic compounds.49 In the thermodynamic cycle, electrostatic free energy, ∆GV,pol , is evaluated by discharging the solute in 1 the first solvent (here in vacuum). The solute, now completely nonpolar, is transferred between the two solvents. This step gives the nonpolar, or repulsive/attractive contributions to the solvation energy, ∆Gnp 2 . The solute is then recharged in the second solvent (here in water), and the associated electrostatic interaction energy, ∆GW,pol , is calculated. The total solvation free energy, 3 ∆Gsol 4 , is obtained from the sum of the free energies of the three separate steps:

7678 J. Phys. Chem. B, Vol. 112, No. 25, 2008 V,pol W,pol ∆Gsol + ∆Gnp 4 ) ∆G1 2 + ∆G3

Tan et al.

(6)

The total electrostatic solvation free energy at steps 1 and 3 is usually termed the polar solvation free energy, and the nonelectrostatic solvation free energy at step 2 is usually termed the nonpolar solvation free energy. Thus, the total solvation free energy of a molecule can be expressed as a sum of polar (electrostatic) and nonpolar (van der Waals) solvation free energies. Both polar and nonpolar components need to be optimized for the continuum polarizable force field. Our previous study of nonpolar solvation free energy, however, shows that the nonpolar solvation free energy model does not contain any compound-specific parameters that can be adjusted for optimal accuracy.45 Thus only the polar solvation free energy needs to be optimized. In doing so the accuracy of the nonpolar solvation free energy model is on the order of ∼1.0 kcal/mol when it is compared with the TIP3P solvent even if an exhaustive optimization was utilized.45 Thus the lack of compound-specific parameters does limit the overall accuracy of the calculated total solvation free energy. Since only the atomic cavity radii are freely adjustable parameters in the polar component, the optimization of the polar component becomes rather straightforward. Here we have used a large and diversified set of generic organic molecules to optimize the atomic cavity radii. Following our previous effort in the calibration of atomic cavity radii for the nonpolarizable Amber force field,44 the atomic cavity radii can first be grouped by Amber atom types and the net charges of functional groups to increase the ratio between parameters and fitting data constraints. To further increase the ratio between parameters and fitting data constraints, we have removed a few rare atom types that cover less than 10 training molecules. Finally, a total of 20 parameters are used for the 39 Amber atom types against 353 experimental data constraints. Computational Details Polar solvation free energy is calculated by the PB solvent with the continuum polarizable solute force field. The boundary between solute and solvent is the solvent excluded surface with a probe radius of 0.6 Å44 and atomic cavity radii to be optimized. The solute and solvent, with associated charges and dielectric boundaries, are mapped onto a finite-difference grid. Electrostatic potentials are calculated by solving the PB equation with the finite-difference solver.64,65 The total polar solvation energy is the sum of steps 1 and 3 of Figure 1

∆Gpol ) 1/2

∑ qi(φiW - φiV)

(7)

i

where qi and φi are the charge and calculated potential at the ith grid point, for the case of transfer from vacuum (V) to water (W). Nonpolar solvation free energy is decomposed into two terms: ∆Grep, repulsive (cavity) free energy, and ∆Gatt, attractive (dispersion) free energy.45 Free energy simulations of spherical cavities have shown that, for small cavities, ∆Grep correlates with the cavity volume, while for large cavities ∆Grep correlates with the cavity surface.67–69 The crossover occurs around the cavity radius of 10 Å.70,71 Interestingly, a surprising conclusion from our previous analysis is that both molecular surfaces and volumes can be used as estimators of repulsive 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.45 A probable reason for the

apparent discrepancy between our previous conclusion45 and the literature is that our simulations were on “realistic” molecules,45 while the literature simulations were performed on ideal spherical cavities. Our additional tests of both molecular surfaces and volumes on larger biomolecules show that molecular volumes transfer better from small molecules to large biomolecules (Tan and Luo, manuscript in preparation). Thus the molecular volumes were used to estimate ∆Grep. Specifically, the solvent accessible volume (SAV, i.e., the molecular volume enclosed by the solvent accessible surface) was used:45

∆Grep ) p(SAV) + c

(8)

where p ) 0.0378 kcal/mol · Å3 is a solvent pressure parameter and c ) -0.5692 kcal/mol.45 These parameters were optimized with respect to the TIP3P solvent, and the SAV was computed with the atomic van der Waals rmin radius definition and an optimized solvent probe of 1.30 Å.45 ∆Gatt is approximated by the van der Waals attractive interaction potential energy between solute (u) and solvent (v) as72–74 uv ∆Gatt ≈ 〈Uatt 〉

(9)

In general, the solute-solvent van der Waals interaction energy can be analytically expressed as the following volume integral: Ns

uv Uatt )

∑ ∫ Faw(raw) Vatt(raw) draw

(10)

a)1

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 the σ decomposition scheme:45

{ {

uLJ(r) r e σ 0 r>σ

(11)

0 reσ Vatt(r) ) u (r) r > σ Lj

(12)

Vrep(r) )

where uLJ(r) ) (A/r12) - (B/r6) is the Lennard-Jones 6-12 potential. Here, A ) 4εσ12 and B ) 4εσ6; ε and σ are the welldepth and radius parameters of the Lennard-Jones potential, respectively. The Amber ff02 nonbonded parameters were used in this study for all molecules. As a first approximation, a uniform distribution (i.e., constant density) can be used in eq 10. Under the constant water density approximation, the volume integral in eq 10 can be transformed into an integral over the solute accessible surface (SAS) with a simple application of the divergence theorem:75 Ns

uv Uatt ) Fw

Ns

∑ ∑ ∫S

a)1 b)1

b

Θatt(ras)ras · n ˆs dσs

(13)

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ˆ is the outward normal vector associated with the SAS element, dσs. In the σ decomposition scheme, the corresponding Θatt in eq 13 can be shown to be


Θatt ) Θatt )

[

1 Aaw Baw - 3 3|ras|3 3σ9 σ Aaw

Baw 9|ras|12 3|ras|6


]

if ras e σ

(14)

if ras > σ

(15)

Here Aaw and Baw are the van der Waal coefficients for interactions between solute atom a and water oxygen atoms (w).45 Empirical parameter optimization was conducted with a genetic algorithm (GA) based optimizer.76 GA is an efficient stochastic optimization method that has been widely applied to minimization problems because it is ideally suited for multipledimensional global search problems where the search space contains multiple local minima and the search variables may or may not be correlated. The GA begins with the generation of an initial population of a given number of solutions. The fitness of each solution, here the root-mean-square deviations between the computed and experimental solvation free energies, is first evaluated. Then a new population is generated via selection, crossover, and mutation based on their fitness. This process is repeated until a desired fitness is reached or the maximum number of generations exceeded. Molecular Systems A set of 353 organic molecules77 was used in the parametrization of the polar component of the continuum polarizable force field. The complete 353-molecule set is further divided into two sets for the optimization and validation of PB cavity radii: (1) the training set, which consists of odd sequence-numbered molecules; and (2) the test set, which consist of even sequencenumbered molecules. Only the solvation free energies in water were used in the parametrization.78–81 All molecules were built with the Sybyl software package of Tripos Inc.82 and optimized using the MMFF94 force field.83–87 For each molecule, only the global minimum conformation was used. Simple conformational searches were performed for some molecules in cases where the global minimum conformations were not apparent. The geometry of each molecule was further optimized in HF/ 6-31G* using the Gaussian 03 program.60 Results and Discussion Quality of Atomic Charges. The agreement between theoretical and experimental solvation free energies strongly depends on the quality of atomic charges for tested molecules. As shown in our previous study, in principle it is possible to use one set of atomic charges to describe different solvent dielectric environments with the continuum polarizable model.54 However, it is impossible to do so when a traditional nonpolarizable model is used.54 This turns out to be the case in the current study with a much larger set of test molecules. Figure 2 shows the correlations of total molecular dipole moments between the continuum polarizable model and B3LYP/cc-pVTZ for the three dielectric environments used, in vacuo, ε ) 4, and ε ) 78.39. It should be pointed out that the total molecular dipole moment from the continuum polarizable model is the sum of the permanent dipole moment from the atomic point charges and the induced dipole moment computed from eq 4 after solving the PB equation. It can be found that the continuum polarizable model is consistent with B3LYP/cc-pVTZ as far as molecular dipole moments are concerned with correlation coefficients higher than 99% for all three different dielectric environments. It should also be emphasized again that it is impossible to use

Figure 2. Dipole moments in vacuo, ε ) 4, and ε ) 78.39 for the 353-molecule set.

a single set of charges in the nonpolarizable model to achieve the same agreement with B3LYP/cc-pVTZ in all three environments simultaneously, as pointed out in our previous study for a smaller set of molecules.54 Comparison with experimental dipole moments also shows that the continuum polarizable model is reasonable as shown in Table 1. The root-mean-square deviation (rmsd), unsigned average error (UAVG), and correlation coefficient are 0.50 D, 0.22 D, and 90%, respectively, between the continuum polarizable model and experiment. This agreement should be viewed in the context that the same set of charges must be used for all three different dielectric environments. Thus, the agreement with experiment in the gas phase alone is not trivial. The dipole moments with B3LYP/cc-pVTZ are also given. Given that B3LYP/cc-pVTZ and the continuum polarizable model are highly correlated, its agreement with experiment should be similar to that of the continuum polarizable model with an rmsd of 0.52 D, a UAVG of 0.20 D, and a correlation coefficient of 90%. It should be pointed out that noticeable deviations exist between the continuum polarizable model and experiment for formates. However, similar deviations also exist for the same compounds with B3LYP/cc-pVTZ. This is not a surprise because the continuum polarizable model is trained to be consistent with B3LYP/cc-pVTZ. Finally the atomic charges for all molecules are given in the Amber prepin format in the Supporting Information. Quality of van der Waals Parameters. The second important component in the continuum polarizable model that must be properly calibrated is the van der Waals interaction. As discussed in Methods, the van der Waals parameters are assumed to be interchangeable between explicit and implicit polarizable force fields if an overall consistent electrostatic treatment can

7680 J. Phys. Chem. B, Vol. 112, No. 25, 2008

Tan et al.

TABLE 1: Gas-Phase Dipole Moments (D) of Selected Organic Compounds in the 353-Molecule Seta compound

MM

QM

expt

propane 2-methyl-1,3-butadiene propyne 1-butyne 1-hexyne toluene ethylbenzene fluoromethane chlorofluoromethane chloromethane dichloromethane bromomethane dibromomethane methanol ethanol phenol diethyl ether diisopropyl ether acetaldehyde propanal acetone acetic acid methyl formate ethyl formate ethyl acetate methylamine ethylamine propylamine diethylamine pyridine 2-methylpyridine 4-methylpyridine acetonitrile nitroethane 1-nitropropane 2-nitropropane N-methylformamide acetamide dimethyl sulfide diethyl sulfide

0.06 0.23 0.76 0.77 0.93 0.32 0.30 1.54 1.72 1.87 1.66 1.89 1.58 1.39 1.35 1.24 0.95 1.04 2.44 2.63 2.68 1.67 3.94 4.03 1.95 1.1 1.1 1.06 0.78 2.14 1.88 2.61 3.81 3.38 3.62 3.44 4.04 3.63 1.51 1.49

0.09 0.26 0.8 0.77 0.93 0.36 0.35 1.71 1.72 1.94 1.67 1.92 1.52 1.57 1.49 1.28 0.99 1.11 2.55 2.69 2.82 1.7 4.08 4.16 2.04 1.26 1.22 1.16 0.81 2.15 1.86 2.63 3.92 3.59 3.76 3.64 4.08 3.77 1.57 1.59

0.08 0.25 0.78 0.78 0.83 0.38 0.59 1.86 1.82 1.9 1.6 1.82 1.43 1.7 1.69 1.22 1.15 1.13 2.75 2.72 2.88 1.7 1.77 1.93 1.78 1.31 1.22 1.17 0.92 2.21 1.85 2.7 3.92 3.23 3.66 3.73 3.83 3.68 1.55 1.54

a

About two per functional groups are chosen. MM, continuum polarizable model; QM, B3LYP/cc-pVTZ; expt, experimental values from the CRC Handbook of Chemistry and Physics, 87th ed.

TABLE 2: Dimer Binding Potential Energies (kcal/mol)a dimers

MM

MP2

ff02

MP2

NMA/NMA ASH/SER AT CG

-3.85 -3.24 -6.41 -12.53

-3.01 -2.64 -6.50 -11.71

-7.16 -8.93 -13.25 -28.65

-7.47 -11.95 -15.29 -28.04

a NMA, N-methylamine; ASH, acetate; SER, methanol; A, adenine; T, thymine; C, cytosine; G, guanine; ff02, explicit polarizable ff02 model; MP2, MP2/cc-pVTZ; MM, continuum polarizable model. Note that ff02 and associated MP2 calculations were performed in the gas phase.

be maintained as in the case between the ff02 explicit polarizable force field and the newly proposed continuum polarizable force field in Amber. We have used dimer binding energies at the MP2/cc-pVTZ level of theory in the PCM solvent to test the compatibility between the continuum polarizable model and the existing van der Waals treatment in the Amber force fields. Table 2 lists the binding energies of four tested dimers in water, the intended environment for using the continuum polarizable model. It can be found that overall the continuum polarizable model agrees

with MP2/cc-pVTZ very well, with deviations less than 0.9 kcal/ mol. This can be compared with the agreement between the ff02 explicit polarizable force field and MP2/cc-pVTZ, which is also shown in Table 2. Note that the deviations between ff02 and MP2/cc-pVTZ are much larger in the gas phase, as high as 3.0 kcal/mol. However, it should be pointed that all binding energies are also uniformly larger in the gas phase. The electrostatic screening in water in general brings down the deviation in the aqueous phase as in the case for the continuum polarizable model. Thus, whether implicit or explicit polarizable treatment is used, reasonable agreement with MP2/cc-pVTZ can be achieved in the molecular mechanics approaches. Calibration of PB Cavity Radii. The third important issue in the development of the continuum polarizable force field is the accuracy of solvation treatment. Note that this is very different from any polarizable or nonpolarizable force field that is intended for explicit solvent simulations, where the solvent interactions have been incorporated within the calibration of the van der Waals parameters of both solute and solvent molecules. As discussed in Methods, a total of 20 parameters are used for the 39 Amber atom types in the optimization against 353 experimental data constraints. The new PB atomic cavity radii are optimized by reproducing the solvation free energies of 177 molecules in the training set, and are tested against the test set of 176 molecules. Genetic algorithm (GA) was used to optimize the cavity radii. The optimized cavity radii along with the corresponding Amber atom types in the continuum polarizable model are given in Table 3. The statistics for the training set, the test set, and the complete 353-molecule set are shown in Table 4. It can be found that the optimized PB atomic cavity radii transfer well from the training set to the test set. rmsd, UAVG, and correlation coefficient are 1.27 kcal/mol, 1.06 kcal/mol, and 93% for the training set, and 1.34 kcal/mol, 1.08 kcal/mol, and 91% for the test set, respectively. Overall an rmsd of 1.30 kcal/mol, a UAVG of 1.07 kcal/mol, and a correlation coefficient of 92% are achieved for the 353-molecule set. Correlation between theoretical and experimental solvation free energies for the 353-molecule set is also given in Figure 3. Table 5 lists the solvation free energies for all molecules by both theory and experiment. The polar components, the nonpolar cavity components, and the nonpolar attractive components of solvation free energies are also given for all 353 molecules. A detailed analysis of agreement between theory and experiment for each functional group is shown in Table 6. It can be found that the agreements for the first five classes of compounds (alkanes, alkenes, alkynes, aromatic hydrocarbons, and fluorides) are very good, with rmsd’s less than 1.0 kcal/mol. The performance for aromatic groups is particularly excellent with rmsd of 0.54 kcal/mol only. However, agreements for nitro compounds and compounds with nitrogen on hetero rings are much worse with rmsd higher than 2.0 kcal/mol. Note that we use a uniform PB cavity radius value for all nitrogen atom types. This may be an oversimplification of highly diversified chemical environments of different nitrogen atoms. A similar discrepancy between theoretical and experimental solvation free energies for amines and amides was also reported before for free energy simulations in explicit solvent in both the Amber and OPLS force fields.88–90 As shown in Table 1, theoretical dipole moments for nitrogen-containing compounds generally agree well with experiment. Therefore, the poor agreements for solvation free energies are mostly likely related to either the continuum solvation approximation or the classical electrostatic approximation of solute-solvent interactions.



TABLE 3: PB Atomic Cavity Radii (Å) and Atom Types atom type

PB radius

description

h1

0.421

h2

0.327

h4

0.577

h5

2.021

ha hc

0.735 0.753

hn ho o oh os c c1 c2 c3 ca cc cd

0.681 0.000 1.117 1.293 1.094 1.200 1.139 1.200 1.181 1.139 1.139 1.139

ce cf

1.139 1.139

cp

1.139

n n1 n2 n3 na nb nc nd

1.114 1.114 1.114 1.114 1.114 1.114 1.114 1.114

nh

1.114

no s sh ss p5

1.114 2.210 2.210 2.210 1.514

py f cl br

1.514 1.612 2.488 2.616

H bonded to aliphatic carbon with 1 electron-withdrawing group H bonded to aliphatic carbon with 2 electron-withdrawing groups H bonded to non-sp3 carbon with 1 electron-withdrawing groups H bonded to non-sp3 carbon with 2 electron-withdrawing groups H bonded to aromatic carbon H bonded to aliphatic carbon without electron-withdrawing group H bonded to nitrogen atoms hydroxyl group oxygen with one connected atom oxygen in hydroxyl group ether and ester oxygen sp2 C carbonyl group sp C sp2 C sp3 C sp2 C in pure aromatic systems sp2 carbons in nonpure aromatic systems sp2 carbons in nonpure aromatic systems, identical to cc inner sp2 carbons in conjugated systems inner sp2 carbons in conjugated systems, identical to ce head sp2 C that connect two rings in biphenyl system sp2 nitrogen in amide groups sp N aliphatic sp2 N with two connected atoms sp3 N with three connected atoms sp2 N with three connected atoms sp2 N in pure aromatic systems sp2 N in nonpure aromatic systems sp2 N in nonpure aromatic systems, identical to nc amine N connected one or more aromatic rings nitro N S with one connected atom sp3 S connected with hydrogen S with one connected atom phosphate with four connected atoms, such as OP(OH)3 special p5 in conjugated system fluorine chlorine bromine

TABLE 4: Statistics for the Training Set (177), the Test Set (176), and the 353-Molecule Seta rmsd UAVG CC

training

test

all

1.27 1.06 93

1.34 1.08 91

1.30 1.07 92

a rmsd, root-mean-square deviation (in kcal/mol); UAVG, unsigned average error (in kcal/mol); CC: correlation coefficient (%).

Conclusions In this work we have developed and tested a complete set of nonbonded parameters for a continuum polarizable force field. Our analysis shows that the new continuum polarizable model is consistent with B3LYP/cc-pVTZ in modeling the electronic response upon variation of dielectric environment. Comparison with experiment also shows that the new continuum polarizable

Figure 3. Calculated and experimental solvation free energies for the 353-molecule set.

model is reasonable, with accuracy similar to that of B3LYP/ cc-pVTZ in reproduction of dipole moments of selected organic molecules in the gas phase. However, it should be pointed out that the proposed continuum polarization scheme is not intended to describe the atomic-detailed polarization within the molecular environment. For example, the response of the continuum polarizable model under an external electric field is left unstudied. Incidentally, such studies have never been conducted for published explicit polarizable models either. It is likely that none of the existing classical polarizable models can properly reproduce the experiments in such studies. Instead, we intend to offer an efficient and self-consistent approach in treating polar interactions in biomolecular simulations more satisfactorily than existing additive force fields with implicit solvents. The modest goal is justified because implicit solvents have been developed for their efficiency. Apparently there are limitations in these simplified solvation treatments, for example, the limitations of atom-centered dielectric models91 and the related issue in the use of solvent excluded surface in molecular dynamics simulations.65 Further improvements in the nonpolar components are also needed.45 Interestingly, efforts to address these limitations have emerged, such as in the treatment of coupling between polar and nonpolar components.92,93 We have further tested the validity of interchanging the Amber van der Waals parameters between the explicit and continuum polarizable force fields with a series of dimers. It can be found that the continuum polarizable model agrees well with MP2/cc-pVTZ, with deviations in dimer binding energies of less than 0.9 kcal/mol in the aqueous dielectric environment. Finally, we have optimized atomic cavity radii with respect to the experimental solvation free energies of 177 training molecules. To validate the optimized cavity radii, we have tested these parameters against 176 test molecules. It is interesting to note that the optimized PB atomic cavity radii transfer well from the training set to the test set. Overall an rmsd of 1.30 kcal/ mol, a UAVG of 1.07 kacl/mol, and a correlation coefficient of 92% are achieved for all 353 molecules in both sets. In addition, a detailed analysis of the agreement between theoretical and experimental solvation free energies for each functional group shows that the new polarizable model works excellently for alkanes, alkenes, alkynes, aromatic hydrocarbons, and fluorides. However, agreements for nitro compounds and compounds with nitrogen on hetero rings are much worse. Because theoretical dipole moments for nitrogen-containing compounds agree very well with experimental values, a likely reason for the discrepancy between theory and experiment may be the oversimplification in using a single PB cavity radius value for all nitrogen atom types in the highly diversified chemical environments. The poor agreement in solvation free energies may also be related to the lack of modeling electron transfer between the solute and


Tan et al.

TABLE 5: Solvation Free Energies (in kcal/mol) of All 353 Organic Compoundsa no.

compound

POL

NPrep

NPatt

total

expt

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70

methane propane 2-methylpropane 2,2-dimethylpropane n-butane 2,2-dimethylbutane cyclopentane n-pentane 2-methylpentane 3-methylpentane 2,4-dimethylpentane 2,2,4-trimethylpentane methylcyclopentane n-hexane cyclohexane methylcyclohexane cis-1,2-dimethylcyclohexane n-heptane n-octane ethylene propylene 2-methylpropene 1-butene 2-methyl-2-butene 3-methyl-1-butene 1-pentene trans-2-pentene 4-methyl-1-pentene 1-hexene cyclohexene trans-2-heptene 1-methylcyclohexene 1-octene 1,3-butadiene 2-methyl-1,3-butadiene 2,3-dimethyl-1,3-butadiene 1,4-pentadiene 1,5-hexadiene acetylene propyne 1-butyne 1-pentyne 1-hexyne 1-heptyne 1-octyne 1-nonyne benzene toluene 1,2,4-trimethylbenzene ethylbenzene m-xylene o-xylene p-xylene propylbenzene 2-propylbenzene butylbenzene 2-butylbenzene tert-butylbenzene tert-amylbenzene naphthalene anthracene phenanthrene acenaphthene p-chlorotoluene fluoromethane 1,1-difluoroethane tetrafluoromethane hexafluoroethane octafluoropropane fluorobenzene

-0.75 -0.39 -1.19 -1.44 -0.61 -1.48 -0.05 -0.33 -1.08 -0.53 -1.89 -2.09 -0.90 -0.57 -0.02 -0.76 -1.24 -0.72 -0.58 -1.99 -2.03 -2.37 -1.92 -1.35 -2.66 -1.99 -1.27 -2.54 -1.92 -1.07 -1.44 -1.41 -2.13 -3.29 -3.23 -2.93 -3.56 -3.38 -4.28 -3.75 -3.71 -3.94 -3.85 -3.89 -4.07 -3.98 -2.56 -3.10 -3.50 -2.57 -3.49 -3.45 -3.37 -2.99 -3.28 -3.36 -3.09 -3.60 -3.59 -4.55 -5.85 -5.48 -4.43 -2.40 -1.16 -1.54 -0.10 -0.08 -0.27 -2.83

6.35 11.18 13.17 14.97 13.44 16.82 13.83 15.68 17.48 17.26 19.22 20.75 15.93 17.94 15.62 17.63 19.36 20.20 22.49 7.72 10.20 12.44 12.54 14.56 14.61 14.79 15.02 16.62 17.07 14.86 19.54 17.05 21.59 11.44 13.44 15.35 13.84 16.18 6.80 9.40 11.85 14.08 16.37 18.62 20.87 23.14 12.95 15.25 19.40 17.49 17.51 17.11 17.55 19.75 19.37 21.99 21.46 20.85 22.73 17.89 22.82 22.55 20.71 16.85 6.73 9.94 8.94 12.27 15.19 13.61

-3.78 -8.46 -10.40 -12.23 -10.36 -13.80 -11.77 -12.50 -14.23 -14.14 -15.99 -17.35 -13.52 -14.66 -13.52 -15.46 -16.80 -16.56 -18.57 -4.95 -7.38 -9.54 -9.52 -11.63 -11.26 -11.61 -11.70 -13.27 -13.47 -12.78 -15.77 -14.57 -17.58 -8.37 -10.58 -12.50 -10.49 -12.51 -4.20 -6.66 -8.78 -11.00 -12.98 -14.97 -17.32 -19.30 -10.93 -12.97 -16.81 -14.91 -14.92 -15.03 -14.96 -16.70 -16.27 -18.59 -18.41 -17.95 -19.72 -16.48 -21.67 -21.61 -19.41 -14.82 -4.32 -7.38 -6.36 -9.53 -12.63 -11.44

1.82 2.33 1.58 1.31 2.47 1.54 2.01 2.85 2.17 2.58 1.34 1.31 1.52 2.71 2.08 1.41 1.32 2.92 3.34 0.77 0.79 0.53 1.10 1.58 0.69 1.18 2.05 0.80 1.68 1.01 2.33 1.07 1.87 -0.22 -0.37 -0.09 -0.21 0.29 -1.68 -1.01 -0.65 -0.86 -0.45 -0.25 -0.52 -0.15 -0.54 -0.82 -0.91 0.01 -0.90 -1.36 -0.78 0.07 -0.18 0.04 -0.04 -0.69 -0.58 -3.14 -4.69 -4.54 -3.13 -0.37 1.25 1.02 2.48 2.67 2.29 -0.66

1.98 1.96 2.32 2.50 2.08 2.59 1.20 2.33 2.52 2.51 2.88 2.85 1.60 2.49 1.23 1.71 1.58 2.62 2.89 1.27 1.27 1.16 1.38 1.31 1.83 1.66 1.34 1.91 1.66 0.37 1.66 0.67 2.17 0.61 0.68 0.40 0.94 1.01 -0.01 -0.48 -0.15 -0.16 0.01 0.60 0.71 1.05 -0.89 -0.76 -0.86 -0.61 -0.80 -0.90 -0.80 -0.53 -0.30 -0.40 -0.45 -0.44 -0.18 -2.41 -4.23 -4.06 -3.4 -1.92 -0.22 -0.11 3.16 3.94 4.28 -0.78



TABLE 5: Continued no.

compound

POL

NPrep

NPatt

total

expt

71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141

2-chloro-1,1,1-trifluoroethane chlorofluoromethane chlorotrifluoromethane dichlorodifluoromethane fluorotrichloromethane 1,1,2-trichloro-1,2,2-trifluoroethane 1,1,2,2-tetrachlorodifluoroethane chloropentafluoroethane 1,1-dichorotetrafluoroethane 1,2-dichlorotetrafluoroethane 1-bromo-1-chloro-2,2,2-trifluoroethane bromotrifluoromethane 1-bromo-1,2,2,2-tetrafluoroethane chloromethane dichloromethane tetrachloromethane chloroethane 1,1-dichloroethane E-1,2-dichloroethane 1,1,1-trichloroethane 1,1,2-trichloroethane 1,1,1,2-tetrachloroethane 1,1,2,2-tetrachloroethane pentachloroethane hexachloroethane 1-chloropropane 2-chloropropane 1,2-dichloropropane 1,3-dichloropropane 1-chlorobutane 2-chlorobutane 1,1-dichlorobutane 1-chloropentane 2-chloropentane 3-chloropentane chloroethylene cis-1,2-dichloroethylene trans-1,2-dichloroethylene trichloroethylene tetrachloroethylene 3-chloropropane chlorobenzene o-chlorotoluene 1,2-dichlorobenzene 1,3-dichlorobenzene 1,4-dichlorobenzene 2,2′-dichlorobiphenyl 2,3-dichlorobiphenyl 2,2,3′-trichlorobiphenyl bromotrichloromethane 1-chloro-2-bromoethane bromomethane dibromomethane bromoethane 1,2-dibromoethane 1-bromopropane 2-bromopropane 1,2-dibromopropane 1,3-dibromopropane 1-bromo-2-methylpropane 1-bromobutane 1-bromo-3-methylbutane 1-bromopentane 3-bromopropene bromobenzene 1,4-dibromobenzne p-bromotoluene 1-bromo-2-ethylbenzene o-bromocumene methanol ethanol

-2.97 -2.19 -0.01 0.00 -0.01 -0.02 -0.03 -0.04 -0.04 -0.02 -2.67 -0.01 -3.08 -1.77 -2.41 -0.03 -0.85 -1.41 -1.68 -0.79 -2.57 -2.43 -3.65 -2.04 -0.05 -1.33 -1.13 -2.87 -2.18 -1.51 -1.25 -1.84 -1.31 -1.10 -1.00 -1.60 -1.44 -1.16 -0.50 0.00 -1.30 -2.02 -2.47 -1.73 -1.40 -1.29 -3.97 -3.64 -3.56 -0.03 -1.79 -2.32 -3.08 -0.96 -1.85 -1.61 -1.13 -3.48 -2.68 -2.65 -1.85 -2.47 -1.57 -3.48 -1.97 -1.14 -2.40 -1.97 -2.49 -5.21 -4.98

11.98 8.35 9.90 10.89 11.85 14.89 15.76 13.17 14.04 14.06 14.37 10.98 13.50 7.74 9.37 12.78 10.29 11.82 11.72 13.39 13.19 14.69 14.60 16.04 17.40 12.61 12.54 13.91 14.00 14.86 14.69 16.26 17.11 16.95 16.80 9.36 10.91 10.99 12.53 14.06 12.61 14.55 16.60 16.02 16.14 16.17 24.52 24.48 25.93 13.75 12.78 8.90 11.57 11.34 13.85 13.66 13.50 15.82 16.13 15.46 15.94 17.73 18.19 12.70 15.60 18.24 17.89 19.57 21.42 6.71 9.30

-9.98 -6.60 -7.84 -9.36 -10.70 -13.59 -14.80 -10.99 -12.16 -12.34 -12.61 -8.85 -11.28 -5.96 -8.14 -12.05 -8.29 -10.28 -10.21 -12.15 -12.16 -13.79 -13.89 -15.50 -17.17 -10.38 -10.36 -12.15 -12.29 -12.55 -12.28 -14.24 -14.49 -14.14 -14.16 -7.17 -9.21 -9.33 -11.23 -13.25 -10.40 -12.87 -14.71 -14.90 -15.03 -14.77 -21.99 -22.57 -23.88 -12.91 -11.06 -6.98 -9.93 -9.17 -11.92 -11.27 -11.25 -13.82 -14.14 -13.04 -13.33 -14.96 -15.41 -10.34 -13.68 -16.64 -15.63 -17.20 -18.85 -5.32 -7.64

-0.96 -0.44 2.05 1.53 1.14 1.28 0.92 2.15 1.84 1.70 -0.91 2.13 -0.86 0.02 -1.19 0.70 1.15 0.14 -0.18 0.46 -1.54 -1.52 -2.93 -1.50 0.17 0.90 1.05 -1.10 -0.47 0.80 1.16 0.19 1.31 1.71 1.64 0.60 0.26 0.49 0.81 0.82 0.91 -0.35 -0.59 -0.61 -0.29 0.11 -1.45 -1.72 -1.52 0.81 -0.06 -0.41 -1.44 1.21 0.09 0.78 1.13 -1.48 -0.70 -0.23 0.76 0.30 1.21 -1.12 -0.04 0.46 -0.14 0.40 0.07 -3.81 -3.33

0.05 -0.77 2.52 1.69 0.82 1.77 0.82 2.86 2.50 2.31 -0.13 1.79 0.52 -0.56 -1.36 0.10 -0.63 -0.85 -1.73 -0.25 -1.95 -1.15 -2.36 -1.36 -1.40 -0.35 -0.24 -1.25 -1.90 -0.14 0.07 -0.70 -0.07 0.07 0.07 0.49 -1.17 -0.76 -0.44 0.05 -0.57 -1.01 -1.15 -1.36 -0.98 -1.01 -2.73 -2.45 -1.99 -0.93 -1.95 -0.82 -2.11 -0.70 -2.10 -0.56 -0.48 -1.94 -1.96 -0.03 -0.41 0.20 -0.08 -0.86 -1.46 -2.30 -1.39 -1.19 -0.85 -5.07 -4.90


Tan et al.

TABLE 5: Continued no. 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212

compound

POL

NPrep

NPatt

total

expt

2,2,2-trifluoroethanol ethylene glycol 1-propanol 2-propanol 1,1,1-trifluoro-2-propanol 2,2,3,3-tetrafluoropropanol 2,2,3,3,3-pentafluoropropanol hexafluoro-2-propanol 2-methyl-1-propanol 1-butanol 2-butanol tert-butyl alcohol 2-methyl-1-butanol 3-methyl butanol 2-methyl-2-butanol 2,3-dimethylbutanol 1-pentanol 2-pentanol 3-pentanol 2-methyl-1-pentanol 2-methyl-2-pentanol 2-methyl-3-pentanol 4-methyl-2-pentanol cyclopentanol 1-hexanol 3-hexanol cyclohexanol 4-heptanol cycloheptanol 1-heptanol 1-octanol allyl alcohol cyclopentanol phenol 4-bromophenol 4-nitrophenol 4-tert-butylphenol 2-cresol 3-cresol 4-cresol 2,2,2-trifluoroethanol p-bromophenol 2-methoxyethanol dimethoxymethane methyl propyl ether methyl isopropyl ether methyl tert-butyl ether diethyl ether ethyl propyl ether dipropyl ether diisopropyl ether di-n-butyl ether tetrahydrofuran 2-methyltetrahydrofuran anisole ethyl phenyl ether 1,1-diethoxyethane 1,2-dimethoxyethane 1,2-diethoxyethane 1,3-dioxolane 1,4-dioxane 2,2,2-trifluoroethyl vinyl ether acetaldehyde propanal butanal pentanal hexanal heptanal octanal nonanal trans-2-butenal

-6.56 -10.63 -5.93 -5.52 -6.18 -7.53 -5.84 -6.68 -6.44 -5.53 -5.49 -6.07 -5.78 -6.02 -5.99 -6.34 -5.58 -5.36 -4.56 -6.02 -6.17 -5.82 -6.07 -5.44 -6.03 -4.73 -5.22 -4.64 -5.39 -5.91 -5.81 -6.83 -5.44 -8.16 -7.69 -12.28 -8.95 -8.57 -8.52 -8.38 -6.56 -7.69 -8.09 -4.02 -2.39 -3.11 -3.51 -2.08 -2.66 -3.51 -4.12 -2.42 -2.23 -3.04 -4.45 -4.68 -4.17 -4.19 -4.60 -3.96 -4.67 -4.47 -3.67 -3.69 -3.82 -3.90 -3.83 -4.00 -4.04 -3.96 -4.42

11.02 9.72 11.60 11.63 13.19 13.53 14.15 14.66 13.54 13.86 13.82 13.67 15.59 15.70 15.64 17.30 16.12 16.09 16.04 17.86 17.89 17.79 17.85 14.32 18.40 18.32 16.10 20.56 17.86 20.66 22.92 10.63 14.32 13.61 16.25 15.95 21.52 15.75 15.88 15.91 11.02 16.25 12.26 12.45 14.15 13.90 15.68 14.45 16.73 19.04 18.53 23.56 12.28 14.64 15.87 18.45 19.83 14.80 19.98 10.69 12.65 15.28 8.77 11.13 13.41 15.67 17.93 20.17 22.44 24.71 12.51

-9.36 -9.01 -9.88 -9.72 -11.28 -11.74 -12.27 -12.75 -11.67 -11.89 -11.81 -11.70 -13.49 -13.67 -13.39 -15.20 -13.94 -13.63 -13.67 -15.42 -15.38 -15.33 -15.37 -13.11 -15.93 -15.68 -14.89 -17.66 -16.52 -17.97 -19.98 -8.70 -13.11 -12.39 -15.06 -16.97 -19.26 -14.24 -14.33 -14.38 -9.36 -15.06 -11.30 -11.35 -11.91 -11.69 -13.35 -12.07 -14.10 -16.22 -15.24 -20.38 -11.21 -13.06 -14.53 -16.24 -16.73 -13.30 -17.65 -10.70 -12.68 -12.56 -6.95 -9.19 -11.24 -13.29 -15.36 -17.27 -19.37 -21.37 -10.47

-4.90 -9.92 -4.21 -3.61 -4.27 -5.74 -3.96 -4.77 -4.57 -3.56 -3.48 -4.09 -3.67 -3.98 -3.74 -4.25 -3.39 -2.90 -2.18 -3.58 -3.66 -3.36 -3.59 -4.23 -3.56 -2.10 -4.01 -1.74 -4.05 -3.22 -2.87 -4.90 -4.23 -6.94 -6.50 -13.29 -6.69 -7.06 -6.97 -6.85 -4.90 -6.50 -7.13 -2.92 -0.15 -0.90 -1.18 0.30 -0.03 -0.69 -0.84 0.77 -1.17 -1.46 -3.11 -2.46 -1.07 -2.70 -2.28 -3.97 -4.70 -1.75 -1.85 -1.74 -1.66 -1.53 -1.26 -1.10 -0.97 -0.62 -2.38

-4.30 -9.30 -4.85 -4.75 -4.16 -4.90 -4.15 -3.76 -4.51 -4.72 -4.61 -4.51 -4.42 -4.42 -4.43 -3.91 -4.49 -4.39 -4.35 -3.93 -3.93 -3.89 -3.74 -5.49 -4.36 -3.68 -4.95 -4.01 -5.49 -4.25 -4.10 -5.03 -5.49 -6.53 -7.10 -10.6 -5.92 -5.86 -5.49 -6.12 -4.31 -7.13 -6.77 -2.93 -1.66 -2.00 -2.21 -1.75 -1.81 -1.16 -0.53 -0.83 -3.12 -3.3 -2.45 -4.28 -3.27 -4.84 -3.53 -4.09 -5.05 -0.12 -3.50 -3.44 -3.18 -3.03 -2.81 -2.67 -2.29 -2.07 -4.23



TABLE 5: Continued no. 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283

compound

POL

NPrep

NPatt

total

expt

trans-2-hexenal trans-2-octenal trans,trans-2,4-hexadienal benzaldehyde m-hydroxybenzaldehyde p-hydroxybenzaldehyde acetone 2-butanone 3-methyl-2-butanone 3,3-dimethylbutanone 2-pentanone 3-pentanone 4-methyl-2-pentanone 2,4-dimethyl-3-pentanone cyclopentanone 2-hexanone 2-heptanone 4-heptanone 2-octanone 2-nonanone 5-nonanone 2-undecanone acetophenone acetic acid propionic acid butyric acid pentanoic acid hexanoic acid 4-amino-3,5,6-trichloropyridine-2-carboxylic acid methyl formate ethyl formate propyl formate methyl acetate isopropyl formate isobutyl formate isoamyl formate ethyl acetate propyl acetate isopropyl acetate butyl acetate isobutyl acetate amyl acetate isoamyl acetate hexyl acetate methyl propionate ethyl propionate propyl propionate isopropyl propionate amyl propionate methyl butyrate ethyl butyrate propyl butyrate methyl pentanoate ethyl pentanoate methyl hexanonate ethyl heptanoate methyl octanoate methyl benzoate methylamine ethylamine propylamine butylamine pentylamine hexylamine dimethylamine diethylamine dipropylamine dibutylamine trimethylamine triethylamine piperazine

-4.61 -4.48 -5.65 -5.06 -10.33 -11.26 -4.10 -4.04 -4.45 -4.68 -4.18 -3.61 -4.61 -5.02 -3.99 -4.35 -4.23 -4.24 -4.31 -4.45 -4.86 -4.46 -5.23 -9.05 -9.33 -9.27 -8.94 -9.40 -12.64 -4.01 -3.95 -4.23 -5.64 -4.28 -4.46 -4.62 -5.77 -6.30 -5.99 -5.84 -6.45 -6.11 -6.47 -6.35 -5.25 -5.82 -6.63 -6.16 -6.38 -5.70 -5.80 -6.47 -6.15 -6.02 -5.68 -5.90 -6.23 -6.54 -5.24 -5.34 -6.04 -5.61 -5.81 -6.18 -4.32 -3.89 -5.08 -4.01 -5.73 -2.87 -7.80

17.11 21.61 16.23 15.05 15.69 15.70 11.15 13.41 15.37 16.88 15.69 15.67 17.53 19.49 13.90 17.96 20.21 20.23 22.46 24.75 24.76 29.27 17.23 9.44 11.72 13.99 16.26 18.52 20.39 9.43 12.04 14.32 11.96 14.21 16.29 18.44 14.56 16.87 16.66 19.13 18.77 21.40 20.96 23.64 14.25 16.85 19.16 18.95 23.69 16.52 19.11 21.44 18.78 21.40 21.05 25.89 25.57 18.15 7.03 9.60 11.89 14.16 16.43 18.68 9.52 14.64 19.26 23.79 11.68 18.58 13.10

-14.54 -18.70 -13.74 -13.73 -15.14 -15.17 -9.25 -11.26 -13.03 -14.44 -13.16 -13.15 -14.96 -16.39 -12.63 -15.15 -17.32 -17.10 -19.20 -21.20 -21.11 -25.22 -15.94 -8.56 -10.57 -12.5 -14.61 -16.60 -21.81 -8.48 -10.68 -12.82 -10.59 -12.60 -14.61 -16.59 -12.69 -14.82 -14.48 -16.90 -16.55 -18.80 -18.47 -20.90 -12.55 -14.60 -16.73 -16.48 -20.81 -14.54 -16.73 -18.75 -16.64 -18.82 -18.70 -22.67 -22.65 -16.84 -5.50 -7.90 -9.98 -11.99 -14.06 -16.00 -8.10 -12.23 -16.36 -20.64 -10.00 -15.74 -13.13

-2.05 -1.57 -3.16 -3.74 -9.78 -10.73 -2.19 -1.89 -2.12 -2.25 -1.66 -1.09 -2.04 -1.92 -2.73 -1.55 -1.34 -1.12 -1.05 -0.89 -1.21 -0.41 -3.95 -8.17 -8.18 -7.77 -7.29 -7.48 -14.07 -3.06 -2.60 -2.73 -4.27 -2.67 -2.78 -2.76 -3.89 -4.24 -3.81 -3.61 -4.23 -3.50 -3.98 -3.61 -3.56 -3.56 -4.20 -3.69 -3.51 -3.71 -3.43 -3.78 -4.01 -3.45 -3.32 -2.68 -3.31 -5.23 -3.70 -3.64 -4.12 -3.44 -3.44 -3.49 -2.90 -1.48 -2.18 -0.86 -4.05 -0.04 -7.83

-3.68 -3.44 -4.64 -4.02 -9.51 -10.47 -3.80 -3.71 -3.24 -2.89 -3.52 -3.41 -3.06 -2.74 -4.68 -3.41 -3.04 -2.93 -2.88 -2.48 -2.67 -2.15 -4.58 -6.70 -6.46 -6.35 -6.16 -6.21 -11.96 -2.78 -2.65 -2.48 -3.31 -2.02 -2.22 -2.13 -3.08 -2.85 -2.65 -2.55 -2.36 -2.45 -2.21 -2.26 -2.97 -2.80 -2.54 -2.22 -1.99 -2.84 -2.50 -2.28 -2.54 -2.52 -2.48 -2.30 -2.05 -4.28 -4.60 -4.61 -4.50 -4.38 -4.09 -4.04 -4.28 -4.06 -3.65 -3.31 -3.23 -3.03 -7.40


Tan et al.

TABLE 5: Continued no. 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353

compound

POL

NPrep

NPatt

total

expt

N,N′-dimethylpiperazine N-methylpiperazine aniline 1,1-dimethyl-3-phenylurea N,N-dimethylaniline ethylenediamine hydrazinehydrazine 2-methoxyethanamine morpholine N-methylmorpholine N-methylpyrrolidine N-methylpiperidine pyrrolidine piperidine pyridine 2-methylpyridine 3-methylpyridine 4-methylpyridine 2-ethylpyridine 3-ethylpyridine 4-ethylpyridine 2,3-dimethylpyridine 2,4-dimethylpyridine 2,5-dimethylpyridine 2,6-dimethylpyridine 3,4-dimethylpyridine 3,5-dimethylpyridine 2-methylpyrazine 2-ethylpyrazine 2-isobutylpyrazine 2-ethyl-3-methoxypyrazine 2-isobutyl-3-methoxypyrazine 9-methyladenine 1-methylthymine methylindole methylimidazole N-propylguanidine acetonitrile propionitrile butyronitrile nitroethane 1-nitropropane 2-nitropropane 1-nitrobutane nitrobenzene 2-nitrotoluene 3-nitrotoluene N,N′-dimethylformamide N-methylformamide acetamide Z-N-methylacetamide propionamide thioanisole dimethyl sulfide diethyl sulfide methyl ethyl sulfide dipropyl sulfide dimethyl disulfide diethyl disulfide trimethyl phosphate triethyl phosphate tripropyl phosphate 2,2-dichloroethenyl dimethyl phosphate dimethyl 4-bromo-2,5-dichlorophenyl thiophosphate dimethyl 2,4,5-trichlorophenyl thiophosphate dimethyl 4-nitrophenyl thiophosphate diethyl 2,4-dichlorophenyl thiophosphate dimethyl 3-methyl-4-thiomethoxyphenyl thiophosphate diethyl 4-nitrophenyl thiophosphonate O-ethyl O′-4-bromo-2-chlorophenyl S-propyl phosphorothioate

-8.69 -8.19 -9.22 -10.22 -4.470 -11.57 -12.42 -8.15 -6.14 -6.47 -3.44 -3.80 -3.12 -3.67 -5.70 -6.06 -5.73 -6.68 -5.09 -5.67 -6.04 -5.64 -7.00 -5.91 -6.31 -6.47 -5.63 -5.46 -4.93 -6.95 -6.38 -7.82 -15.43 -12.83 -8.17 -9.11 -14.84 -4.95 -4.94 -5.06 -5.09 -5.57 -5.10 -5.70 -6.13 -6.45 -6.66 -5.58 -6.55 -9.72 -7.78 -9.29 -3.73 -3.09 -1.02 -2.43 -2.12 -2.66 -1.29 -8.96 -8.83 -9.60 -7.96 -10.51 -4.37 -9.79 -4.91 -6.88 -9.16 -7.77

17.43 15.25 13.86 21.83 18.00 10.34 5.90 12.57 12.86 15.03 14.74 16.54 12.49 14.36 12.29 14.70 14.56 14.57 16.95 16.83 16.83 16.61 16.97 17.03 17.11 16.49 16.87 14.03 16.29 20.31 19.22 23.25 17.37 17.31 18.48 12.68 15.92 8.75 11.18 13.44 11.24 13.55 13.40 15.80 15.31 17.30 17.60 11.92 9.72 9.72 12.22 12.00 16.97 10.50 15.53 13.04 20.15 12.80 17.82 17.05 24.77 31.69 22.01 28.50 28.97 27.12 32.47 31.24 30.51 35.75

-17.09 -15.11 -12.53 -20.90 -16.78 -9.41 -4.96 -11.45 -12.93 -14.83 -13.49 -15.18 -11.42 -13.17 -11.11 -13.03 -13.07 -13.04 -14.88 -14.96 -14.98 -14.95 -15.02 -15.09 -15.00 -14.83 -15.03 -13.23 -15.04 -18.73 -18.03 -21.36 -18.74 -17.84 -17.13 -11.84 -14.86 -6.70 -8.95 -11.09 -10.98 -12.91 -12.79 -15.00 -15.63 -17.12 -17.53 -11.02 -8.87 -8.71 -10.95 -10.79 -15.22 -8.64 -12.92 -10.80 -17.15 -10.89 -15.16 -16.31 -22.38 -28.16 -21.29 -27.93 -27.50 -27.24 -29.28 -28.94 -29.04 -32.49

-8.35 -8.05 -7.89 -9.30 -3.25 -10.64 -11.47 -7.03 -6.20 -6.26 -2.20 -2.45 -2.05 -2.49 -4.52 -4.39 -4.24 -5.16 -3.03 -3.80 -4.20 -3.98 -5.05 -3.97 -4.20 -4.81 -3.79 -4.67 -3.68 -5.37 -5.19 -5.93 -16.80 -13.36 -6.83 -8.26 -13.79 -2.91 -2.71 -2.71 -4.83 -4.94 -4.48 -4.90 -6.45 -6.27 -6.59 -4.68 -5.70 -8.70 -6.51 -8.08 -1.98 -1.23 1.58 -0.19 0.87 -0.75 1.38 -8.22 -6.43 -6.06 -7.24 -9.94 -2.89 -9.91 -1.72 -4.58 -7.68 -4.50

-7.58 -7.77 -5.49 -11.87 -2.90 -9.75 -9.30 -6.55 -7.17 -6.34 -3.97 -3.89 -5.47 -5.10 -4.69 -4.62 -4.77 -4.92 -4.32 -4.60 -4.72 -4.81 -4.85 -4.70 -4.60 -5.21 -4.84 -5.51 -5.45 -5.05 -4.39 -3.68 -13.60 -10.40 -5.91 -10.25 -10.92 -3.89 -3.85 -3.64 -3.71 -3.34 -3.14 -3.08 -4.12 -3.59 -3.45 -4.90 -10.00 -9.72 -10.00 -9.42 -2.73 -1.54 -1.43 -1.49 -1.27 -1.83 -1.63 -8.70 -7.80 -6.10 -6.61 -5.70 -5.06 -7.62 -3.86 -6.92 -6.27 -4.09

a POL, the polar/electrostatic component; NPrep, the repulsive nonpolar component; NPatt, the attractive nonpolar component; total, total theoretical values; expt, experimental values.

solvent molecules in the continuum solvation treatments, which has been argued for the particular difficult modeling of solvation free energies of amines by classical force fields.88–90

Based on the development documented here, our next natural step is the construction of a full protein/nucleic acid force field within the new continuum polarization framework. Given



TABLE 6: Performance by Solute Functional Groupa solute group

no. of solute

rmsd

UAVG

alkanes alkenes alkynes aromatic hydrocarbons fluorides chlorides bromides alcohols ethers aldehydes ketones acids esters amines amides nitriles nitro compounds compounds with N in hetero rings compounds with S compounds with P

19 (1-19) 20 (20-39) 7 (40-46) 18 (47-64) 19 (65-83) 38 (84-121) 18 (122-139) 45 (140-184) 19 (185-203) 15 (204-218) 17 (219-235) 6 (236-241) 29 (242-270) 27 (271-297) 23 (298-320) 4 (321-324) 6 (325-330) 5 (331-335) 7 (336-342) 11 (343-353)

0.74 0.78 0.84 0.54 0.89 1.10 1.35 1.07 1.44 1.44 1.59 1.55 1.10 1.65 1.43 1.05 2.24 2.59 1.94 1.97

0.59 0.68 0.78 0.41 0.74 0.98 1.16 0.90 1.25 1.33 1.52 1.52 1.00 1.36 1.10 1.05 2.15 2.07 1.66 1.59

a No. of solute: number of solute molecules. The ranges of compound sequence numbers as in Table 5 are also shown for the listed solute groups. rmsd and UAVG are in kcal/mol.

implicit solvents that are based on the Poisson-Boltzmann theory, for example, either numerical Poisson-Boltzmann approaches50,51 or generalized Born approaches52,53 that were developed to deal with nonvacuum solute interiors, the dielectric forces can be computed on the fly in molecular dynamics simulations. Therefore, such a polarizable force field can be applied in implicit-solvent molecular dynamics simulations with virtually no additional cost over a nonpolarizable force field. Nevertheless, more detailed refinement efforts in both the covalent terms and the van der Waals terms in the Amber polarizable force fields, which have been left untouched, should also be addressed in future publications. Acknowledgment. We thank the Amber developer team for exciting discussion in the development of the new continuum polarizable force field. This work is supported by NIH (GM069620 and GM079383). Supporting Information Available: Additional theoretical information. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Weiner, P. K.; Kollman, P. A. J. Comput. Chem. 1981, 2, 287– 303. (2) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, S.; Weiner, P. J. Am. Chem. Soc. 1984, 106, 765– 784. (3) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179–5197. (4) Wang, J. M.; Cieplak, P.; Kollman, P. A. J. Comput. Chem. 2000, 21, 1049–1074. (5) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G. M.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; Caldwell, J.; Wang, J. M.; Kollman, P. J. Comput. Chem. 2003, 24, 1999–2012. (6) Yang, L. J.; Tan, C. H.; Hsieh, M. J.; Wang, J. M.; Duan, Y.; Cieplak, P.; Caldwell, J.; Kollman, P. A.; Luo, R. J. Phys. Chem. B 2006, 110, 13166–13176. (7) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187–217. (8) Mackerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; JosephMccarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick,

S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. B 1998, 102, 3586–3616. (9) Oostenbrink, C. J. Comput.Chem. 2004, 25, 1656–1676. (10) Jorgensen, W. L.; Tiradorives, J. J. Am. Chem. Soc. 1988, 110, 1657–1666. (11) Jorgensen, W. L.; Maxwell, D. S.; Tiradorives, J. J. Am. Chem. Soc. 1996, 118, 11225–11236. (12) Kaminski, G. A.; Friesner, R. A.; Tirado-Rives, J.; Jorgensen, W. L. J. Phys. Chem. B 2001, 105, 6474–6487. (13) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926–935. (14) Hermans, J.; Berendsen, H. J. C.; Vangunsteren, W. F.; Postma, J. P. M. Biopolymers 1984, 23, 1513–1518. (15) Beachy, M. D.; Chasman, D.; Murphy, R. B.; Halgren, T. A.; Friesner, R. A. J. Am. Chem. Soc. 1997, 119, 5908–5920. (16) Van Der Vaart, A.; Bursulaya, B. D.; Brooks, C. L.; Merz, K. M. J. Phys. Chem. B 2000, 104, 9554–9563. (17) Roux, B.; Berneche, S. Biophys. J. 2002, 82, 1681–1684. (18) Rappe, A. K.; Goddard, W. A. J. Phys. Chem. 1991, 95, 3358– 3363. (19) Rick, S. W.; Stuart, S. J.; Berne, B. J. J. Chem. Phys. 1994, 101, 6141–6156. (20) Liu, Y. P.; Kim, K.; Berne, B. J.; Friesner, R. A.; Rick, S. W. J. Chem. Phys. 1998, 108, 4739–4755. (21) Kitao, O.; Ogawa, T. Mol. Phys. 2003, 101, 3–17. (22) Mitchell, P. J.; Fincham, D. J. Phys: Condens. Matter 1993, 5, 1031–1038. (23) Dykstra, C. E. J. Am. Chem. Soc. 1989, 111, 6168–6174. (24) Cieplak, P.; Caldwell, J.; Kollman, P. J. Comput. Chem. 2001, 22, 1048–1057. (25) Lybrand, T. P.; Kollman, P. A. J. Chem. Phys. 1985, 83, 2923– 2933. (26) Ren, P. Y.; Ponder, J. W. J. Comput. Chem. 2002, 23, 1497–1506. (27) Ren, P. Y.; Ponder, J. W. J. Phys. Chem. B 2003, 107, 5933–5947. (28) Patel, S.; Brooks, C. L. J. Comput. Chem. 2004, 25, 1–15. (29) Patel, S.; Mackerell, A. D.; Brooks, C. L. J. Comput. Chem. 2004, 25, 1504–1514. (30) Lamoureux, G.; Mackerell, A. D.; Roux, B. J. Chem. Phys. 2003, 119, 5185–5197. (31) Lamoureux, G.; Harder, E.; Vorobyov, I. V.; Roux, B.; Mackerell, A. D. Chem. Phys. Lett. 2006, 418, 245–249. (32) Noskov, S. Y.; Lamoureux, G.; Roux, B. J. Phys. Chem. B 2005, 109, 6705–6713. (33) Vorobyov, I. V.; Anisimov, V. M.; Mackerell, A. D. J. Phys. Chem. B 2005, 109, 18988–18999. (34) Lopes, P. E. M.; Lamoureux, G.; Roux, B.; Mackerell, A. D. J. Phys. Chem. B 2007, 111, 2873–2885. (35) Banks, J. L.; Kaminski, G. A.; Zhou, R. H.; Mainz, D. T.; Berne, B. J.; Friesner, R. A. J. Chem. Phys. 1999, 110, 741–754. (36) Stern, H. A.; Kaminski, G. A.; Banks, J. L.; Zhou, R. H.; Berne, B. J.; Friesner, R. A. J. Phys. Chem. B 1999, 103, 4730–4737. (37) Gresh, N. J. Chim. Phys. Phys.-Chim. Biol. 1997, 94, 1365–1416. (38) Gresh, N. J. Phys. Chem. A 1997, 101, 8680–8694. (39) Masella, M.; Gresh, N.; Flament, J. P. J. Chem. Soc., Faraday Trans. 1998, 94, 2745–2753. (40) Hermida-Ramon, J. M.; Brdarski, S.; Karlstrom, G.; Berg, U. J. Comput. Chem. 2003, 24, 161–176. (41) Palmo, K.; Mannfors, B.; Mirkin, N. G.; Krimm, S. Biopolymers 2003, 68, 383–394. (42) Nina, M.; Beglov, D.; Roux, B. J. Phys. Chem. B 1997, 101, 5239– 5248. (43) Swanson, J. M. J.; Adcock, S. A.; Mccammon, J. A. J. Chem. Theory Comput. 2005, 1, 484–493. (44) Tan, C. H.; Yang, L. J.; Luo, R. J. Phys. Chem. B 2006, 110, 18680– 18687. (45) Tan, C. H.; Tan, Y. H.; Luo, R. J. Phys. Chem. B 2007, 111, 12263– 12274. (46) Maple, J. R.; Cao, Y. X.; Damm, W. G.; Halgren, T. A.; Kaminski, G. A.; Zhang, L. Y.; Friesner, R. A. J. Chem. Theory Comput. 2005, 1, 694–715. (47) Schnieders, M. J.; Baker, N. A.; Ren, P. Y.; Ponder, J. W. J. Chem. Phys. 2007, 126, 124114. (48) Sharp, K.; Jeancharles, A.; Honig, B. J. Phys. Chem. 1992, 96, 3822–3828. (49) Sitkoff, D.; Sharp, K. A.; Honig, B. J. Phys. Chem. 1994, 98, 1978– 1988. (50) Gilson, M. K.; Honig, B. Proteins: Struct. Funct. Genet. 1988, 4, 7–18. (51) Honig, B.; Nicholls, A. Science 1995, 268, 1144–1149. (52) Sigalov, G.; Scheffel, P.; Onufriev, A. J. Chem. Phys. 2005, 122, 094511.

7688 J. Phys. Chem. B, Vol. 112, No. 25, 2008 (53) Feig, M.; Im, W.; Brooks, C. L. J. Chem. Phys. 2004, 120, 903– 911. (54) Tan, Y. H.; Luo, R. J. Chem. Phys. 2007, 126, 094103. (55) Schutz, C.; Warshel, A. Proteins: Struct. Funct. Genet. 2001, 44, 400–417. (56) Honig, B.; Sharp, K.; Yang, A. S. J. Phys. Chem. 1993, 97, 1101– 1109. (57) Gilson, M. K. Curr. Opin. Struct. Biol. 1995, 5, 216–223. (58) Woods, R. J.; Khalil, M.; Pell, W.; Moffat, S. H.; Smith, V. H. J. Comput. Chem. 1990, 11, 297–310. (59) Dupradeau, F. Y.; Pigache, A.; Zaffran, T.; Cieplak, P. Submitted for publication. (60) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao; O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03; Gaussian, Inc.: Wallingford, CT, 2004. (61) Tomasi, J.; Persico, M. Chem. ReV. 1994, 94, 2027–2094. (62) Tomasi, J.; Mennucci, B.; Cammi, R. Chem. ReV. 2005, 105, 2999– 3093. (63) Case, D. A.; Darden, T. A.; Cheatham, T. E., III; Simmerling, C. L.; Wang, C. L.; Duke, R. E.; Luo, R.; Merz, K. M.; Pearlman, D. A.; Crowley, M.; Walker, R. C.; Zhang, W.; Wang, B.; Hayik, S.; Roitberg, A.; Seabra, G.; Wong, K. F.; Paesani, F.; Wu, X.; Brozell, S.; Tsui, V.; Gohlke, H.; Yang, L.; Tan, C.; Mongan, J.; Hornak, V.; Cui, G.; Beroza, P.; Mathews, D. H.; Schafmeister, C.; Ross, W. S.; Kollman, P. A. Amber 9; University of California, San Francisco, 2006. (64) Luo, R.; David, L.; Gilson, M. K. J. Comput. Chem. 2002, 23, 1244–1253. (65) Lu, Q.; Luo, R. J. Chem. Phys. 2003, 119, 11035–11047. (66) Kaminski, G.; Duffy, E. M.; Matsui, T.; Jorgensen, W. L. J. Phys. Chem. 1994, 98, 13077–13082. (67) Pierotti, R. A. Chem. ReV. 1976, 76, 717–726.

Tan et al. (68) Irisa, M.; Nagayama, K.; Hirata, F. Chem. Phys. Lett. 1993, 207, 430–435. (69) Chandler, D. Nature 2005, 437, 640–647. (70) Lum, K.; Chandler, D.; Weeks, J. D. J. Phys. Chem. B 1999, 103, 4570–4577. (71) Huang, D. M.; Chandler, D. Proc. Natl. Acad. Sci. U.S.A. 2000, 97, 8324–8327. (72) Pratt, L. R.; Chandler, D. J. Chem. Phys. 1980, 73, 3434–3441. (73) Gallicchio, E.; Kubo, M. M.; Levy, R. M. J. Phys. Chem. B 2000, 104, 6271–6285. (74) Huang, D. M.; Chandler, D. J. Phys. Chem. B 2002, 106, 2047– 2053. (75) Floris, F.; Tomasi, J. J. Comput. Chem. 1989, 10, 616–627. (76) Holland, J. H. SIAM J. Comput. 1973, 2, 88–105. (77) Wang, J. M.; Wang, W.; Huo, S. H.; Lee, M.; Kollman, P. A. J. Phys. Chem. B 2001, 105, 5055–5067. (78) Hine, J.; Mookerjee, P. K. J. Org. Chem. 1975, 40, 292–298. (79) Pearson, R. G. J. Am. Chem. Soc. 1986, 108, 6109–6114. (80) Hawkins, G. D.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. B 1997, 101, 7147–7157. (81) Li, J. B.; Zhu, T. H.; Hawkins, G. D.; Winget, P.; Liotard, D. A.; Cramer, C. J.; Truhlar, D. G. Theor. Chem. Acc. 1999, 103, 9–63. (82) Sybyl Manual; Tripos Inc., 1699 South Hanley Rd., St. Louis, MO 63144-2913, 1998. (83) Halgren, T. A. J. Comput. Chem. 1996, 17, 490–519. (84) Halgren, T. A. J. Comput. Chem. 1996, 17, 520–552. (85) Halgren, T. A. J. Comput. Chem. 1996, 17, 553–586. (86) Halgren, T. A.; Nachbar, R. B. J. Comput. Chem. 1996, 17, 587– 615. (87) Halgren, T. A. J. Comput. Chem. 1996, 17, 616–641. (88) Rizzo, R. C.; Jorgensen, W. L. J. Am. Chem. Soc. 1999, 121, 4827– 4836. (89) Morgantini, P. Y.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 6057–6063. (90) Ding, Y. B.; Bernardo, D. N.; Kroghjespersen, K.; Levy, R. M. J. Phys. Chem. 1995, 99, 11575–11583. (91) Swanson, J. M. J.; Mongan, J.; Mccammon, J. A. J. Phys. Chem. B 2005, 109, 14769–14772. (92) Dzubiella, J.; Swanson, J. M. J.; Mccammon, J. A. J. Chem. Phys. 2006, 124, 084905. (93) Dzubiella, J.; Swanson, J. M. J.; Mccammon, J. A. Phys. ReV. Lett. 2006, 96, 087802.

JP7110988

Continuum Polarizable Force Field within the ... - ACS Publications

Recommend Documents