Untitled - Jay Ponder Lab

17 Pages·2012·0.71 MB·English

J. Phys. Chem. B 2010, 114, 2549–2564 2549
Current Status of the AMOEBA Polarizable Force Field
Jay W. Ponder and Chuanjie Wu
Department of Biochemistry and Molecular Biophysics, Washington UniVersity, St. Louis, Missouri 63110
Pengyu Ren
Department of Biomedical Engineering, UniVersity of Texas, Austin, Texas 78712-1062
Vijay S. Pande,†,‡ John D. Chodera,† Michael J. Schnieders,† and Imran Haque‡
Departments of Chemistry, Computer Science, Stanford UniVersity, Stanford, California 94305
David L. Mobley
Department of Chemistry, UniVersity of New Orleans, New Orleans, Louisiana 70148
Daniel S. Lambrecht, Robert A. DiStasio, Jr., and Martin Head-Gordon
Department of Chemistry, UniVersity of California, Berkeley, California 94720
Gary N. I. Clark, Margaret E. Johnson, and Teresa Head-Gordon*
Department of Bioengineering, UniVersity of California, Berkeley, California 94720
ReceiVed: NoVember 9, 2009; ReVised Manuscript ReceiVed: December 17, 2009
Molecular force fields have been approaching a generational transition over the past several years, moving
awayfromwell-establishedandwell-tuned,butintrinsicallylimited,fixedpointchargemodelstowardmore
intricateandexpensivepolarizablemodelsthatshouldallowmoreaccuratedescriptionofmolecularproperties.
TherecentlyintroducedAMOEBAforcefieldisaleadingpubliclyavailableexampleofthisnextgeneration
of theoretical model, but to date, it has only received relatively limited validation, which we address here.
We show that the AMOEBA force field is in fact a significant improvement over fixed charge models for
smallmoleculestructuralandthermodynamicobservablesinparticular,althoughfurtherfine-tuningisnecessary
to describe solvation free energies of drug-like small molecules, dynamical properties away from ambient
conditions,andpossibleimprovementsinaromaticinteractions.Stateoftheartelectronicstructurecalculations
revealgenerallyverygoodagreementwithAMOEBAfordemandingproblemssuchasrelativeconformational
energiesofthealaninetetrapeptideandisomersofwatersulfatecomplexes.AMOEBAisshowntobeespecially
successfulonprotein-ligandbindingandcomputationalX-raycrystallographywherepolarizationandaccurate
electrostatics are critical.
Introduction relyonapproximationsandempiricalinputinordertoformulate
tractabledescriptionsofthe(bio)materialinarealisticchemical
Molecular simulation is now an accepted and integral part
environment.
ofcontemporarychemistry,biology,andmaterialscience.The
Nonpolarizable(fixedcharge)modelsprovideaninexpensive
allure of molecular simulation is that most if not all relevant
description or “effective” potential with approximations that
structural,kinetic,andthermodynamicobservablesofachemical
cannot fully capture many-body effects such as electronic
system can be calculated at one time, in the context of a
polarization. Fixed charge protein and water models went
molecularmodelthatcanprovideinsightandnewhypotheses.
through an extensive period of validation for several decades
The predictive quality of these observables depends on the
after they were introduced.1 The consensus of a number of
accuracy of the potential energy surface and the ability to
validation studies is that, while fixed charge models offer
characterize it through effective sampling of configurations or
tractabledescriptionsandarerobustforequilibriumproperties
phase space. Over the last two decades, the field of molecular
forhomogeneoussystems,evidentdiscrepancieswereidentified
simulation has been dominated by research problems such as
between simulations and experiments away from ambient
proteinfoldingwheredynamicaltimescalesorconfigurational
conditions, for dynamical properties, and for heterogeneous
sampling are the biggest bottlenecks to reaching testable
chemicalsystemsingeneral.1,2Polarizableempiricalforcefields,
hypotheses or comparisons to experimental results. Given the
which offer a clear and systematic improvement in functional
demands of sampling over so many degrees of freedom to
formbyincludingmany-bodyeffects,havebeenintroducedinto
convergence,potentialenergysurfacesformolecularsimulations
the chemical and biochemical simulation community over the
pasttwodecades,andonlyrecentlyforbiomolecularsimulation.1,3
†DepartmentofChemistry,StanfordUniversity.
‡DepartmentofComputerScience,StanfordUniversity. Thequestioninmolecularcomputationcurrentlyiswhethernew
10.1021/jp910674d  2010 American Chemical Society
Published on Web02/05/2010
2550 J. Phys. Chem. B, Vol. 114, No. 8, 2010 Ponder et al.
Jay Ponder (Ph.D. 1984, Harvard University, 1985-1990, Postdoctoral sincesomanyacademicandindustryresearchersusemolecular
Fellow,YaleUniversity)isonthefacultyatWashingtonUniversity,where mechanics and molecular dynamics methodology to tackle
hisgrouphasalongstandinginterestindevelopmentofsoftwaretoolsfor
biological,chemical,andmaterialscienceproblemsofinterest.
molecularmodeling,withparticularemphasisonaccurateconformational
analysisandcalculationofintermolecularinteractions. For example, the TINKER software program for molecular
mechanicsanddynamicssimulation4hasbeendownloadedby
ChuanjieWu(Ph.D.2006,TianjinUniversity)didhisdoctorateinforce
close to 60000 separate external users, including essentially
fielddevelopmentandapplicationandnowisapostdoctoralresearcherin
AMOEBAforcefieldparametrizationandrelatedmethodologydevelopment every major research university and many biotech and phar-
atWashingtonUniversity,St.Louis. maceuticalcompanies,andevengreaternumbersareexpected
for other simulation software packages such as Amber,5
PengyuRen(Ph.D.1999,UniversityofCincinnati;2000-2005Postdoctoral
CHARMM,6GROMACS,7andNAMD.8Forcefieldvalidation
Researcher,WashingtonUniversity)isanAssistantProfessorofBiomedical
Engineering at University of Texas at Austin, developing a range of and subsequent improvement has been the admirable history
computationaltoolstostudyproteinsandnucleicacids,withemphasison of the large community effort on fixed charge force fields led
protein-ligandbindingandcomputationaldrugdiscovery. by developers of the Amber,9 CHARMM,10 GROMOS,11 and
VijayPande(Ph.D.1995,MIT;1996-1999MillerFellow,UCBerkeley) OPLS12 potential energy models over many decades. In this
isaProfessorofChemistry,StructuralBiology,andComputerScienceat Feature Article, we hope to continue that tradition by sum-
StanfordUniversity,workingontheoreticalandcomputationalmethodsfor marizingsomeimportantearlyvalidationtestsbyaconsortium
biomolecularsimulationandapplicationstothebiophysicsofproteinfolding
of research groups at Washington University St. Louis, Uni-
andmisfoldingdiseases.
versityofTexasatAustin,UCBerkeley,andStanfordUniver-
JohnD.Chodera(Ph.D.2006,UCSanFrancisco;2006-2008Postdoctoral sity conducted on the general purpose polarizable force field
Researcher, Stanford University) is a Distinguished Postdoctoral Fellow AMOEBA (atomic multipole optimized energetics for biomo-
with the California Institute of Quantitative Biosciences (QB3) at the
lecular applications) developed by Ponder and co-workers.13
UniversityofCalifornia,Berkeley,workingonthestatisticalmechanicsof
biomolecularfunction,withemphasisonconformationaldynamics,single- The first level of comprehensive testing of any force field
moleculeexperiments,anddrugdiscovery. willincludepredictionsmadebythatpotentialagainstthebest
experiments and theoretical calculations available on a wide
MichaelJ.Schnieders(Ph.D.2007,WashingtonUniversity,St.Louis)is
a Postdoctoral Fellow at Stanford University working on computational arrayofsmallmoleculedatainbothgasphaseandcondensed
X-raycrystallography. phaseenvironments.Infact,AMOEBAbelongstotheclassof
molecular mechanics force fields that aims for high fidelity to
ImranHaque(B.S.2006,UCBerkeley)ispursuingaPh.D.atStanford
abinitiocalculationsbutatacomputationalcostthatmakesit
Universityworkingoncomputer-aideddrugdesignandhigh-performance
methodsformolecularsimulation. suitedforbothsmallmoleculeandbiomoleculecondensedphase
studies where statistical mechanical sampling is necessary. In
David L. Mobley (Ph.D. 2004, UC Davis; 2005-2008 Postdoctoral
practicalterms,AMOEBAisintermediateincomputationalcost
Researcher, UC San Francisco; 2008 Chief Science Officer, Simprota
Corporation)isanAssistantProfessoratUniversityofNewOrleansapplying between other transferable polarizable force fields such as
computational and theoretical methods to understand and quantitatively SIBFA (sum of interactions between fragments ab initio),14
predictprotein-ligandbinding,solvation,andsolubility. NEMO (non-empirical molecular orbital),15 and QM/MM ap-
proaches such as DRF16 and inexpensive polarizable biomo-
DanielSebastianLambrecht(Dipl.-Chem.2003,U.Du¨sseldorf,Dr.rer.
nat. 2007, U. Tu¨bingen) completed his doctorate on the development of lecularforcefieldsfromtheAmber,3iCHARMM,3e,f,handOPLS/
linear-scalingapproachesinquantumchemistry.Heisnowapostdoctoral PFF consortiums.3d,g In this paper, we review the AMOEBA
scholaratUCBerkeley.
modelanditsperformanceinseveralareasincludinggasphase
propertiesagainststate-of-the-artquantummechanicalcalcula-
RobertDistasio,Jr.(Ph.D.2009,UCBerkeley)completedhisdoctorate
ondevelopmentandapplicationoffastelectronicstructuremethodsthat tions, aqueous peptide solvation, structure and dynamics,
includeelectroncorrelation.HeisnowaPostdoctoralResearchAssociate solvationfreeenergiesofsmallmoleculeproteinanaloguesand
atPrincetonUniversity.
drug-likemoleculeswithhighprecision,earlystructuralstability
studies of aqueous solvated proteins, computational X-ray
MartinHead-Gordon(Ph.D.1989,CarnegieMellon,1989-1992,Post- crystallography, and protein-ligand binding.
doctoralMemberofTechnicalStaff,AT&TBellLaboratories)leadsagroup
atUCBerkeleythatdevelopsandapplieselectronicstructuretheory.
The AMOEBA Force Field
Gary N. I. Clark (B.Sc. 2002, U. Loughborough; Ph.D. 2008, Imperial The AMOEBA force field has the following general func-
CollegeLondon)isapostdoctoralresearcheratUCBerkeleyworkingon tional form for the interactions among atoms
water structure, clustering, and dynamics at hydrophobic/hydrophilic
interfaces.
U ) U + U + U + U + U + U +
bond angle bθ oop torsion vdW
MargaretE.Johnson(B.S.2004ColumbiaU.Ph.D.2009,UCBerkeley) Uperm + Uind (1)
completedherdoctorateatUCBerkeleyonstudiesofbulkandhydration ele ele
waterstructureanddynamicsandisnowanNIHpostdoctoralresearcher.
where the first five terms describe the short-range valence
TeresaHead-Gordon(Ph.D.1989,CarnegieMellon,1990-1992,Post-
interactions(bondstretching,anglebending,bond-anglecross
doctoralMemberofTechnicalStaff,AT&TBellLaboratories)leadsagroup
at UC Berkeley that develops theoretical/experimental methods to study term,out-of-planebending,andtorsionalrotation)andthelast
biomaterialsassemblyandbulkandhydrationwaterproperties. threetermsarethenonbondedvdWandelectrostaticcontribu-
tions. AMOEBA contains a number of differences from
“traditional”biomolecularpotentialssuchasthecurrentAmber
polarizableforcefieldparametrizationshavesuccessfullyreached ff99SB,9c CHARMM27,10 OPLS-AA,12b,c and GROMOS
a new level of predictive power over their nonpolarizable 53A611b in the use of bond-angle cross terms, a formal
predecessors. Wilson-Decius-Crossdecompositionofanglebendingintoin-
Demonstrable testing of empirical biomolecular and water plane and out-of-plane components, and a “softer” buffered
forcefieldsissomethingthatthesimulationcommunityrequires, 14-7vdWform.However,themajordifferenceisreplacement
Feature Article J. Phys. Chem. B, Vol. 114, No. 8, 2010 2551
of the fixed partial charge model with polarizable atomic isusedtoaidinmergingtheshort-range“valence”termswith
multipolesthroughthequadrupolemoments.Oneadvantageof the long-range “nonbonded” interactions. For dihedral angles
theAMOEBAmodelisitsemphasisonreplicationofmolecular involving two joined trigonal centers, such as the amide bond
polarizabilities and electrostatic potentials, instead of just of the protein backbone, a Bell torsion20 functional is applied
interactionenergies.Theuseofpermanentdipolesandquadru- inadditiontotheregulartorsionalterms,whereφusedineq6
poles allows accurate reproduction of molecular electrostatic isthedihedralanglecomputedfromthep-orbitaldirectionsat
potentials, and fine-tuning of subtle directional effects in thetwotrigonalcenters,ratherthanfromtheusualbondvectors.
hydrogen bonding and other interactions. The inclusion of The rotational barrier around the amide bond is much higher
explicit dipole polarization allows the AMOEBA model to thanforasinglecovalentbond,andthebiggerbarrierislargely
respondtochangingorheterogeneousmolecularenvironments, duetothedoublebondnatureoriginatingintheoverlapofthe
andallowsdirectparametrizationagainstgasphaseexperimental adjacentp-orbitals.UseoftheBelltorsionallowsappropriately
dataandhigh-levelquantummechanicalresults.TheAMOEBA increased flexibility of atoms bonded to trigonal centers (e.g.,
model also presents a consistent treatment of intra- and aromatichydrogenatoms).21Thetorsionalparametersarerefined
intermolecularpolarizationthatisachievedthroughaphysically after the nonbonded parameters are determined with the hope
motivated damping scheme for local polarization effects.17 A thattheimprovedAMOEBAintramolecularelectrostaticmodel
further attractive aspect of AMOEBA is its use of multipole willleadtoamore“physical”balancebetweenthelocal(vdW
moments derived directly from ab initio quantum mechanical + electrostatic + torsional) and long-range (vdW + electro-
electrondensitiesforsmallmoleculesandmolecularfragments. static) interactions in the conformational energy.
ThedesigngoalforAMOEBAhasbeentoachievea“chemical vanderWaalsInteractions.Thepairwiseadditivevander
accuracy” of 0.5 kcal/mol or better for small molecule and Waals(vdW)interactioninAMOEBAadoptsthebuffered14-7
protein-ligand interactions. We describe the functional form functional form22
of the AMOEBA force field below and provide the current
( ) ( )
standardparametersetforsmallmoleculesandproteinsinthe
1.07 7 1.12
StiuopnpaorretinggivIennfoirnmraetfio1n3,.whilefurtherdetailsofitsparametriza- Uvdw(ij) ) εij Fij + 0.07 Fi7j + 0.12 - 2 (7)
Short-RangedValenceInteractions.TheAMOEBAmodel
includesfullintramolecularflexibility.Foratomsdirectlybonded where ε in kcal/mol is the potential well depth and F ) R /
(1-2)andseparatedbytwobonds(1-3),thecovalentenergy R0, wheirje R in angstroms is the actual separation beitjweeniji
ij ij
isrepresentedbyempiricalfunctionsofbondlengthsandangles. andjandR0istheminimumenergydistance.Forheterogeneous
ij
Thefunctionalformsforbondstretching(eq2),anglebending atom pairs, the combination rules are given by
(eq 3), and the coupling between the stretching and bending
(eq 4) are those of the MM3 force field,18 and include an
accounting of anharmonicity through the use of higher-order ε ) 4εiiεjj and R0 ) (Ri0i)3 + (Rj0j)3 (8)
deviations from ideal bond lengths (b0) and angles (θ0): ij (ε1/2 + ε1/2)2 ij (R0)2 + (R0)2
ii jj ii jj
Ubond ) Kb(b - b0)2[1 - 2.55(b - b0) + Thebuffered14-7functionyieldsaslightly“softer”repulsive
(7/12)2.55(b - b )2] (2) region than the Lennard-Jones 6-12 function but achieves a
0
steeperrepulsionatveryshortrangethantypicalBuckingham
exp-6 formulations. The buffered 14-7 form was considered
U ) K (θ - θ )2[1 - 0.014(θ - θ ) + superior,asitprovidesabetterfittogasphaseabinitioresults
angle θ 0 0 andliquidpropertiesofnoblegases.22TheAMOEBAvander
5.6 × 10-5(θ - θ0)2 - 7.0 × 10-7(θ - θ0)3 + Waals parameters are derived by fitting to both gas phase and
2.2 × 10-8(θ - θ )4] (3) bulk phase experimental properties.
0 Each atom in AMOEBA possesses a vdW site. For non-
hydrogenatoms,thesiteislocatedatthepositionoftheatomic
U ) K [(b - b ) + (b′-b ′)](θ - θ ) (4) nucleus. For a hydrogen atom connected to an atom X, it is
bθ bθ 0 0 0 placedalongtheH-Xbondsuchthatthedistancebetweenthe
atomXandthevdWsiteofHisapercentageofthefullbond
U ) K (cid:2)2 (5) length,namely,the“reductionfactor”.Applicationofreduction
oop (cid:2) factors to shift hydrogen sites off of the nuclear centers dates
fromearlyworkbyStewartetal.,23andX-raystructuralanalyses
wherethebondlength,borb′,andbondangle,θ,andenergies ofglycylglycineandsulfamicacidalsosupportthisview.24A
areinunitsofÅ,degrees,andkcal/mol,withtheforceconstants, similar approach is used in MM3 and other force fields from
K,givenincorrespondingunits.Inaddition,aWilson-Decius-Cross the Allinger group.18 The use of a reduction factor was found
function is used at sp2-hybridized trigonal centers to restrain tosimultaneouslyimprovethefittoaccurateQMwaterdimer
theout-of-planebending(eq5),19whereforsequentiallybonded structures and energies for several configurations.
centersi,j,k,andl,(cid:2)referstotheanglebetweenthejlvector Permanent Electrostatic Interactions. The electrostatic
and the ijk plane. energy in AMOEBA includes contributions from both perma-
A traditional Fourier expansion (a 1-fold through 6-fold nentandinducedmultipoles.Thepermanentatomicmultipoles
trigonometric form) torsional functional (PAM) at each atomic center include the monopole (charge),
dipole, and quadrupole moments
∑
U ) K [1 + cos(nφ ( δ)] (6)
torsion nφ M ) [q,µ ,µ ,µ ,Q ,Q ,Q ,...,Q ]t (9)
n i i ix iy iz ixx ixy ixz izz
2552 J. Phys. Chem. B, Vol. 114, No. 8, 2010 Ponder et al.
where R is the rotation matrix transforming the local into the
global reference frame.27
ElectronicPolarization.Electronicpolarizationreferstothe
distortionofelectrondensityundertheinfluenceofanexternal
field. It represents a major contribution to the overall many-
bodyenergeticdescriptionofmolecularclustersandcondensed
phases,eventhoughtherearesituationswhereothercontribu-
tionsrelatedtodispersionandrepulsionarenotnegligible.28In
Figure1. Anexampleofthe“z-then-bisector”localframedefinition. AMOEBA,aclassicalpointdipolemomentisinducedateach
Shownformethylaminethatensuresthattheatomicmultipolesremain
polarizable atomic site according to the electric field felt by
constantwithtimewithinthislocalreferenceframe.
that site. Molecular polarization is achieved via an interactive
induction model with distributed atomic polarizabilities based
whereqiisthepointchargelocatedatthecenterofatomi,µis on Thole’s damped interaction method.17a This interactive or
thedipole,andQisthequadrupole,allinCartesianrepresenta- mutual induction scheme requires that an induced dipole
tion, and t is[ the transpose. In the Cartesia]n polytensor produced at any site i will further polarize all other sites, and
formalism,25 the interaction energy between atoms i and j such mutual induction will continue until the induced dipoles
separated by rji is represented as Upeleercm(rij) ) MiTTijMj where at each site reach convergence. One key aspect of Thole’s
approachisdampingofthepolarizationinteractionatveryshort
∂ ∂ ∂ range to avoid the so-called polarization catastrophe, a well-
1 ∂x ∂y ∂z L known artifact of point polarizability models. The damping is
j j j
effectively achieved by smearing one of the atomic multipole
∂ ∂2 ∂2 ∂2
L moments in each pair of interaction sites (the result is
∂xi ∂xi∂xj ∂xi∂yj ∂xi∂zj ( ) independentofwhichoneissmeared).29Thesmearingfunction
1
T ) ∂ ∂2 ∂2 ∂2 (10) for charges adopted by AMOEBA has the functional form
ij L r
∂y ∂y∂x ∂y∂y ∂y∂z ji
i i j i j i j 3a
∂ ∂2 ∂2 ∂2 F ) exp(-au3) (12)
L 4π
∂z ∂z∂x ∂z∂y ∂z∂z
i i j i j i j
M M M M O whereu)r/(RR)1/6istheeffectivedistanceasafunctionof
linear separaijtioni rj and atomic polarizabilities of sites i (R)
There are typically five independent quadrupole components and j (Rj). The facitjor “a” is a dimensionless width parameteir
duetosymmetry(QR(cid:5))Q(cid:5)R)andtheuseoftracelessmoments ofthesmearedchargedistribution,andeffectivelycontrolsthe
(ΣQRR)0).Furthermore,theµy,Qxy,andQyzcomponentsare dampingstrength.Correspondingdampingfunctionsforcharge,
zeroexceptforchiralatomssuchasthebackboneCRinamino dipole,andquadrupoleinteractionswerederivedthroughtheir
acids. Therefore, most nonchiral atoms will carry six unique, chain rule relationships.13b
permanent electrostatic multipole parameters. The Thole model has the advantages of simplicity and
AspreviouslydescribedfortheAMOEBAwatermodel,the transferability, as evidenced by the fact that it reasonably
dipole and quadrupole are defined with respect to a local reproducesthemolecularpolarizabilitytensorofnumeroussmall
reference frame defined by neighboring atoms.13b A new “z- moleculesusingjustoneisotropicatomicpolarizabilityforeach
then-bisector” local frame definition has been developed for element,plusauniversaldampingfactor.17aHowever,therehas
atoms with single lone pairs such as the sp3 nitrogen. An beencontroversyastowhetherpolarizabilitydecreases,andif
example of this new frame is given for the N atom in so to what extent, when a molecule moves from gas to
methylamineinFigure1.Inprinciple,thechoiceofframeshould condensed phase. Morita recently estimated that water polar-
respect local symmetry such that axes are placed along major izabilitydecreasesby7-9%,30reducedfrom13-18%reported
chemical determinants. As the molecules vibrate, rotate, and inanearlierpublication.31Mennucietal.showedthattheeffect
diffuse over the course of a dynamic simulation, the atomic of Pauli exclusion is to reduce the dipole polarizability of a
multipoles remain constant with respect to the local frame soluteby2%.Incontrast,GubskayaandKusiliksuggestedan
definition. increase of the polarizability of water in condensed phases.32
Atomic multipole moments in this study are derived from In light of uncertainty in theoretical estimates of liquid polar-
ab initio calculations of the small molecules using Stone’s izability,wehavechosentousethesameatomicpolarizability
distributedmultipoleanalysis(DMA).26Convergencetoreason- valuesforbothgasandcondensedphase.Theresultingaverage
ablechemicalaccuracygoalsof0.5kcal/molrequiresinclusion dipolemomentofAMOEBAliquidwater,usingaconstantgas
oftermsthroughquadrupolemoments.Alternativeapproaches, phase polarizability value, is 2.8 D, only slightly lower than
such as electrostatic potential fitting and electron density recentquantummechanicalestimatesof2.95D.33Furthermore,
partitioning,havealsobeenexplored.13aFormoleculessuchas in the AMOEBA polarization model, the damping factor
alanine dipeptide that possess conformational degrees of free- providesanothercontrolovertheabilityofanatomtopolarize;
dom, an extra step is necessary to obtain the conformation- theuniversaldampingfactoradoptedbyAMOEBAisa)0.39,
independent permanent atomic multipoles (PAM), as will be which effectively leads to a stronger damping and less short-
discussedintheintramolecularpolarizationsectionbelow.The range polarization than the original value of 0.572 suggested
original DMA-derived multipoles26a are converted to the final by Thole. We have kept the same atomic polarizabilities (Å3)
electrostatic parameters in the corresponding local frame for givenbyThole,i.e.,1.334forcarbon,0.496forhydrogen,1.073
each atom type: for nitrogen, and 0.837 for oxygen. The only exception is for
carbon and hydrogen in aromatic rings, where we found the
useofsomewhatlargervaluesgreatlyimprovesthemolecular
Mlocal ) R-1M (11) polarizability tensor of benzene and polycyclic aromatics.
Feature Article J. Phys. Chem. B, Vol. 114, No. 8, 2010 2553
Validation Studies against Electronic Structure
Calculations
Anecessaryifnotsufficientconditionforrobustperformance
of a force field is its ability to reproduce or predict relative
conformationalenergiesofmodelsystemscomplexenoughto
contain realistic features but simple enough to be treated by
electronicstructuremethods34thatcanyieldreliablebenchmark
results. First principles electronic structure calculations can
provide benchmarks of uncompromising accuracy for relative
energiesofmoleculesindifferentconformations,sincequantum
mechanics provides an essentially exact description of the
behaviorofelectronsinsmallmolecules.35However,inpractice,
approximate electronic structure calculations for larger mol-
eculessufferfromerrorsassociatedwiththeimperfecttreatment
ofelectroncorrelationsandtheuseofincompleteatomicorbital
basis sets, which can render the results too inaccurate to be
useful, or even worse, potentially misleading. Incomplete
treatments of electron correlation such as commonly used
density functional theory (DFT) methods omit dispersion
interactions that are very important in biological macromol-
ecules, while incomplete basis sets give rise to, for example,
intramolecular basis set superposition error (BSSE), which
favors compact relative to extended conformations.
In the work summarized here, we cannot claim to have
completely eliminated either of these problems, but we have
certainlyreducedsomelimitationsofearliercalculations,using
newalgorithmsandfastercomputers.Withrespecttoelectron
correlation, we have used the second order Møller-Plesset
(MP2)method,whichincludeslong-rangeelectroncorrelation
effects in a reasonably accurate manner, and then tested for
remaining errors by using local coupled cluster theory with a
smaller basis set. With respect to basis set errors, we have
performedcalculationswiththeDunningaugmentedcorrelation Figure2. EnergycorrelationsbetweenAMOEBAandMP2energies
consistent basis sets up to the aug-cc-pVQZ level, which we for(a)AMOEBAminimizedwater-sulfateanionclustersand(b)MP2
used with aug-cc-pVTZ results to perform an extrapolation to minimizedwater-sulfateanionclusters.Shownfor(H2O)nSO42-,where
n)3(X),4(9),and5(4).Correlationcoefficientsare0.88(n)3),
the complete basis set limit (TQ extrapolation). Comparisons
0.77(n)4),and0.79(n)5)forAMOEBAgeometriesandcorrelation
against smaller basis sets show that this level of theory is coefficientsare0.92(n)3),0.92(n)4),and0.90(n)5)forMP2
required to obtain reasonable convergence of MP2 relative geometries.
energies.Onthebasisofthesequantummechanicalbenchmarks,
weevaluatetheAMOEBAperformanceonnanosolvationand configurationsforthesamplingofanexpensiveabinitioenergy
conformational energetics of alanine tetrapeptide. function.Theprimaryprobleminusinghybridenergyschemes
Conformational Searching for Global and Low-Lying isthatwehavenoknowledgeorguaranteethatthedistribution
EnergyMinimaforWaterClusterSystems.Nanodropletsand of configurations generated with the lower quality energy
nanosolvation are interesting and demanding test cases for function overlaps sufficiently with the higher quality energy
modelingwaterviapolarizableforcefieldsbecausetheycontain function. However, we have found AMOEBA to be a reliable
watermoleculesindifferentextremesofenvironment,ranging generator of viable minima with sound energy ordering when
fromsurfacemoleculesexposedtovacuumtoburiedmolecules benchmarked against a reliable ab initio theoretical model. In
thatexperiencebulk-likeenvironment.Beyonddirectuseasa arecentstudyonn)3,4,and5water-sulfateanionclusters
“stress test” for polarizable force fields, the science of nano- (H O) SO 2-, we used replica exchange simulations over the
2 n 4
dropletsisinterestinginitself,sincethesespeciesareintermedi- temperature range from 140 to 500 K using the AMOEBA
ate between small clusters (20-30 molecules and below) that model, and all samples collected every 0.5 ps at every
are currently intensively studied by high-accuracy electronic temperature were energy minimized using the BFGS local
structuretheoryaswellasbeamexperiments,andsolvationin optimizationalgorithm.Samplingforallclustersizesconsidered
thebulkliquidsolvent.Theywillhavesomeofthefeaturesof appeared to be exhaustive, since all of the 10000 structures
water in confined regions, and may well exhibit interesting collectedforeachclustersizereducedtoasmallersetofupto
structuralmotifsthatlieinbetweensmallclusterbuildingblocks 200 local minima.
and the hydrogen-bonding patterns of the bulk. The behavior Figure 2a shows that the quantitative correlation between
of a solute in these clusters in terms of whether it appears on AMOEBA and the ab initio theory is very good (correlation
the surface or in the interior may have similarities to the coefficient,r2∼0.9)whilethequalitativecomparisonisexcellent
partitioning of solutes at interfaces. giventheagreementontheglobalminimumstructureforn)
Studies of the properties of the nanodroplets require an 3 and 4 that will likely dominate the nanosolvation properties
effectivesamplingtechniqueduetotheexponentiallyfastrise ofthissystemsize,andverycompetitivelowlyingminimafor
in the number of minima with cluster size. In hybrid energy n)5.Thelowestminimumenergystructuresdeterminedfrom
approaches, a cheap energy function is used to provide theempiricalpolarizablemodelwereinturnenergyminimized
2554 J. Phys. Chem. B, Vol. 114, No. 8, 2010 Ponder et al.
TABLE1: ComparisonofRelativeEnergies(kcal/mol)for to within 1 kcal/mol of the larger TQ results. One must also
Sulfate-WaterClusters(H2O)3SO42-a assumethattherewouldbeafurthershiftontheorderofperhaps
RIMP2/ RIMP2/ up to 0.1 kcal/mol upon further improvement of the basis set
isomer benchmark aug-cc-pVQZ aug-cc-pVDZ AMOEBA beyond the TQ extrapolation.
The second comparison of importance in Table 2 is with
1 0.00 0.00 0.00 0.00
standardDFTandthewidelyusedB3LYPfunctional,usinga
2 0.29 0.29 0.25 0.72
3 0.57 0.60 0.33 0.37 very large cc-pVQZ basis set. While DFT calculations cannot
4 0.65 0.54 0.32 1.31 be soundly extrapolated to the complete basis set limit, they
5 0.71 0.59 0.29 1.74 alsoconvergemorerapidlywithbasissetsizethanMP2theory,
6 2.38 2.68 2.08 2.63 so we can consider these results to be quite well-converged.
7 2.66 2.80 3.04 2.04
However, B3LYP is known to perform fairly poorly for
8 3.62 3.54 3.27 2.27
intermolecularinteractions(andhenceconformationalenergies),
aThe geometries of each cluster isomer were optimized at the as a result of limitations in its exchange and correlation
RIMP2/aug-cc-pVTZ level, and single point quantum mechanical functionals (for instance, it neglects dispersion interactions).
energies were calculated at a benchmark level (RIMP2/
Indeed,thiscausesseriousdiscrepanciesrelativetothebestMP2
aug-cc-pVQZ+∆CCSD(T)/6-31+G*), as well as RIMP2 using
results.Theoverallenergyrankingofconformersisquitepoor
aug-cc-pVQZ and aug-cc-pVDZ basis sets. AMOEBA results are
reportedfortheAMOEBAminimizedstructures. using this conventional DFT method, emphasizing its lack of
suitability as sources of benchmark conformational energies.
Thedevelopmentofnewfunctionalsthatimproveexchange
by the RI-MP2 level of theory using an augmented cc-pVDZ
functionalsandincludeempiricalvanderWaalscorrectionsis
basis set. Figure 2b shows that for minimized MP2 structures
likely to yield significantly improved performance. We assess
the quantitative correlation between AMOEBA single point
theroleofimprovedexchangewiththerange-separatedωB97
energies and the ab initio theory is still very good (r2 ∼0.8),
andωB97Xfunctionals,38andtheadditionaleffectofdispersion
showingthatAMOEBAgeometriesareinverygoodagreement
with the recently proposed ωB97X-D functional.39 Table 2
with the benchmark calculation.
showsconformerenergiesfortheωB97,ωB97X,andωB97X-D
For the smaller n ) 3 clusters, we can benchmark against
long-rangecorrectedfunctionals.Allofthemyieldasignificant
high-levelQMresults.Table1showstherelativeenergiesfor
improvementoverB3LYPandexhibitanexcellentagreement
theeightlowest-lyingconfigurations.TheRIMP2+∆CCSD(T)
withthebenchmarkenergies(correlationcoefficientsof0.910,
referenceenergieswereobtainedattheRI-MP2/aug-cc-pVQZ
0.932, and 0.908, respectively). Interestingly, however, the
leveloftheory,whichwerecorrectedattheCCSD(T)/6-31+G*
dispersion correction in ωB97X-D does not improve the
levelforhigher-ordercorrelationeffects.ComparingtheRIMP2/
performanceofthefunctionalinthepresenttestcase,suggesting
aug-cc-pVDZ and RIMP2/aug-cc-pVQZ results, we find basis
that intramolecular dispersion effects may be adequately cap-
seteffectsofupto0.4kcal/mol.Higher-ordercorrelationeffects
tured by other parts of the functional.
are on the order of up to 0.3 kcal/mol, as seen by comparing
RIMP2 and RIMP2+∆CCSD(T) results. This emphasizes the Finally,wereporttheoriginalLMP2resultsbutusingmore
tightlyconvergedgeometriesthanreportedoriginally,36aswell
importance of both effects, given that the energy differences
as the AMOEBA results which used LMP2 conformational
between most low-lying isomers are on the same order of
energiesofthealaninedipeptideaspartoftheparametrization
magnitude.WeseethattheAMOEBAforcefieldhasquitegood
of the AMOEBA protein model. It is interesting to see that
energy ordering of the isomers relative to the highest level of
AMOEBA (using AMOEBA relaxed geometries) gives a
theory,showingitsvalidityoutsidethequantumchemistrylevels
competitiveenergyrankingoveralltheconformationscompared
of theory used in the parametrization scheme reported in ref
to the RI-MP2 benchmark (r2 ∼0.88), comparable to that
24.
exhibitedbytheLMP2leveloftheory(r2∼0.95),andfarbetter
Electronic Structure Calculations of Conformational
Energies of the Alanine Tetrapeptide. Alanine tetrapeptide than conventional DFT (r2 ∼0.46). AMEOBA is essentially
competitivewiththenewgenerationdensityfunctionals,ωB97,
is a system which has at least several dozen low-lying
conformationalminimarangingfromglobulartoextended,and ωB97X, and ωB97X-D, which illustrates that it is very well
includeshydrogenbondingandpackinginteractionsthatmake balanced for polypeptide conformation energies.
it quite a rich biochemical system even in the gas phase. For From a biophysical viewpoint, one of the most important
thisreason,benchmarkcalculationsonalaninetetrapeptidefirst comparisonsisbetweentheextendedconformation(conformer
appeared roughly a dozen years ago.36 In this section, we 1) and a compact globular conformation with a tight hairpin
summarize,andinsomeinstancesextend,recentcalculations37 turn(conformer3).Thistypeofenergydifferenceisparticularly
that significantly improve the accuracy and reliability of the sensitivetobasissetconvergenceproblemsbecauselimitations
earlierbenchmarks.Table2containsrelativeenergiescalculated of the basis set will favor the globular conformation, where
with different popular electronic structure methods all using atomsinnonbondedcontactcanartificiallylowertheirenergy
geometriesoptimizedatthesameHartree-Fockleveloftheory by making fractional use of the functions on their nonbonded
with the 6-31G** basis set, for 27 conformations of alanine neighbors. This intramolecular basis set superposition error is
tetrapeptide using the labels reported in ref 36. Our best essentially absent in the extended conformation. As a result,
(benchmark) level of theory (MP2 with TQ extrapolation) as thebenchmarkextended-globularenergygap,E )3.56kcal/
gap
well as the MP2 theory with a less complete basis set (DT mol,isoverestimatedby∼1.3kcal/molattheDTextrapolated
extrapolation) are beyond originally published benchmarks at level. This emphasizes the importance of carrying out the
the lower double-(cid:7) level of quality.36 The comparison of the calculations to the largest feasible basis set size. Errors
first two columns of Table 2 indicates how troublesome associated with neglect of dispersion interactions, which are
obtainingfullyconvergedresultsisusinganelectroncorrelation relativelynonspecific,cansometimesapproximatelycancelout
methodsuchasMP2.TheDTleveloftheoryisalreadybeyond forconformationsofapproximatelysimilarcompactness.How-
mostliteraturecalculationsyetinsomecasesisnotconverged ever, the energy difference between extended and globular
Feature Article J. Phys. Chem. B, Vol. 114, No. 8, 2010 2555
TABLE2: ComparisonofBenchmarkRI-MP2CalculationsApproachingtheBasisSetLimitagainstOtherElectronic
StructureMethodsandAMOEBAfor27AlanineTetrapeptideConformationsa
conf. MP2/TQ MP2/DT ωB97/LP ωB97X/LP ωB97X-D/LP B3LYP/Q LMP2/cc-pVTZ(-f) AMOEBA
11 0.000 0.000 0.000 0.000 0.000 2.184 0.000 0.090
12 0.290 0.346 1.099 1.187 0.902 3.266 0.699 0.372
3 0.571 0.693 0.723 0.425 0.523 1.251 0.195 0.000
26 0.674 1.223 2.367 2.187 2.398 1.806 0.373 1.509
20 1.755 2.335 3.032 2.666 3.190 1.270 1.061 2.432
18 1.913 2.468 1.944 1.577 2.668 0.000 0.718 1.938
15 2.194 1.707 2.136 2.347 2.591 5.291 2.261 1.164
25 2.495 3.118 3.575 3.389 4.030 2.442 1.784 2.935
6 2.895 3.148 2.601 2.433 3.196 3.228 2.383 2.422
21 2.918 3.009 2.474 2.547 2.749 3.336 2.300 2.828
17 3.418 3.414 2.474 2.547 2.749 4.638 2.980 2.520
16 3.549 3.784 4.713 4.292 4.438 3.925 3.021 2.575
13 3.655 4.538 4.034 3.474 4.509 1.225 1.965 3.519
19 3.816 4.319 3.950 3.682 4.648 2.033 3.029 3.610
24 3.976 4.115 4.280 4.241 5.131 4.521 3.171 3.424
27 4.020 4.513 5.197 4.989 5.423 4.155 3.378 4.355
1 4.130 5.553 4.745 4.088 5.576 0.742 2.690 4.162
2 4.190 5.390 4.892 4.358 5.650 1.056 2.780 4.001
8 4.640 4.477 5.024 5.050 5.390 6.333 4.364 4.258
14 4.679 5.395 5.811 5.336 5.758 3.862 3.877 4.029
5 5.261 6.353 6.431 5.835 5.758 3.263 4.074 4.152
4 5.730 6.884 6.907 6.219 7.317 3.350 4.062 4.831
23 5.815 5.979 5.944 5.809 6.383 6.335 5.018 5.618
22 5.824 5.899 5.667 5.520 6.126 6.318 5.019 4.295
7 6.665 6.931 6.648 6.614 6.126 7.730 5.927 4.385
10 7.791 7.766 7.637 7.707 8.286 8.367 7.189 5.613
9 7.923 8.197 7.932 7.631 8.033 6.264 7.129 8.066
aAllrelativeenergiesareinkcal/molandgeometriesoptimizedattheHF/6-31G**level.AMOEBAresultsusedminimizedstructuresbased
ontheAMOEBAforcefieldforeachconformation.
TABLE3: EffectoftheLevelofTheoryUsedforGeometry
conformationalenergyatthehighestleveloftheory(RI-MP2/
OptimizationontheEnergyDifference(inkcal/mol)between
theExtendedandGlobularConformationsofAlanine TQ) depending upon the geometry that is used, with a new
Tetrapeptidea benchmarkvalueofEgap)4.994kcal/mol.Thereisashiftof
over2kcal/molbetweenHFandMP2geometries,withtheDFT
leveloftheoryfor
geometry in much closer agreement (0.65 kcal/mol) with RI-
geometryoptimization
energy MP2.EvenbetweensmallbasisHF(Table1)versusthelarger
evaluation RI-MP2/T B3LYP/T HF/T HF/6-31G**
basis results shown in Table 2, there is a shift of roughly 0.7
RI-MP2/TQ 4.994 4.414 2.884 3.559 kcal/mol. Against the new MP2 geometry benchmark, E is
gap
RI-MP2/DT 6.720 5.582 3.942 4.860 overestimatedby∼1.7kcal/molattheRI-MP2/DTextrapolated
B3LYP/Q -1.320 -0.093 -0.560 -0.590 level,whileDFTunderestimatesthegapbynowroughly6kcal/
aThe benchmark value of 4.994 kcal/mol for the energy gap is mol. By contrast, AMOEBA now undershoots the benchmark
highlightedinbold. result by ∼0.9 kcal/mol, showing that AMOEBA geometries
arethemostrobustwhencomparedtotheRI-MP2geometries
conformations is quite sensitive to the neglect of dispersion, and energies.
resulting in a calculated B3LYP E ) -0.51, a large error
gap Validation against Hydration Free Energies
thatunderestimatesthebenchmarkcalculationbyroughly4kcal/
mol. The corresponding gap measured by LMP2 is underesti- Theevaluationofhydrationfreeenergiesisanaturaltestof
matedby∼1.1kcal/mol,likelyduetobasissetsizelimitations anyforcefield,sinceitincorporatesmanychallengingaspects
andthelocalapproximationofthemodel.AMOEBAperforms ofaheterogeneouschemicalenvironmentthatarenotinvolved
thebestonthisbenchmark,overshootingtheRI-MP2/TQresult in the parametrization of the protein fragments or water force
by only ∼0.6 kcal/mol (close to the AMOEBA chemical fields by themselves.40 The AMOEBA solvation free energies
accuracygoalof0.5kcal/mol),althoughagainitisbasedona werecomputedusingafreeenergyperturbationprocedurebased
comparison using AMOEBA relaxed geometries and not the on three thermocycle steps and processed with the Bennett
HF/6-31G** geometries. acceptanceratio(BAR)method.41Foreachsmallmolecule,the
Theeffectofgeometryoptimizationontheextended-globular thermodynamic cycle corresponded to first solute discharging
gap is probed further with the calculations shown in Table 3, in a vacuum over 7 windows, followed by a soft core
wherelargebasissetgeometryoptimizationsattheMP2,DFT, modification of eq 7 to introduce the solute-solvent van der
andHFlevelsarecomparedviasinglepointenergycalculations Waalscouplingover16windows,andfinallysoluterecharging
using the three different sets of structures. While it has been in water over 7 windows. The statistical samples of the first
shownthattheMP2geometriesaresuperiortoHFgeometries, thermocyclestepinavacuumwerecollectedevery0.5psfrom
itiscommonlyassumed(generallyforgoodreason)thatDFT a10nsstochasticdynamicssimulationwithanintegrationtime
or HF structures are adequate, because errors in electron step of 0.1 fs, while the thermocycle steps in the condensed
correlation treatment cancel for small displacements of the phasewererunfor1nsintheNVTensemblewiththedensity
geometry. However, there are significant shifts in relative fixedat1.000gcm-3.Induceddipoleswereconvergedto10-5
2556 J. Phys. Chem. B, Vol. 114, No. 8, 2010 Ponder et al.
TABLE4: AccuracyofAMOEBASolvationFreeEnergiesforSmallMoleculesa,b
compound AMOEBA experiment compound AMOEBA experiment
isopropanol -4.21(0.34 -4.74 propane 1.69(0.17 1.96
methylether -2.22(0.38 -1.92 methane 1.73(0.13 1.98
HS -0.41(0.17 -0.44 methanol -4.79(0.23 -5.10
2
p-cresol -5.60(0.23 -6.61 n-propanol -4.85(0.27 -4.85
ethylsulfide -1.74(0.24 -1.14 toluene -1.53(0.25 -0.89
dimethylsulfide -1.85(0.21 -1.83 ethylbenzene -0.80(0.28 -0.79
phenol -5.05(0.28 -6.62 N-methylacetamide -8.66(0.30 -10.00
benzene -1.23(0.23 -0.90 water -5.86(0.19 -6.32
ethanol -4.69(0.25 -4.96 aceticacid -5.63(0.20 -6.69
ethane 1.73(0.15 1.81 methylsulfide -1.44(0.27 -1.24
n-butane 1.11(0.21 2.07 methylethylsulfide -1.98(0.32 -1.50
dinitrogen 2.26(0.12 2.49 imidazole -10.25(0.30 -9.63
methylamine -5.46(0.25 -4.55 acetamide -9.30(0.27 -9.71
dimethylamine -3.04(0.26 -4.29 ethylamine -4.33(0.24 -4.50
trimethylamine -2.09(0.24 -3.20 pyrrolidine -4.88(0.29 -5.48
aTheuncertaintyisthestatisticaluncertaintyintheBARfreeenergycalculation.Allunitsarekcal/mol.Whencomparedtotheexperimental
results, the RMS error for the 30 AMOEBA solvation free energies is 0.68 kcal/mol and the mean signed error is +0.14 kcal/mol.
bExperimentalvaluesarereportedinref42a,exceptforimidazoltakenfromref42b,N-methylacetamidetakenfromref42c,andmethylsulfide
andwatertakenfromref42d.
Dperstepperatomforsimulationsinavacuumand10-2Din moleculesforthe2009SAMPLexercise(http://sampl.eyesopen.
the liquid during the trajectory, and the energies of the com/).Alchemicalhydrationfreeenergycalculationstocompute
condensedphasesnapshots(savedevery0.5ps)werereevalu- the transfer free energy from 1 M gas phase to 1 M aqueous
ated with the induced dipole converged to 10-5 D. BAR was solutionwerecarriedoutinamannersimilartothatdescribed
thenusedtoestimatethefreeenergybetweentheneighboring in ref 44 using a preview release of Tinker 5 modified to add
steps, and the final free energy was taken as the sum over all a numerically computed analytical long-range dispersion cor-
windows. rection,45soft-coreformsoftheHalgrenpotential,andtheability
Table 4 reports the AMOEBA solvation free energies of to periodically evaluate potential energies at all alchemical
common small molecules found in biochemistry, including intermediates. In a vacuum and solvent, seven discharging
commonaminoacidsidechainanalogues,withcorresponding intermediatestateswereusedtoscalecharges,multipoles,and
statisticaluncertaintiesobtainedviaablockaveragingapplied polarizabilitiesbyfactorslambda,crudelyoptimizedtoreflect
to each simulation step, and the final statistical error bar is a the quadratic dependence of charging self-energies, while
sumoftheuncertaintiesoverallsteps.Whencomparedtothe torsionalbarrierswerecorrespondinglyscaledbylinearfactors.
experimental results,42 the rms error for AMOEBA solvation In solvent, a decoupling parameter λ was used to modify the
h
freeenergiesis0.68kcal/mol,withameansignederrorof+0.14 Halgrenpotentialshiftconstantsandwelldepthtomimicasoft-
kcal/mol. Calculated solvation free energies using traditional corepotentialatintermediatevaluesofλ .Vacuumsimulations
h
fixedchargeforcefieldstypicallyhaveanaveragermserrorof (dischargingonly)ateachalchemicalintermediatewererunfor
1.0-1.25 kcal/mol compared to available experiments for 5nsusingLangevindynamicswithacollisionrateof5/ps,with
similar sets of molecules and a general shift in solvation free energies at all alchemical states written every 10 ps. Solvated
energy with a mean error of approximately 1 kcal/mol,40 systems consisted of the solute molecule immersed in an
demonstrating that for chemical spaces similar to proteins isotropicboxcontaining850AMOEBAwater.Solvatedsimula-
AMOEBA offers significant improvement over corresponding tions (discharging and decoupling) for each alchemical inter-
fixed charge force fields. mediate were run for 300-600 ps using the Berendsen weak-
PredictionofSolvationFreeEnergiesfor2009OpenEye coupling algorithm46 for both thermal (coupling time 0.1 ps)
SAMPL Competition. The Statistical Assessment of the andvolume(couplingtime2ps)controlwhichareavailablein
ModelingofProteinsandLigands(SAMPL)blindchallengeis Tinker,thoughthedistributiongeneratedbyBerendsenshould
anassessmentofforcefieldsandsamplingmethodsforprotein approachthecorrectNPTensembleinthethermodynamiclimit.
andligandmodeling.Onepredictionaspecthighlightedinthe Potentialenergiesfromsolvatedsimulationswerecomputedat
first SAMPL contest in 2008 consisted of predicting 63 all alchemical intermediates and stored every 0.5 ps. Particle
vacuum-watertransferenergies.Anumberofresearchgroups mesh Ewald was employed with a real-space cutoff of 7 Å,
using fixed charge force field models with water represented interpolation order of 5, and a grid of 42 × 42 × 42 points.
explicitlycalculatedsolvationfreeenergiesusingstandardfree Dynamicswereintegratedusingthe“betterBeeman”algorithm
energyperturbationMDcalculations,andultimatelytheirblind with a time step of 1 fs.
prediction results were compared to available experimental Correlation times were computed for the potential energy
literaturenumbers.Theoverallperformanceoftheseapproaches history and the trajectories subsampled to produce a set of
gaveanrmserrorofover3kcal/molcomparedtotheSAMPL uncorrelatedsamples.Allrecordedsampleswereprocessedwith
reportedexperimentaldata.43Thegoaloftheseblindassessment the multistate Bennett acceptance ratio (MBAR)47 to estimate
approaches is to not criticize the underperformance of fixed free energies and uncertainties for each leg of the thermody-
chargeforcefieldsbuttobetterunderstandwhentheydowell, namic cycle corresponding to transfer from 1 M gas to 1 M
and when additional physics of the computational model is aqueoussolution:discharginginvacuum,decouplinginwater,
neededforpredictingmorechallengingclassesofcompounds. anddischarginginwater.Thefirst1nsofvacuumsimulations
The AMOEBA force field was used to predict vacuum-to- and50psofsolvatedsimulationswerediscardedtoequilibra-
watersolvationfreeenergiesof43drug-likeandotherorganic tion, and the remainder (up to 4 ns for vacuum simulations)
Feature Article J. Phys. Chem. B, Vol. 114, No. 8, 2010 2557
dynamics near chemically heterogeneous protein surfaces.49
ThesesimulationsaretightlycoupledtoX-raydiffractionand
quasi-elasticneutronscattering(QENS)perfomedonthesesame
systems at the same concentrations. Unlike a majority of
macromolecularsimulationsthatmodelasinglesolvatedprotein,
these studies included ∼30-50 individual peptides that can
interact with one another as well as the water molecules. The
abilitytoaccuratelymodeltheinteractionsofindividualpeptide
fragmentsinacrowdedsolutionisimportantforeventualstudies
of protein-ligand binding and protein-protein interactions,
wherein the proteins can form temporary and reversible
complexes. Hence, these peptide simulation studies represent
animportantbiologicalenvironmentwithwhichtotestanyforce
field.
For the fixed charge case, we used the AMBER ff039b all-
atomproteinforcefieldandpotentialparameterstomodelthe
NALMA and NAGMA solutes, and the rigid, nonpolarizable
TIP4P-Ewmodel50forthewater.Wehavechosenanonstandard
protein-water model combination because we know that
Figure 3. Comparison of the AMOEBA solvation free energies vs
reportedvaluesfromSAMPL2009.SeeTable5fordetails.Allunits transport properties of TIP4P-Ew are excellent over a large
arekcal/mol. temperaturerange,unlikethedefaultTIP3Pmodeltypicallyused
with biomolecular solutes. Unfortunately, we found that the
simulated solution structure with nonpolarizable force fields
analyzed with MBAR; each leg of the thermodynamic cycle
predicts too much aggregation of both the hydrophobic and
was processed individually, but all simulations within the leg
hydrophilicpeptidesolutes(Figure4),indisagreementwithour
wereusedtogethertoobtainthemostaccurateestimatesoffree
liquid diffraction experiments. This in turn frees up too much
energies and their uncertainties.
bulk-like water, so as to yield water diffusion constants that
Figure3showstheAMOEBApredictionagainsttheOpenEye
are faster and with an Arrhenius temperature dependence,
reported experimental literature values for a few classes of
compounds, while Table 5 reports all of the results submitted contradictingourquasi-elasticneutronscatteringexperiments.
toSAMPL2009.ItisevidentfromTable5thatAMOEBAdid However,whenwefixthesolutestoremainsolvent-separated
especially well in areas where traditional force fields failed, asthatdeterminedfromthestructuralexperiments,wefindthat
especially for very soluble molecules such as D-xylose and thesimulatedhydrationdynamicswiththenonpolarizableforce
D-glucose. In addition, the 2009 SAMPL experimental values fieldsareclosetoquantitativewithrespecttotheexperimental
forglycerolandcyanuricacidwerelaterrevised(aftersubmis- dynamicaltrendswithtemperatureforNAGMA(Figure5a)and
sion of this manuscript), bringing AMOEBA into far better NALMA. It is clear that reparameterization of a biomolecular
agreement with experiment. Other compounds such as the forcefieldsuchasAmberff03(orotherfixedchargeforcefields)
uracils, parabens, and NSAIDs have a range of reported toimprovesolvationpropertiesusingTIP4P-Ewisanimportant
experimental values. For example, the AMOEBA predictions direction for future nonpolarizable force field efforts.
for the uracils are between the SAMPL experimental values Due to the unphysical perturbation introduced by fixing the
(taken from Cabani et al.48) and other more recently reported solutes, we also performed the same simulations with the
experimental values, suggesting that experimental uncertainty AMOEBA polarizable force field.13b In contrast to the fixed-
is much greater than the SAMPL error bars that are typically chargesimulations,thepolarizableforcefieldnicelyreproduces
reported to be below 1 kcal/mol. anonaggregated,uniformdistributionofsolutesthroughoutthe
It is noteworthy that AMOEBA tended to do poorly on the volume(Figure4).Itappearsfromtheseresultsthattheability
polyhalogenatedcompounds,whichtypicallyhavelargeatomic of the peptides to respond dynamically to their electrostatic
polarizabilities on the halogen atoms, values that were not environment via polarization is important for reproducing a
derived in the original work by Thole. The AMOEBA force correct uniform mixture of peptides in water. Given the
field derived the atomic polarizabilities for the halogens by qualitativeimprovementinsolutionstructureusingtheAMOEBA
fittingtojustacoupleofmonohalogenatedorganicliquids.The model, we also compared the changes in water dynamics as a
results suggest that the reason AMOEBA underestimates the function of temperature against our experimental data. On the
solvationfreeenergyisthattheatomicpolarizabilitiesneedto basis of quasi-elastic neutron scattering (QENS) experiments,
increase.Thenitrocompoundswerealsoachallenge,withsome the amphiphilic NALMA peptide solution exhibits two trans-
evidenceoflargebondlengthchangesbetweengasphaseand lationalrelaxationsatlowtemperatures,whilethehydrophilic
liquid (as there are for amides!), as well as more complicated peptideshowsonlyasingletranslationalprocess,withtransport
“push-pull”polarizationthatisnotfullycapturedbythecurrent properties of water near both peptide chemistries being very
“simple” polarization model. suppressedwithrespecttobulkdynamics.49d,eThisisarealstress
testforanyforcefieldgiventherangeofdynamicaltrendsthat
Condensed Phase Structure and Dynamics dependonaminoacidchemistryandtemperature.Wenotethat
We have completed molecular dynamics simulations using weconvergedtheinduceddipolesverytightlyinordertoensure
nonpolarizable and polarizable protein force fields to contrast energy conservation in the NVE ensemble under which we
thewaterdynamicsnearhydrophilic,N-acetyl-glycine-methy- collectedtimecorrelationfunctionsforcalculatingthediffusion
lamide (NAGMA), and amphiphilic, N-acetyl-leucine-methy- coefficients.
lamide (NALMA), peptides as a function of temperature, as AMOEBA provides reasonable agreement with the experi-
models for understanding temperature dependent hydration mental temperature trends in regards to translational diffusion