akin 
 to developments spurred by the rapid expansion of computer 
 power around 2000 , the burgeoning capacity provided by programmable 
 graphics processing units ( gpus ) has extended the utility of molecular 
 simulations as a practical tool for assessing biophysical processes . notably , gpu - accelerated computing has enabled routine simulations 
 on the microsecond ( s ) time scale , a critical regime on which 
 biological processes including protein recognition , ligand binding , 
 and protein conformational changes occur . access to these longer time scales may reveal flaws in the simulation 
 models that were not previously apparent , driving refinements and 
 leading toward improved predictive power of these models . historically , 
 efforts in force field development have been largely 
 focused on selecting only a subset of the parameters in a complex 
 model for reoptimization , e.g. , reoptimizing certain torsion parameters 
 while keeping the set of atomic charges fixed to improve the accuracy 
 in modeling particular behaviors , while retaining what is already 
 successful . however , such efforts are limited by the accuracy of the 
 unoptimized parameters . however , these refinements may be compensating for deficiencies 
 in the modeling of electrostatic and nonbonded interactions that limit 
 the maximum attainable accuracy of the model . in addition , contemporary 
 force fields have a tendency to borrow from a similar set of values 
 for bond lengths , angles , and atomic radii that were fit many years 
 ago , leading to interdependencies that may be difficult to untangle 
 when optimizing only a portion of the parameters . more recently , 
 semiautomated schemes have been developed to simultaneously 
 optimize hundreds of parameters , thereby enabling the rapid development 
 of new force fields . in particular , the force balance and implicitly 
 polarized charge ( ipolq ) methods have yielded the amber 
 ff15fb and ff14ipq force fields , respectively . these methods rely largely on automated tools for parameter optimization , 
 but still require some amount of manual intervention in the form of 
 fitting set composition or user - specified settings of the fitting 
 algorithm . the engine behind the ipolq workflow is the mdgx module 
 of the amber software package . this module 
 combines its molecular dynamics ( md ) facility with linear algebra 
 routines for solving least - squares problems , manages extensive bookkeeping 
 to organize parameters , provides user control over the fitting process , 
 and interprets statistics to aid in further refinements . the mdgx 
 module contains charge and torsion fitting routines that were built 
 throughout the development of amber ff14ipq , the first complete protein 
 force field based on the ipolq scheme . here , we have developed the new amber ff15ipq force field 
 using 
 the ipolq workflow . the original motivation for the development of 
 ff15ipq was to tackle concerns that its predecessor , ff14ipq , overestimates 
 the stability of salt - bridge interactions  a limitation shared 
 with many other contemporary force fields . however , in contrast to recently developed variants of force fields 
 that address such concerns , ff15ipq is far from a limited adjustment of its predecessor . rather , 
 ff15ipq is a complete rederivation , comprising new atomic charges , 
 a greatly expanded torsion parameter set with several new atom types 
 to decouple distinct amino acids , and new backbone angle bending terms . in addition , whereas ff14ipq employed the tip4p - ew model for the solvent 
 in the ipolq scheme , ff15ipq uses spc / eb , a recently developed three - point water model that yields 
 more - accurate rotational diffusion for proteins in solution . the use of a three - point water model instead 
 of a four - point water model , leaving out a virtual site , affords a 
 modest improvement in the speed of cpu - based simulations and a larger 
 acceleration in computing under the amber gpu engine . in addition , the more - accurate rotational diffusion afforded 
 by spc / eb opens new avenues for validating ff15ipq through 
 direct calculation of nmr relaxation parameters . with the aid of gpus 
 and the amber gpu engine , we have extensively validated the force 
 field by running md simulations of peptide and protein systems on 
 the s time scale , yielding over 200 s of aggregate simulation 
 time . we expect ff15ipq , or a close relative , to be valuable 
 for a long 
 time , even as we explore more expensive alternatives by adding virtual 
 sites to both the protein and the standard water model . in addition , 
 we are working to apply the mdgx workflow to other classes of biopolymers 
 such as carbohydrates and nucleic acids , and to small organic molecules . we hope that the sweeping reoptimization made possible by mdgx and 
 tools similar to it will inspire initiatives with other force fields 
 and create complete chemical representations with predictive power 
 in biomolecular simulations . of future interest will be comparisons 
 to contemporary force fields that have been developed in the traditional 
 manner such as amber ff14sb , opls - aa / m , and charmm36 , as well as those developed using 
 alternative sweeping reoptimization schemes , such as amber ff15fb . the implicitly polarized 
 charge ( ipolq ) method is a protocol for 
 parametrizing fixed - charge force fields for solution - phase simulations 
 that is comprised of two main components , implemented in the mdgx 
 program of ambertools . the first component 
 is a protocol for deriving nonpolarizable atomic charges that implicitly 
 represent the energy of polarization by the presence of a solvent 
 such as water . the ipolq charge derivation 
 draws on approximations of dipole interactions in an external electrostatic 
 field to arrive at the optimal nonpolarizable representation of a 
 solute s atomic charges in the presence of a solvent such as 
 water : precisely halfway between the charges that would reproduce 
 the solute s electrostatic field in the gas phase and those 
 that would reproduce the solute s electrostatic field after 
 solvent - induced polarization . this averaging 
 comes about from the fact that the energy of a set of polarizable 
 dipoles in an external field is identical to the energy of a set of 
 fixed dipoles whose polarizations are halfway between the field - polarized 
 dipoles and their gas - phase counterparts . ipolq fits such fixed charges 
 by applying the restrained electrostatic potential ( resp ) method , using a pair of representations of the solute s 
 electrostatic field corresponding to the vacuum and solution phases . 
 while the former representation is straightforward to obtain from 
 qm calculations , the latter is computationally unfeasible using a 
 pure qm representation . instead , the ipolq method represents the polarizing 
 solvent reaction field potential ( srfp ) in its qm calculation using 
 a field of point charges , derived from an md simulation in which water , 
 represented by the model with which the solute will ultimately be 
 simulated , moves in equilibrium around the fixed solute . atomic charges are subsequently fit to reproduce 
 the qm electrostatic potential at a set of grid points surrounding 
 the molecule . as described previously , grid points are selected within the first and second solvation shells , 
 with the inner boundary defined by excluding points for which the 
 energy of the lennard - jones interaction between the solute and a probe 
 representing the water model exceeds a selected maximum cutoff . while , 
 in this work , we have applied equal weights to all of the selected 
 grid points , others have found that more consistent charges may be 
 obtained by applying a weighting function to de - emphasize points close 
 to or distant from the solute . such improvements 
 will be investigated as the ipolq method is applied to other classes 
 of molecules beyond peptides and proteins . the second component 
 of the ipolq method is an extension for the fitting of bonded parameters 
 that accounts for the discrepancy between the desired solution - phase 
 conformational preferences and vacuum - phase qm calculations . this 
 is accomplished by fitting a pair of solute charge sets : one appropriate 
 for the vacuum phase ( qvac ) , and the other 
 for the solution phase ( qsolv ) . in the 
 presence of the qvac charge set , the force 
 field s bonded parameters are fit to reproduce the relative 
 vacuum - phase qm energies of a diverse set of solute conformations . 
 in subsequent simulations in the solution phase , these same bonded 
 parameters are paired with the polarized qsolv charge set with the intention that the difference in the charge 
 sets would account for the difference in solute conformational preferences 
 between the vacuum and solution phases . in contrast to the standard resp method 
 of fitting atomic charges for amber force fields , the ipolq method explicitly considers the influence of 
 the water model on the solute s charge distribution . while the atomic charges of the ff14ipq force 
 field were fit using the tip4p - ew water model , we have elected to fit the charges of ff15ipq using the 
 spc / eb water model . this recently 
 developed water model offers two advantages : in addition to the reduction 
 in computational cost that is obtained by switching from a four - point 
 water model to a three - point water model , this water model has been 
 parametrized to yield accurate rotational diffusion of solvated proteins . a key advantage to performing simulations with accurate rotational 
 diffusion is the ability to directly calculate the nmr relaxation 
 parameters n r1 and r2 and n  h heteronuclear 
 noe . these parameters provide information about fast dynamics ( picosecond 
 ( ps ) or nanosecond ( ns ) scale ) of individual backbone n  h bond 
 vectors within a protein , potentially offering a powerful means with 
 which to validate md simulations . in 
 principle , these nmr relaxation parameters may be calculated directly 
 from an md trajectory from the autocorrelation functions of the backbone 
 n  h vectors . in practice , however , the poor reproduction of 
 protein rotational diffusion in md simulations using popular water 
 models such as tip3p limits the utility of such calculations . this limitation has historically been addressed using approaches 
 such as model - free analysis that attempt to separate the global rotational 
 diffusion of the protein from its residue - specific internal dynamics , 
 comparing only the internal dynamics between experiment and simulation . however , these approaches require extensive fitting of models to the 
 experimental data , and further require that the global and local dynamics 
 occur on separable time scales , limiting their applicability to highly 
 flexible systems such as disordered peptides and proteins . therefore , 
 it would be preferable for md simulations to yield accurate rotational 
 diffusion , such that the simulated and experimental relaxation parameters 
 can be compared directly . our decision to fit the charges of 
 ff15ipq to the spc / eb water model was based on preliminary 
 tests in which we ran a series 
 of 24 simulations of the proteins gb3 , ubiquitin , and binase using 
 two different force fields ( amber ff99sb - ildn and charmm22 * ) paired with four different water models ( tip3p , tip4p - ew , tip4p - d , 
 and spc / eb ) ( see figure s1 in the supporting information ) . consistent with prior published 
 work , tip3p and tip4p - ew yielded rotational diffusion 
 significantly faster than experiment , and spc / eb yielded 
 the most accurate result . earlier versions of mdgx were capable of fitting 
 angle stiffnesses 
 alongside torsions , but recent advances proposed by vanommeslaeghe 
 et al . permit the calculation of both the optimal stiffness constant 
 and equilibrium value from the same linear least - squares problem . the strategy generalizes work by hopkins and 
 roitberg , representing the parabolic 
 angle as the sum of two parabolic basis functions and solving for 
 both scaling coefficients to interpolate the optimal parameters . basis 
 functions were chosen such that their minima lay at 0.2 radians 
 from the equilibria of the original angles , which had been inherited 
 from the amber ff94 force field . as explained 
 by vanommeslaeghe et al . , the optimized angle s stiffness is 
 given by the sum of coefficients solved for each basis function , and 
 its equilibrium is given by the average of the two basis functions 
 minima weighted by their coefficients . restraints on the optimized 
 angle parameters follow from these definitions : a restraint equation 
 setting the sum of the two coefficients ci,1 and ci,2 to a target value ki , such as the stiffness 
 of the original angle in the input force field , will harmonically 
 penalize solutions which depart from the original stiffness value : a similar restraint on the ratio of the two 
 coefficients can be used to penalize solutions that depart from the 
 original equilibrium value ti , if the 
 minima of the basis functions scaled by ci,1 and ci,2 are bi,1 and bi,2 , respectively : as explained in the previous study , these 
 restraint equations will have a more pronounced 
 effect if the dataset contains only 10 data points rather than 1000 
 data points . therefore , both sides of each restraint equation were 
 scaled by a user - defined constant ( scl ) times the 
 number of instances in which each optimizable angle appeared in the 
 dataset ni , analogous to the scaling constant 
 applied to torsion restraints . the scaling constants may appear to 
 have no effect on the solution to the equations when either these 
 constants are present on both sides of the equation or one side of 
 the equation is zero . however , since the least - squares fit finds an 
 approximate solution to each equation , the scaling constants do , in 
 fact , influence the relative importance of each restraint . in addition , 
 because of the fact that these restraints penalize numerical deviations 
 from the target values but angle stiffnesses and equilibria are expressed 
 in different units by numbers of different scale , a separate scaling 
 factor ( cpl ) was introduced to control the way in 
 which restraints on the equilibria scale , relative to restraints on 
 the stiffness . for example , an cpl of 57 
 applied to restraints on equilibria would penalize 1 deviations 
 from the original value by the same amount as a 1 kcal/(mol rad ) deviation from the original stiffness constant . after some 
 experimentation , however , we found that , in our very large and heterogeneous 
 datasets , much smaller values of cpl ( 0.51.0 ) 
 result in the best fits to the data while partitioning the changes 
 between equilibria and stiffness constants . alongside 
 the fitting of torsions and angles , several new atom types were added 
 to ff15ipq in order to more accurately capture residue - specific conformational 
 preferences . most protein force fields use gly and ala as templates 
 to develop backbone torsion parameters that are then inherited by 
 other residues . the amber ipolq force fields adopt an unconventional , 
 concerted approach in which  ,  ,  ,  , 
 and all other torsions are simultaneously fit to the conformational 
 preferences of all residues in which they appear . while the resulting 
 backbone torsions therefore consider the conformational preferences 
 of residues other than gly and ala , their overall accuracy may decrease 
 as different residues pull the parameters in different directions . 
 within the context of the ipolq fitting method , a set of backbone 
 torsion parameters that more accurately capture the conformational 
 preferences of different residues may be obtained by introducing new 
 atom types , creating decoupled classes of backbone torsions , which 
 are applied to subsets of residues . during the development of 
 ff14ipq , three such classes were introduced : one for gly , one for 
 pro , and one for all other residues . in order to decouple the  
 and  torsions of gly , which lacks c and therefore has 
 no  and  torsions , a unique c 
 atom type was assigned . similarly , the 
 backbone torsions of pro were decoupled by assigning a unique atom 
 type for the backbone n. this additional atom type not only created 
 unique  ,  , and  terms for pro , but also 
 a set of separate  and  torsions for residues 
 preceding pro , thereby enabling the force field to capture the unique 
 conformational preferences of these contexts . the remaining residues were further divided into three subclasses 
 based on their c types , which determine the applied  
 and  torsions . the c types of ff14ipq yielded 
 four subclasses : ( i ) flexible positively charged residues ( arg , lys ) , 
 ( ii ) residues whose c atoms are bonded to two heavy atoms and 
 whose side - chains are not aromatic ( asn , asp , cys , gln , glu , leu , 
 met , ser ) , ( iii ) residues whose c atoms are bonded to three 
 heavy atoms ( ile , thr , val ) , and ( iv ) all other residues ( ala , his , 
 phe , trp , tyr ) . notably , this last subclass yielded  
 and  torsions shared between ala and the bulky aromatic 
 residues . this unusual coupling was a consequence of the force field s 
 lineage from ff12sb , where refitting of x1 torsions alongside 
 fixed  ,  ,  , and  did not 
 require a unique c type for the aromatic residues , which already 
 have unique c types that yield unique x1 . in order to further improve the accuracy 
 of residue - specific conformational 
 preferences in ff15ipq , several new atom types were added to further 
 decouple the backbone torsion parameters of different residues , leading 
 to a total of five backbone classes . in order to restrict each class 
 of backbone torsions to a single set of scaled 14 electrostatic 
 terms , negatively charged ( asp , glu ) and positively charged ( arg , 
 lys ) residues have been given unique c types , decoupling their 
  ,  ,  , and  from those 
 of the neutral residues . while the backbone n of pro was decoupled 
 in ff14ipq , it retained a shared  torsion ; to break 
 this dependency , pro has now been assigned a new c atom type . finally , the coupling between ala and the bulky aromatic residues 
 has been removed by assigning his , phe , trp , and tyr a unique c 
 type , decoupling their  and  terms from 
 ala . this decoupling divides the neutral residues into four subclasses : 
 ( i ) ala , ( ii ) residues whose c atoms are bonded to two heavy 
 atoms ( asn , cys , gln , leu , met , ser ) , ( iii ) residues whose c 
 atoms are bonded to three heavy atoms ( ile , thr , val ) , and ( iv ) bulky 
 aromatic residues ( his , phe , trp , tyr ) . the backbone torsion classes 
 of the 28 residue forms supported by ff15ipq are listed in table s1 in the supporting information . to compare the accuracy of ff15ipq in modeling the stability of protein 
 salt bridges to its predecessor ff14ipq and contemporary force fields , 
 we simulated the association of three pairs of oppositely charged 
 amino acid side - chain analogues : guanidinium cation / acetate anion 
 ( arg / asp ) , butylammonium cation / acetate anion ( lys / asp ) , and imidazolium 
 cation / acetate anion ( his(+)/asp ) . for comparison , such simulations 
 were carried out using the polarizable force fields charmm drude-2013 
 and amoeba , in addition to other fixed - charge 
 force fields that we had previously tested using these three model 
 systems . simulations with ff15ipq , 
 ff14ipq , and amoeba were carried out using the amber 15 software package , while those with charmm drude-2013 were run 
 with namd 2.10.0 , following a protocol analogous 
 to that used for the previously evaluated fixed - charge force fields 
 ( full details are provided in the supporting information ) . systems were constructed to be consistent 
 with the experimental conditions under which the association constants 
 ( ka ) of guanidinium acetate and butylammonium 
 acetate have been measured , i.e. , each 
 system consisted of 100 molecules of cation ( guanidinium , butylammonium , 
 or imidazolium ) , 2 molecules of acetate , and 98 chloride counterions 
 solvated by 18  000 water molecules . for the fixed - charge 
 force fields , parameters of the side - chain analogues were based on 
 those of the complete amino acids . for the charmm drude-2013 polarizable 
 force field , parameters of guanidinium , imidazolium , and acetate were 
 those distributed alongside the force field . since methylammonium rather than butylammonium was used as the analogue 
 of lys during the development of drude-2013 , the butylammonium acetate system was not tested with this force 
 field . for the amoeba force field , parameters of guanidinium , imidazolium , 
 and acetate were generated using the poltype derivation protocol ( details 
 are provided in the supporting information ) . as it was with charmm drude-2013 , 
 the butylammonium acetate system was not tested with amoeba . for all of the simulations mentioned above , the probability that 
 an acetate molecule was bound to one or more cation molecules was 
 calculated by assigning each pair to either the bound or the unbound 
 state . for each force field and pair of side - chain analogues , definitions 
 of the unbound and bound states were based on the potential of mean 
 force ( pmf ) as a function of the minimum distance between nitrogen 
 atom(s ) of the cation and the oxygen atoms of acetate . in particular , 
 the cutoff between the bound and unbound states was defined as the 
 point of inflection between the free energy minimum of the bound state 
 ( 2.53  ) and the free - energy maximum , which corresponds 
 to the desolvation barrier ( 33.5  ) . pairs whose 
 minimum n  o distances were below this cutoff were assigned 
 to the bound state , while those beyond were assigned to the unbound 
 state . in addition to species in which a single acetate molecule was 
 bound to a single cation molecule , forming a 1:1 complex ( e.g. , the 
 guanidinium / acetate complex ) , species in which acetate was bound to 
 two or more cation molecules ( e.g. , the 2:1 diguanidinium / acetate 
 complex ) were observed and counted separately . the atomic charges of ff15ipq 
 were fit using the ipolq module of mdgx , as described previously for 
 ff14ipq . during charge fitting , each 
 amino acid was represented by a blocked dipeptide including acetyl 
 ( ace ) and n - methylamide ( nme ) caps ; terminal forms 
 were represented by omitting one of the blocking groups , while the 
 disulfide form of cysteine ( cyx ) was represented by a pair of dipeptides 
 linked by a disulfide bond . to expand on the set of amino acids and 
 protonation states that were supported by ff15ipq , atomic charges 
 were also derived for the following : the n- and c - terminal forms of 
 protonated aspartate ( ash ) and glutamate ( glh ) , the c - terminal form 
 of neutral lysine ( lyn ) , the terminal and nonterminal forms of deprotonated 
 cysteine ( cym ) , and the noncanonical amino acid norleucine ( nle ) . each solute of interest was solvated in a cubic box of spc / eb water with a clearance of 10  between the solute and 
 the edge of the box , and subjected to a high - temperature md simulation 
 at 450 k , from which were collected a set of 20 conformations . each 
 conformation was subsequently re - equilibrated at 298 k before being 
 input to the ipolq module of mdgx . this module was used to run an 
 md simulation with the solute fixed , during which the coordinates 
 of surrounding solvent molecules were collected and used to generate 
 a collection of point charges representing the solvent reaction field 
 potential . this collection consists of an inner cloud of point charges 
 taken directly from the coordinates of solvent molecules within 5 
  of the solute , and three outer shells of point charges fit 
 to reproduce contributions to the solvent reaction field potential 
 from the infinite periodic system beyond 5  . a pair of 
 qm calculations for the solute were then run at the mp2/cc - pvtz 
 level of theory : one in vacuum and the other including the solvent reaction field 
 potential . as modeled by the collection of point charges . these calculations 
 were run using the orca 3.0.3 software package for each conformation 
 of each residue , requiring over 3000 
 density calculations . the resulting densities were then input to mdgx s 
 fitq module , yielding a pair of charge sets : one valid for simulation 
 under vacuum ( qvac ) and the other for 
 simulation in solution ( qsolv ) . the bonded parameters of ff15ipq 
 were fit to reproduce the relative vacuum - phase qm mp2/cc - pvtz potential 
 energies of a set of diverse conformations of short peptides using 
 an iterative cycle of refinement , similar to that used for its predecessor , 
 ff14ipq . this cycle involved the following 
 steps : ( i ) md simulations were carried out to generate a set of peptide 
 conformations , ( ii ) these conformations were subjected to energy minimization 
 in vacuum using the molecular mechanics ( mm ) energy function with qvac and the current generation of bonded parameters , 
 ( iii ) qm energies of the energy - minimized conformations were calculated , 
 ( iv ) the conformations and energies were used to fit an improved set 
 of bonded parameters , and ( v ) steps ( i ) through ( iv ) were repeated 
 to fit the next generation of bonded parameters . in this way , subsequent 
 generations of the force field  learned  from the biases 
 of their ancestors , provided those biases were captured in the qm 
 energies of the additional conformations that resulted from step ( i ) 
 of the iterative cycle . during the development of ff14ipq , selected 
 conformations from an initial fitting set of 28  000 
 were subjected to energy minimization with each new generation of 
 bonded parameters to yield new conformations , accumulating a total 
 of 65  000 structures and single - point energies . the first generation of ff15ipq fitting data 
 was created by pairing ff14ipq with generalized born implicit solvent 
 md simulations of amino - acid dipeptides 
 at 450 k , followed by vacuum energy minimization of many snapshots 
 from each simulation . while we have not tested how well ff14ipq behaves 
 with implicit solvent , the purpose was to capture any spurious conformational 
 preferences that might remain in the original force field . in addition , 
 we included 1400 conformations of the ace - ala - pro - ala - nme 
 tetrapeptide , while the second generation added numerous tripeptides 
 containing gly , and conformations of the disulfide - bridged cyscys 
 system ( among the largest of all the systems used in qm single - point - energy 
 calculations ) . the next three generations of refinement 
 were designed to cover 
 sampling of the multiple classes of backbone parameters applied to 
 different residues , as described in section 2.4 . in order to ensure sampling of diverse 
 backbone conformations , conformations were generated by progressively 
 restraining  and  at 20 intervals , using a 16 
 kcal / molrad harmonic restraint over the course of 
 the md simulation , yielding 324 conformations of each . since the unique 
 backbone nitrogen type of pro creates unique  and  
 terms for preceding residues , the third generation of conformations 
 consisted of 51 pro - containing tripeptides for which the non - pro residue s 
  and  were restrained . at this point in development , 
 it was decided to branch the positively- and negatively charged residues 
 into unique backbone classes , and as such the fourth and fifth generations 
 of refinement consisted of 57 tripeptides containing the charged residues 
 asp , glu , cym ( deprotonated cys ) , arg , lys , and hip ( doubly protonated 
 his ) , in which  and  of the charged residues were restrained . in order to cover the unique backbone parameters applied to the terminal 
 forms of each residue , additional conformations were added for a set 
 of 78 terminal nxaa - nme and ace - cxaa monopeptides . for these terminal 
 systems , scans of either the unique  of the n - terminal forms 
 or the unique  of the c - terminal forms since the unique backbone 
 nitrogen types of the n - termini and pro in tandem yield an additional 
 set of  terms for nxaa - pro , scans of  were run for an 
 additional set of nxaa - pro - nme dipeptides . finally , in order to cover 
 the unique backbone  and  terms of the amide 
 blocking group ( nhe ) , scans of  and  were run for 17 
 ace - xaa - nhe dipeptides , yielding 324 conformations of each . during 
 these three generations of refinement , 60  000 conformations 
 were added to the ff15ipq fitting set . after the fifth generation 
 of refinement , support for the fitting 
 of angle equilibria and force constants alongside torsions was implemented 
 in mdgx , and subsequent generations emphasized comprehensive sampling 
 of backbone angles . the sixth , seventh , and eighth generations of 
 refinement consisted of perturbations of the angles around n , c , 
 and c. starting from an initial conformation , a selected angle of 
 interest was subjected to a random perturbation within a range of 
 20 of its original equilibrium value ( as inherited from 
 the ff94 force field and retained in contemporaries such as ff14sb ) . target values for the other angles around the same central atom were 
 then chosen by taking their initial values and adjusting them such 
 that the total sum of angles around the central atom was appropriate 
 for the known geometry ; target sums of 360 for planar geometry 
 around n and c and 660 for tetrahedral geometry around c 
 were used . during subsequent mm minimization , the target values for 
 these angles were restrained using 256 kcal / mol rad harmonic 
 restraints . during the eighth and final generation of refinement , angle perturbations were resampled in the context of new scans of 
  and  backbone torsions at 10 intervals for each 
 ace - xaa - nme dipeptide , yielding 1296 conformations of each , alongside 
 additional sampling of terminal monopeptides . during these three generations , 
 125  000 conformations were added , yielding a final 
 fitting set of ff15ipq consisting of > 250  000 single - point 
 qm energies , which is over four times larger than that used for ff14ipq . as done 
 previously for ff14ipq , the torsion parameters of ff15ipq 
 were fit using a linear least - squares fit implemented in the param 
 module of mdgx ; extensions to the module 
 for angle fitting are described in section 2.2 . this module selects a set of torsional 
 barrier heights , angle equilibria , and angle stiffnesses that best 
 reproduce the relative conformational energies of the systems included 
 in the fitting set . during the fitting process , the fourier series 
 lengths and phase angles of the torsional terms were not optimized , 
 and phase angles were set to either 0 or 180 to enable 
 the development of parameters that are transferable to alternative 
 chiralities . all backbone  ,  ,  ,  , 
 and side - chain x torsions of nonterminal forms of the amino acid residues 
 were allocated four terms in their fourier series . torsions unique 
 to the terminal forms of residues and residues preceding pro were 
 restricted to only three terms since these terms were less exhaustively 
 sampled in the fitting set . this restriction was applied to limit 
 the risk of overfitting . while all torsion parameters of residues 
 in the fitting set were fit , only angles in which the central atom 
 was the n , c n type of pro introduces a large number of parameters that are dependent 
 on the preceding residue . in order to avoid overfitting torsional 
 barrier heights , torsions were restrained toward 0 with a force 
 constant of 2  10 kcal / mol . similarly , angles 
 were restrained to their original values , inherited from ff94 , with 
 the equilibria and stiffness force constants set to 5  10 kcal / mol and 2  10 kcal / mol , 
 respectively . to characterize the backbone conformational preferences of ff15ipq 
 in explicit spc / eb water , we carried out umbrella sampling 
 simulations of blocked tetrapeptides ace - ala - xaa - ala - nme , calculating 
 the potential of mean force as a function of the backbone  
 and  torsions of the central residue xaa . in order to identify 
 differences in the conformational preferences between ff15ipq , its 
 predecessor ff14ipq , and contemporary force fields , simulations of 
 ace - ala - ala - ala - nme were carried out using the amber force fields 
 ff15ipq , ff14ipq , and ff14sb ; the opls force field opls - aa / m ; and the charmm force fields charmm36 and drude-2013 . analogous simulations were carried out to compare the conformational 
 preferences of other central amino acid residues using the amber ff15ipq , 
 ff14ipq , ff14sb , and charmm36 force fields . the backbone  and 
  torsions of the central residue were restrained in a series 
 of 1296 windows spaced at 10 intervals , using a harmonic penalty 
 function with a force constant of 8 kcal / mol rad . each 
 window was seeded from a continuous , incrementally restrained simulation , 
 and sampled for 2.0 ns , following a 0.2 ns equilibration . from each 
 set of 1296 windows were reconstructed the unbiased potentials of 
 mean force using the weighted histogram analysis method ( wham ) . to validate 
 ff15ipq as a general force field for peptides and proteins , 
 extensive md simulations on the s time scale were carried out 
 for a variety of benchmark systems consisting of both structured and 
 disordered peptides and proteins . for each system , the amino - acid 
 sequence or pdb code , sources of initial coordinates , and temperatures 
 maintained throughout the simulations are listed in table 1 . hp35 double - norleucine 
 mutant mutant 
 ( lys24nle , asn27his , and lys29nle ) . the p53 peptide used contained residues 
 1729 of the full - length protein and included an n - terminal 
 acetyl ( ace ) and c - terminal amide ( nhe ) blocking group . mdm2 included 
 residues 25109 of the full - length protein , omitting a mobile 
 n - terminal region unresolved in the crystal structure . the n- and 
 c - termini of mdm2 were blocked with acetyl ( ace ) and n - methylamide ( nme ) blocking groups , respectively . in order to accurately match the 
 amino acid sequences used in nmr experiments , residues not resolved in the crystal structure were built using 
 avogadro . residues stsaa were appended 
 to the c - terminus of the s - peptide , and sss to the n - terminus of the 
 s - protein . in addition , ga residues were appended to the n - terminus 
 of the s - peptide , representing a cloning artifact present in the nmr 
 experiments . these residues were not restrained during equilibration 
 of the system . given that the nmr experiments on the s - peptide / s - protein 
 complex were conducted at ph 3.7 , which is close to the average pka values of asp ( 3.7 ) and glu ( 4.1 ) , 
 we elected to run our simulations with negatively charged asp and 
 neutral glu . as done for ff14ipq , we validated ff15ipq 
 by simulating penta - alanine ( ala5 ) , the -helical 
 k19 peptide , the gb1 -hairpin from the c - terminal fragment 
 of protein g , the designed -hairpin chignolin , trp - cage , gb3 , 
 and lysozyme . we also carried out simulations of the -helical 
 ( aaqaa)3 peptide , the cln025 mutant of the chignolin -hairpin , 
 the double - norleucine variant of the villin headpiece subdomain , bovine 
 pancreatic trypsin inhibitor ( bpti ) , ubiquitin , and binase . in order to 
 evaluate the ability of ff15ipq to model disordered proteins , we simulated 
 two classic systems for studying the binding processes of disordered 
 peptides that fold only upon binding their partner proteins : ( a ) the 
 n - terminal p53 peptide and mdm2 oncoprotein , and ( b ) the s - peptide 
 and s - protein cleavage products of the rnase a protein . both of these 
 peptides fold into -helical conformations only upon binding 
 their partner proteins . simulations were performed with these peptides 
 both in isolation and in complex with their protein binding partners . simulations of benchmark systems were carried out using the gpu 
 implementation of the pmemd module in the amber 15 software package . each system was solvated in a truncated octahedral box of spc / eb explicit water with a 12  buffer for the disordered 
 peptide / protein systems and 10  buffer for all other systems . prior to production simulation , each system was subjected to energy 
 minimization , followed by a three - stage equilibration . in the first 
 stage , a 20 ps simulation of the energy - minimized system was carried 
 out at constant temperature while restraining the solute heavy atoms 
 to their initial positions using a harmonic potential with a force 
 constant of 1 kcal / mol  . in the second stage , a 
 1 ns simulation was carried out at constant pressure with the same 
 harmonic position restraints . finally , an additional 1 ns unrestrained 
 simulation was carried out at constant temperature and pressure . temperatures 
 were maintained at selected values ( between 270 k and 370 k ) , using 
 a langevin thermostat ( frictional constant of 1 ps ) , while pressure was maintained at 1 atm using a monte carlo barostat 
 ( 200 fs between attempts to change the system volume ) . van der waals and short - range electrostatic 
 interactions were truncated at 10  ; long - range electrostatic 
 interactions were calculated using the particle mesh ewald method . to enable at least a 2 fs time step , bonds to 
 hydrogen were constrained to their equilibrium values using the shake 
 and settle algorithms . for the k19 , ( aaqaa)3 , gb1 hairpin , chignolin , and cln025 systems , hydrogen mass repartitioning was used to enable the use of longer timesteps . in particular , the masses of solute hydrogen atoms were 
 increased by a factor of 3 , and that of their attached heavy atoms 
 decreased by a corresponding amount such that the total mass remained 
 constant ; the masses of water molecules were not repartitioned . this 
 mass repartitioning scheme enables a 4 fs time step for simulations 
 at 300 k ; for simulations at > 300 k , shorter timesteps 
 were 
 used and set to be equal to 1200 k fs , divided by the set temperature . diagnostics included dssp , rotational diffusion , and nmr relaxation calculated by the ired method . to evaluate the accuracy of ff15ipq 
 in modeling the strengths of 
 protein salt bridges , we simulated the association of three pairs 
 of oppositely charged amino acid side - chain analogues : guanidinium 
 cation / acetate anion ( arg / asp ) , butylammonium cation / acetate anion 
 ( lys / asp ) , and imidazolium cation / acetate anion ( his(+)/asp ) . for 
 each salt bridge , the resulting probability of an anion binding to 
 one or more cation molecules ( pbound ) 
 was compared to experiment ( if available ) , as well as those of six 
 other fixed - charge force fields , including ff14ipq , ff14sb , ff03 , 
 charmm22 * , charmm36 , and opls_2005 , and two polarizable force fields 
 ( charmm drude-2013 and amoeba ) . the original motivation for 
 the development of ff15ipq was to correct for the overstabilization 
 of protein salt bridges by its predecessor , ff14ipq . during the development 
 of ff14ipq , the lennard - jones radii of several polar heavy atoms were 
 refit to reproduce the experimental solvation free energies of side - chain 
 analogues . although the resulting set 
 of radii were initially intended to be applied globally , several of 
 the larger radii resulted in increased 14 repulsion during 
 torsion fitting , which made the torsion parameters more difficult 
 to fit . to overcome this difficulty , mixed lennard - jones combining 
 rules ( called ljedit within amber software or nbfix within charmm ) 
 were applied to these polar groups , assigning different radii for 
 their solute  solvent interactions from those used for their 
 solute  solute interactions . for example , for the carboxylate 
 oxygen atoms of the side - chains of asp and glu larger lennard - jones 
 radii for interactions with water were used than for interactions 
 with solute atoms . an undesirable effect 
 of this strategy , however , was the overstabilization of salt bridges  to 
 the point that in our simulations each acetate molecule was bound 
 to three or more cation molecules ( e.g. , guanidinium ) for most of 
 the simulation ( see figure 1 , as well as figure s2 in the supporting 
 information ) . essentially , the larger radii used for solute solvent 
 interactions forced the carboxylate group out of solution and into 
 interactions with available solute atoms . probability of binding 
 ( pbound ) between 
 acetate and one or more molecules of three cationic side - chain analogues 
 using seven fixed - charge and two polarizable biomolecular force fields , 
 each paired with either the water model with which it was derived 
 or that with which it is most - commonly used . the pbound values corresponding to the experimentally determined ka values of guanidinium acetate and butylammonium 
 acetate are depicted as horizontal gray bars ; no experimental value is available for the imidazolium acetate system . results for the charmm36 , 
 charmm22 * , opls_2005 , amber ff14sb , and amber ff03 force fields are 
 taken from previous simulation studies . for ff15ipq , we addressed the 
 problem of overstabilized salt bridges 
 by discarding the mixed lennard - jones radii of ff14ipq and instead 
 applied empirical corrections to the radii of polar hydrogen atoms 
 bonded to nitrogen ( atom type h  ) in both the protein 
 backbone and side - chains ( note that the original  value of 
 1.07  for this atom type may equivalently be expressed as an 
 r * of 0.6000 , and the details of its fitting appear to have been lost 
 to history ) . these corrections were determined 
 from simulations of the three oppositely charged side - chain analogue 
 systems with h  ranging from the ff94 value of 1.07  
 up to 1.5  , calculating the probability of salt bridge formation 
 ( pbound ) , and comparing this probability 
 to that from experiments . based on the results ( figure s2 ) , we selected a  value of 1.3  for 
 nitrogen - attached hydrogens in both the protein backbone and side - chains . 
 for guanidinium acetate , we found that a further increase in  
 to 1.5  for the side chain of arg was necessary to achieve satisfactory 
 agreement with the experimental value of this system . as shown in figure 1 , ff15ipq yields pbound values that are 
 among the most reasonable , relative to the other fixed - charge force 
 fields that were tested ( ff03 , ff14sb , ff14ipq , opls_2005 , charmm22 * , 
 and charmm36 ) . for guanidinium acetate , our results with ff15ipq are 
 roughly consistent with charmm22 * and ff03 , which we had previously 
 found to provide the most accurate modeling of this system ( note that 
 the corresponding figure of our previous work omitted complexes including 
 multiple cation molecules bound to a single anion molecule and therefore 
 under - represented the extent to which the force fields overstabilize 
 salt bridges ) . for butylammonium acetate , ff15ipq yields a pbound value that is 
 slightly higher than that of amber ff03 , but similar to those of charmm22 * 
 and ff14sb . for imidazolium acetate , ff15ipq yields a pbound that is higher than those of charmm22 * and ff03 , 
 but similar to that of ff14sb . of particular note is that charmm22 * 
 was also parametrized to reproduce the experimental association of 
 guanidinium acetate , but via adjustments to the atomic charges of 
 only the side - chains of arg , asp , and glu ; as a result , these adjusted parameters are inconsistent with the 
 rest of the force field , whose charges had been fit years earlier , 
 using a different method . along the same 
 lines , a recently developed variant of the amber ff99sb - ildn force 
 field has involved the application of mixed lennard - jones combining 
 rules exclusively to interactions between the side chains of arg , 
 asp , and glu . in contrast to the posthoc adjustments 
 of these two other force fields , our approach involves first adjusting 
 the lennard - jones radii , followed by refitting of atomic charges and 
 bonded parameters . this approach  which has been an onerous 
 one in the past  has been significantly streamlined by the mdgx 
 software . we note that all of the fixed - charge force fields 
 are outperformed 
 by the more expensive , polarizable charmm drude-2013 and amoeba force 
 fields . while it is likely that much of their superior performance 
 results from the more complex charge model for the solutes , it is 
 also possible that the solute  solvent interactions  which 
 compete with solute  solute interactions are more accurately 
 represented by the use of polarizable water models . in particular , 
 fixed - charge water models similar to the spc / eb model used 
 here have recently been found to generally underestimate the strength 
 of solute  solvent interactions , and it is possible that this limitation of the water models restricts 
 the accuracy with which the solute models may represent salt bridges . a key metric for assessing the accuracy of the torsion and angle 
 parameters of ff15ipq was the ability to reproduce the target qm potential 
 energy surface . figure 2 shows the distribution and root - mean - square error ( rmse ) of ff15ipq 
 energies , with respect to their target qm potential energies for the 
 20 canonical amino acids . the rmse values for all neutral residues 
 are < 1.3 kcal / mol , while those of the negatively charged residues 
 asp and glu are < 1.9 kcal / mol , and those of the positively charged 
 residues arg and lys are < 2.4 kcal / mol . as shown in the supporting 
 information ( figure s3 ) , the neutral forms 
 asp , gln , and lys ( ash , glh , and lyn , respectively ) have rmse values 
 that are consistent with the other neutral residues , suggesting that 
 the increased rmse values of lys and arg , relative to uncharged residues , 
 are related to their net charge , rather than to their additional flexible 
 x torsions . optimization of the backbone angle parameters introduced 
 a 5%15% improvement in rmse and enabled expansion of the fitting 
 set by more than 4-fold without sacrificing the ability to reproduce 
 those parts of the qm potential energy surface represented in the 
 original dataset ( see the section entitled  timeline of development  
 in the supporting information ) . distributions 
 of residuals of relative molecular mechanical energies , 
 with respect to their qm target potential energies , for 18 ace - xaa - nme 
 dipeptides and the ace - ala - ala - ala - nme and ace - gly - gly - gly - nme tetrapeptides . the 25th , 50th , and 75th percentiles are represented by horizontal 
 lines , and root - mean - square values are represented by white circles . each dataset is colored based on its corresponding backbone torsion 
 class . to assess the accuracy of ff15ipq 
 in modeling the backbone conformational preferences of proteins within 
 computationally tractable systems , we carried out a series of simulations 
 of short peptides that may be affordably simulated to convergence . initially , we focused on simulations of the ace - ala - ala - ala - nme tetrapeptide , 
 for which we calculated the pmf as a function of the backbone  
 and  torsions of the central residue . figure 3 shows the results for ff15ipq , its predecessor 
 ff14ipq , and several contemporary force fields . relative to ff14ipq , 
 ff15ipq has larger free energy barriers ( by 1 kcal mol ) between the  well (   70 , 
   20 ) and  well (  
  80 ,   60 ) and between 
 the  well (   150 ,   
 150 ) and ppii well (   70 ,  
  140 ) . in addition , ff15ipq has a more clearly defined 
  well (   140 ,   
 50 ) . on the left half of the ramachandran plot , 40 ) 
 has decreased slightly , and that of the  well (   
 70 , 40 ) has decreased by 1 
 kcal / mol , while the pii well (   60 , relative to the ff14sb 
 and charmm36 force fields , ff15ipq shows similar  and ppii 
 well depths , although ff14sb and charmm36 do not exhibit  
 or  wells and the precise positions of the various wells differ 
 between the force fields . larger differences are observed relative 
 to the opls - aa / m and polarizable charmm drude-2013 force fields , which 
 have shallower and deeper  wells , respectively . potentials of mean force 
 for the central residue of blocked alanine 
 tetrapeptides , as a function of backbone  and  torsions 
 of the central residue using five fixed - charge force fields and one 
 polarizable force field . each force field was paired with either the 
 water model with which it was derived or that with which it is most 
 commonly used . next , we extended our 
 validation of ff15ipq by examining residue - specific 
 backbone conformational preferences . in particular , we carried out 
 simulations of ace - ala - xaa - ala - nme tetrapeptides containing each of 
 the 20 canonical residues at the central position , including the 25 
 protonation states of these residues that are supported by the force 
 field . for comparison , analogous simulations were carried out using 
 the ff14ipq , ff14sb , and charmm36 force fields . the resulting / 
 backbone torsional preferences of the central residues were then compared 
 to those of ala - xaa - ala obtained from the neighbor - dependent ramachandran 
 distribution ( ndrd ) dataset , derived from conformations observed in 
 the loop regions of proteins ( non--helix/-sheet secondary 
 structures ) . the ndrd dataset is drawn 
 from a collection of 3000 high - resolution crystal structures 
 in the protein data bank , and accounts 
 for the influence of preceding and following residues on the / 
 backbone torsional preferences of the central residue . given the considerable 
 differences in the contexts of our simulations and the ndrd experimental 
 data set , i.e. , solution vs crystal environment , we focused solely 
 on qualitative differences between the simulated and experimental 
 conformational preferences of each peptide . in particular , we compared 
 the conformational preferences of peptides containing a nonalanine 
 central residue , relative to that of the reference ace - ala - ala - ala - nme 
 peptide . generally , both ff15ipq and ff14ipq show greater variation 
 between 
 amino acids than ff14sb and charmm36 , which apply the same backbone 
 torsions to all residues ( figure s4 in 
 the supporting information ) . several differences between ff15ipq and 
 ff14ipq are apparent . for the neutral residues whose c atoms are bound to two heavy atoms , the clearest difference is the decreased 
 favorability of the 180 <  < 90 
 region for asn , gln , leu , and met ; ff15ipq is more consistent with 
 ndrd distributions in which such conformations are rare , because of 
 the broader sampling of such uncommon backbone conformations in the 
 ff15ipq fitting set . an exception is ser , which retains this region 
 and exhibits overall broader sampling , in contrast to the ndrd distribution , 
 in which conformations are restricted largely to the canonical wells . 
 for the neutral residues whose c atoms , the ndrd dataset shows increased conformational 
 preferences in the  region and in the region adjacent to the 
  well , centered at   120 ,  
  60. these preferences are captured by both 
 ff14ipq and ff15ipq , but the lower region is erroneously disfavored 
 by both ff14sb and charmm36 . while ff15ipq is improved relative to 
 ff14ipq by disfavoring the l well of thr , the near - absence 
 of sampling of this well in the ndrd distribution suggests that it 
 may still be too favorable , relative to the experiment . differences 
 between ff15ipq and ff14ipq for the bulky aromatic residues are less 
 pronounced ; the conformational preferences of these residues may be 
 more dependent on sterics as modeled by the lennard - jones parameters , 
 which have not been changed in ff15ipq from those of ff14ipq . the 
 greatest differences between ff15ipq and ff14ipq are observed for 
 the negatively charged residues asp and glu , which have been granted 
 their own  ,  ,  , and  torsions 
 in ff15ipq . these residues largely restrict sampling to the of ppii , 
  , and  wells , lacking clearly defined  
 wells and any wells on the right side of the ramachandran plot . the 
 differences for the positively charged residues arg and lys are much 
 smaller , likely because these residues were already assigned unique 
  and  torsions in ff14ipq . to 
 complement the above qualitative comparisons , we obtained quantitative 
 measures of the accuracy of ff15ipq s backbone conformational 
 preferences by calculating j - coupling constants for the ala5 peptide and comparing these values to experiment . this peptide was 
 the focus of a study by best et al . in 
 which multiple force fields were compared , in terms of their ability 
 to reproduce experimental j - coupling constants using the karplus equation 
 and three different sets of karplus coefficients : the original coefficients , 
 as used by graf et al . , and two sets of dft - based coefficients ( dft-1 and 
 dft-2 ) by case et al . in this study , 
 a suggested criterion for a high - quality force field is that the  value between calculated and experimental j - coupling constants 
 should be 2.25 for all three sets of karplus coefficients . three useful points of reference are the recently developed ff14sb , 
 charmm36 , and ff03w force fields , which were empirically corrected 
 to improve reproduction of experimental ala5 j - coupling 
 constants . the ff14sb 
 force field yielded  values of 0.9 and 1.2 with 
 the original and dft-2 coefficients , respectively , but a higher  of 2.7 with the dft-1 coefficients ; charmm36 and ff03w were 
 tested only with the dft-2 coefficients , yielding  values of 1.16 and 0.9 , respectively . as shown in table 2 , ff15ipq yields values 
 of  = 0.53 with the original 
 coefficients ,  = 1.08 with the dft-2 coefficients , 
 and  = 0.67 with an additional set of karplus coefficients 
 from lindorff - larsen et al . however , 
 similar to ff14sb , we obtained a higher  value 
 of 2.91 with the dft-1 coefficients ; in our case , the higher  value is driven primarily by a single outlier that deviates 
 greatly from the experiment , jhnc. based 
 on results from preliminary versions of ff15ipq ( see the  timeline 
 of development  section in the supporting information ) , it appears that lower  values 
 with the original karplus coefficients may come at the expense of 
 higher  values with the dft-1 coefficients . notably , 
 ff15ipq  which employs a general parametrization to reproduce 
 qm potential energies  performs at least as well as ff14sb , 
 charmm36 , and ff03w , which have been parametrized specifically to 
 reproduce ala5 j - coupling constants . all 
 four force fields yield results that are much improved , relative to 
 force fields developed only a few years ago . uncertainties on  values represent one standard error of the mean , calculated using 
 a block averaging method . uncertainties 
 on individual j - coupling constants are omitted for the sake of clarity . note that the calculated j - couplings 
 are dependent on the ratios 
 of various backbone conformations , between which transitions may be 
 relatively rare . in our studies of ala5 , we found that 
 1.5 s of aggregate simulation time , which we had previously 
 used in our validation of ff14ipq , did 
 not yield sufficiently precise calculations of the j - couplings . although 
 the j - couplings may appear to be converged based on their relatively 
 small statistical variances ( evaluated using block averaging ) , these 
 variances may be misleading . for example , we observed what appeared 
 to be small , but statistically significant differences in the j - couplings 
 between simulations run with and without hydrogen mass repartitioning 
 after 1.5 s of simulation , but these differences ultimately 
 disappeared after 6 s ( see the  timeline of development  
 section in the supporting information ) . the extensive sampling needed to obtain converged j - couplings illustrates 
 the challenge of mapping the conformations of just a few residues 
 using brute - force md simulation . to assess the propensity 
 of ff15ipq to form -helices , 
 we studied the temperature - dependent behavior of two model -helical 
 peptides : k19 and ( aaqaa)3 . both peptides 
 are variants of the motif ( ala - ala - xaa - ala - ala)n , in which xaa is lys in k19 and gln in ( aaqaa)3 ; their sequences are listed in table 1 . for each peptide , we carried out six 4-s simulations 
 at different temperatures and monitored the formation of various types 
 of secondary structure . as shown in figure 4 , both peptides undergo multiple folding 
 and unfolding events , although our simulations are not sufficiently 
 long to obtain converged estimates of the probability of adopting 
 -helical conformations . qualitatively , k19 adopts -helical 
 conformations for a greater proportion of the simulation than ( aaqaa)3 , which is consistent with the experimental observation that 
 k19 and ( aaqaa)3 are 40% and 20% -helical , 
 respectively , at 300 k. both peptides transiently 
 form -sheet contacts , which do not appear to be stable for 
 more than 100 ns , indicating that ff15ipq correctly identifies the 
 favored secondary structures of these peptides . secondary structure of 
 model -helical peptides ( a ) k19 and 
 ( b ) ( aaqaa)3 at various simulated temperatures over the 
 course of 4-s simulations . while the two peptides differ in length , the observed difference 
 in -helical stability is likely due to parameters of gln and 
 lys residues at the central positions of the peptides . in our umbrella 
 sampling simulations of tetrapeptides , which are too small to form 
 an -helix ( figure s4 ) , we observed 
 a broader , deeper  well for lys than for gln , suggesting that 
 the observed difference in -helical stability between the two 
 peptides has already been  built - in  to ff15ipq at the 
 residue level . in addition , the two residues have different backbone 
 charges : gln shares its n , h , c , and o charges with the other neutral 
 residues , and lys shares its charges with the positively charged residues . while the backbone h charges for neutral and positively charged residues 
 are similar , the backbone o of the positively charged residues is 
 0.05 e more negative than that of the neutral residues , which 
 may result in more stable hydrogen bonding . in order to assess the stability of -sheet 
 structures in ff15ipq , we simulated three model -hairpin systems : 
 the gb1 hairpin , the designed peptide chignolin , and its hyper - stable 
 variant cln025 . we simulated the gb1 hairpin at 6 temperatures , ranging 
 from 275 to 325 k , and we simulated cln025 at 10 temperatures , ranging 
 from 280 k to 370 k ; we simulated chignolin only at 298 k. figure 5 shows the secondary 
 structures observed during 4-s simulations of these systems . as with our simulations of the -helical peptides , our -hairpin 
 simulations are not sufficiently long to precisely quantify secondary 
 structure stability , although qualitative trends may be identified . 
 as shown in figure 5a , the gb1 hairpin is metastable over the tested temperature range 
 of 275325 k , and in two of our simulations unfolds and refolds . our simulations at 285 k are in qualitative agreement with 
 the experiment , which have indicated that the gb1 hairpin is 85% 
 folded at 275 k ,  50% folded at 295 k , and 20% folded 
 at 325 k. however , an anomaly is observed 
 in our 275 k simulation , in which the gb1 hairpin unfolds after 200 
 ns and does not refold . this unfolding event may simply be an artifact 
 of our limited sampling that would disappear if the simulations were 
 to run to convergence . alternatively , it may reflect limitations of 
 the spc / eb water model at temperatures distant from those 
 at which it was parametrized ; while the temperature - dependent behavior 
 of spc / eb has not been characterized , to our knowledge , 
 three - point water models including the parent spc / e water model are 
 known to poorly reproduce the temperature dependence of properties 
 such as density . secondary 
 structures of model -hairpin peptides ( ( a ) gb1 
 hairpin , ( b ) chignolin , and ( c ) cln025 ) at various simulated temperatures 
 over the course of 4-s simulations . in contrast , our simulations of chignolin and cln025 suggest 
 that 
 these -hairpin systems may be more stable than observed experimentally . 
 as shown in figure 5b , chignolin maintains its -hairpin configuration throughout 
 our 4-s simulation at 298 k , including two hydrogen bonds in 
 an antiparallel sheet configuration , while , experimentally , the peptide 
 is only 60% folded at this temperature . chignolin s hyper - stable variant cln025 has an experimental 
 melting temperature of 343 k. as shown 
 in figure 5c , in our 
 simulations at temperatures ranging from 280 k to 370 k , we observe 
 unfolding and refolding events at several temperatures , although the 
 overall folded population is larger than that measured experimentally . in particular , cln025 is > 80% folded in our simulation at 370 k , 
 while , 
 experimentally , the peptide is only 25% folded at this temperature . as with -helices , we expect ff15ipq 
 to yield residue - specific 
 propensities in -sheet stability , although the large difference 
 in sequence between the two tested types of model -hairpins 
 makes comparing them difficult . the aforementioned lack of a clear 
  well in asp may destabilize the gb1 hairpin , which contains 
 two adjacent asp residues , one of which forms part of the antiparallel 
 -sheet and the other , the turn . the observed stabilities of 
 chignolin and cln025 generally preclude the notion that ff15ipq is 
 biased against the -sheet structure . further studies including 
 additional hairpin sequences and parallel -sheet structures 
 will be necessary to quantify and mitigate residue - specific biases 
 for future ipolq force fields . in order 
 to assess the stability of proteins with ff15ipq , we simulated the 
 trp - cage miniprotein and a series of six globular proteins : bpti , 
 villin , gb3 , ubiquitin , binase , and lysozyme . extensive experimental 
 data are available for all of these model systems , providing excellent 
 opportunities for validation of our simulations . we carried out 24 
 s of aggregate equilibrium simulations of trp - cage at temperatures 
 ranging from 275 k to 325 k and simulations 210 s in 
 duration of the six globular proteins at temperatures between 298 
 k and 303 k. details of our simulations are listed in table 1 . the designed miniprotein 
 trp - cage is central to a long - running computational success story 
 for amber force fields . folding simulations of this miniprotein using 
 the amber ff99sb force field have successfully recovered the folded 
 structure , yielded multiple folding and unfolding events , and provided 
 a melting temperature ( tm ) of 318 k , which 
 is in reasonable agreement with the experimental tm value of 315 k. as shown in figure 6 , trp - cage remained stable in our simulations between 
 275 k and 295 k , with an average backbone rmsd from the experimental 
 nmr structure of < 1  , and unfolded between 305 k and 325 
 k. while our simulations are not extensive enough to obtain precise 
 estimates of the tm , these results suggest 
 that the tm value is somewhere between 
 295 k and 305 k , which is slightly lower than the experiment . each 
 unfolding event is marked by an initial shift of the backbone / 
 of pro 12 from the  well to the ppii well , followed by the 
 loss of the n - terminal -helical component of the polypeptide . notably , in our simulation that was run at 325 k , the protein refolded 
 for 500 ns , indicating that the folded state is a stable free - energy 
 minimum . thus , despite the extensive reoptimization of the parameter 
 set , the important success of the amber force fields in modeling the 
 stability of trp - cage has been maintained . stability of trp - cage 
 over the temperature range of 275325 
 k over the course of 4-s simulations , as monitored by the backbone 
 root - mean - square deviation ( rmsd ) , relative to the experimental nmr 
 structure . note that the pro - rich sequence 
 of trp - cage allows the opportunity 
 to validate the unique pre - pro  and  terms of 
 ff15ipq , because it contains gly , arg , and pro residues that precede 
 pro . whenever trp - cage was folded in our simulations , gly 11 remained 
 stably in its ppii well , while arg 16 was sampled broadly 
 across the  and ppii regions without a clear barrier between 
 them . these results are in good agreement with the experimental nmr 
 ensemble , within which the  and  backbone torsions of 
 arg 16 are distributed in a line across these two regions , and with the pre - pro distributions observed 
 for these residues in the ndrd dataset . also consistent with both the nmr ensemble and ndrd dataset are 
 the observed distributions of pro 17 and pro 18 , which strictly maintained 
 their positions in the ppii well , even as the protein unfolded . as shown in figure 7 , the overall structures of all six globular proteins bpti , 
 villin , gb3 , ubiquitin , binase , and lysozyme  remained stable 
 over their entire simulations . the bpti protein remained closest to 
 its crystal structure , yielding an average backbone rmsd of 0.7  , 
 which may be a result of its three disulfide bonds among a total of 
 58 residues . all -helical and -sheet regions of this 
 protein were retained for the entire simulation . the only notable 
 deviation from the crystal structure of bpti was observed for ala 
 16 and arg 17 , which form the end of the loop preceding the first 
 -sheet . these two residues , which occupy the  and  
 wells , respectively , in the crystal structure , made temporary excursions 
 of several hundred nanoseconds to alternative conformations before 
 returning to the crystal conformation ( figures s5 and s6 in the supporting information ) . the villin headpiece 
 subdomain also remained close to its crystal structure , yielding an 
 average backbone rmsd of 1.1  . in order to test our parameters 
 for the noncanonical amino acid norleucine ( nle ) , which was included 
 alongside the canonical residues during development of ff15ipq , we 
 have simulated the fast - folding double - nle mutant of villin . the nle residues , both of which are located 
 in the third -helix , strictly sampled the -helical well , 
 suggesting that our parameters for nle are appropriate . the only significant 
 deviation from the crystal structure of this mutant of villin is a 
 0.5 s excursion made by residues 1012 , which 
 form the end of helix 1 and the loop linking helices 1 and 2 , to an 
 alternative conformation different from that observed in the crystal 
 structure ( figures s7 and s8 in the supporting 
 information ) . stability of folded proteins over the course of 10-s 
 simulations 
 as monitored by the backbone rmsd , relative to the experimental structures . 
 for binase , the mean rmsd relative to the ensemble of 20 nmr structures 
 is shown , along with the range between the minimum and maximum values 
 ( light blue shaded region ) . our simulation of gb3 yielded an average backbone rmsd of 
 1.0  
 from the nmr structure . three of the residues exhibit significant 
 deviations from the crystal structure : leu 12 , asp 40 , and thr 55 
 ( figures s9 and s10 in the supporting information ) . the conformation of leu 12 , which is located in the turn linking the 
 first and second -strand , falls precisely between the  
 and ppii wells in the nmr structure , while both wells were almost 
 equally sampled in our simulation . as a result of this increased conformational 
 flexibility at leu 12 , adjacent residues also occasionally sampled 
 conformations outside their nmr structure . the presence of a free - energy 
 barrier between the  and ppii wells is a necessity for maintaining 
 stable conformations within these wells ; it is likely that the forces 
 contributing to the stabilization of leu 12 s unusual conformation 
 in the nmr structure are simply not captured by the functional form 
 of ff15ipq . consistent with this hypothesis , almost - identical deviations 
 were observed for the ff14sb force field that shares this functional 
 form . the conformation of asp 40 , which 
 is located in the loop between the -helix and third -strand , 
 occupies the  well in the crystal structure while other conformations , 
 predominantly  , were sampled in our simulations . indeed , based 
 on this result and others presented below , the negatively charged 
 residues of ff15ipq completely lack an  well . as with leu 12 , the increased conformational flexibility of asp 40 led to broader 
 sampling by adjacent residues . it is worth noting that gb3 contains 
 two negatively charged residues , glu 15 and asp 22 , which remained 
 stable in -sheets , indicating that the limited sampling of 
 this region observed in our umbrella sampling simulations may be overcome 
 within the context of a folded protein . finally , a notable deviation 
 from the nmr structure occurs for thr 55 during the last microsecond 
 of our simulation when the antiparallel -sheet hydrogen bonds 
 between this residue and val 42 are broken , although the remainder 
 of the -sheet remains in place . while this deviation may be 
 a transient event , it could also be a consequence of the conformational 
 deviations of the nearby asp 40 . similar to gb3 , ubiquitin in 
 our simulations exhibited a low overall 
 average backbone rmsd of 1.2  from the crystal structure , but 
 significant deviations in certain regions . in particular , transient 
 deviations from the crystal structure were observed for residues 811 , 
 which form the turn connecting the first and second -strands 
 ( see figures s12 and s13 in the supporting 
 information ) . while these residues sampled conformations that were 
 different from the crystal structure , their turn conformation was 
 retained throughout 80% of the simulation , and experimental nmr relaxation 
 data suggests that the region is truly flexible . a more significant deviation was observed for asp 52 and 
 gly 53 , which are located in a loop region ; in the crystal structure , 
 these residues both occupy the  well , whereas , in our simulation , 
 asp 52 and gly 53 shift to the ppii and  wells , respectively . 
 unlike the deviations observed in gb3 , this shift does not lead to 
 broader sampling by adjacent residues or appear to otherwise destabilize 
 the protein , suggesting that the observed alternative conformation 
 may simply be erroneously modeled by ff15ipq to be lower in energy 
 than the conformation found in the crystal structure . finally , after 
 8.5 s of simulation , glu 34 , which is the last helical 
 residue in the central -helix , shifts to the ppii and  
 regions , leading to shifts in residues 3336 . combined with 
 the observations made for asp 52 , this shift suggests that , for negatively 
 charged residues , ff15ipq may overstabilize ppii conformations , relative 
 to . among the simulated globular proteins , the greatest 
 deviations 
 from the initial structure were observed for binase with an average 
 backbone rmsd of 3.4  from the nmr ensemble of 20 models . the same average rmsd was also obtained , with 
 respect to the crystal structure of wild - type binase , which differs 
 in the amino acid sequence at six positions ( pdb code : 1gou ) . these larger differences are primarily caused by variability 
 in loop regions , as noted for the experimental structures , and the core structure of binase remained relatively close to the 
 experimental structures with an average backbone rmsd of 1.9  . 
 the first loop , which consists of residues 3439 , adopted multiple 
 conformations in our simulation , which is consistent with the nmr 
 ensemble ( figures s15 and s16 in the supporting 
 information ) . moreover , this loop adopts the conformation observed 
 in the crystal structure for 75% of the simulation . the second loop 
 consists of residues 5662 , which also sampled broadly in our 
 simulations , consistent with diverse conformations in the nmr ensemble 
 and poorly defined electron density in the crystal structure . notably , in both our simulation and the nmr 
 ensemble , flexibility in this region extends to gly 67 . this difference 
 may relate to the difference in sequence between the nmr and crystal 
 structure proteins , in which ser 66 is replaced by gly and gly 67 
 by ser . the third loop , comprised of residues 7683 , is also 
 broadly sampled in our simulation , and is the source of greatest difference , 
 relative to the experimental structures . the final loop is comprised 
 of residues 99104 and , in our simulation , we observe several 
 residues sampling two different conformations , consistent with the 
 observation of two states in crystal structures determined under different 
 conditions . the largest protein 
 system , lysosome , was simulated for 2 s , 
 over which it exhibited an average backbone rmsd of 1.2  , with 
 respect to the crystal structure . as shown in figures s17 and s18 in the supporting information , the largest 
 deviations were found in the loop comprised of residues 100104 . residues 101104 adopted an alternative conformation , with 
 the flanking residue val 99 no longer part of the preceding  
 helix and gly 104 part of the following helix . similar to our observations 
 for ubiquitin , this difference appears to be related to the shift 
 of a negatively charged residue , asp 101 , from the  well to 
 the ppii well . bpti , 
 gb3 , ubiquitin , 
 and lysozyme  have also been used for validation of previously 
 developed force fields , including ff14ipq , ff14sb , and charmm36 . similar to these force fields , 
 ff15ipq yielded low average backbone rmsd values for these proteins , 
 relative to their initial structures ( i.e. , 1.2  ) . a 
 key point to be considered while comparing our results to those of 
 previous force fields is that advancements in gpu computing over the 
 last several years have enabled us to 
 validate ff15ipq using simulations up to 10 s , which is up 
 to 10 times longer than those used for the other force fields . in 
 particular , ff14sb was validated using sets of four 1-s simulations , while charmm36 was validated using 200 ns simulations . many of the key deviations that we observe do 
 not occur for several microseconds ; for example , we observe changes 
 in the c - terminus of gb3 after 9 s , and in the loops of ubiquitin 
 after 8.5 s . these deviations are informative and will guide 
 development of successors to ff15ipq , illustrating the utility of 
 long - time scale simulations for force field development . a particularly 
 appealing feature of the spc / eb water 
 model with which we have developed ff15ipq is its ability to more 
 accurately reproduce the rotational diffusion of solvated proteins 
 relative to water models , such as tip3p and tip4p - ew . in order to measure how accurately the combination of ff15ipq 
 and spc / eb are able to reproduce rotational diffusion , 
 we branched off sets of ten 200-ns simulations in the microcanonical 
 ensemble ( nve ) from our 10-s simulations in the canonical ensemble 
 ( nvt ) for gb3 , ubiquitin , and binase , thereby avoiding perturbation 
 of the dynamics by the use of a thermostat . as shown in table 3 , ubiquitin , gb3 , and binase 
 diffused 14% , 15% , and 22% more slowly than 
 that measured experimentally by nmr , respectively . note that the experimental 
 values were corrected for differences in temperature , isotopic labeling , 
 and solvent d2o content , potentially introducing error 
 into the comparison of simulated and experimental values . interestingly , while the errors that we obtained 
 were consistent with our test simulations with the amber ff99sb - ildn 
 force field and spc / eb ( see figure s1 ) , we found charmm22 * and spc / eb to yield lower 
 errors of 7% , 6% , and 16% , illustrating the coupling of solute and 
 solvent parameters on the motions of proteins through solution . since 
 the spc / eb water model had been empirically optimized for 
 proteins with the amber ff99sb force field , this is no fault of ff15ipq 
 ( or charmm22 * , for that matter ) , but suggests that improved performance 
 might be obtained by optimizing the protein and solvent models in 
 tandem . finally , it is worth noting that our use of the langevin thermostat 
 in the nvt simulations results in 42%52% longer rotational 
 diffusion times , relative to those of the nve simulations , demonstrating , 
 for the first time ( to our knowledge ) , that the use of a thermostat 
 can significantly perturb dynamical properties for proteins and not 
 just for small molecules and polymer chains . uncertainties represent one standard 
 error of the mean . experimental 
 rotational diffusion 
 measured using nmr relaxation and corrected for differences in temperature and d2o 
 content between simulation and experiment . uncertainties represent 
 one standard 
 error of the mean calculated from 50 consecutive 200-ns blocks from 
 a single 10-s simulation . uncertainties represent one standard 
 error of the mean calculated from ten independent 200-ns simulations . a major advantage of performing 
 simulations with accurate rotational 
 diffusion is that one can directly calculate nmr relaxation parameters n r1 and r2 , and the n h heteronuclear 
 noe , which report on the dynamics of individual residues , and compare 
 these values with experiment . therefore , we calculated relaxation 
 parameters for gb3 and ubiquitin , for which experimental data are 
 available at five and four magnetic field strengths , respectively 
 ( see figures s11 and s14 , respectively , 
 in the supporting information ) . overall , our calculated r2 values are in excellent agreement with the experiment , 
 with average mean absolute percent error ( mape ) values of 8% for gb3 
 and 9% for ubiquitin . our r1 values are also in good agreement , 
 with average mape values of 10% and 12% , with a consistent offset 
 observed across all residues . however , our heteronuclear noe values 
 are somewhat poorer agreement , with average mape values of 22% and 
 30% for the two systems . residues for which our calculated r1 and r2 differ 
 by > 20% from the experimental 
 values may be due to limitations of our force field . notably , we find 
 that several such residues are those for which we also observed deviation 
 in sampled backbone conformations , relative to the experimental structures . 
 for gb3 , leu 12 and asp 40 yielded above - average errors in r1 ( 25% ) , while gly 41 had a larger error in 
 both r1 and r2 ( 47% and 33% ) . similarly for ubiquitin , lys 11 and asp 52 yielded 
 above average errors in r1 ( 20% ) . also 
 in ubiquitin , asn 25 yielded an error of 22% in r2 ; however , this residue is found experimentally to undergo 
 chemical exchange , fluctuating at time 
 scales beyond those captured by our simulations . among residues with 
 errors below 20% , one particular trend is apparent : ile 7 , thr 16 , 
 thr 49 , and thr 51 of gb3 and thr 12 and ile 36 of ubiquitin all yielded 
 errors in r2 of 15% , despite having 
 tightly restricted / sampling consistent with their 
 experimental structures . this trend suggests that ff15ipq may have 
 some discrepancy with the three branched residues that is not apparent 
 when examining only backbone / preferences , and this 
 will be the subject of further study . in order to test the suitability 
 of ff15ipq for simulating disordered 
 proteins , we focused on two model peptides : the n - terminal , 13-residue 
 peptide fragment of the tumor suppressor p53 , and the 22-residue s - peptide 
 fragment of rnase a. both of these disordered peptides ( p53 peptide 
 and s - peptide ) only adopt -helical conformations when bound 
 to their structured partner proteins ( mdm2 and s - protein , respectively ) . for each of these peptides , we carried out two 10-s simulations : 
 one of the isolated peptide and the other of the native peptide both peptides in their isolated states adopted conformations 
 distant from their partner - bound conformations , sampling a diverse 
 set of conformations with average backbone rmsds of 5  from 
 their corresponding bound conformations in the crystal structures 
 of the peptide  protein complexes and maintaining only 25% 
 of their native intrapeptide contacts . furthermore , the p53 peptide 
 only transiently adopted -helical conformations that resembled 
 its partner - bound conformations ( backbone rmsd value of <3  ) 
 with these conformations unfolding within 200 ns ( see figure s19 in the supporting information ) . in 
 contrast , the s - peptide did not even transiently sample -helical 
 conformations resembling its bound state ; instead , the peptide formed 
 -hairpins that persisted for periods as long as 4 s 
 ( see figure s22 in the supporting information ) . stability 
 of p53 and s - peptide alone and in complex with binding 
 partners mdm2 and s - protein over the course of 10-s simulations 
 as measured by backbone rmsd relative to ( a , b ) the crystal structures 
 and ( c , d ) the percent of native contacts formed . in our simulation of the p53/mdm2 complex , the peptide remained 
 stably bound to its partner protein for the entire 10 s of 
 simulation , with an average backbone rmsd of 1  from its bound 
 conformation in the crystal structure ( figure 8) . curiously , while most of the intramolecular 
 native contacts of the p53 peptide and mdm2 were retained ( 95% 
 and 85% , respectively ) , only 60% of the intermolecular 
 p53/mdm2 native contacts persisted . examination of the structure showed 
 that many of these contacts lay just below the threshold distance 
 of 5.5  between their heavy atoms in the crystal structure . thus , these contacts were no longer  formed  when slightly 
 different conformations were adopted in our simulations . throughout 
 our simulation , the p53 peptide retained the -helical structure 
 of residues 1925 ( recall figure s19 ) . in contrast , our simulation of the s - peptide / s - protein complex 
 sampled conformations more distant from its crystal structure , with 
 average backbone rmsd values of 3  for both s - peptide and s - protein . 
 although the s - peptide partially lost its -helical structure 
 near the n - terminus from 1.5 s to 5 s 
 in our simulation ( see figure 8 and figure s22 ) , the structure 
 reformed and persisted for the remainder of the simulation during 
 which both the s - peptide and s - protein retained most of their intramolecular 
 native contacts ( 90% and 80% , respectively ) . similar 
 to the p53/mdm2 complex , the s - peptide / s - protein complex retained 
 60% of its peptide  protein contacts throughout our 
 simulation . our simulations of these flexible , disordered peptides 
 provide 
 an additional opportunity to validate the backbone conformational 
 preferences of ff15ipq . from our four simulations , we have calculated 
 the backbone / sampling of the p53 peptide ( figure 9 ) and s - peptide ( figure s23 in the supporting information ) . for 
 comparison , we have provided the distributions for the peptide sequences 
 obtained from the ndrd dataset , which accounts for the influence of 
 adjacent residues on each distribution . in our simulation of the p53/mdm2 complex , residues 1827 
 of p53 occupied exclusively the wells observed in the crystal structure , 
 aligning with the observed low rmsd and high percentage of native 
 contacts . our simulation of the isolated p53 peptide exhibits similarities 
 with the ndrd distribution for most residues , but several inconsistencies 
 are informative : in particular , the neutral residues exhibit reasonable 
 overall agreement and the bulky aromatic residues phe 19 and trp 23 
 closely resemble the ndrd distributions , while leu 22 and leu 25 sampled 
 of the -region more extensively than suggested by the ndrd . the negatively charged residues glu 17 and glu 28 sampled consistently with the ndrd , while 
 asp 21 missed sampling in the  region . the sole positively 
 charged residue , lys 24 , sampled the -well more backbone conformational sampling of the disordered 
 p53 peptide 
 observed in 10-s simulations ( a ) in complex with the mdm2 protein 
 and ( b ) alone . for comparison , panel ( c ) shows distributions for the 
 p53 sequence obtained from the neighbor - dependent ramachandran distribution 
 ( ndrd ) dataset , derived from conformations 
 observed in the loops of crystal structures . in our simulation of the s - peptide / s - protein 
 complex , the crystal 
 conformations were retained for most of the simulation ( see figure s23 ) . as described above , several residues 
 near the n - terminus left the -helical well , and the adjacent 
 thr 3 , which is not helical in the crystal structure , eventually joined 
 the helix as it reforms . the clearest difference from the crystal 
 structure is found for asp 14 ; this residue occupies the  well 
 in the crystal structure , which is not present for asp in ff15ipq , 
 causing it to adopt a ppii conformation . within the s - protein , conformation 
 (   100 ,   130 ) , which was retained throughout our simulation ( figure s25 in the supporting information ) . in our simulation 
 of the isolated , unbound s - peptide , the formation of long - lived -hairpin 
 structures prevented us from obtaining converged conformational preferences 
 for comparison with the ndrd distributions , despite the long duration 
 of the simulations ( 10 s ) . taken together , the above 
 results indicate that ff15ipq can reliably 
 predict disorder as well as order for peptides that fold upon binding 
 their partner proteins . these encouraging results are worth pointing 
 out since ff15ipq was not specifically parametrized for disordered 
 peptides / proteins , as is the case for contemporary force fields such 
 as ff03w and its subsequent variants . as shown in figure 1 , both ff03 ( whose 
 atomic charges and radii are shared by ff03w ) and ff15ipq are able 
 to reliably model propensities of salt - bridge formation , which can 
 be critical for such systems that are rich in polar and/or charged 
 residues . thus , ff15ipq is a reasonable alternative to ff03w for the 
 simulation of disordered peptides / proteins . since the establishment of the  one atom , 
 one site  
 model of all - atom fixed - charge force fields over 20 years ago , several 
 major lineages of protein force fields and countless branches have 
 been developed through cycles of validation and refinement . in this 
 work , we present the ff15ipq force field , the latest in the amber 
 ipolq lineage , and validate its accuracy by > 200 s of aggregate 
 md simulation . the distinguishing features of ff15ipq are ( i ) a charge 
 set that accounts directly for water induced polarization , ( ii ) the 
 incorporation of two related charge sets for creating new force fields 
 based purely on ab initio calculations , ( iii ) the 
 scope of the parameter optimization , including backbone angles alongside 
 torsions , and ( iv ) the degree of automation and transferability of 
 the methods to other regions of chemical space . our simulations suggest 
 that ff15ipq yields reasonable salt - bridge propensities , maintains 
 secondary structures and globular protein folds on the s time 
 scale , predicts order as well as disorder in protein structures , and 
 yields strong agreement with nmr j - couplings and relaxation rates . however , even with this extensive amount of validation , several unconverged 
 results remain and will be explored further using enhanced sampling 
 techniques such as replica exchange or recent 
 variants of the weighted ensemble path 
 sampling strategy . a major motivation for creating ff15ipq was to address concerns 
 about the overstabilization of salt - bridge interactions by ff14ipq , 
 which is a limitation that is shared by other contemporary force fields . we addressed these concerns by abandoning the 
 mixed lennard - jones radii of ff14ipq and instead increasing the radii 
 of polar hydrogen atoms bonded to nitrogens in both the protein backbone 
 and side - chains to yield more accurate salt - bridge propensities . the 
 atomic charges of the force field were then rederived by applying 
 the automated machinery , which had created ff14ipq , exchanging the 
 tip4p - ew water model for spc / eb . finally , all torsion and 
 selected backbone angle parameters were refit to a qm dataset over 
 four times as large as that used for ff14ipq . in doing so , we found 
 that angle optimization was essential for recovering reliable results 
 with our more - extensive dataset ( see the  timeline of development while the angle optimization feature and other particulars of the 
 torsion fitting are subjects of ongoing development in our force field 
 engine mdgx , our results with ff15ipq indicate that the workflow illustrated 
 in figure 10 is a 
 viable approach to creating new force fields . each step entails additional 
 layers of details in order to address practical considerations such 
 as infinite electrostatics or the forms of molecular mechanics basis 
 functions and the shape of the target energy surface . we have addressed 
 each of these issues in the theory and methods sections as well as in prior publications , but the synthesis of all these details reflects the physical arguments 
 behind the ipolq charge derivation . starting from an existing 
 model , selected global changes are optionally first applied to obtain 
 an initial model for optimization . the ipolq charge deviation protocol 
 is then used to fit a pair of atomic charge sets for the vacuum ( qvac ) and solution ( qsolv ) phases . the vacuum - phase charges are used to fit parameters for 
 bonded terms to vacuum - phase qm targets , and these parameters are 
 subsequently paired with the solution - phase charges to yield a complete 
 force field for solution - phase simulations . the force field is then 
 validated through extensive md simulation , informing future development . this approach should be viewed 
 in context with the contemporary 
 force balance approach , which also performs sweeping optimization 
 of hundreds of parameters simultaneously . unlike our mm - minimized conformations , force balance typically considers 
 conformations qm - minimized at the same level of theory at which the 
 target qm energies are calculated , although this is not a strict requirement . going beyond the capabilities of mdgx , force balance includes numerous 
 nonlinear optimization methods and offers the capability to incorporate 
 results from sources beyond qm single - point energies , including in vitro experiments , directly into the parameter optimization . 
 in the future , such diverse targets might be paired with the ipolq 
 method , for example , by using the vacuum charge set ( qvac ) for comparison with vacuum - phase qm data , but using 
 the polarized charge set ( qsolv ) in simulations 
 for comparison to experimental results in solution . it is rather 
 remarkable that a viable protein force field can be 
 produced within months , almost entirely from qm data . also noteworthy is the fact that features such as angle optimization and generational 
 refinement , which had incremental but definite effects on the accuracy 
 of data fitting , could be so influential in the final result . one 
 way of considering the remaining error in our mm model is to partition 
 it between two sources : bonded and nonbonded interactions . the improvements 
 in ff15ipq that were obtained relative to ff14ipq , whose nonbonded 
 parameters are of similar accuracy , resulted from optimization of 
 angles and the branching of bonded parameters . while the inclusion 
 of anharmonicity in bond and angle stretching or a spline - based treatment 
 of torsion cross terms ( cmap ) may reduce errors further , greater improvements might be accomplished in the nonbonded interactions . for this reason , the next planned advance of the amber ipolq force 
 field lineage is to improve the accuracy of electrostatic interactions 
 by making liberal use of virtual charge sites . paradigm used for ff15ipq , which was established 
 several decades ago , appears sufficiently accurate for most purposes , 
 economical by construction , and thoroughly optimized in existing md 
 engines . models with significant numbers of virtual charge sites are 
 presently in the process of becoming established . these models offer 
 improved accuracy for various chemistries with clear physical motivations , 
 accompanied by a more modest increase in computational cost than that 
 afforded by polarizable functional forms . future amber ipolq development at the one - site - per - atom level will 
 continue to explore the applicability of our methodology to the chemical 
 space of other biologically important molecules , including nucleic 
 acids , carbohydrates , and small molecules . further ahead along the 
 path lie multisite ipolq models , which will push the mimicry of qm 
 potential energy surfaces  within the confines of a nonpolarizable 
 model  to new levels .
we present the amber 
 ff15ipq force field for proteins , the second - generation 
 force field developed using the implicitly polarized q ( ipolq ) scheme 
 for deriving implicitly polarized atomic charges in the presence of 
 explicit solvent . 
 the ff15ipq force field is a complete rederivation 
 including more than 300 unique atomic charges , 900 unique torsion 
 terms , 60 new angle parameters , and new atomic radii for polar hydrogens . 
 the atomic charges were derived in the context of the spc / eb water model , which yields more - accurate rotational diffusion of 
 proteins and enables direct calculation of nuclear magnetic resonance 
 ( nmr ) relaxation parameters from molecular dynamics simulations . 
 the 
 atomic radii improve the accuracy of modeling salt bridge interactions 
 relative to contemporary fixed - charge force fields , rectifying a limitation 
 of ff14ipq that resulted from its use of pair - specific lennard - jones 
 radii . 
 in addition , ff15ipq reproduces penta - alanine j - coupling constants 
 exceptionally well , gives reasonable agreement with nmr relaxation 
 rates , and maintains the expected conformational propensities of structured 
 proteins / peptides , as well as disordered peptides  all on the 
 microsecond ( s ) time scale , which is a critical regime for 
 drug design applications . 
 these encouraging results demonstrate the 
 power and robustness of our automated methods for deriving new force 
 fields . 
 all parameters described here and the mdgx program used to 
 fit them are included in the ambertools16 distribution .