molecular simulations of biomolecules are expressions of their 
 underlying force fields , sampling the interactions of chemically bonded 
 atoms through newtonian approximations of the quantum systems . the 
 typical model , given in eq 1 , is well - known . it approximates the energy of a molecular system as a thoroughly decomposable , 
 easily differentiable sum of terms involving harmonic bonds , harmonic 
 angles , lennard  jones repulsion / dispersion terms , electrostatic 
 interactions , and additional parameters guiding the dihedral preferences 
 of the model . the  additional parameters  
 are both the least physically grounded and the most frequently edited 
 part of most force fields ; numerous forms of the expressions modifying 
 the dihedral potential energy include screening factors applied between 
 the electrostatic and lennard  jones interactions of atoms connected 
 by chains of three bonds , fourier series in the dihedral angles made 
 by such atoms , and coupling terms between fourier series of consecutive 
 dihedral atom chains . changes in these parameters that 
 have made considerable improvements in numerous force fields should be viewed in light of the fact that their adjusments to the 
 system energy are small . the terms describe a regime between the high - frequency 
 motions of bonded atoms and the low - frequency rearrangements of nonbonded 
 chemical groups and appear to be a sort of keystone in numerous molecular 
 models.1many molecular force fields can be classified 
 into several lineages . while some force 
 fields have been developed to reproduce bulk liquid properties , new 
 models today tend to emerge as new levels of quantum mechanical theory 
 become accessible and also as updates based on more extensive fitting 
 data or inconsistencies with known biochemical data . within each lineage , the models tend to evolve while retaining a significant portion of 
 their parameters from previous work . even comparatively minor changes 
 can take years to gain acceptance , however , and any jump to parameter 
 development based on a new level of quantum theory also requires significant 
 effort to validate . the lineages of force fields are a natural consequence 
 of the economics of force field development , particularly in the academic 
 community where novelty and publications are essential . this environment 
 has also driven the evolution of composite models where the various 
 energy terms are shaped by different levels of quantum theory . , incorporated first into the amber ff94 force 
 field , has remained in service for nearly 20 years . these charges 
 were developed according to the kollmann resp method , using hf/6 - 31 g * quantum calculations to provide the target 
 electrostatic potentials . numerous other force fields have adopted 
 the charge model , most of them distinguished by new dihedral fourier 
 series , typically based on quantum calculations 
 at the mp2 level with either a cc - pvtz or 6 - 31++g basis set . while the mismatch in quantum targets implies that the torsion terms 
 are effectively modifying the short - ranged behavior of the charge 
 set to fit a different potential energy surface , the force fields 
 derived in this way have led to impressively stable protein behavior . 
 in 2003 , duan et al . derived an alternative charge model , based on 
 b3lyp / cc - pvtz quantum calculations in a polarizable continuum ( pcm ) 
 solvent intended to mimic the interior of a protein . backbone torsion fourier series were derived specifically 
 for this new charge set at the mp2/cc - pvtz level of theory , also in 
 the context of pcm solvent , to complete the amber ff03 force field . while the charges and single - point energies are derived with similar 
 styles of qm theory , the standard torsion fitting procedure does not 
 incorporate its own pcm solvent , and it is not trivial to extract 
 the energy of the pcm - polarized wave functions in vacuum . molecular 
 mechanics gas - phase energies computed with charges derived in the 
 context of pcm solvent have accordingly been shown to double - count 
 polarization effects , as the statically polarized charges are likely 
 to be further affected by direct interaction with the solvent . it is unclear how this imbalance may have affected 
 the quality of simulations performed with ff03 , but the force field 
 has not become as widely used or further refined as ff94 and its derivatives . in recent work , we devised the implicitly polarized charge model 
 ( ipolq ) , based on a new quantum mechanical 
 method that integrates the condensed - phase environment due to explicit 
 water molecules into quantum calculations at the mp2 level of theory . 
 in this method , the appropriate atomic partial charges of a nonpolarizable 
 model are estimated to be halfway between the charges of the system 
 in vacuum and those of a system fully polarized in the presence of 
 a condensed - phase reaction field potential . the averaging accounts for the polarization energy and instantaneous 
 rearrangements of dipoles in otherwise nonpolarizable models . the 
 choice of tip4p - ew to represent the solution - phase 
 environment while developing the protein force field supports the 
 use of this water model for later simulations . furthermore , the ipolq 
 method provides a basis for understanding why essentially all fixed - charge 
 water models carry a dipole moment of roughly 2.3 d : the dipole is 
 halfway between the true condensed phase dipole of water , estimated 
 to be 2.6 to 2.9 d , and the 1.85 d dipole 
 of water in the gas phase . charges for solutes can therefore be derived 
 under physical approximations that also obtain very similar results 
 to existing water models . while the ipolq approach offers numerous 
 benefits for creating a self - consistent force field , the process of 
 deriving charges is laborious , and our previous work did not assign 
 torsion fourier terms to complete the model . in this article , we will 
 extend the ipolq method to facilitate derivation of complete models , 
 and we will also describe automation of the procedure , which greatly 
 reduces the human effort needed to create new force fields for proteins 
 and drug - like molecules . the amber 
 ff14ipq force field is intended to be a direct alternative 
 to other nonpolarizable amber force fields and contains no new functional 
 forms . however , ff14ipq is also designed to consistently adhere to 
 a specific level of quantum theory , mp2/cc - pvtz . over the course of 
 this study , we found two major obstacles to fitting a molecular mechanics 
 ( mm ) expression to a quantum mechanical ( qm ) potential energy surface 
 ( pes ) and devised our own solutions to each of them . the first concerns 
 the means for adapting the charge set and torsion fourier series terms 
 to work in concert . the second focuses on whether to make corrections 
 for the inability of common mm models to reproduce certain high - energy 
 features of the qm pes . the implicitly polarized charge method produces charges 
 that , given certain assumptions detailed earlier , reflect the appropriate partial atomic charges of a nonpolarizable 
 model intended for simulations in aqueous solution . the atomic partial 
 charges of our ipolq amino acid charge set approximate electrostatic 
 potentials around dipeptides ; these potentials are halfway between 
 the potential calculated for an unpolarized dipeptide ( in numerous 
 conformations ) in vacuo and the potential calculated for the dipeptide 
 ( in the same conformations ) after its electron density has been polarized 
 by a solvent reaction field potential due to the time - averaged tip4p - ew 
 water density sampled around each conformation . the electrostatic 
 potentials needed for fitting charges in the condensed phase are laborious 
 but straightforward to compute . it is less obvious how to extract 
 the internal potential energy of the dipeptides in the condensed phase 
 or make an equivalent averaging to derive an appropriate pes for fitting 
 torsion parameters . zgarbov and colleagues solved this problem by computing quantum mechanical energies 
 of solutes at the hartree fock level of theory in the presence 
 of the cosmo continuum solvent model , removing the solute solvent 
 interaction energy and then comparing to molecular mechanics potential 
 energies computed with a pb solvent model . in effect , they compared 
 the internal energies of the qm and mm systems in the presence of 
 equivalent continuum solvents . in our case , it is not as straightforward 
 to decompose single - point energies at the mp2 level of theory , nor 
 is it as certain whether a quantum - mechanical implicit solvent model 
 could be substituted for tip4p - ew . we instead chose to extend the 
 ipolq methodology and offer an alternative method for deriving torsion 
 fourier series from gas - phase pess when the accompanying charge sets 
 approximate polarization in a condensed - phase environment . it 
 would be simple to derive torsion parameters for a set of charges 
 derived strictly from the electrostatic potentials for solutes in 
 vacuum : compute the single - point energies of many additional solute 
 conformations at the same level of quantum theory , again in vacuum , 
 and fit the mm energies of each conformation to match the qm results . furthermore , the role of the torsion fourier series is to artificially 
 correct errors in the nonbonded electrostatic and van der waals interactions 
 between atoms on either end of a dihedral ; the pes of a rotatable 
 bond is mostly captured by the nonbonded parameters . we therefore 
 sought a way to express the set of ipolq charges for a solute of interest , q , as a perturbation q of charges q derived for the solutes 
 using only the vacuum quantum calculations:2it is simple to compute torsion 
 parameters 
 to match an mm pes q to a qm pes computed 
 in vacuum . changes in the pes of a molecule s rotatable bond 
 due to immersing it in water could then be assumed to be adequately 
 represented by changes in local dipoles as expressed by q . this appears to be a safe assumption , given that q  q and that the 
 van der waals parameters , which also strongly influence the potential 
 energy changes due to rotation about a bond , are constant across both 
 phases . however , the linear least - squares fit from which our ipolq 
 charges q and our proposed vacuum 
 charges q are derived is able to obtain 
 many different solutions with similar levels of accuracy , particularly 
 for buried atoms whose charges are not well determined . in our earlier 
 work , we used large numbers of conformations of each solute in order 
 to guard against this behavior . however , because residual indeterminacy 
 could exist in either q or q , the difference q could still be amplified . we therefore extended the ipolq least - squares 
 fitting procedure to fit both q and q simultaneously , with additional restraint 
 equations to keep q small . the results are 
 solutions for each of q and q related by a minimal , smooth perturbation 
 q . the original matrix equation constructed 
 in the ipolq fitting procedure 
 may be expressed in shorthand:3here , a is the fitting matrix 
 whose elements are derived from the molecular coordinates and the 
 kernel of coulomb s law . each row of a describes 
 one linear equation through which the partial charges of a molecule 
 in some conformation create an electrostatic potential at a nearby 
 point in space . the electrostatic potential found by the target quantum 
 method is stored in the vector u. each column of a solves for an independent charge variable ; if multiple 
 atoms of a molecule are constrained to have the same charges , then 
 they contribute to the same column . appended to a , the matrix v contains additional restraint equations 
 that penalize the least - squares fit when particular charges stray 
 from target values encoded in the vector vt appended to u. a more detailed statement of the 
 linear least - squares problem is as follows:4here , the elements of the fitting 
 matrix a(p , i)j pertain to the influence of all instances 
 of the jth adjustable charge on the ith fitting point around 
 the pth molecular conformation . the solution vector u on the right - hand side is filled with values of the electrostatic 
 potential measured at many points rp , i . as was mentioned , there may be multiple 
 instances of the same adjustable charge in each molecule . a general 
 definition of the elements of a is then5where the summation runs over all atoms c bearing 
 the fitted charge described by variable j. the various 
 molecular conformations may be different 
 poses of the same molecule or even a collection of poses of many different 
 molecules . in the latter case ( which best describes the original ipolq 
 fit we performed ) , the matrix a becomes a sparse 
 matrix due to the fact that not all systems contain all charge variables , 
 but it is not often sparse enough to warrant a special storage format . this can require a large amount of memory , but simultaneous solution 
 of all charge variables permits multiple systems to share the same 
 charges . in our ipolq fit , for instance , we chose to follow the cornell 
 charge set convention of giving similar charges to backbone n , h , 
 c , and o atoms involved in the peptide bonds for neutral , positively 
 charged , and negatively charged amino acids , compressing what could 
 have been more than 80 charge variables into just 12 , in the hope 
 that this would expedite development of a common set of backbone torsion 
 potentials . the restraint matrix v is appended to 
 the fitting matrix a in order to keep the charges 
 of certain atoms small and also to restrain the overall charges of 
 certain groups of atoms : the fit is heavily penalized by vnc = 1.0  10 kcal / mole if the sum of atomic partial charges of neutral or ionic 
 residues differs from the appropriate net charge ( nc ) , qk for the kth system . 
 in eq 5 , p , j is 1 if the pth system contains 
 the jth charge variable . to put all restraint constants vsa , j for specific atomic charges 
 on the same scale , these constants were set proportional to the number 
 of times nj the jth charge variable appeared in a6 in our ipolq fit , we set vsa0 to 1.0 
  10 kcal / mole . the extended least - squares fit for deriving ipolq 
 charges as a 
 perturbation of vacuum charges solves for twice as many variables 
 based on the same molecular conformations . the extended problem to 
 solve is simply7 in this system 
 of equations , the original 
 fitting matrix is replicated into several blocks of the extended system , 
 and the solution vector contains electrostatic potentials computed 
 from each molecular conformation in vacuum as well as the vacuum potential 
 averaged with the electrostatic potential computed in the condensed 
 phase u. the original restraint equations 
 are still present , but they apply only to the vacuum charges q ; additional restraint equations are added 
 to strongly force the sum of all perturbation charges in each system 
 to zero . an additional block matrix loosely restrains perturbation 
 charges individually to zero using stiffness constants vgp scaled by the number of instances of each ( perturbation ) 
 charge variable nj ; numerical 
 values of vgp are discussed below . in eq 7 , i is the identity matrix ; 
 every perturbation charge is individually restrained . the extended 
 system , which requires slightly more than four times the memory of 
 the original problem and eight times the computation cost , yields q by applying eq 2 . the traditional least - squares approach for fitting torsion parameters 
 is composed as shown in eq 8 . the fitting matrix 
 is composed of the kernels t(p , i)j of the jth torsion term to be fitted for the ith conformation 
 of the pth system , while other columns store energy 
 adjustment constants cp applied to the pth chemical system present in the 
 fitting data . the solution vector u in this equation 
 contains the qm single - point energies of each conformation , less the 
 average single - point energy of all systems with the same chemical 
 composition , less the mm energy arising from other terms in the force 
 field . the torsion terms are fitted so that the sum of their contributions 
 compensates for inaccuracies in the molecular mechanics model to bring 
 the relative energies estimated for different conformations of the 
 same molecule into agreement with the qm pes.8 in principle , 
 it would be feasible to leave 
 the average single - point energy of each chemical system in u , but this would lead to a poorly conditioned matrix , as 
 the energy adjustment constants cadj , 
 here c1 ... ck , would need to take on very large values . it is 
 also possible to compose u based on system conformational 
 energies relative to the structure with the lowest quantum mechanical 
 energy and to omit the constants cadj for 
 each system . however , no quantum method is perfect , and omission of cadj would amplify the influence of a particular 
 conformation over the fitted parameters . it is well - known to force field developers 
 that the bond and angle terms that approximate high - frequency motions 
 in molecular simulations are not consistent with the qm target models ; 
 aside from the breakdown of the harmonic approximation with increasing 
 strain , the equilibrium bond length or angle is dependent on the chemical 
 context . the inconsistency elevates the mm energies for structures 
 optimized by qm methods and is therefore a sort of contaminant when 
 fitting lower - frequency degrees of freedom , in particular torsion 
 fourier series terms . concern about introducing error by this method 
 has led numerous groups to fit torsion fourier series by positing 
 that a  conformation  of a molecule is defined strictly 
 by its torsional degrees of freedom and making two slightly different 
 variants of that conformation : relaxing all bond , angle , and nonbonded 
 degrees of freedom by mm and qm approximations , respectively . the 
 objective then becomes to fit torsion parameters such that the coordinates 
 optimized by mm produce an energy most like the single - point energy 
 found for the nuclear coordinates optimized according to qm . this 
 method , hereafter the tandem optimization approach , relaxes much of 
 the inconsistency arising from high - frequency degrees of freedom in 
 the two pess , but it also introduces a new source of error in the 
 nonbonded interactions between the different coordinate sets . while 
 many investigators have accepted the trade - off , we performed an independent 
 analysis , which led us to reject the approach in favor of a direct , 
 one - to - one mapping between coordinates and energies . we analyzed 
 the trade - off between error removed by relaxing high - frequency degrees 
 of freedom and error introduced by changes in nonbonded interactions 
 by computing mm energies for mm and qm optimized variants of 648 conformations 
 in each of 15 amino acid dipeptides ( over 9000 pairs of energies ) . the conformations sampled 1 and 2 at 20 intervals while the backbone was held in an -helical 
 or -sheet conformation ; 1 , 2 ,  , and  dihedrals were held fixed during both optimizations . mm optimizations and all mm energy evaluations were performed with 
 a new variant of the amber ff99sb force field provided to us courtesy 
 of james maier and professor carlos simmerling s research group . qm optimizations , also performed for us by james maier , were performed 
 at the mp2 level with the 6 - 311++g basis set . the mm and qm approximations are similar to those we have chosen 
 for development of amber ff14ipq . first , we focused on the difference 
 in the molecular mechanics 
 energy computed for each variant of a given conformation ; by construction , 
 the mm optimized variant was always lower in energy , as scored by 
 amber ff99sb . table 1 gives the differences 
 computed for 628 conformations of each dipeptide . by decomposing the 
 mm energy , it is apparent that the mm and qm models do disagree about 
 the optimal bond lengths and angles , and the energy differences in 
 there terms are greater , sometimes much greater , than the energy differences 
 arising from nonbonded terms . at face value , this result would appear 
 to support the tandem optimization approach . each difference 
 is given as an 
 average  standard deviation , in kcal / mol . includes lennard  jones ( lj ) 
 and electrostatic ( elec ) contributions . when fitting torsion fourier series to match a molecular 
 mechanics 
 pes to a qm pes , the objective is correct relative energies ; the mm 
 pes is otherwise very far removed from the qm pes , which includes 
 factors such as electron  nuclei interactions and nuclear repulsion . the mean energy of the qm pes is therefore subtracted , and an additional 
 constant cadj is included in the fit to 
 arbitrarily adjust the energy of a particular molecule , regardless 
 of conformation , up or down , to bring the two pess into agreement . because of this , the mean energy differences between the mm and qm 
 optimized variants of each conformation will  fall through  
 the fit , absorbed into cadj . if the molecular 
 mechanics and quantum models optimize bonds to different lengths , 
 but the disagreement is consistent across all conformations of the 
 molecule , then the mismatch in these high - frequency terms will have 
 no effect on the fitted torsion parameters . standard deviations of 
 the difference , the variability of the disagreement between the qm 
 and mm pes , is a much better indicator of possible contamination in 
 a parameter fit . table 1 shows that bond energy 
 differences have a small deviation but angle energy differences carry 
 a much higher deviation and therefore might contaminate the torsion 
 parameter fit . comparing the standard deviations in the energy differences 
 arising from the sum of bond and angle terms to those arising from 
 all nonbonded terms suggests that the tandem optimization approach 
 offers a marginal benefit in most cases but is detrimental in cases 
 such as ser and thr . the total nonbonded energy differences 
 appear to be smaller than 
 those of bonds , but the nonbonded energy is a sum of many interactions . the major amber molecular dynamics engines conveniently break nonbonded 
 interactions into  14  contributions between 
 atoms at either end of a dihedral group and  all other  
 contributions . these two parts of the lennard  jones interactions 
 are uncorrelated , but the sum of electrostatic 14 interactions 
 tends to be strongly anticorrelated with the sum of all other electrostatic 
 contributions , as shown in table 2 . as a consequence , 
 the standard deviations of the differences in the electrostatic 14 nonbonded terms are as large or larger than the deviations in total 
 nonbonded energy differences . jones short - ranged interactions . ) these 14 
 interactions form the base of the torsional pes around each rotatable 
 bond and , of all nonbonded interactions , are the most strongly connected 
 to the torsion potentials . this suggests that the contaminant introduced 
 by the tandem optimization method , a mismatch in the nonbonded interactions , 
 is actually worse than the contaminant being removed , the mismatch 
 in energies due to high - frequency modes . like table 1 , this table 
 compares mm energies computed for each of two optimized 
 variants of a dipeptide conformation . we conclude that it is at least as sane to accept 
 errors that might 
 be introduced by a mismatch between molecular mechanics bond and angle 
 geometries and those of the quantum target as to accept errors arising 
 from a mismatch in electrostatic or steric interactions . for this 
 reason , we chose to submit the mm optimized coordinates , which reflect 
 the states that the molecular mechanics model will actually explore 
 in simulations , to single - point quantum calculations and demand that 
 our fitted molecular mechanics model reproduce the qm single - point 
 energies calculated for precisely the same set of nuclear positions ; 
 details of our procedure for generating the actual fitting data for 
 torsion fourier series can be found in the methods . due to the one - to - one mapping of coordinates and energies , the 
 fitting data suffers some contamination from angle terms in the mm 
 approximation being inadequate to describe the qm pes . these terms 
 are roughly 1 order of magnitude higher in energy than the torsion 
 terms or nonbonded interactions in the mm model ; as such , they could 
 make erroneous contributions to the energy with the risk of torsion 
 parameters becoming fitted against noise . however , the majority of 
 angle strain is orthogonal to the torsional subspace of molecular 
 motions . with adequate sampling of the torsional degrees of freedom 
 in the fitting set , and with a consistent optimization of each set 
 of coordinates relative to one model or the other we designed the amber ff14ipq force 
 field to be similar to previous 
 amber nonpolarizable force fields . bond and angle parameters were 
 taken from the existing ff99sb force field , along with most lennard  jones parameters . atomic partial 
 charges were refitted according to the updated ipolq fitting procedure 
 described in the theory section , using the 
 same data produced in our previous study , but for future force field development , we automated the ipolq fitting 
 cycle in the amber mdgx program . in our 
 derivation of ipolq charges for amino acid side - chain analogues , we 
 adjusted the lennard jones  parameters of polar atoms 
 in order to bring the computed hydration free energies into agreement 
 with experiment . in this respect , the development philosophy of ff14ipq 
 resembles that of the new gromos 54a8 force field . however , the interactions of these atoms with water were 
 essentially the only changes that were of consequence during our hydration 
 free energy calculations , and we found that larger  radii fitted 
 for most of the atom types intensified 14 nonbonded repulsion , 
 making torsion parameters much more difficult to fit . to bring the 
 lennard  jones changes into ff14ipq , we made them applicable 
 only between the polar atom types and the tip4p - ew oxygen type . the 
 lennard  jones lorentz  berthelot mixing rule is broken 
 for the interactions of these atom types with water , in the manner 
 that charmm36 makes use of nbfix terms . with the bonded and nonbonded parameters established , the majority 
 of the development in ff14 focused on torsion fourier series terms . a new module 
 was added to the mdgx program to aid users in computing 
 ipolq charges for arbitrary molecules . for a standard resp procedure , 
 researchers must have a set of conformations of their molecule of 
 interest ; the conformations serve as inputs to tools such as the red 
 server , which manages the necessary quantum calculations and performs 
 the restrained charge fit . with the new ipolq module , researchers must have a set of conformations 
 of their molecule immersed in the solvent of interest . the mdgx program 
 will read the solvated conformations as restart files along with an 
 appropriate topology and begin dynamics with the solute molecule held 
 in a fixed position . the mdgx program automates the process of collecting 
 the solvent charge density , computing the solvent reaction field potential 
 ( srfp ) , preparing inputs to a quantum program , and launching the calculations 
 to ultimately produce grids of electrostatic potential computed for 
 the molecule in vacuum and in the influence of the srfp . the mdgx 
 ipolq module improves on the original protocol , always applying a 
 shell of point charges around the average solvent charge density taken 
 from the simulation . the charges in the shell are fitted to reproduce 
 the srfp due to infinite electrostatics present in the simulation 
 in the context of a quantum calculation on an isolated system . furthermore , 
 users can specify up to three concentric charge shells to increase 
 the accuracy of the srfp , evaluate the srfp at additional sites throughout 
 the solute volume , and even specify an interior shell of charges to 
 reproduce the srfp in and around the solute without including point 
 charges nearer than an arbitrary distance from solute atoms . the final 
 feature could be useful for researchers who are concerned about qm 
 basis functions adversely interacting with solvent charges , although 
 we did not notice any effects on the ipolq results for some systems 
 we tested with mp2/cc - pvtz calculations used in the original ipolq 
 protocol ( data not shown ) . in summary , mdgx manages a cycle of the 
 ipolq procedure for a conformation of the molecule of interest ; the 
 output of independent runs on multiple solute conformations can be 
 pooled and sent to the mdgx charge fitting module to derive a new 
 charge set , which can then be used to update the solvated system s 
 topology and start another round of calculations with the ipolq module . 
 in this manner , mdgx manages nearly all aspects of the ipolq procedure 
 until the solute charge model converges . the module supports both 
 gaussian and orca quantum chemistry packages , performs dynamics with modest 
 parallelism and efficient cpu execution , and launches quantum chemistry 
 programs for parallel execution in accord with the molecular dynamics 
 run to avoid wasting cpu cycles . data for the fitting matrix a was derived from the same quantum data as the original 
 ipolq charge set . the same protocol was followed for selecting fitting 
 points from the electrostatic potentials evaluated around each conformation , 
 but , due to memory constraints in the much larger matrix problem , 
 only 3750 points per conformation were selected . ( our earlier work 
 showed that anywhere from 3000 to 5000 points per conformation yielded 
 convergent results in the fitted charges . ) this method implies two 
 charge sets , one valid in vacuo and the other in solution . while the 
 complete release version of ff14ipq takes the charges valid in solution , q + q in eq 7 , we also considered a variant based on q in some studies of systems in vacuo , hereafter named 
 v - ff14 and not distributed for general use . the fitting set for amber ff14ipq was created 
 by molecular simulations and energy optimizations performed with the 
 amber ff99sb force field and new variants under development in the 
 simmerling group . all of these mm models are based on the cornell 
 charge set and have been edited over more than a decade . the most 
 significant feature of the force fields , to us , was the similarity 
 in form and some parameters to the proposed amber ff14ipq force field . 
 while it is impossible to completely sample the available conformational 
 space in the fitting data , we chose the most recent amber force fields 
 to simulate each molecule in the hope that the these models would 
 sample the unrestrained degrees of freedom in a manner that reveals 
 their biases so that those biases could be eliminated in ff14ipq . nevertheless , critical dihedral angles near the center of each molecule 
 were restrained to various values as we populated the fitting data 
 set ; the exact details of the simulations , such as cutoff , time step , 
 and solvent model , are therefore of secondary importance to the choice 
 of protein - like molecules , the array of restraints , and the number 
 of snapshots collected for single - point energy calculations . the most common structures included in our fitting data set were 
 blocked dipeptides , ace - xx - nme , where xx is an amino acid , although 
 we also included ace - ala - ala - ala - nme and ace - gly - gly - gly - nme tetrapeptides 
 to sample backbone  and  angles and ace - xx - xx - nme tripeptides 
 to further reduce the influence of the blocking groups on these backbone 
 dihedrals . in all , the torsion fitting data set consisted of nearly 
 28  000 structures whose single - point energies were computed 
 in vacuum by mp2/cc - pvtz calculations . to sample side chain 
 conformations , we obtained a set of some 17  000 
 dipeptide conformations for neutral asp , asn , deprotonated asp , cys , 
 -protonated his , -protonated his , ionic his , ile , leu , 
 phe , ser , thr , trp , tyr , and val from james maier and the simmerling 
 group . each amino acid was sampled in 628 conformations , exploring 
 1 and 2 angles at 20 intervals 
 while holding the backbone in either an -helical or -sheet 
 conformation . these amino acid conformations were among those used 
 in our analysis of the effects of unrelaxed high - frequency degrees 
 of freedom described in the theory section . 
 to obtain a better sampling of backbone conformations and to reduce 
 the risk of side chain dihedrals becoming coupled to particular backbone 
 conformations , we generated another 180 conformations of each dipeptide 
 by sampling  and  individually at 20 intervals 
 with no restraints on any other degrees of freedom . to decorrelate 
 the unrestrained degrees of freedom , we ran molecular dynamics simulations 
 of each dipeptide while progressively advancing the  and  
 restraints . dynamics were run in baths of tip4p - ew water at 300 k , 
 with 100 ps intervals between each snapshot . or  were 
 incremented with every snapshot , but each angle was completely rotated 
 five times ; as a result , separate snapshots at the same  or 
  coordinate were spaced by 1.8 ns of dynamics . to sample backbone conformations , reduce the dependence of the 
 fitted parameters on ace and nme blocking groups , and increase side - chain 
 sampling , we included a set of over 7000 tripeptides . we sampled all 
 20 amino acids with ala , gly , or ser adjacent in either direction 
 in the peptide chain , restraining  and  individually 
 at 15 intervals with no other restraints and running dynamics 
 on each system at 450 k to generate multiple conformations at each 
 value of  or . because there were so many of these systems , 
 a generalized born solvent was used to 
 generate the conformations . to further sample 
 backbone conformations of two key amino acids , 
 alanine and glycine tetrapeptides were sampled in 1296 conformations 
 each , sampling the central residue s  and  space 
 at 10 intervals on a two - dimensional grid . aside from simultaneous 
 restraints on  and  , no other restraints were used as 
 dynamics were performed for 100 ps between snapshots to decorrelate 
 other degrees of freedom . snapshots of molecular dynamics trajectories 
 are not suitable for 
 direct incorporation into quantum calculations because of the potential 
 for bond and angle high - frequency energy terms to contaminate the 
 fitting data , as explained in the theory section . 
 however , if the contributions from these terms are relaxed out with 
 respect to either a mm or qm approximation , then the single - point 
 qm energies calculated for each conformation then constitute an acceptable 
 fitting set . we optimized each conformation with restraints on dihedrals 
 involving four heavy atoms ( but not hydrogens ) in order to preserve 
 the diversity of conformations created in our simulations while relaxing 
 as much angle and bond strain as possible . we chose to optimize conformations 
 with respect to the mm approximation because the calculations were 
 cheap , allowing us to devote more time to evaluating single - point 
 energies . also , we felt that it was more instructive to evaluate the 
 local minima of an mm model with respect to the mp2/cc - pvtz target , 
 to examine states that an mm model might have a propensity to populate 
 and to verify their energies with qm in the fitting data . all torsion parameters 
 for the amber ff14ipq force field were fitted 
 simultaneously . while the length of each fourier series and also the 
 phase angles were taken from ff99sb , new atom types added by the simmerling 
 group as well as new atom types required by the ipolq charge set for 
 amino acids were included , and the glycine c and 
 proline backbone nitrogen atoms were given their own unique atom types 
 distinct from other amino acids . including a new atom type implied 
 replicating all bond , angle , and fourier series terms pertaining to 
 the original type , thereby creating additional parameters to fit . all fourier series amplitudes were reoptimized by the standard linear 
 least - squares approach in a single matrix equation to 
 obtain the best overall fit for terms that appear in different contexts 
 across multiple residues , as described in the theory section . another new module of the mdgx program was created to perform 
 the atom type branching , set up the matrix equation , and perform the 
 linear least - squares fit . while it was an advantage to have 
 the extensive fitting set described in the preceding section , all 
 torsion fourier terms were still restrained loosely toward zero by 
 a penalty of 2.0  10 kcal / mol ( multiplied 
 by the number of times each parameter appeared in the fitting matrix ) 
 to keep the amplitudes small and avoid overfitting . torsion amplitudes 
 were optimized to make molecular mechanics energies of the di- , tri- , 
 and tetrapeptide systems computed in vacuo with charges from q agree with mp2/cc - pvtz single - point energies 
 of the peptides also computed in vacuo . as explained in the theory section , we assume these torsion parameters 
 to be transferable to describe the behavior of solvated peptides when 
 paired with q. as a preliminary 
 test of force field performance , we evaluated fidelity to the underlying 
 mp2/cc - pvtz benchmark over the course of molecular mechanics energy 
 minimization . the molecular conformations in the training set were 
 already optimized with respect to a very similar mm model , but in 
 the presence of one or more restraints on dihedral angles . energy 
 minimization of each conformation of all systems found in the fitting 
 set was performed in vacuo with v - ff14 , the variant of ff14ipq substituting q for the implicitly polarized charges found 
 in the release version . new qm energy calculations were then performed 
 on the resulting structures optimized by v - ff14 . these new conformations 
 and energies were added to the training set , and new energy optimizations 
 were performed in an iterative fashion until v - ff14 could score structures 
 created by its own optimization consistently with those in its training 
 set . building on the fitting data of small peptides , 
 we computed potentials of mean force ( pmfs ) for pairs of dihedral 
 angles in blocked dipeptides at 298 k with a standard two - dimensional 
 umbrella sampling technique . we collected data in 1296 windows 
 spaced by 10 in  and  or in 1 and 2 , depending on the amino acid and pmf . each 
 window was seeded from a continuous , incrementally restrained simulation 
 similar to that used to generate conformations for the corresponding 
 tetrapeptides prior to mp2/cc - pvtz calculations . after seeding , each 
 window dihedral 
 angles were initially restrained by an 8 kcal / molrad harmonic penalty function , but if some of the 10 bins were 
 left undersampled , then more windows were added on a 5 grid 
 with 16.0 kcal / molrad dihedral restraints to completely 
 fill out the pmf . we began with simulations of 
 penta - alanine ( hereafter , ala(5 ) ) , the -helix k19 , and -hairpins 
 chignolin ( starting structure pdb entry 1uao(32 ) ) and the 
  gb1 hairpin from the c - terminal fragment of protein 
 g. furthermore , we performed microsecond length simulations of the 
 globular proteins gb3 ( starting structure pdb entry 1p7e , similar to protein g and containing a motif homologous 
 to the gb1 hairpin ) and lysozyme ( 4lzt ) in solution . each system was equilibrated with protein backbone atoms held under 
 progressively decreasing restraints for up to 7 ns , depending on the 
 size of the system . all simulations were performed with tip4p - ew water 
 ( the water model used to develop our solution phase charge set ) in 
 sufficient quantity to enclose the peptide and solvated protein systems 
 by at least 10.0  within octahedral boxes after equilibration 
 under constant pressure dynamics . nonbonded interactions were calculated 
 with a 10.0  cutoff on lennard  jones interactions , a 
 homogeneity approximation for long - ranged van der waals interactions , 
 and smooth particle - mesh ewald electrostatics with direct and reciprocal 
 space accuracies near 1 part in 200  000 ( this corresponds to 
 a direct sum tolerance of 5.0  10 and a 
 direct space cutoff of 9.0  with the default mesh grid spacing 
 of up to 1  ) . a 2 fs time step was used in all simulations , 
 along with the shake and settle algorithms to constrain the lengths of bonds 
 to hydrogen . most of the small systems were simulated at 277 k , by 
 a langevin thermostat with collision 
 frequency 3/ps , to replicated nmr conditions . the gb1 hairpin system 
 was simulated at 298 k to investigate -sheet stability at room 
 temperature . the updated ipolq procedure described in the theory section uses the same fitting data as the original method to derived q , but it also derives q , a set of charges appropriate for modeling electrostatics 
 in vacuo . these two charge sets would be independent , and the updated q would be the same as the original , but 
 for the fact that the restraint equations used to temper the fit of 
 the original q instead temper values 
 in q , and the perturbation q that relates q to q by eq 2 is subject 
 to its own set of restraints . in our protocol , the size of q is controlled by setting the parameter vgp in eq 7 . larger values of q could make it difficult to transfer torsion parameters 
 developed for systems in vacuum to simulations in solution , but too 
 small a q would imply q  q , and neither 
 set of charges would describe its intended environment well . we chose 
 to set vgp to 0.005 kcal/(mole - instance ) , where e is the 
 charge of a proton . this is half the stiffness of restraints holding 
 charges of underdetermined sites to zero in this and earlier ipolq 
 derivations . the fitting matrix is not well conditioned in either 
 case : for main - chain amino acids , we found the ( 2-norm ) condition 
 number of the matrix a in eq 3 to be 1.8  10 and that of the extended matrix 
 in eq 7 to be 2.2  10 . the 
 possibility of a worst - case loss of precision compels us to use double 
 precision arithmetic for what is already a very large matrix , but 
 the results are nonetheless consistent . the earlier ipolq protocol , 
 for instance , judged convergence of the charge set by the point at 
 which successive generations changed the fitted charges by less than 
 the inclusion of one of the fitting conformations . ( these changes 
 were on the order of 0.01e . ) as before , we are more 
 concerned with the effect that a different set of constraints might 
 affect the outcome . the value of vgp we chose seems to 
 have minor effects on the overall accuracy of q and q. averaged over all 
 systems and conformations , the original q reproduces the target electrostatic potential with 1.82 kcal / mole root - mean - squared error ( rmse ) , whereas q + q reproduces the same potentials 
 with 1.89 kcal / mole rmse . the effects on selected 
 systems are shown in the lower panel of figure 1 . over the range we tested vgp , the effects 
 on electrostatic potential rmse appear to begin to approach an asymptotic 
 limit ( half the difference between the electrostatic potentials computed 
 in vacuo and in the condensed phase ) . if vgp is tuned even lower , then the rmse of q + q actually improves over the original q due to the fact that the perturbation charges 
 are then much less restrained than the original charges . however , 
 tuning vgp too low leads to undesirable 
 effects such as ala c taking a charge of 0.21 
 at vgp = 0.0025 , as opposed to 0.04 
 in the original q. accuracies and 
 charge perturbations for all amino acids at vgp = 0.005 are shown in table 3 . effects of 
 the coupling constant vgp on q and the accuracy of electrostatic potentials 
 for condensed - phase systems . vgp determines 
 the strength of the harmonic penalty restraining all charges q toward zero . values of the root - mean - squared error ( rmse ) 
 describe electrostatic potentials projected by each dipeptide s 
 molecular mechanics charge set relative to the qm target . vgp = 0 corresponds to partial charges by the 
 original ipolq method before the extension , 
 which allows us to express them as a perturbation to charges appropriate 
 for simulations in vacuo . errors and charge variations 
 are expressed for a restraint of vgp = 
 0.005 kcal/(mole - instance ) applied 
 to all perturbation charges . root - mean - squared error ( rmse ) of 
 fitted mm charges in replicating the qm target electrostatic potential . the target is the average electrostatic 
 potential of the solute s mp2/cc - pvtz wave function in vacuum 
 and in the solvent reaction field potential due to tip4p - ew water . the target is the electrostatic 
 potential of the solute s mp2/cc - pvtz wave function in vacuum . maximum absolute deviation 
 in partial 
 charges unique to this residue ; backbone atoms frequently showed q of 0.060.10 , as shown in the supporting information . at vgp = 0.005 , the majority of values 
 in q are smaller than 0.05e , as shown by the top panel in figure 1 ; the 
 larger values tend to be in buried methyl carbons , whose charges are 
 small to begin with , and the atoms of polar head groups . neither of 
 these changes are likely to make torsion parameters developed with q less transferable because the electrostatic 
 potential energy surface will change only if atoms with large values 
 of q can rotate around nearby atoms possessing 
 large of charges of their own . the perturbation does not , in fact , bring q + q to the polarity 
 of the original q : the carbonyl carbon 
 atom becomes less positively charged by as much as 0.05e , the oxygen becomes less negatively charged by 0.03e , and the polarity of the n however , the perturbation does make the backbone significantly more 
 polar than q alone would describe it . 
 as a consequence , the ipolq protocol continues to model protein backbones 
 with more polarity than the cornell and duan charge sets , maintaining 
 a principal finding of our earlier study . torsion fourier series 
 terms for all amino acids were fitted by 
 a linear least - squares approach as described in the methods . the preliminary set of 28  000 structures and 
 mp2/cc - pvtz single - point energies is , to our knowledge , the most extensive 
 potential energy surface ever employed for this type of molecular 
 mechanics parameter development . in addition to the atom types present 
 in amber ff99sb and new atom types for c atoms introduced 
 by the simmerling group , we included new types for each of the lennard jones 
 modifications introduced to adjust hydration free energies of amino 
 acid side - chain analogues in the initial ipolq derivation . we also 
 included new types for the glycine c atom and proline 
 backbone nitrogen because of the unique chemical and bonding structures 
 around these atoms . when the electrostatics of each di- , tri- , 
 and tetrapeptide are described by q , the fitted torsion fourier series terms complete a force field 
 and describe the molecular mechanics energy of the system in vacuo . the rmse of these energy estimates relative to the mp2/cc - pvtz single - point 
 energies for di- and tetrapeptides is given in table 4 . also given in this table are the molecular mechanics ( mm ) 
 energies obtained by reoptimizing a much smaller set of torsion fourier 
 series terms , the parameters found in amber ff99sb , in conjunction with our newly derived q and the mm energy estimates that would have been obtained 
 if no torsion terms were used . the overall contributions from the 
 torsion terms are often small , only reducing the rmse of mm energy 
 estimates by 1.3 to 1.5 kcal / mol relative to a model that has no such 
 terms . the overall size of the torsion terms contributions 
 may understate their importance , given the number of published force 
 field improvements based on changes in these terms . in all cases , 
 the charge set q fitted to reproduce 
 the electrostatic potentials 
 of blocked dipeptides in vacuo was used to estimate the molecular 
 mechanics energy of each blocked dipeptide in vacuum . all energies 
 are given in kcal / mol . the 
 v - ff14 force field : q has been substituted 
 for the implicitly 
 polarized charge set in the release version . the ff99 force field , with q as derived for v - ff14 ( identical to the 
 force field in the first column , but with a smaller torsion parameter 
 space ) . the new ff14ipq force 
 field has many more torsion fourier series 
 terms than the ff99sb force field : 427 to 67 . most of the new parameters 
 were added by including the simmerling group s unique c atom types . while the total number of torsion parameters 
 increases nearly 7-fold from ff99sb to ff14ipq , table 4 shows that the number of parameters expressed in any particular 
 system doubles at most . new atom types c8 , 3c , and 2c for c have added new parameters for 1 and 2 , which distinguish the side - chain rotamer energetics for 
 protonated his , arg , and lys , the amino acids ile , thr , and val , and 
 other amino acids . this partitioning is just one of many possible 
 approaches , but it seems to have achieved a similar effect to an earlier 
 extension of the ff99sb force field , which improved the rotamer propensities 
 of residues ile , leu , asp , and asn by adding new atom types and fitting 
 the newly minted torsion parameters to mp2/(aug)cc - pvtz energy profiles . the ff14ipq parameter space makes notable improvements 
 over the ff99sb parameter space for ile , asp , and asn ; it appears 
 that providing only a single atom type for c forces 
 the same set of torsion parameters to average between disparate potential 
 energy surfaces and causes difficulty fitting the data . tg  atom type for glycine c likewise decouples two distinct potential energy surfaces : the dihedral 
 angles between a backbone -hydrogen and the backbone polar 
 hydrogen , or the backbone carbonyl oxygen , are sampled twice ( roughly 
 120 apart ) in glycine and only once in all other amino acids . 
 these dihedrals contribute directly to the protein  and  
 propensities , and the new tg atom type improves the fit for glycine 
 tetrapeptide by roughly 0.5 kcal / mol rmse and the fit for all other 
 dipeptides by up to 0.1 kcal / mol rmse ( data not shown ) . we also tried 
 adding distinct atom types to the glycine -hydrogens , but this 
 was not as effective as adding the tg atom type . in contrast , distinguishing 
 the proline backbone nitrogen as tn has a negligible effect on the 
 fit for any residues but the proline itself . given the inherent difficulty 
 in sampling dihedrals related to this atom type , we considered whether 
 to include it at all , but it does decouple the proline  profile 
 from other amino acids and serves as a handle for future parameter 
 development . with the expanded parameter set , the torsion potentials 
 contribute 
 more and more to the total molecular mechanics energy ; their total 
 contributions in ff14ipq are frequently double their contributions 
 in the refitted ff99 parameter set . however , as shown by the final 
 column of table 4 , a model with no torsion 
 potentials is only perhaps 1 to 2 kcal / mol less accurate . the larger 
 torsion potential contributions in ff14ipq are being counterbalanced 
 by larger energy adjustment constants cadj , which never appear in the force field . the number of adjustable 
 terms and the size of cadj are only weakly 
 correlated ( pearson coefficient 0.40 ) , perhaps because many copies 
 of the torsion terms can appear in the total mm energy and also because 
 the amplitudes of the terms themselves are so variable . despite this , 
 the sizes of cadj may be a useful indicator 
 of whether a model is overfitted . the increases in cadj led us to track 
 the sampling of each torsion parameter throughout the fitting set , 
 as shown in the supporting information . 
 with our exceptionally thorough data set , the torsion parameters related 
 to rotatable bonds seem to be well - sampled , on the whole and in each 
 system - dependent context in which they appear . even with this degree 
 of coverage , however , these fitted torsion parameters are not the 
 final settings distributed as ff14ipq , whether in combination with 
 the partial charge set q or q ; they are , rather , a first draft . the following 
 section presents results obtained when the draft model was allowed 
 to guide geometry optimizations and thereby expose the ways in which 
 its parameters could conspire to accumulate errors over the course 
 of a simulation . energy optimization served as a preliminary test 
 of our fitting program and also of the robustness of our fitting data . energy minimization of the fitting set s structures using the 
 amber pmemd program and the preliminary v - ff14 force field confirmed 
 that the mm energy computed for each initial structure matched that 
 produced by the mdgx fitting module after solving eq 8 . after energy minimization in vacuo , some initial structures 
 converged to the same final configuration , but a number of local minima 
 were still produced for each di- , tri- , and tetrapeptide system . we 
 computed mp2/cc - pvtz single - point energies for new structures in three 
 tripeptide systems shown in figure 2 and compared 
 them to the mm energy estimates to test whether the new model s 
 427 parameters were prone to overfitting . energies of dipeptide 
 and tripeptide conformations before and after 
 energy minimization with the preliminary v - ff14 force field . all molecular 
 mechanical energies are adjusted according to the adjustment constants 
 found while fitting torsion parameters ; quantum mechanical energies 
 are normalized to a mean of zero . hence , the energies of conformations 
 found in the fitting data ( black diamonds ) lie directly on the trendlines , 
 and the energies of system conformations optimized according to the 
 preliminary v - ff14 ( red , open squares ) may not track it . the results in figure 2 show 
 that v - ff14 
 is often able to guide small peptides into configurations that mp2/cc - pvtz 
 calculations agree have lower potential energy . however , after energy 
 minimization , v - ff14 tends to estimate the energy of each configuration 
 as being more favorable than the mp2/cc - pvtz single - point energy . the degree to which v - ff14 departs from the 1:1 trendline with its 
 benchmark may be small , as shown in table 5 . however , geometry optimizations in systems with his , arg , and lys 
 led to mm energies that departed severely from the qm benchmark , despite 
 the agreement maintained in the training set . by adding the structures 
 freely optimized by v - ff14 back into the training set , the second generation 
 v - ff14 nearly eliminated departures from the qm benchmark in arg and 
 his and also showed minor improvements in its ability to optimize 
 the structures of most other residues . v - ff14 
 is able to guide unrestrained 
 energy minimizations of structures in its own training set and reduce 
 the internal potential energy by up to tens of kcal / mol . ( this may 
 be realistic , as most training set structures were restrained in one 
 or more torsional degrees of freedom . ) however , when re - evaluated 
 at the mp2/cc - pvtz level , the resulting structures were often not 
 as optimal as molecular mechanics depicted . negative numbers in the 
 table indicate that v - ff14 strayed from its mp2 benchmark and estimated 
 its optimizations to be too favorable . gaps in the table indicate 
 that a system was omitted from one generation , due to compute cluster 
 downtime or sufficiently low error in the previous generation . lysine presented more of a challenge : 
 as shown in figure 3 , the second generation 
 v - ff14 apparently contained 
 a new artificial minimum that was subsequently found in all optimizations 
 of the dipeptide . we examined the molecular mechanics energies of 
 lysine conformations in the original training set , the second - generation 
 training set , and the final training set with respect to each generation 
 of torsion parameters . torsion parameters describing two of the rotatable 
 bonds generate the severe departures from the benchmark seen in the 
 first generation of the force field . first , a wildcard parameter with 
 3-fold periodicity describing rotation of the lysine side chain amino 
 terminus takes on an amplitude of 4.99 . in the original training set , every conformation sampled the orientation of the terminus in 
 a staggered conformation of the hydrogens at the crest of the cosine 
 wave . the eclipsed conformation , favored by nearly 10 kcal / mol by 
 this spurious parameter , was not sampled in the training set , and 
 the result was six permutations of amino and aliphatic hydrogen interactions 
 creating a nearly 60 kcal / mol fictitious energy release by adopting 
 the unnatural , eclipsed conformation in every structure optimized 
 by the first generation of the force field . a more serious problem 
 occurred in the four - term fourier series describing interactions between 
 the atom type c8 , coined by the simmerling group for lysine and arginine 
 c atoms , and the backbone carbonyl carbons . because of its 
 unique appearance in lysine and arginine , this fourier series strongly 
 influences the backbone  angles adopted by the residues . the 
 original training set did not sample a 100 arc in  for 
 either residue , however , and the fourier series became fitted to produce 
 an unnaturally low energy in the first generation , approximately 14 
 kcal / mol more favorable than any conformation of the  angle 
 should have allowed . only some of the lysine structures fell into 
 this trap when optimized by steepest descent energy minimization under 
 the first generation of the force field , and the result was the two 
 striations , which can be seen in figure 3 . 
 in the second - generation training set , sampling in these torsion parameters 
 improved considerably by including the structures trapped in the first 
 generation s spurious minima , as shown in table 6 . rotation of the lysine amino terminus needed to be sampled 
 completely in order to get a transferable model , but this was accomplished 
 by the third generation . energies of lysine dipeptide conformations estimated 
 by three generations 
 of the v - ff14 force field . each 
 generation s fitting set contained all of the initial fitting 
 set data plus conformations created by all previous generations . the 
 torsion fourier series terms 
 describing dihedral interactions of atom types c  n  cx c8 
 ( backbone carbonyl carbon of any residue n - terminal to lysine or arginine , 
 backbone nitrogen , c , and c of lysine or arginine ) and 
 x  c8nl x ( generic torsion affecting amino terminal 
 hydrogens ) evolve rapidly over three generations as sampling of the 
 backbone  angles and amino headgroup orientations becomes more 
 complete . the set of torsion 
 parameters optimized in the third generation 
 of training was chosen for pairing with q and distributed as ff14ipq in ambertools14 ( see http://ambermd.org ) . blocked dipeptides , the simplest systems exhibiting protein - like 
 backbone and side - chain dynamics , were studied extensively to characterize 
 the effects of slight alterations in the charge set , q , which differentiate ff14ipq , our force field for simulations 
 of proteins in water , from v - ff14 , a force field used as a parameter 
 fitting apparatus that is otherwise appropriate only for simulations 
 of proteins in vacuum . furthermore , while we found that torsion parameter 
 refinement could eliminate catastrophic traps in the potential energy 
 surfaces of some amino acids , the approach made modest improvements 
 in the behavior of every other amino acid as well . we computed multiple 
 two - dimensional potentials of mean force ( pmfs ) for alanine , glycine , 
 and serine dipeptides to assess backbone propensities with either 
 charge set and with alanine dipeptide to examine the impact of torsion 
 refinement . as shown in figure 4 , the 
 pmfs of simple , nonpolar amino acids are nearly identical for each 
 charge set , particularly in the populated regions of the ramachandran 
 plot . by inspection , the polarization of charges slightly increases 
 the alanine propensity toward ( right - handed ) -helices , decreases 
 its propensity toward -sheets , and leaves the model s 
 strongest tendency , toward poly proline ii backbone conformations , 
 unchanged . larger differences in all dipeptides appear near ( , ) 
 = ( 0,0 ) , as the backbone n  h and c = o dipoles become 
 favorably aligned to create a more favorable free energy under the 
 ff14ipq model , which strengthens both dipoles considerably . however , 
 the ( , ) = ( 0,0 ) arrangement remains strongly disfavored 
 by steric clashes . the pmfs for alanine and glycine amino acids are 
 also very close to those of the ff99sb force field . the lowest row of panels in figure 4 suggests 
 that solvent effects that make the side - chain charges more polar also 
 facilitate transitions between poly proline ii and left - handed -helicies , 
 but the relative weights of each major conformation are unchanged . while it is not certain what differences in protein folding or dynamics 
 a g of 0.25 kcal / mol in a low - energy 
 region of ( , ) space could lead to , these plots suggest 
 that q changes protein folding transition 
 states more than equilibria . potentials of mean force for blocked alanine , 
 glycine , and serine 
 dipeptides in water . the color scale for the preliminary versions 
 of ff14ipq ( leftmost panels ) and v - ff14 ( middle panels ) measures g , the energy difference between any point in ( , ) 
 space and the minimum free energy attainable in each model at 298 
 k. differences between models are shown on the rightmost panels in 
 a separate color scheme . the importance of fitting torsion parameters in the context 
 of 
 an appropriate charge model is shown in figure 5 . in contrast to the changes introduced by fitting torsion parameters 
 for v - ff14 to a gas - phase qm pes and then transferring them to work 
 in ff14ipq , fitting torsion parameters in the context of ipolq charges , 
 retroactively forcing this charge model to mimic a gas - phase qm pes , 
 makes more substantial changes to the alanine pmf . the difference 
 between this pmf and the ff14ipq pmf is much more frustrated than 
 the difference between the ff14ipq and v - ff14 pmfs . the strong gradients 
 evident in the difference map drive the -helical minimum approximately 
 15 to the south if the torsion parameters are fitted in the 
 context of q and gas - phase quantum 
 data . when working against an already polarized 
 electrostatic model , the torsion parameters were also less effective 
 at reproducing the gas - phase qm pes ( data not shown ) . difference plot of the alanine dipeptide pmf 
 with torsion parameters 
 derived for q rather than q. the color scheme is similar to difference 
 plots in figure 4 : here , solid red implies 
 that a hypothetical ( and incorrect ) model fitting gas - phase quantum 
 data in the context of charges appropriate to the solution - phase estimates 
 a point in / space more than 1 kcal / mol more favorably 
 than a properly tuned model ; solid blue would imply that the incorrect 
 model disfavors the conformation . refining the torsion parameters has its own effects on the 
 alanine 
 dipeptide pmf , as shown in figure 6 . even though 
 the quantum mechanics target is unchanged and the initial set of torsion 
 parameters had no serious problems optimizing the geometries of alanine 
 ( tetra)peptide , the third - generation model produces noticeable differences 
 in the pmf . notably , what was a saddle point between -helical 
 and poly proline ii conformations in the first - generation ff14ipq 
 becomes a local minimum in the third . transitions between poly proline 
 ii , -helical , and l--helical conformations 
 also tranverse lower barriers according to the third generation of 
 torsion parameters , and in some of the high - energy regions of the 
 pmf , the differences between the first and final generations are as 
 great as those shown in figure 5 . while the 
 magnitudes of the differences are comparable to those in figure 5 , however , the locations of all the major backbone 
 conformational minima remain in the same places , and the most significant 
 differences reside in high - energy regions of the pmf , which stands 
 in contrast to the consequences of fitting torsion parameters with 
 the wrong charge set . as before , the influence of these changes on 
 protein folding is not discernible from the pmf alone , but the differences 
 again appear mostly in highly strained configurations , suggesting 
 that the refined torsion parameters depict more frequent transitions 
 between major backbone conformations . potential of mean force for blocked alanine 
 dipeptide in the release 
 version of ff14ipq and comparison to the initial model . both plots 
 refer to combinations of condensed - phase charges with torsion parameters 
 fitted to reproduce gas - phase quantum data . peformance of the release version 
 of ff14ipq , specifically the third generation of torsion parameters 
 paired with the condensed - phase appropriate q , was evaluated on ala(5 ) , -hairpins chignolin and 
 protein g c - terminal fragment , the -helical miniprotein k19 , 
 two variants of the miniprotein trp cage , and globular proteins gb3 
 and lysozyme . sources of each protein structure , as well as sequences 
 of the miniproteins , are given in table 7 . each of the miniproteins simulated in this study was selected to evaluate 
 ff14ipq s performance on a particular secondary structure element . simulations combined third - generation 
 torsion parameters with q. protein sequence ; blocking groups 
 are indicated by ace- and -nme . the ala(5 ) system has recently become a standard diagnostic of 
 a force field s ability to balance three major backbone configurations 
 and reproduce nmr j - coupling results . best and colleagues 
 made a comprehensive assessment of modern force fields with respect to nmr data from graf and co - workers , calculating mean  values 
 for each model s reproduction of 11 order parameters . the karplus 
 relations used to calculate order parameters from the md simulations 
 are sensitive to their own coefficients , but best and colleagues took 
 three different sets of karplus coefficients and posited that a  value of 2.25 or less under all three karplus relations indicated 
 a high - quality force field . the results for ff14ipq using the same 
 karplus coefficients and quadruplicate 375 ns simulations of the unblocked 
 peptide with protonated c - terminus are shown in table 8 . the details of the simulated system were intended to match 
 the acidic conditions of the nmr experiments ; new charges were derived 
 specifically for the protonated c - terminal alanine and are given in 
 the supporting information . the mean  values obtained with the original karplus coefficients used 
 by graf and two sets of dft - based karplus coefficients from case and 
 colleagues were 1.3  0.0 , 2.6  
 0.1 , and 1.5  0.0 , respectively . ( the error bars are standard 
 deviations of the  from each of the four 375 ns 
 simulations . ) the major drivers of the scores are disagreement with 
 the j(c , c ) coupling for the second alanine 
 residue and the j(hn , c ) coupling for the third residue . while ff14ipq does 
 not meet best s  standard by one of the 
 dft results , it does score very well by the other two , and the j(hn , c ) coupling 
 is known to be difficult for the karplus relation itself . by best s definitions of each secondary 
 structure element , ff14ipq models the central residue of ala(5 ) in 
 poly proline ii , -helical , and -sheet conformations 
 for 56 , 18 , and 14% of the pooled trajectories , respectively . dft - based karplus coefficients from 
 case . ff14ipq also stabilizes larger -sheet structures , 
 as indicated 
 by the plots in figure 7 . the gb1 hairpin from 
 protein g appears to be challenging for other fixed - charge 
 force fields to stabilize ( emilio gallicchio , personal communication ) 
 but maintains its secondary structure throughout the 250 ns simulations 
 whether simulated with blocking groups or without . while it is reassuring 
 that ff14ipq stabilizes a -sheet structure , the hairpin should 
 be only approximately 30% folded at 298 k. one element that appears 
 to stabilize the system considerably is the ionic interaction between 
 the termini : it remains intact for virtually the entire simulation 
 with a mean and modal distance between the termini of 3.5  0.7 
  . we investigated the stability of the blocked peptide with 
 quadruplicate 600 ns runs and found no significant unfolding of the 
 antiparallel -sheet arrangement ( data not shown ; the 250 ns 
 simulation reported in figure 7 is representative ) . given that the melting curves for this peptide and other -hairpins are well - established , it may 
 be possible to chip away at the stability of these structures in ff14ipq , 
 and this is probably easier than to try and obtain the right structure 
 from a model that can not stabilize the peptide at all . backbone stability of 
 the -hairpin from protein g over 250 
 ns . blocked and unblocked forms of the peptide were simulated , and 
 backbone rmsd is plotted for both variants . the dssp chart below the 
 rmsd plots refers to the blocked peptide system and indicates that 
 the antiparallel -sheet is maintained throughout the simulation . the chignolin system could be 
 simulated economically on equally 
 long time scales and appears to fluctuate between two major backbone 
 states . neither of them is very far from the structure obtained after 
 equilibration ; figure 8 shows backbone rmsd 
 relative to the first structure in the nmr ensemble , but if it is 
 calculated relative to the first frame of the simulation after restrained 
 equilibration , then the rmsd would fluctuate between 0.7 and 1.4  . in addition to calculating the standard , overall backbone positional 
 root - mean - squared deviation ( rmsd ) relative the first nmr model , we 
 computed the way in which individual residues deviate from the nmr 
 ensemble as a whole . this was accomplished by making optimal alignments 
 of the simulated backbones to each of the 18 chignolin nmr models 
 found in the pdb and then computing the rmsd of each particular residue 
 without further position alignment . the minimum residue rmsd to any 
 of the models in the nmr ensemble was recorded for each frame of the 
 trajectory to indicate how far from any of the plausible structures 
 our simulations might have departed . the lower panel of figure 8 shows histograms of these per - residue backbone 
 position rmsds for each residue over the course of the simulation . this confirms that the backbone adopts two major conformations , with 
 larger fluctuations in the zwitterionic , unblocked terminal residues . the arrangement of the termini follows the gb1 hairpin system : the 
 zwitterionic termini begin the simulation in contact with one another . 
 in chignolin , the termini fluctuate in concert most of the time , although 
 for 7% of the simulation the two ionic groups separate by more than 
 4.5  . the -hairpin fold is maintained throughout the 
 simulation , but the state adopted in the middle of the simulation 
 appears to be slightly outside the nmr structure ensemble . further 
 simulations may be able to fold the peptide , test the pathway against 
 the known mechanism , and suggest whether 
 changes to the backbone  backbone interactions would tighten 
 up the equilibrium structural results . backbone rmsd in 
 the top plot is calculated relative to the first nmr structure ; per - residue 
 backbone rmsd reflects the deviation of each residue s backbone 
 from the closest possible match out of the entire nmr ensemble . over the course of 1 s , residues 112 
 exhibit consistent helicity , while the c - terminus transiently explores 
 alternative coil conformations and sometimes even forms a packed double - helical 
 structure , as shown in figures 9 and 10 . the experimental data suggest that the protein s 
 first 18 residues should consistently maintain an -helix , but the ff14ipq results are otherwise in qualitative 
 agreement with experiment and previous simulations of this system . 
 of particular interest was the degree to which highly solvent - exposed 
 lysine residues in the c - terminal region might make hydrogen bonds 
 with the backbone . jones radii in both the 
 lysine amino headgroup and the backbone oxygen were altered when interacting 
 with water for agreement with hydration free energies but allowed 
 to maintain their original interaction with one another in the interest 
 of keeping repulsion between nearby atoms from disrupting the internal 
 potential energy surface . interactions between lysine and the backbone 
 oxygen may not be well tuned : much has been done to reduce the backbone 
 oxygen s affinity for water , but its affinity for amino groups 
 is only balanced , if at all , by the fact that amino groups 
 lennard  jones radii were reduced to provide greater affinity 
 for water . we calculated the minimum distance between each lysine 
 headgroup on k19 with any backbone oxygen atom throughout the simulation . as shown in figure 11 , lysines 4 and 9 rarely 
 make contact with backbone carbonyl groups , but lysines 14 and 19 
 are much more likely to form hydrogen bonds to the peptide backbone . 
 defining an amino group to be hydrogen - bonded to a carbonyl if the 
 nitrogen and oxygen atoms come within 3.2  of each other , lysines 
 4 , 9 , 14 , and 19 are bonded to the backbone in 4 , 7 , 17 , and 27% of 
 the snapshots , respectively . while the interaction between some groups 
 may be stronger than is realistic , their interactions remain transient . the nitrogen atom types in the lysine headgroup and the charged amino 
 teminus of our -hairpin systems are the same , and the oxygen 
 atom types in carbonyl and carboxylate groups were altered in a similar 
 manner with respect to their precursors in ff99 . the artificial behaviors 
 seen in each system are likely to have a common origin . the time axis 
 applies to both the rmsd plot in the top panel and the dssp plot in 
 the lower panel . the k19 peptide is predominantly -helical , 
 with some instability at the c - terminus . residues 1619 begin 
 to adopt a metastable 310 helical conformation near the middle 
 of the simulation . all conformations 
 have been aligned relative to the stable backbone of residues 214 . radial distributions of lysine head groups 
 and backbone oxygen 
 atoms in the k19 system .. as shown in figure 12 , both trp cage 
 simulations 
 showed very stable backbone configurations over the 500 ns simulations . the baseline rmsd of approximately 1.2  may indicate different 
 preferences of the ipolq charge set compared to the parameters used 
 in nmr refinement , but larger departures from the backbone configuration 
 depicted in the nmr models were merely transient . figure 12 shows histograms of per - residue backbone position 
 rmsds as were calculated for chignolin . most residues display low 
 rmsds under this test , and most importantly the pro - pro - pro sequence 
 near the c - terminus remains stable . some residues , in particular the 
 asp - gly - gly - pro sequence in simulations of neidigh s trp cage 
 ( pdb code 1l2y(45 ) ) , show a weakly bimodal distribution , 
 suggesting that the simulations explore an alternate conformation 
 not seen in the nmr ensemble . when simulating a hyperstable trp cage 
 mutant ( pdb code 1rij(46 ) ) , the asp - gly - gly - pro sequence again 
 departs from the nmr ensemble more significantly than other regions 
 of the protein . the conformations of these residues are good candidates 
 for future analysis and refinement of ff14ipq . backbone positional 
 root - mean - squared deviations ( rmsds ) for trp 
 cage miniprotein simulations . simulations of each of two trp cage 
 proteins are indicated by their respective pdb codes . the top panel 
 shows overall backbone rmsd to the first published nmr model for each 
 system . lower panels show histograms of per - residue backbone rmsd 
 to the closest possible match out of all published nmr models , darkened 
 to indicate increasing occupancy at a particular deviation . even though the ipolq charge model 
 is directed toward solvated 
 molecules , ff14ipq is expected to have application to globular proteins . most residues on proteins of even a few hundred residues have some 
 degree of solvent exposure , and the majority of buried residues will 
 be nonpolar and therefore little different when represented by ipolq 
 as opposed to other charge models . simulations 
 of both gb3 and lysozyme showed that , overall , ff14ipq stabilizes 
 the crystallographic backbone configurations of both proteins ( figure 13 ) . gb3 and 
 lysozyme 
 ( 56 and 129 residues , respectively ) were simulated with ff14ipq in 
 baths of tip4p - ew water for 1 s . multiple simulations of gb3 
 are shown , each performed with a different generation of the ff14ipq 
 torsion parameters . see figure 3 for the accuracy 
 of each generation in predicting energetics of lysine dipeptide . because we began simulations of 
 gb3 before realizing the importance 
 of iterative torsion parameter refinement , the system provides an 
 indication of the symptoms of underfitted parameters when simulating 
 a complex biomolecule . the first generation of torsion parameters 
 caused one of the lysine - containing loops ( residues 916 ) to 
 take on alternative conformations that drove the overal rmsd significantly 
 higher . application of the second - generation torsion parameters , which 
 fixed an artificial minimum in the lysine backbone  angle , 
 greatly diminished these excursions , as shown in figure 14 . however , another very short loop of the protein , 
 residues 3941 connecting -sheet to -helical 
 structures , makes excursions from the x - ray backbone structure in 
 any generation . it may be significant that this loop contains an aspartate 
 residue : in the torsion fitting results ( table 4 ) , aspartate presented one of the most difficult potential energy 
 surfaces for our torsion parameters to capture . while ff14ipq seems 
 to provide a better fit than would be possible with ff99 , the greatly 
 expanded parameter space that is almost certainly the basis for this 
 improvement may have allowed some overfitting that we have not yet 
 been able to eliminate with our iterative scheme . alternatively , aspartate 
 polarization in solution may not be well captured by our approximations . ( a third possibility remains , whereby in solution this loop really 
 does take on conformations not seen in the crystal lattice . ) we intend 
 that future releases of ff14ipq will build on all of the existing 
 quantum data with further refinements to improve the description of 
 larger protein systems . per - residue backbone rmsd for gb3 , with three 
 generations of the 
 ff14ipq model . per - residue rmsd was calculated in the same manner 
 as was done for trp cage , but the reference ensemble comprised only 
 the one x - ray structure . loops that make significant departures from 
 the x - ray structure in solution - phase simulations are emphasized on 
 the x axis . the simulations performed to date with ff14ipq suggest that 
 the 
 force field is ready to guide much longer time scale experiments with 
 proteins and solvated peptides . while some results indicate possible 
 weaknesses in ff14ipq , it is straightforward to develop hypotheses 
 for their origins in the parameter development . it is encouraging 
 to get this level of performance from a new force field and to see 
 the model improve over generations of new fitting data . before completing our molecular mechanics 
 model based on the ipolq 
 charge set , we needed a set of charges compatible with the gas - phase 
 quantum calculations available to guide our torsion parameter optimization . this , in turn , required us to revisit the way in which we derived 
 charges , decomposing the ipolq charges into foundational values describing 
 the electrostatics of solutes in vacuo and perturbations describing 
 the manner in which atoms polarize in aqueous solution . because 
 of the iterative nature of the ipolq procedure , we must consider whether 
 this latest change will require additional iterations of molecular 
 simulations with fixed solutes followed by charge fitting . first , although the 
 polarities of certain bonds have changed somewhat , the individual 
 charges of polar side - chain atoms and overall dipoles of the molecules 
 did not change very much . restraints on q seem 
 to have the strongest effects on buried atoms or groups : increasing vgp reveals indeterminacy in the model fitting 
 but is not likely to affect hydration free energies . indeed , we chose 
 this approach to exploit the fact that there exist many partial charge 
 models that reproduce a molecule s electrostatic potential 
 with nearly the same accuracy : increasing vgp chooses two models that are most similar in terms of the mean - squared 
 value of q . second , the objective of the ipolq 
 procedure is a solvent reaction field potential ( srfp ) that is consistent 
 with the way a solute polarizes ; the result is a set of charges that 
 is consistent with the srfp . restraints on q can make minor changes to the result , but they should not change 
 the objective . it would therefore be appropriate to apply the original 
 ipolq procedure until convergence and then to make other stipulations 
 about how the model s partial charges reproduce the fitting 
 data . one surprising result of this analysis is that the electrostatic 
 potentials of solutes in vacuo are harder to fit than their potentials 
 in the condensed phase when the wave function has been computed in 
 the presence of a polarizing charge density . this occurs in spite 
 of the fact that the charges for polar groups , including the backbone , 
 become significantly larger in the condensed phase to represent stronger 
 electrostatic fields . while an analysis of why this occurs is beyond 
 the scope of this work , it is a worthy subject for future studies 
 and may inform the design of polarizable charge models . in another 
 study , extension 
 of the standard resp protocol , which 
 is similar in spirit to our extended ipolq scheme if one takes the 
 vacuum charges to be dresp s baseline charge model . dresp employs 
 variable restraint stiffnesses on polar and nonpolar atoms to accommodate 
 different polarizabilities , something that our arbitrary constant vgp does not support . we had thought that , with a complete charge parameter scheme and 
 no obvious reasons to make further changes to bonded parameters or 
 lennard  jones terms , producing viable torsion parameters would 
 be the simplest step in developing ff14ipq . beyond the aforementioned 
 issues with double - counting solvent contributions , there were many 
 hurdles in deriving a transferable torsion parameter set , and this 
 aspect of the model development required as much effort as the charge 
 derivation . we hope that our experience and the programs we designed 
 to meet each challenge will benefit future ab initio force field development . we found inclusion of new atom types and the new fourier series 
 terms they generated to be very helpful in fitting the vacuum - phase 
 potential energy data . the new types introduced to accommodate the 
 ipolq charge set did not introduce many new parameters , as the atom 
 types they replaced were already found in unique bonded arrangements . however , introducing atom types specific to individual residues introduced 
 many new terms that could be fitted to the unique potential energy 
 surface of the residue and relieve the burden that would have been 
 placed on generic terms . in this manner , including a new atom type 
 for c made improvements in the 
 rmse for nearly all dipeptide systems , as the  and  
 torsions of the uniquely achiral center of glycine could then be decoupled 
 from those of all other , chiral amino acids . while we had wanted to 
 include separate atom types for the c atoms of positively and negatively charged residues to match the 
 unique charges derived for all backbone atoms in these residues , we 
 were not confident in our ability to generate enough fitting data 
 to sample the  and  torsion space for these individual 
 residues . in particular , sampling  and  for consecutive 
 arginine or lysine residues would require mp2/cc - pvtz calculations 
 on large tripeptide systems at very high computational cost . expansion 
 of the data set in this manner will likely be part of the next version 
 of ff14ipq . our study also highlights the perils of adding new 
 parameters without 
 adequate sampling . while we believe it is correct to assign unique 
 parameters to unique chemical environments , if the degrees of freedom 
 governed by each parameter might take on values in simulations that 
 are not present in the training set , then catastrophic excursions 
 into artificial minima can occur . data collected with the first - generation 
 ff14ipq showed unnatural backbone  angles adopted by lysine 
 residues in the gb3 protein , which disappeared in subsequent generations 
 trained with complete sampling of this rotatable bond . while such 
 cases might be apparent to an especially careful investigator , we 
 programmed the mdgx fitting routines to produce a great deal of information 
 on the sampling of each parameter , and we did not connect higher backbone 
 fluctuations in the gb3 protein loops to undersampling in particular 
 parameters until significant improvements were seen in the second 
 generation of the force field and a complete breakdown of the molecular 
 mechanics energy was determined for each of the hundreds of lysine 
 structures in each training set . perhaps ore important is our finding 
 that the descriptions of residues that do not exhibit such catastrophic 
 artificial minima can still be improved by the simple approach of 
 reintroducing the products of molecular simulations back into the 
 model s training data . whatever artificial minima were removed 
 from the other residues by the third generation of ff14ipq must have 
 been shallow and probably involved combinations of several parameters , 
 implying a very large conformational space to sample in order to obtain 
 transferable parameters on the first attempt . if allowed to 
 drive simulations , then any molecular model should 
 be expected to overpopulate regions of the conformational space in 
 which it scores too favorably . however , because all of the torsion 
 parameters discussed in this work are based on cosines , an unrefined 
 model s propensity to overstate the favorable nature of some 
 conformations may conversely create a tendency to overestimate barrier 
 heights . if this is true , then it may not be a coincidence that the 
 third generation of ff14ipq models transition states in alanine dipeptide 
 with lower barriers than the first generation . in our charge parameter 
 development , we were forced to alter 
 the radii of certain polar atoms , mostly by making them larger , to 
 adjust hydration free energies of amino acid side - chain analogues 
 into agreement with experiment . in this work , we were forced to break 
 the combining rule when integrating these atoms into a complete force 
 field because the larger radii of nearby atoms would clash and make 
 the internal potential energy surfaces of polar amino acids very hard 
 to model . we have covered the possible implications of this decision 
 in the k19 system , where lysine head groups appear to interact too 
 strongly with the backbone . the special combining rules for these 
 polar atoms would also affect the stability of backbone backbone 
 contacts , particularly -sheets and -helicies , and are 
 likely the source of the overstabilization reported for the gb1 -hairpins . goetz and colleagues reported strong performance 
 for the ipolq charge set when studying aggregation of ionic amino 
 acids , but in that system the correct 
 result was for the amino acids to aggregate . while a prototype version 
 of ff14ipq was able to form these aggregates with the correct radial 
 distributions , their study did not test whether the contacts formed 
 might , in fact , be too strong . more recently , debiec , gronenborn , 
 and chong reported that the ipolq charge set produced the correct 
 association constants for charged amino acid side chains when the 
 full amino acids were considered : positively charged side chains showed 
 affinity for both negatively charged side chains and the backbone 
 carbonyl group , and the competition drove the equilibrium of side 
 chains contacting one another in the right direction . however , they did not rule out the possibility that the 
 side chain to backbone contacts were themselves overstated , as our 
 results on k19 lead us to believe . debiec and co - workers also did 
 not have available to them the special combining rules we adopted 
 for ff14ipq . when these rules are applied , they found that ff14ipq 
 overstabilized salt bridge formation between charged amino acid side - chain 
 analogues by about 1 kcal / mol ( personal communication ) . we are 
 collaborating with debiec and others to find a pair - specific 
 lennard  jones matrix that will not create excessive 1:4 strain 
 but will properly balance the strength of polar interactions . we also 
 intend to add new residues to the palette , including nucleic acids 
 and phospholipids . the folding and aggregation of these residues is 
 dominated by electrostatic interactions , which will demand the proper 
 balance between charged groups interacting with one another and also 
 with water . for proteins jones 
 interactions will still contain special combining rules but that the 
 departures from a standard rule will be smaller than they are in ff14ipq . some of the lennard  jones interactions in ipolq itself can 
 be traced to limitations of the nuclear - centered charge model . for 
 instance , the nuclear - centered charges can not portray the hydroxyl 
 group with as strong a polarity as it should have without making larger 
 errors elsewhere in the electrostatic field . the hydroxyl oxygen  radius was reduced to correct alcohol 
 hydration free energies , but this probably would not have been necessary 
 with a virtual site scheme that could model the correct polarity . 
 for lipids and nucleic acids , any new departures from a standard lennard jones 
 combining rule may hold clues to the inherent limitations of a nonpolarizable , 
 nuclear - centered charge model .
we 
 present the ff14ipq force field , implementing the previously 
 published ipolq charge set for simulations of complete proteins . 
 minor 
 modifications to the charge derivation scheme and van der waals interactions 
 between polar atoms are introduced . 
 torsion parameters are developed 
 through a generational learning approach , based on gas - phase mp2/cc - pvtz 
 single - point energies computed of structures optimized by the force 
 field itself rather than the quantum benchmark . in this manner , we 
 sacrifice information about the true quantum minima in order to ensure 
 that the force field maintains optimal agreement with the mp2/cc - pvtz 
 benchmark for the ensembles it will actually produce in simulations . 
 a means of making the gas - phase torsion parameters compatible with 
 solution - phase ipolq charges is presented . 
 the ff14ipq model is an 
 alternative to ff99sb and other amber force fields for protein simulations 
 in programs that accommodate pair - specific lennard 
 jones combining 
 rules . 
 the force field gives strong performance on -helical 
 and -sheet oligopeptides as well as globular proteins over 
 microsecond time scale simulations , although it has not yet been tested 
 in conjunction with lipid and nucleic acid models . 
 we show how our 
 choices in parameter development influence the resulting force field 
 and how other choices that may have appeared reasonable would actually 
 have led to poorer results . 
 the tools we developed may also aid in 
 the development of future fixed - charge and even polarizable biomolecular 
 force fields .