molecular 
 mechanics ( mm ) force fields for biomolecular simulations 
 are undergoing increased scrutiny as access to improved algorithms 
 and hardware allow for validation of , for example , protein dynamics , protein  ligand binding free energies , and liquid properties on an unprecedented scale . their success likely stems from a functional form 
 that is based on physical laws , and very meticulous parametrization 
 by many research groups on increasingly large databases of ab initio and experimental data . force fields that model the nonbonded 
 component of the total energy by coulomb interactions between fixed 
 atom - centered point charges and lennard - jones ( lj ) interactions are , 
 by far , the most widely used functional form today , and form the basis 
 of the opls , amber , charmm , gromos , and many other force fields . the force field parameters have undergone extensive fitting against ab initio binding energies of complexes and properties such 
 as free energies of hydration and liquid densities and heats of vaporization 
 for small molecules . furthermore , it has been shown that , when small molecules are assembled 
 into macromolecules such as proteins , the same parameters perform 
 reasonably well in the prediction of , for example , temperature - dependent 
 structural changes and protein ligand 
 interactions . with the wealth of validation 
 tests available in the literature , a user can be confident that dynamical 
 simulations of simple proteins with a state - of - the - art transferable 
 force field will provide an ensemble of structures that is representative 
 of the experimental reality . given the widespread use of transferable 
 force fields , it is important 
 to consider the long - term routes to their improvement . as computational 
 chemistry advances , one can continue to fit parameters to more - accurate 
 quantum mechanical electrostatic potentials and energy surfaces within 
 the confines of a fixed functional form . one can even consider amending 
 the functional form , and important steps in this direction are being 
 taken in the development of polarizable force fields . each time that the method for the derivation 
 of point charges is updated , for example , the library of parameters 
 that describe the lj interactions must be refit to experimental data . 
 for a typical force field , this library may consist of many hundreds 
 of parameters , although automated parameter optimization programs , 
 such as forcebalance , promise to bring 
 down the cost of these fitting efforts in the future . yet , even 
 if force field functional forms and parameters are found 
 that perfectly reproduce the empirical and ab initio properties of small molecules , the question of transferability remains . it is well - established that adding electron - withdrawing or electron - donating 
 functional groups to a molecule can significantly alter its chemistry 
 through polarization effects . for example , trends in hydrogen bond 
 strength between a para - substituted phenol and a water molecule were 
 investigated using two force fields : one nonpolarizable ( opls - aa ) 
 and one in which the charge distribution of the phenyl ring responds 
 to the presence of the substituents ( opls / cm1a ) . the opls / cm1a variant reproduced ab initio binding energies more accurately than the nonpolarizable version 
 of the force field , which was attributed to the ability of the charge 
 derivation scheme to respond to the chemical environment . in this regard , it is common in small - molecule force field parametrization 
 for the user to assign atom - centered charges based on quantum mechanical 
 calculations . care must be taken that the obtained charges are compatible 
 with both the lj parameters and the force field used to describe the 
 surrounding environment , which could contain water , protein , lipids , 
 and so forth . finally , the effects of environmental polarization are 
 not limited only to the force field charges . it has been shown that 
 the dispersion ( or van der waals ( vdw ) ) coefficients also are sensitively 
 dependent on both the local environment of the atom and long - ranged electrodynamic screening . hence , the accuracy of the force field may also be improved if the 
 lj parameters , as well as the charges , were able to respond to the 
 chemical environment . interestingly , a fundamentally different 
 approach to force field 
 parametrization was used recently by grimme in the development of 
 the quantum mechanically derived force field ( qmdff ) . here , instead of relying on a small number of research 
 groups to fit parameters and package them into a force field with 
 the implicit assumption that the parameters are transferable and consistent 
 with small molecule force fields , the user derives all of the nonbonded 
 parameters that are specific to their system in an automated fashion 
 from quantum mechanical calculations . such a scheme 
 might be termed an environment - specific , rather than transferable , 
 force field . more generally , we envisage future environment - specific 
 force fields 
 being derived in the same manner as today s small molecule 
 force fields , except that it would be possible to also ( i ) derive 
 lj parameters and ( ii ) derive parameters for both small molecules 
 and macromolecules such as proteins . environment - specific force fields 
 have the potential to be more accurate than their transferable counterparts , 
 since they would be derived explicitly for the system under study . however , of course , it would be necessary to extensively benchmark 
 their accuracy against a wide range of biomolecular data , such as 
 nuclear magnetic resonance ( nmr ) structural data , protein  ligand binding free energies , and so forth . a second advantage is that there 
 would be no need for the user to mix force fields . protein and small - molecule 
 force fields would be computed by the user simultaneously , thus ensuring 
 their consistency . furthermore , charges and lj parameters would be 
 computed from the same quantum mechanical data and 
 so the parameters would be inherently consistent . if the user wishes 
 to use a more advanced quantum chemistry method to perform the quantum 
 mechanics ( qm ) calculation , then both the charges and vdw parameters 
 would respond to the change in qm electron density by construction . thus , the user would be able to experiment with more advanced chemistry 
 or with alternative functional forms of the force field without relying 
 on another research group to fit the library of force field parameters 
 to experimental data . a disadvantage of such a scheme is that the 
 user would require sufficient computational resources to perform a 
 qm calculation on their system of interest . however , with ongoing 
 advances in large - scale density functional theory ( ls - dft ) software and increasing access to improvements in hardware , 
 qm calculations in excess of 1000 atoms are becoming routine . as an example , ufimtsev et al . have computed the electronic structure 
 of a 2634-atom protein water complex at the b3lyp/6 - 31 g level 
 using the terachem software and a single desktop workstation with 
 eight gpus within <3 h. in recent 
 years , we have started to develop the techniques necessary 
 to derive environment - specific force fields for large systems . the 
 force field is based on the concept of atoms - in - molecule ( aim ) electron 
 density partitioning , whereby the total qm electron density of the 
 system under study is computed and then partitioned into approximately 
 spherical atomic basins . although there is no unique way to perform 
 this partitioning , it has been shown , by careful choice of the weighting 
 functions , that aim densities can be optimized to simultaneously produce 
 chemically meaningful partial atomic charges , yield an efficiently 
 converging multipole expansion of the electrostatic potential , and 
 be insensitive to small conformational changes or to buried atoms . furthermore , by combining aim partitioning with ls - dft calculations , 
 atomic partial charges may be computed for systems that are comprised 
 of many thousands of atoms , including proteins . using this combination of aim electron density partitioning 
 and 
 ls - dft , we have previously derived partial charges for many small 
 proteins and shown that they reproduce very well the electrostatics 
 of the underlying dft calculation , as well as experimental nmr order 
 parameters when used in mm molecular dynamics simulations . however , a more sensitive test of intermolecular 
 interactions is to compare hydration free energies and liquid properties 
 of small organic molecules with the experiment . these properties are 
 commonly used to validate new force fields since there 
 exists a wealth of experimental data , which is directly comparable 
 with simulation . furthermore , any difference between the experiment 
 and the theory may be unambiguously attributed to errors in the force 
 field , rather than finite sampling errors and so forth , which may 
 affect more - complex validation datasets . initial tests revealed that 
 the aim charges are not as accurate as we would like when coupled 
 with the lj parameters of commonly used force fields . for example , 
 the combination of aim charges with opls lj parameters gave a computed 
 free energy of hydration of 12.8 kcal / mol for acetamide , which 
 may be compared with the experimental result of 9.7 kcal / mol . one solution might be to reparameterize libraries of lj parameters 
 for compatibility with aim charges in the knowledge that this procedure 
 would need to be repeated periodically whenever an adjustment to the 
 charge model is made . instead , we show here that the lj parameters 
 may be obtained , in conjunction with the atomic partial charges , directly 
 from aim electron densities . the lj parameters are thereby specific 
 to the system under study and respond automatically to the chemical 
 environment of the atom . similar to that observed with the charges , 
 they may be derived for both small and large molecules , including 
 proteins . the lj parameters are computed from the same qm electron 
 density as the partial charges and , hence , the entire parameter set 
 is consistent and there is no additional computational cost to the 
 user . perhaps most importantly , there are a very small number of fitting 
 parameters in the model . as we will describe , our nonbonded force 
 field contains only seven fitting parameters , which 
 is sufficient to describe the chemistry of a wide range of small organic 
 molecules and proteins . in contrast , typical transferable force fields 
 require many tens or hundreds of fitting parameters to model nonbonded 
 interactions in the same systems . in what follows , the possible applications of biomolecular 
 force fields are extremely wide and varied , ranging from protein ligand 
 binding free energies for drug discovery , to pka prediction , to dynamics for studying protein function , to 
 hybrid qm / mm simulations for enzymology . as such , it is not possible 
 to validate a new force field approach in all of these contexts in 
 a single study . instead , following the example of many other development 
 efforts , we begin by validating our proposed nonbonded force 
 field against experimental free energies of hydration , liquid densities , 
 and heats of vaporization of small organic molecules . as we will show , 
 aim techniques are extremely competitive with existing transferable 
 force fields for this dataset , which will motivate future validation 
 tests beyond those presented in this study . furthermore , since our 
 goal is to eventually derive force fields for large systems , such 
 as protein  ligand complexes , we present here a proof - of - principle 
 derivation of a force field for an 1600 atom model of the 
 l99a mutant of t4 lysozyme . we compare and contrast the derived nonbonded 
 parameters with a standard transferable force field and compute the 
 relative binding free energy of indole and benzofuran to the protein . the quantum 
 mechanical interaction 
 energy is typically decomposed into four components : ( i ) electrostatics , 
 ( ii ) induction , ( iii ) dispersion , and ( iv ) exchange - repulsion . here , we shall discuss how each of these terms , in turn , is incorporated 
 into an effective force field , which , for two atoms i and j separated by a distance rij , is typically written in the form of eq 1:1 the first term of 
 the nonbonded expression describes the classical electrostatic interactions 
 between the charge distributions of the two atoms . it is usually described 
 by fixed atom - centered point charges ( qi ) , although off - site charges to model , for example , lone - pairs and -holes are relatively common . typically , the charges are fit to reproduce ab initio electrostatic properties of small molecules , either by fitting to 
 the electrostatic potential or by distributed multipole analysis . however , neither of these methods has been applied 
 to very large systems , such as proteins , due to the difficulties of 
 fitting charges to buried atoms or the expense of the qm calculation . 
 here , we instead use an aim approach whereby the total qm electron 
 density ( n(r ) ) is partitioned into overlapping 
 atomic densities ( ni(r)):2 the 
 atomic partial charges are then 
 computed by integrating the atomic electron densities over all space:3where ni is the 
 number of electrons assigned to atom i and zi is its effective nuclear charge . similarly , 
 higher - order atomic multipoles may be computed as first - order , second - order ,  ... 
 moments of the atomic electron densities . for example , in hirshfeld partitioning , they are set equal to neutral 
 gas - phase atomic densities , although 
 these weights are known to result in atomic populations that are too 
 close to zero . instead , we employ density - derived 
 electrostatic and chemical ( ddec ) electron density partitioning . the form of the ddec weighting function is described in detail elsewhere , but in brief the atomic weights are simultaneously optimized to 
 resemble the spherical average of ni(r ) and the density of a reference ion of the same element 
 with the same atomic population ni . in 
 this way , the assigned atomic densities yield a rapidly converging 
 multipole expansion of the qm electrostatic potential and the computed 
 populations are chemically reasonable . the induction term results 
 from the distortion of a molecule s electron density in response 
 to its environment . for example , moving a molecule from the gas phase 
 into water enhances the molecule s dipole moment . in fixed - charge 
 mm force fields for condensed - phase simulations , which are the focus 
 of this paper , induction this is generally achieved either by using an 
 empirical bond - order correction , a scaling 
 factor , or by computing 
 the charges using an artificially polarized quantum chemistry method 
 or an implicit solvent model . here , the 
 latter approach is used whereby the qm electron density is computed 
 via direct solution of the inhomogeneous poisson equation in a medium 
 with a dielectric constant . to determine 
 which value of  should be used , let us consider the case of 
 the transfer of a molecule from the gas phase to water . use of  
 = 1 would neglect the attractive inductive interactions between the 
 polarized charge distribution of the molecule and water and lead to 
 computed free energies of hydration that are too positive . on the 
 other hand , the use of  = 80 ( to model the dielectric constant 
 of water ) neglects the energetic cost of distorting the electronic 
 wave function , which has been shown to amount to several kcal / mol , and would lead to computed free energies of 
 hydration that are too negative . in fact , the magnitudes of both of 
 these energetic contributions may be estimated using the implicit 
 solvent model , and it can be shown that 
 they approximately counterbalance each other for a wide range of molecules 
 if the electron density is computed in a dielectric medium of  
 = 4 ( section s2 in the supporting information ) . this is also the dielectric constant that is frequently used to 
 model protein interiors , and , as we will 
 show , produces a net polarization effect that is consistent with that 
 recommended by karamertzanis et al . and 
 the developers of the latest generation of amber force fields . therefore , a background dielectric constant value of  = 4 
 is used to compute all condensed - phase qm electron densities in this 
 study . dispersion ( or vdw ) 
 interactions are typically modeled in mm force fields by the attractive bijrij term in eq 1 . the bij parameters are atom - dependent but are almost always determined 
 empirically and are assigned to new molecules from outside the fitting 
 set , using a library of atom types . despite their empirical treatment 
 in mm force fields , dispersion interactions are inherently a quantum 
 mechanical phenomenon , which arise even in nonpolar molecules , because 
 of interactions between spontaneous electron density fluctuations . standard dft exchange - correlation functionals fail to capture the 
 correct long - ranged behavior of the dispersion interaction , but much 
 progress has been made in recent years in computing the strength of 
 the dipole dipole pairwise dispersion interactions ( that is , 
 the bij parameters ) from the ground - state 
 electron density of a molecule . in this paper , the 
 tkatchenko  scheffler ( ts ) scheme is employed to compute bij coefficients by rescaling accurate free atom 
 reference data by the effective volume of the atom in the molecule . specifically , the ratio of the aim volume to 
 the volume of the atom in free space is computed using the ddec atomic 
 densities from eq 2:4where r is the distance from 
 the nucleus of atom i , and nifree(r ) is the electron density of the free atom i in vacuum . vifree was computed for 
 each of the seven elements used in the current study , using the mp4(sdq)/aug - cc - pvqz 
 method in gaussian 09 , and the chargemol 
 code and the data are provided in table s1 in the supporting information . the environment - specific 
 atomic dispersion coefficient is then computed by5 the 
 bifree parameters are obtained from 
 time - dependent density functional theory ( tddft ) calculations from 
 the literature and are tabulated in table s1 . however , 
 it has been shown that computed dispersion coefficients of ionic systems 
 are substantially improved by allowing the weights in eq 2 to update self - consistently with 
 the aim densities . hence , in what follows , we use 
 ddec electron density partitioning to compute both the atomic charges 
 and volumes . finally , in the opls force field , heteronuclear bij parameters are computed via a geometric combining 
 rule , and the same approach is adopted here:6 figure 1 demonstrates 
 the potential utility of this approach for classical force field design . 
 for benzene , the ddec aim population and dispersion coefficient are 
 essentially identical to those used in the opls / cm5 force field , which 
 has been extensively parametrized to reproduce the liquid properties 
 of the molecule . for a druglike molecule , with 
 more - complex bonding environments , both methods agree that the n atom 
 in the pyridine ring ( blue ) has a high electron population , and therefore 
 is expected to be more polarizable and is assigned a high dispersion 
 coefficient . in contrast , the n atom highlighted in red on the triazole 
 ring has a lower electron population . the aim approach responds to 
 the local environment of the n atom and correspondingly reduces the 
 dispersion coefficient to account for electron depletion . opls , like 
 most force fields , assigns identical dispersion parameters to n atoms 
 in pyridine and triazole rings and is thus expected to overestimate 
 the dispersion coefficient in the latter case . fundamentally , there 
 should be coupling between the electron density on an atom and its 
 capacity for dispersion interactions . comparison between the ddec aim parametrization 
 approach used in 
 this work and the opls / cm5 force field for ( left ) benzene and ( right ) 
 a larger druglike aryl-1,2,3-triazole molecule . exchange - repulsion 
 is primarily intended to account for repulsion between overlapping 
 electron clouds due to the pauli exclusion principle , but it also 
 effectively describes other terms that are dependent on electron density 
 overlap , such as electrostatic penetration energies . it is typically modeled in mm force fields by a repulsive lj term 
 of the form aijrij . this term is problematic because there 
 is no clear method to parametrize the atomic ai values directly from qm calculations . here , one approach to derive the strength of the lj repulsion 
 is hypothesized and justified a posteriori by comparison 
 between computed and empirical condensed - phase properties of small 
 molecules . in the opls force field , in order to ensure that the lj 
 potential has a minimum that is coincident with the van der waals 
 ( vdw ) radii of the atoms ( riaim ) , the repulsive parameters 
 of the lennard - jones expression are related to the dispersion coefficients 
 via7 here , the scaling 
 relationship suggested 
 by tkatchenko and scheffler is used to 
 estimate the vdw radius of the atom in the molecule from the radius 
 of the free atom in a vacuum:8 thus , by substituting eq 8 into eq 7 , environment - specific ai terms are estimated from free atom reference 
 data and the aim - partitioned volume of the atom in the molecule ( eq 4 ) . heteronuclear aij parameters are derived via the geometric 
 combining rule:9 the atomic 
 radius of an atom in a vacuum can not be measured directly , 
 and so the rifree are treated as variable fitting 
 parameters . since computed liquid densities are very sensitive to 
 the ai coefficients used in mm force fields , 
 we have fit the rifree parameters by hand to reproduce 
 empirical liquid densities . the parameters for h , c , n , o , s , f , and 
 cl are provided in table s1 . these seven 
 radii are the only fitting parameters used in the construction of 
 our nonbonded force field . this may be contrasted with typical transferable 
 force fields which require many tens , or even hundreds , of fitting 
 parameters to describe charge and lj interactions in biomolecular 
 systems . to summarize , figure 2 shows the suggested workflow for derivation of nonbonded 
 force field 
 parameters from qm calculations . we use the onetep linear - scaling 
 dft code to perform the ground - state electronic structure calculations , and the ddec atomic population analysis and 
 atomic volume calculations have been implemented here as post - processing 
 routines . apart from the qm calculation , the only 
 input data required is summarized in the green box in the figure ( and table s1 ) . these are the atomic volumes and dispersion 
 coefficients of free atoms in vacuo , and they have 
 been computed using qm calculations as described earlier . the free 
 atom radii ( rifree ) have been fit to reproduce 
 experimental liquid densities . these parameters are optimized for 
 the qm method described below and they do not need to be refit by 
 the force field user . a simple script is used to rescale the free - atom 
 parameters according to eqs 5 , 7 , and 8 to 
 provide environment - specific force field parameters that are unique 
 to each atom in the system . the force field 
 input parameters are tabulated in the green box ( and table s1 in the supporting information ) . in the yellow box , 
 a script computes the lj parameters and exports them into a suitable 
 format . the acetonitrile example shown is designed for boss input , 
 in the standard opls format q ( e ) ,  (  ) , 
  ( kcal / mol ) . for the testing 
 on liquid - state properties , 41 small molecules ( listed in table s4 in the supporting information ) were 
 optimized under vacuum at the mp2/cc - pvtz level using the gaussian 
 09 package . the ground - state electron 
 density was computed using the onetep linear - scaling dft code with the pbe exchange - correlation functional . onetep combines computational efficiency with 
 accuracy equivalent to traditional plane - wave dft codes by optimizing 
 a minimal set of spatially truncated nonorthogonal generalized wannier 
 functions ( ngwfs ) on each atom . one ngwf 
 was used for hydrogen and four for all other elements used in the 
 current study . the ngwfs were expanded in a periodic sinc ( psinc ) 
 basis with an equivalent plane - wave cutoff energy of 1020 ev , and 
 were localized in real space with 10 bohr radii . figure s3 in the 
 supporting information shows good agreement between the gas - phase 
 dipole moments of the set of small molecules computed with onetep 
 as described here and those computed at the mp2/cc - pvtz level in gaussian 
 09 , although the dft approach overestimates 
 the dipole moment of more polar molecules by up to 0.4 d in some cases . to account for induction in an effective manner , dft calculations 
 were performed using a smeared ion representation under open - boundary 
 conditions with a relative dielectric constant of  = 4 . the dielectric cavity is defined by an isosurface of the electronic 
 density under vacuum , with parameters controlling the smoothness and 
 density threshold taken from the literature . partitioning of the polarized ground - state electron density was 
 performed using the ddec scheme in onetep as described in detail elsewhere , with the mixing parameter (  ) set to 0.02 , which appeared 
 in preliminary tests to give reasonable agreement between ddec and 
 opls force field charges for a wide range of functional groups . lj 
 parameters were computed as described in figure 2 . to avoid nonphysical asymmetries 
 in force field simulations , lj parameters for polar h atoms ( that is , h 
 atoms bonded to o , n , or s ) were set to zero . however , this neglects 
 important dispersive interactions in the current format , and so the bi coefficient of the neighboring heavy atom 
 is correspondingly increased according to eq 10 , 10where bx and bx are the old and new dispersion parameters for atom x (= o , n , or 
 s ) , nh is the number of neighboring h 
 atoms , and bh represents their dispersion 
 coefficients . bonded interactions were represented by the opls - aa 
 force field and , for simulations involving water , the rigid tip4p 
 model was used . the environment - specific 
 force field for a 1646-atom cluster extracted 
 from the l99a mutant of t4 lysozyme ( pdb : 185l ) was also 
 derived in an identical manner from a single point qm calculation , 
 which is a routine computation using ls - dft . the cluster 
 included the 105 amino acid residues nearest to the ligand . for computational 
 feasibility , the cluster did not undergo prior structural optimization 
 and the nonbonded parameters were not symmetrized . setup of the protein 
 complexes , liquid simulations , and free - energy calculations were performed 
 using standard procedures with the boss 
 and mcpro codes , as detailed in the supporting information ( section s1 ) . figure 3a compares 
 the computed qm dipole moments 
 in implicit solvent (  = 4 ) with the same quantity under vacuum 
 (  = 1 ) for the 41 neutral organic molecules . as expected , the 
 dielectric medium increases the polarity of the molecules . interestingly , 
 the dipole moment is increased by an approximately constant factor 
 of 1.19 , independent of the chemical nature of the molecule . this 
 finding helps to justify the common use of scaled gas - phase charges 
 in condensed - phase modeling to account for polarization effects . in fact , the scaling relationship between the gas - phase and condensed - phase 
 dipole moments found here is very similar to the cm5 scaling factor 
 ( 1.20 ) , which was optimized to reproduce a range of experimental data . we have also plotted the line of best fit through 
 the dipole moments of the molecules computed with a dielectric constant 
 of  = 80 to model water ( dashed line in figure 3a ) . the polarity of the molecules is now 
 increased further , by a factor of 1.41 , relative to that under vacuum . it has been argued that charges for condensed - phase simulation should 
 be polarized halfway between those charges computed under vacuum and 
 those computed in a condensed - phase reaction field potential in order 
 to account for induction in an effective manner.figure 3a reveals that aim charges 
 computed in a dielectric medium (  = 4 ) should be fully consistent 
 with this charge - fitting philosophy . ( a ) comparison between condensed phase 
 (  = 4 ) and vacuum 
 dipole moments of 41 small molecules ( red crosses ) . dashed lines indicate 
 the dipole moments of molecules computed using the implicit solvent 
 model with  = 80 and  = 1 . ( b ) comparison between point - charge 
 and qm dipole moments of the same database , after the addition of 
 extra point sites to 12 molecules . an important measure of the likely success of aim charges 
 in condensed - phase 
 modeling is the ability of the point - charge model to reproduce the 
 electrostatic properties of the underlying qm calculation . the ddec 
 atomic densities are optimized to be close to spherically symmetric , 
 and hence the atomic multipole moments should be small and the electrostatic 
 potential surrounding the molecule well - approximated using an atom - centered 
 point charge model . however , there are certain well - documented cases , 
 such as atoms that possess lone pairs or -holes , in which the 
 atomic electron density is anisotropic and the aim atom - centered charges 
 are expected to produce a poor approximation to the qm electrostatic 
 potential . in section 
 s3 in the supporting information , the anisotropy in the atomic 
 electron densities is analyzed up to quadrupole order for every atom 
 in the 41 molecule test set . fifteen ( 15 ) molecules were identified 
 as containing at least one atom with significant anisotropy . all amines 
 ( 6 ) and molecules containing sulfur ( 5 ) or chlorine ( 2 ) were identified 
 by this analysis . the remaining two molecules were dimethyl ether 
 and methanol , although these were borderline cases . a common 
 approach in the design of mm force fields is to move charge 
 from the atoms to off - center sites to model anisotropic electron density , where the positions and magnitudes of the charges are parameters 
 that are fit to empirical data . in keeping with our intentions in 
 this work of minimizing the number of fitting parameters , the positions 
 and charges of the extra sites are automatically optimized to reproduce 
 the dipole and quadrupole moments of the aim atomic electron densities . the optimization procedure is described in section 
 s3 in the supporting information , and the errors in the atomic 
 dipole and quadrupole moments before and after the addition of extra 
 sites are given in table s3 in the supporting 
 information . figure 3b compares the dipole moments of all 41 molecules , using the aim 
 point charge model with the qm dipole moments , after the addition 
 of extra sites on 12 molecules . the mean unsigned error ( mue ) across 
 the entire set ( excluding two molecules whose dipole moment is zero 
 by symmetry ) is 0.10 d , compared to that without the extra sites ( mue 
 = 0.17 d ) . the positions and signs of the derived charges are 
 for the most 
 part chemically intuitive ( figure 4 ) . chlorine carries an off - site positive charge in 
 the position of the -hole , which has been shown to be important 
 for the accurate description of halogen bonding in mm force fields . the majority of the other molecules carry off - site 
 negative charges on lone - pair sites . it is striking that the number 
 of extra point sites required to reproduce the qm multipole moments 
 is larger than that typically used in mm force fields . this is perhaps 
 best exemplified by chloromethane to which the optimization procedure 
 has assigned three extra sites , rather than one extra site . the dipole 
 moment of chloromethane in the aim atom - centered point charge model 
 is 2.55 d , which is larger than the qm dipole moment in implicit solvent 
 ( 2.17 d ) . the -hole extra point site has a charge of + 0.20 
 e at 1.40  from the cl atom , which may be compared with the 
 extra site used in the opls - aax force field of + 0.075 e at 1.60  . however , with this single extra site , the dipole 
 moment of chloromethane decreases to 1.19 d. by allowing a second 
 site on the cl atom and one site on the neighboring c atom , the dipole 
 and quadrupole moments of the qm atomic electron densities are reproduced 
 by the extra sites , as is the net dipole moment of the molecule ( 2.14 
 d ) . it should be emphasized that the fitting procedure is automated 
 and parameter - free , thus all of the chemistry of the intermolecular 
 interaction is derived directly from the electron density , rather 
 than relying on chemical intuition . positions of extra point charge sites 
 on 12 molecules from the 
 benchmark set . it has 
 been shown in the previous section that , to a good approximation , 
 aim charges are capable of reproducing the qm electrostatics of small 
 molecules , especially when off - center charges are used to model anisotropy 
 in the electron density . our long - term goal is to investigate the 
 accuracy of aim nonbonded interactions in condensed - phase modeling , 
 particularly in the context of protein  ligand binding free 
 energies where the ligand is transferred from bulk solvent to the 
 binding site . pure liquid and hydration properties are commonly used 
 as a surrogate for this process ; liquid densities , heats of vaporization , 
 and free energies of hydration can be precisely measured experimentally 
 and the sampling requirements are low , so any discrepancies between 
 simulation and experiment can be directly attributed to the force 
 field . as such , these properties are widely used as a starting point 
 for the validation of new force fields . to investigate the accuracy 
 of aim charges for pure liquid simulations and free - energy calculations , 
 they were combined with the aim lj parameters described in section 2.1 and opls - aa 
 bonded parameters to form a complete mm 
 force field ( eq 1 ) . the 
 extra charge sites do not have lj parameters in the current model . the full list of nonbonded parameters for the 41 molecules studied 
 here are provided in the supporting information . figure 5 shows 
 the liquid densities and heats of vaporization , as well as the free 
 energies of hydration , for the 41 organic molecules computed using 
 the derived force field for the cases where experimental data are 
 available . the mue values , 
 compared with those of the experiment , were , respectively , 0.014 g / cm , 0.65 , and 1.03 kcal / mol . liquid densities are particularly 
 well - reproduced using the current method , which supports the derivation 
 of the aij parameters from the aim vdw 
 radii ( eq 7 ) . a modest 
 outlier is acetic acid , which is predicted to have a higher density 
 ( 1.11 g / cm ) than observed experimentally ( 1.04 g / cm ) . a notable feature of the computed energetics is that nonpolar 
 molecules , particularly those containing benzene rings , have a tendency 
 to be more weakly bound ( lower heat of vaporization and less exoergic 
 free energy of hydration ) than suggested by the experiment . possible 
 reasons for this discrepancy include the neglect of higher - order dispersion 
 terms beyond the dipole  dipole interaction , under - estimation of c ho hydrogen 
 bonding interactions , and a lack of explicit 
 treatment of polarization , which has been shown to induce the long - ranged 
 hydrophobic effect . properties of organic 
 molecules : ( a ) liquid densities , ( b ) heats 
 of vaporization , and ( c ) free energies of hydration . mean unsigned 
 error ( mue ) values , relative to the experiment , are also shown . the 
 raw data and experimental references are provided in tables s4 and s5 in the supporting information . nevertheless , the accuracy of the aim force field 
 is extremely 
 competitive with existing force fields . for example , the mue values 
 for a similar dataset computed using the opls / cm5 force field are , 
 respectively , 0.024 g / cm , 1.06 , and 0.94 kcal / mol . the mue values for free energies of hydration 
 for a large database of 239 small molecules was reported to be 1.93 
 kcal / mol ( charmm ) , 1.17 kcal / mol ( gaff ) , and 0.73 kcal / mol ( opls2.1 ) . the reported root - mean - square deviations for 
 the densities and heats of vaporization are 0.083 g / cm and 2.53 kcal / mol for gaff , and 0.026 g / cm and 1.12 
 kcal / mol for cgenff , which were computed for databases consisting 
 of more than 100 molecules . it should 
 be emphasized that the number of empirical fitting parameters in the 
 current study ( seven in total ) is substantially lower than in the 
 generation of these other force fields . all of the nonbonded parameters 
 in the present case arise directly from the dft calculations . a notable 
 feature of our aim - based nonbonded force field is that , 
 in principle , the methods scale to arbitrarily large system sizes , 
 potentially allowing the derivation of virtually parameter - free environment - specific 
 force fields for systems comprised of many thousands of atoms . while 
 our ultimate goal is to benchmark the accuracy of this force field 
 for the study of protein ligand binding free energies , this 
 is beyond the scope of the current study . instead , here , we demonstrate 
 the feasibility and potential benefits of using the aim approach in 
 computer - aided drug design efforts and compare the derived parameters 
 with a standard transferable force field . for our proof - of - principle 
 example , we have chosen to study the binding of two small molecules , 
 indole and benzofuran , to the l99a mutant of t4 lysozyme . this target 
 is an engineered hydrophobic binding pocket for which extensive binding 
 data and x - ray crystal structures are available . importantly , there are no water molecules in the binding pocket , 
 the protein adopts an almost - identical conformation in the two crystal 
 structures , and the ligands are small , relatively rigid , and isosteric . thus , the transformation 
 can be studied via a simple one - step free - energy calculation , sampling 
 requirements are low , and any discrepancy between experiment and theory 
 may be attributed directly to force field error . intriguingly , despite 
 the structural similarities between the two bound crystal structures , 
 free - energy methods are unable to reproduce the experimental binding 
 free energy of benzofuran , relative to indole ( 0.57 kcal / mol ) , with estimates of 2.3 kcal / mol ( am1-bcc 
 charges ) and 1.9 kcal / mol ( resp 
 charges ) . first , an aim force field 
 was derived for the two small molecules , 
 as in the previous section . the relative free energy of hydration of 
 the two molecules was computed and is shown in table 1 . no experimental data are available , to 
 the best of our knowledge , but the observation that indole has a much 
 lower free energy of hydration than benzofuran is consistent with 
 calculations performed using the opls / cm1a force field , the qm dipole 
 moments of the two molecules ( 2.19 d vs 0.67 d under vacuum ) , and 
 their observed octanol / water partition coefficients ( 2.14 and 2.67 ) . thus , it is not surprising that benzofuran should 
 be more strongly bound , inside the hydrophobic lysozyme pocket , than 
 the more polar , isosteric indole . all simulations 
 have been run in 
 triplicate , and the standard error in the mean is shown . x - sites 
 denotes the use of extra 
 off - center point charges on residue met102 . data taken from ref ( 73 ) . to 
 investigate whether the inclusion of native state polarization 
 in the force field derivation procedure affects the description of 
 protein  ligand binding , a protein - specific force field was 
 built for a significant portion of the t4 lysozyme protein ( 1646 atoms ) . figure 6a shows the correlation between the atom - centered 
 point charges derived using the aim approach , and opls - aa force field 
 charges for lysozyme . generally , there is good correlation , although , 
 as expected , the aim charges have a wider spread , because they can 
 respond to the macromolecular environment . the contrast in the computed 
 dispersion parameters is much more pronounced ( figure 6b ) . similar to all common force fields , opls 
 uses a library of atom types , while the aim dispersion coefficients 
 are derived from the atomic electron densities via eq 5 . as an example , both opls and aim 
 assign charges of 0.50 e and 0.14 e to the backbone 
 n atoms of residues leu84 and pro86 , respectively . however , opls assigns 
 the same b parameter to the two atoms ( 58 a.u . ) , 
 while the aim force field gives 51 a.u . the cluster of atoms highlighted by an arrow in figure 6b are carbonyl carbon atoms , which are electron - deficient 
 and thus intuitively expected to interact weakly through dispersive 
 effects . interestingly , of the 1646 atoms in the protein cluster , 
 only five atoms required off - center point sites to model anisotropy . four sites occurred on s atoms of met residues , and one on a neutral 
 lys residue , which is used by mcpro to neutralize the system . this 
 unexpectedly low number may , in part , be a selection effect , since 
 halogens and neutral amines are not found in the structure ( there 
 are no his residues in the modeled cluster ) . however , table s3 in the supporting information shows 
 that the anisotropy for met102 is much lower than on , for example , 
 dimethyl sulfide . in what follows , extra point sites were added to 
 met102 , since it is the only anisotropic site in contact with the 
 ligands ( figure 7 ) . comparison 
 between aim and opls ( a ) charge and ( b ) dispersion parameters 
 for the l99a mutant of t4 lysozyme . final snapshots of ( left ) benzofuran and ( right ) indole bound to 
 the hydrophobic pocket of the l99a mutant of t4 lysozyme from monte 
 carlo free - energy calculations . the relative free energy of binding of benzofuran and indole 
 to 
 the l99a mutant of t4 lysozyme was computed with free - energy perturbation 
 theory , using the mcpro software ( see table 1 ) . the standard opls - aa 
 force field , with scaled cm1a charges for the ligands , overestimates 
 the binding free energy of benzofuran , relative to indole , by 2 
 kcal / mol , in agreement with previous force field studies . replacing the nonbonded parameters by aim - derived charge and lj 
 coefficients reduces the relative binding free energy by 1 
 kcal / mol , which may be attributed to the lower ghyd in table 1 , as well as the altered protein ligand nonbonded 
 interactions . finally , adding the extra point sites to met102 ( figure 7 ) results in a relative 
 binding free energy of 0.37 kcal / mol , which is in excellent 
 agreement with the experiment ( 0.57 kcal / mol ) . thus , it appears 
 that the accurate treatment of the electrostatic potential around 
 the nh group of indole is an important determinant of binding . in 
 this regard , note that several experimental and theoretical studies 
 have highlighted the importance of n hydrogen 
 bonds in protein biochemistry . while the studied substitution 
 demonstrates the feasibility and benefits of the developed methods , 
 clearly many more comparisons with benchmark data must be performed 
 to ascertain the accuracy for drug discovery applications . in this paper , the feasibility of deriving 
 the nonbonded parameters 
 of molecular mechanics ( mm ) force fields from large - scale density 
 functional theory ( dft ) calculations has been demonstrated . atomic 
 charges and lj coefficients are all computed directly from partitioned 
 atoms - in - molecule electron densities , thus incorporating local and 
 long - ranged polarization that is specific to the environment under 
 study , while requiring only a minimal number of fitting parameters . this latter feature is a key part of the new force field . in the current 
 study , all of the nonbonded interactions of 41 organic molecules and 
 the lysozyme protein were described using just seven fitting parameters 
 ( the van der waals ( vdw ) radii of the free atoms in a vacuum ( r ) , presented in table 
 s1 in the supporting information ) . users of the force field 
 will be able to derive parameters for their own environment - specific 
 force fields using these simple input parameters and the procedure 
 outlined in figure 2 . in contrast , typical transferable force fields require a 
 large 
 number of fitting parameters to describe the wide ranging chemistry 
 of intermolecular interactions in organic molecules . in these standard 
 mm force fields , if one wishes to investigate the effect of computing 
 atomic charges with a method different from that was used to parametrize 
 the lennard - jones ( lj ) interactions ( for example , using a larger quantum 
 mechanics ( qm ) basis set ) , then , for the sake of consistency , the 
 entire library of lj parameters should be refit , which is a formidable 
 undertaking and one that discourages advances in force field accuracy . 
 in our proposed atoms - in - molecule ( aim ) force field , both charge and 
 lj parameters update consistently with changes in the qm atomic electron 
 densities . since liquid densities are extremely sensitive to the atomic 
 size , the r parameters are straightforward 
 to fit to the experiment and , thus , a new force field could be developed 
 within less than a week . this simplicity of design will allow the 
 straightforward investigation of new qm methods and potentially even 
 new functional forms of the force field . in this regard , a goal of 
 this investigation has been to investigate the accuracy of aim force 
 field parameters within the fixed functional form of eq 1 , such that the methods are consistent 
 with the majority of widely used mm codes . future studies will investigate 
 alternative functional forms of the repulsive vdw potential , such 
 as an exponential born - mayer form , and 
 the inclusion of explicit atomic polarizability . in this regard , combinations 
 of density - derived electrostatic and chemical ( ddec ) electron density 
 partitioning with standard tkatchenko scheffler ( ts ) methods and also with a new self - consistent screening method to compute 
 atomic polarizabilities and dispersion coefficients are currently 
 under development . note that only 
 the nonbonded parameters of the force field have 
 been derived in this study , while the bonded parameters have been 
 taken directly from the opls - aa force field . the compatibility of 
 the aim nonbonded parameters with the bonded parameters of standard 
 force fields  in particular , the torsional parameters  will 
 need to undergo substantial benchmarking against experimental properties 
 such as nuclear magnetic resonance ( nmr ) measurements and temperature - dependent structural propensities , which is beyond the scope of the current study . the accuracy 
 of our results for the small molecule benchmark set implies that the 
 use of opls - bonded parameters is a reasonable approximation , although 
 these molecules are relatively rigid . possible long - term strategies 
 for improving the compatibility between bonded and nonbonded parameters 
 involve on - the - fly parametrization of torsional parameters for an 
 aim protein and small molecule force field using existing and automatically generated qm 
 torsional scans , or machine learning of bonded force fields . the data presented in figure 5 reveal that the free energies 
 of hydration appear 
 to be difficult to predict computationally . part of the reason for 
 this may be that we have used the standard tip4p water model for these 
 calculations . while this model has been extensively validated against 
 experimental data , other combinations of nonbonded parameters also accurately describe many physical properties of water . in future research , it would be interesting 
 to develop a water model based on our aim force field parametrization 
 strategy to investigate whether it improves the correlation between 
 computed free energies of hydration and the experiment . one 
 of the advantages of the proposed methods , namely , the specificity 
 of the parameters to their environment , is also one of the potential 
 disadvantages . the derived force field is not expected to be transferable 
 between different systems , and , hence , a new qm calculation is required 
 to parametrize each new molecule under investigation . however , this 
 is not expected to add a substantial time cost to the investigation . 
 to draw an example from our own interests in computer - aided drug design , 
 it is already commonplace to perform a series of inexpensive qm calculations 
 to parametrize a set of small molecule drug candidates , while the 
 only additional cost in the new scheme would be a single large - scale 
 dft calculation to parametrize the target protein . in this investigation , 
 we have used ddec electron density partitioning 
 as implemented in the onetep linear - scaling dft code to compute the 
 atomic electron densities , and a variant of the tkatchenko  scheffler 
 method to compute vdw coefficients . ddec partitioning has been shown 
 to yield an efficiently converging multipole expansion of the electrostatic 
 potential , which is important for minimizing anisotropies in the atomic densities . 
 here , it has been shown that residual anisotropy may be accurately 
 represented by a small number of off - center charges . a further advantage 
 of ddec partitioning for flexible force field design is that the computed 
 atomic electron densities are relatively insensitive to small changes 
 in conformation , thus dynamical simulations of 
 a protein fluctuating about its native state are expected to be accurate , 
 while the applicability of the derived force fields to protein folding 
 may be questionable and will require extensive validation . the use 
 of the ts method for the computation of environment - specific dispersion 
 coefficients is , to the authors knowledge , its first use in 
 mm force field design , although grimme has recently used the d3 dispersion 
 scheme as part of a quantum mechanically derived force field ( qmdff ) 
 with extremely encouraging results . the 
 agreement between computed and experimental liquid properties and 
 free energies of hydration in the current study further supports the 
 use of nonempirical vdw parameters in force field design . in future 
 work , it will be interesting to investigate the effects of long - ranged 
 electrodynamic screening effects that are thought to dampen dispersion 
 interactions in the condensed phase . in this 
 paper , we have benchmarked the accuracy of aim nonbonded 
 parameters in the description of 120 experimental liquid densities , 
 heats of vaporization , and free energies of hydration . the obtained 
 results are extremely competitive with existing transferable force 
 fields that have been specifically parametrized against similar datasets . furthermore , all of the methods described here scale linearly with 
 system size and are accurate even for buried atoms , thus allowing the derivation 
 of environment - specific force fields for systems of unprecedented 
 size . recent advances in ls - dft algorithms and access to gpu computing 
 will allow rapid and routine electronic structure analysis of systems 
 comprising tens of thousands of atoms in the future . the largest molecule studied 
 in this paper is a substantial portion of the l99a mutant of t4 lysozyme , 
 consisting of 1646 atoms . the computed relative free energies of binding 
 of indole and benzofuran to the protein are encouraging and demonstrate 
 the importance of accurately modeling anisotropic electron density . while we would encourage users to experiment with deriving their own 
 environment - specific force fields by following the scheme suggested 
 in figure 2 , we emphasize 
 that further validation is required before their widespread adoption . our own priorities are ( i ) to test the compatibility between aim nonbonded 
 parameters and standard bonded force fields , and to develop alternative 
 bonded force fields if necessary , and ( ii ) to run hundreds of free - energy 
 calculations across multiple protein targets to extensively compare 
 the accuracy of environment - specific force fields with state - of - the - art 
 transferable force fields in the description of protein
molecular mechanics 
 force fields , which are commonly used in biomolecular 
 modeling and computer - aided drug design , typically treat nonbonded 
 interactions using a limited library of empirical parameters that 
 are developed for small molecules . 
 this approach does not account 
 for polarization in larger molecules or proteins , and the parametrization 
 process is labor - intensive . using linear - scaling density functional 
 theory and atoms - in - molecule electron density partitioning , environment - specific 
 charges and lennard - jones parameters are derived directly from quantum 
 mechanical calculations for use in biomolecular modeling of organic 
 and biomolecular systems . the proposed methods significantly reduce 
 the number of empirical parameters needed to construct molecular mechanics 
 force fields , naturally include polarization effects in charge and 
 lennard - jones parameters , and scale well to systems comprised of thousands 
 of atoms , including proteins . 
 the feasibility and benefits of this 
 approach are demonstrated by computing free energies of hydration , 
 properties of pure liquids , and the relative binding free energies 
 of indole and benzofuran to the l99a mutant of t4 lysozyme .