proteins 
 are essential components of all organisms , carrying out 
 tasks defined by the information encoded within genes , for example 
 catalyzing biochemical reactions , mediating cell signaling , or providing 
 structural rigidity .
simulations range from elucidation of enzymatic 
 reaction mechanisms , to the study of folding pathways , to design of 
 therapeutic molecules against disease . in biomolecular simulations such as these
, molecular mechanics ( mm ) 
 force fields are often used in which electrostatic interactions are 
 described by atom - centered point charges . however , there is no unique 
 method for partitioning the rigorously calculated quantum mechanical 
 ( qm ) electron density among the individual atoms and different charge 
 derivation schemes often lead to very different results . in 

commonly used force fields such as amber , the mm partial charges of protein molecules are optimized 
 by fitting them to reproduce the qm electrostatic potential ( esp ) 
 of small molecules .
these esp charges are well - suited 
 for mm force fields , as they reproduce ab initio multipole 
 moments and electrostatic interactions between molecular fragments .
a disadvantage of such techniques is the neglect 
 of polarization by the environment  indeed , a recent density 
 functional theory ( dft ) natural population analysis of an entire protein 
 in water found that net charges of residues can vary by up to 0.5 e from their putative integer values . while mean field approaches for charge fitting are the 
 most appropriate for deriving transferable force field parameters , 
 often , as in the example of the study of protein  ligand binding , 
 we are only interested in sampling in the vicinity of the protein s 
 native state . in these cases
, it would be ideal to incorporate electrostatic 
 polarization that is specific to that native state into the charge 
 fitting procedure .
recent studies have calculated atom - centered 
 charges for entire 
 proteins , accounting for native state polarization by including background 
 point charges in a series of iterative fragment - based esp fits .
the 
 resulting polarized protein - specific charges perform better than standard 
 amber charges in determining free energies of ligand binding , in pka calculations , and in comparisons with nmr data .
these results point to the potential improvements that can be made 
 by using polarized protein - specific point charges .
however , such an 
 approach potentially requires a large number of qm calculations to 
 iterate all charges in the system to self - consistency and requires restraints or 
 conformational averaging to treat buried atoms 
 and to address the sensitivity of the charges to small conformational 
 changes .
an ideal charge derivation scheme should efficiently 
 account for 
 the surrounding environment , while the resulting charges should be 
 chemically intuitive , reproduce ab initio electrostatic 
 properties , be robust with respect to conformational changes , and 
 be insensitive to buried atoms . the charges 
 should be derived from first principles , with no empirical parameters , 
 applicable to a wide range of systems without requiring specialized 
 treatments based on specific chemical knowledge of a particular molecule , 
 and preferably computable from a single qm calculation of the whole 
 system .
recently , there has been renewed interest in electronic density - based 
 atoms - in - molecule ( aim ) charge partitioning based on the hirshfeld 
 approach .
such methods differ conceptually from esp in that the net atomic 
 charges are assigned by dividing a converged , qm electronic density 
 into a union of overlapping atomic basins .
the density derived electrostatic 
 and chemical charges ( ddec ) method , developed by manz and sholl , combines two aim approaches , iterative hirshfeld 
 ( ih ) and iterated stockholder atoms ( isa ) , to assign atomic charges 
 from the electron density .
the resulting charges have already been 
 shown to be suitable for force field development .
they adapt to the atom s environment , 
 reproduce the ab initio electrostatic potential , 
 and , where applicable , correlate well with x - ray diffraction and x - ray 
 absorption near - edge spectroscopy data .
the method 
 can be applied with no adjustable parameters to buried atoms and to 
 either periodic or nonperiodic systems .
ddec charges have already 
 been used to develop force fields for molecular adsorption inside 
 metal  organic frameworks .
the ddec method 
 is implemented in a freely available code ( http://ddec.sourceforge.net/ ) , which is interfaced with codes such as vasp and gaussian among 
 others .
thus , ddec charges are suitable for environment - specific , 
 flexible 
 force field development for biomolecular simulations but are limited 
 by the unfavorable computational scaling of the underlying qm calculation 
 to systems of a few hundred atoms . in this paper , we overcome this 
 limitation by implementing the ddec scheme in the onetep linear - scaling 
 dft code .
onetep combines high basis set accuracy , comparable to that 
 of plane - wave dft methods , with a computational cost that scales linearly 
 with the number of atoms , which allows for an accurate , fully qm description 
 of systems of thousands of atoms , including 
 entire proteins .
we begin by outlining the various underlying aim 
 schemes , followed by a brief description of the linear - scaling dft 
 code onetep .
we validate our internal implementation of the ddec methodology 
 against quantum chemistry calculations in gaussian09 for a benchmark set of 25 representative small molecules 
 and show that charges are derived with linear - scaling computational 
 cost , allowing analysis of proteins comprising thousands of atoms 
 from a single dft calculation .
we demonstrate that , for these large 
 systems , the features of the ddec scheme that make the charges desirable 
 for flexible force field development are maintained .
namely , the charges 
 ( i ) respond to their environment in a chemically intuitive manner ; 
 ( ii ) reproduce electrostatic properties of the dft calculation ; and 
 ( iii ) are not overly sensitive to small conformational changes or 
 the presence of buried atoms .
finally , we construct a mm force field 
 based on the ddec charges for three proteins and compare the results 
 of our mm simulations with experimentally measured nmr dynamic observables 
 and a standard biomolecular force field .
the ddec method is based upon two recently developed 
 extensions to the original hirshfeld aim scheme . in the original formulation , an electronic density n(r ) is divided into overlapping atomic basins na(r ) for each atom a according 
 to the weighting formula:1where na0(r ) is the atomic 
 reference density , whose overlapping sum over all atoms in the system 
 bnb0(r ) is termed the promolecular 
 density .
this form of stockholder partitioning , where the electronic 
 density at each point r is distributed based on the proportional 
 contribution from the reference densities of each atom at that point , 
 has been shown to minimize the information distance faim ( kullback  liebler entropy ) between the real 
 and promolecular density , maximizing 
 the information retained in the reference atomic states
na0(r ) when transferring to the molecular environment:2subject to 
 the constraint that ana(r ) = n(r ) ; that is , the 
 electronic density is completely partitioned . 
 in other words , the stockholder partitioning exhaustively divides 
 the real electronic density into a set of overlapping atomic densities na(r ) in such a way as to maximize 
 the density distribution similarity to their respective reference 
 density counterparts na0(r ) .
shortcomings 
 of the original hirshfeld method included an arbitrariness in the 
 choice of reference atomic densities na0(r ) .
neutral , gas - phase atomic densities were often chosen as references , 
 although these commonly led to atomic populations that were too close 
 to zero .
such problems are addressed 
 in recently proposed iterative extensions to the hirshfeld method , in which reference densities are successively improved until self - consistency 
 is achieved . in the iterative hirshfeld ( ih )
, the ih reference densities na0(i , r ) are derived from the partitioned 
 atomic density at iteration i , na(i , r ) , by the following 
 procedure . first
, the ( noninteger ) electronic populations of each 
 atom qa(i ) are computed 
 in the usual way from the partitioned atomic densities:3instead of choosing neutral reference states 
 as in the original hirshfeld method , new reference states are generated 
 by linear interpolation between densities of free atoms or ions with 
 the next lowest integer (  = lint(qa(i ) ) and the next highest integer (  + 1 ) 
 number of electrons:4the 
 purpose of the interpolation is to obtain 
 a suitable gas - phase reference density for a hypothetical ion comprising qa(i ) electrons .
the partitioned 
 atomic densities for the next iteration i + 1 are 
 then derived from the reference states generated at iteration i:5the procedure 
 is iterated until the changes 
 in the set of ih charges fall below a specified threshold .
the resulting 
 charges have been shown to reproduce ab initio electrostatic 
 properties of small polypeptides and to be relatively insensitive 
 to small conformational changes .
an alternative approach by lillestolen and wheatley , meanwhile , named iterated stockholder atoms 
 ( isa ) , takes the spherical average of the partitioned atomic density na(i , r ) at iteration i as the ( isa ) reference density that enters into eq 5 , where a denotes spherical averaging about the center of atom a:6 in practice , the averaging is performed on 
 a set of discrete radial shells up to a maximum radius rmax .
the isa scheme is argued to be less empirical than 
 ih , as the latter still relies on a library of externally generated 
 ionic densities .
isa also produces a better fit to the esp due to 
 the low - order multipoles possessed by the converged na(r ) resulting from the spherical - averaging 
 procedure used to generate the reference densities .
the ddec scheme by manz and sholl combines the 
 ih and isa methods by minimizing a combined information entropy functional 
 ( eq 2):7where fih / isa are 
 constructed with ih / isa reference densities and weighted by an adjustable 
 parameter . minimizing eq 7 with respect 
 to na(r ) with the same constraint 
 as eq 2 leads to a partitioning of the form:8where naih / isa(i , r ) are the respective reference 
 densities given by eqs 4 and 6 . by allowing a fraction 
 of fih to contribute toward curvature 
 in regions that otherwise have shallow optimization landscapes , for 
 example
, buried atoms , this technique alleviates the slow convergence 
 of the isa method for such regions while retaining the appealing attributes 
 of the isa scheme .
an important addition to the ddec scheme is the 
 employment of charge - compensated reference ih densities , referred 
 to as the ddec / c2 scheme .
these densities 
 are generated from dft ground - state ionic calculations in the presence 
 of a charge compensation sphere , akin to a conductor - like polarizable 
 continuum model ( cpcm ) treatment with a solvent of infinite dielectric 
 constant .
the justification for using 
 such compensated densities instead of free ionic states , as in the 
 original ih scheme , stems from the dielectric screening experienced 
 by an ion embedded within a molecule , whereby its density profile 
 can be modified by the effective dielectric constant .
the charge compensation 
 sphere acts to expand ( contract ) the reference density in the case 
 of cations ( anions ) .
onetep is a linear - scaling dft package based on a reformulation 
 of conventional kohn 
sham dft in terms of the single - particle 
 density matrix:9where { (r ) } 
 are nonorthogonal generalized wannier functions ( ngwfs ) that are 
 localized in real space , and k is a representation of the density 
 matrix in the biorthogonal duals of the ngwfs .
onetep achieves linear - scaling 
 by exploiting the  nearsightedness  of the single - particle 
 density matrix in nonmetallic systems . in practice ,
linear - scaling arises from enforcing strict localization 
 of the ngwfs onto atomic regions and through the optimization of the 
 density kernel and ngwfs , subject to localization constraints .
optimizing 
 the ngwfs in situ allows for a minimal number of 
 atom - centered orbitals to be used while maintaining plane - wave accuracy . 

the ngwfs are represented in a basis of highly localized periodic 
 cardinal sine ( psinc ) functions ( otherwise known as fourier 
the psinc functions 
 are related to plane waves via a fourier transform , meaning that systematic 
 improvement is possible through adjustment of the psinc grid spacing , 
 analogous to converging the kinetic energy cutoff in traditional ( n ) plane - wave 
 dft codes . in order to expedite optimization
, the ngwfs can be initialized 
 closer to their ground states by using an in - built pseudoatomic solver , 
 which self - consistently solves the dft kohn  sham equations 
 for isolated atoms using the same pseudopotentials and exchange - correlation 
 functional as the full calculation . implicit 
 solvation , whose inclusion is essential both for an accurate description 
 of the protein s aqueous environment and to aid optimization 
 of the density kernel , is implemented 
 within onetep .
this is a minimal parameter , self - consistent model 
 based on direct solution of the inhomogeneous poisson equation for 
 a solute cavity defined by the isosurface of the electron density .
we 
 have implemented the ddec method within onetep 
 as a postprocessing module to take advantage of its parallel and efficient 
 algorithms for sparse matrix algebra and operations with localized 
 orbitals and electron distributions .
a single dft calculation is performed 
 on the system to obtain the ground - state electronic density , which 
 is then processed to extract the ddec net atomic charges .
reference 
 densities for the ih part of the calculation are generated internally 
 at run time using the same exchange - correlation functional , pseudopotentials , 
 and ngwf cutoff radius as the full dft calculation , with the c2 charge 
 compensation scheme .
specifically , using 
 the pseudoatomic solver module , as described 
 in section 2.2 , the kohn  sham equations 
 for a set of atomic orbitals of an isolated atom are solved self - consistently 
 on a radial grid with an additional spherical surface of compensation 
 charge .
the effect of this charged surface 
 is expressed equivalently via gauss law as an additional radial 
 electrostatic potential term vcomp in 
 the hamiltonian:10where the 
 first three operators are the kinetic 
 energy , local and nonlocal potentials , and vcomp = qcomp/|r| for 
 |r| > rcomp , with qcomp and rcomp the 
 charge and radius of the compensation sphere . following manz 
 and sholl , for cations , the magnitude 
 of the compensation charge ( qcomp ) is 
 chosen to neutralize the ion , with the compensation sphere radius 
 ( rcomp ) set to the average of nlm|r|nlm 
 expectation values of all occupied orbitals { |nlm } in the neutral species that are vacant in the cation . for 
 anions ,
the compensation charge is chosen to make the electrostatic 
 potential at an infinitesimal distance outside the compensation sphere 
 zero .
the radius of the compensation sphere is incrementally adjusted 
 in 0.1 bohr steps until the total energy of the system is minimized
. whole - molecule densities ( total and promolecular ) are stored on 
 regular cartesian grids , while spherically averaged promolecular density 
 profiles for individual atoms are computed and stored on atom - centered 
 sets of equally spaced radial shells up to a predefined maximum cutoff 
 radius .
linear - scaling computation with respect to the number of atoms 
 ( for each iteration ) is achieved at a modest cost of computational 
 memory by computing and storing the total promolecular density ( bnb0(i , r ) ) once per 
 ddec iteration , on the same grid as the molecular density .
figure 1a shows the computational scaling of our ddec implementation 
 for periodic supercells of bulk water of increasing size , up to around 
 2500 atoms ( supporting information ) .
even 
 for the largest system , the computation time of the ddec analysis 
 is under 20 min .
the dft calculation itself requires around 12 h on 
 160 cores , making the calculation of ddec charges for systems comprising 
 thousands of atoms feasible on a near - routine basis .
( a ) computation time 
 required for ddec postprocessing calculations 
 of bulk water .
( b ) ddec atomic and molecular charge distributions 
 of bulk water for three system sizes .
importantly , the number of iterations required for charge 
 convergence 
 remains virtually constant for increasing numbers of atoms , ensuring 
 that the total computational effort remains linear - scaling with respect 
 to system size .
figure 1b reveals that the 
 distribution of charges obtained for bulk water is independent of 
 system size .
the atomic charges distributions are centered on values 
 close to those typically used in mm force fields ( tip3p ) , while they show some spread as one might expect 
 for atoms in a liquid with differing local environments .
recently , 
 manz and sholl introduced an update to the ddec method , 
 named ddec / c3 .
the ddec / c3 method improves 
 convergence of charges for nonporous solids with short bond lengths 
 between diffuse atoms by enforcing exponential decay of atomic electron 
 densities and conditioning the reference ih densities naih(i , r ) that enter eq 8 to 
 the material of interest . in our implementation
, we use the conditioned 
 reference densities but do not enforce radial decay of the partitioned 
 atomic densities na(r ) .
our 
 reasoning is that our studies are aimed at biomolecular systems , which 
 are porous materials lacking compacted regions , and we anticipate 
 that the partial weighting by ih reference densities during optimization 
 will be sufficient to ensure the stability of atomic charges for embedded 
 atoms .
in addition , in our implementation , core charge handling has 
 been completely excluded , as we employ norm - conserving pseudopotentials . 

literature studies have found little dependence of ih charges on the 
 treatment of core electrons .
nevertheless , 
 in the following section , we test both of these assumptions by validating 
 onetep - calculated ddec charges against the full ddec / c3 method employing 
 the all - electron code gaussian09 .
the validity of our ddec 
 implementation was tested using a set of 25 diverse , neutral , small , 
 organic molecules .
benchmarking was performed 
 against the standard ddec method as implemented in the chargemol package 
 ( version 2.1 beta , obtained from http://ddec.sourceforge.net/ ) , employing 100 radial shells ( nrad ) 
 up to a maximum cutoff radius ( rmax ) of 
 5  , together with the supplied c3 charge - compensated reference 
 density library ( ddec / chargemol ) .
molecular 
 geometries were optimized in vacuo with a 6 - 311g(d , p ) 
 basis set , and electronic densities were generated with an aug - cc - pvqz 
 all - electron basis set in gaussian09 .
the mixing parameter 
  was set to /14 , which has been shown 
 to give the optimal balance between minimizing the atomic multipoles 
 and maximizing chemical accuracy .
the 
 set of net atomic charges was considered converged when the maximum 
 absolute change for any atom for three successive iterations was less 
 than 2  10e .
onetep 
 calculations were performed in vacuum using a cubic simulation cell 
 of 30  , with the spherical cutoff coulomb approach to avoid 
 electrostatic interactions between molecules and their periodic images .
ngwfs were initialized 
 as orbitals obtained from solving the kohn  sham equation for 
 free atoms , with a 1s configuration for h , a 2s2p configuration 
 for c , n , o , and f , and a 3s3p configuration 
 for s and cl .
the ngwfs were expanded in a psinc basis with an energy 
 cutoff of 1000 ev , corresponding to a grid spacing of 0.45 bohr , and 
 were localized in real space with radii of 10 bohr .
convergence of 
 the atomic charges generated by our implementation of the ddec method 
 in onetep ( referred to as ddec / onetep ) with respect to the onetep 
 psinc grid spacing and ngwf cutoff radii , as well as the ddec parameters rmax and nrad was 
 investigated for nitroethane and discussed in supporting information figure s1 .
figure 2 shows the correlation between the 
 ddec / onetep and ddec / chargemol charges for every atom in the benchmark 
 set . despite the different approaches used in obtaining the ground - state 
 electronic densities and the subsequent charge analysis ,
the difference 
 between the two charge sets is very small with a mean absolute deviation 
 ( mad ) of less than 0.02 e. correlation between ddec / chargemol 
 and ddec / onetep charges for 
 all atoms in the 25 molecule benchmark set .
also shown is the mean 
 absolute deviation ( mad ) between the two sets .
atomic charges for use in force fields should approximately 
 reproduce 
 the ab initio electrostatic potential outside the 
 molecule s electron density .
we can measure the error v in the coulombic potential of the ddec charges , compared 
 with dft , as11where the sum is performed for all 
 points i lying within 1.4 and 2.0 times the van der 
 waals radii 
 of the nuclei , on the same grid mesh 
 used to calculate the electron density . in order to remove the arbitrary 
 vacuum level of the dft potential v(r ) , the potentials are displaced by the averaged difference 
 over the included grid points .
our calculations of v are performed for a test charge of 1 e in vacuum .
typical errors in , for example , interaction energies 
 of neutral ligands with proteins in solution will tend to be lower 
 than this error estimate . nevertheless , 
 it is a good indicator of the ability of a method to reproduce the ab initio esp .
figure 3 ( bottom ) 
 shows that ddec charges 
 reliably reproduce the electrostatic potential around the molecule , 
 with v errors of approximately 1 kcal / mol 
 for both ddec / chargemol and ddec / onetep . for neutral molecules ,
the 
 dipole is the leading order term in the multipole expansion and , thus , 
 is the most significant in determining intermolecular electrostatic 
 interactions . in figure 3 ( top ) , we compare 
 the ddec vector dipole moments with those computed directly from their 
 respective ab initio electronic densities . again , 
 onetep performs as well as chargemol . for comparison
, we have also 
 plotted the error in the dipole moment and v ( eq 11 ) for esp charges . in agreement with 
 previous charge comparisons ,
the esp charges reproduce the ab initio electrostatic 
 potential and dipole moment more closely . despite this better performance 
 of esp charges for small molecule electrostatics
, ddec charges will 
 be suitable for polarized protein - specific charges if , as we shall 
 investigate in the next section , they are able to reproduce the electrostatic 
 properties of much larger molecules containing buried atoms better 
 than standard mm force fields .
we conclude that , at least for organic systems comprising 
 light 
 atoms , our use of ddec charge analysis within the pseudopotential 
 approximation and without enforcing radial decay of the partitioned 
 atomic densities is valid .
it should be emphasized that the ddec analysis 
 can only be as accurate as the electron density output by the underlying 
 qm calculation . in the current paper , we have used the pbe exchange - correlation 
 functional in our onetep calculations , although future advances will 
 allow alternative treatments of electron  electron interactions , 
 such as b3lyp and dynamical mean field 
 theory . for comparison , supporting information figure s2 compares gaussian09 ddec 
 charges calculated using pbe and b3lyp exchange - correlation functionals . 

the mad is less than 0.02 e , implying that the pbe 
 functional gives a reliable electron density for the calculation of 
 ddec charges .
in order 
 to aid in the interpretation of qm simulations , atomic 
 point charges should respond in a chemically intuitive manner to their 
 environment . as an example , figure 4 , b and 
 c , shows onetep / ddec charges for the phenol molecule in two different 
 environments : first in water and , second , in a small model cluster 
 representing the negatively charged binding pocket of the l99a / m102e 
 mutant of t4 lysozyme ( supporting information ) .
the two environments are typical 
 of those simulated in , for example , the optimization of small molecule 
 inhibitors for drug design . for comparison , 

figure 4a shows the resp charges for phenol , 
 calculated in vacuum at the hf/6 - 31 g * level to approximate aqueous 
 polarization , the same level of theory 
 that has been used to parametrize the amber force field that we employ as a benchmark in later sections . 

the ddec charges in water agree with the resp charges with a root - mean - square 
 ( rms ) deviation of 0.04 e. similar results are obtained 
 by comparing with resp charges derived using a polarizable continuum 
 model ( pcm ) to model solvation ( supporting information figure s3 ) .
there are no large changes in the ddec charges of the 
 phenol molecule on moving to the protein binding pocket .
charges on 
 the aromatic ring change by up to 0.05 e in the more 
 hydrophobic environment , while the oh group , which is hydrogen - bonded 
 to residue e102 , becomes more strongly polarized .
( a ) resp charges for 
 phenol calculated in vacuo at the hf/6 - 31 g * level . 

( b ) ddec / onetep charges calculated in a 
 20 molecule water cluster using with identical simulation parameters 
 as section 2.4 .
( c ) same as in part b , but 
 within a 86 atom cluster representing the t4 lysozyme l99a / m102e binding 
 pocket .
standard error of the mean is less than 0.01 e for the oh group , and the maximum on any atom 
 is 0.018 e. we have utilized the linear - scaling nature of our ddec implementation 
 to calculate ddec charges of three proteins : ubiquitin , the smn tudor 
 domain , and hen egg lysozyme ( pdb : 1ubq , 1mhn , 6lyt ) .
initial structures were prepared by 
 protonation using the amber11 tleap module or by the molprobity software .
dft calculations were performed using the minimal parameter 
 continuum solvation model of onetep with a relative dielectric permittivity 
 of 80 to represent water solvent .
a postprocessing 
 calculation was performed on the converged electronic density to extract 
 ddec charges ( supporting information ) . 

it should be emphasized that the entire charge set is derived from 
 a single dft calculation , on systems up to 1960 atoms , and thus naturally 
 incorporates native state polarization . as we found for the 
 small test systems , ddec charges for the proteins 
 are similar to the mean field amber ff99sb charges , with a correlation coefficient of 0.96 ( figure 5a ) .
however , the rms deviation between the charge 
 sets is 0.11 e , as expected since the aims were optimized 
 based on the ab initio electronic density for atoms 
 embedded in their local environment .
comparisons with the amber ff03 force field charges yield similar conclusions 
 ( supporting information figure s4 ) . for 
 chemically identical atoms that share the same charges in amber ,
the 
 ddec charges span a wider range , which is indicated by the vertical 
 spread of points .
for example , ddec charges for backbone o atoms of 
 leucine vary between 0.72 and 0.56 e , with an average of 0.63 e , compared to 
 the constant amber value of 0.57 e. the prominent 
 outliers in figure 5a ( with absolute charge 
 disagreement larger than 0.25 e ) are mostly nonbackbone sp carbon atoms , most of which are close to 
 neutral in amber , but adopt a range of values in ddec , depending on 
 the residue type and their particular local environments .
for example , 
 the c atom in glutamine has a charge of 0.06 e in amber , but ranges from 0.41 to 0.37 e in ddec .
similar deviations in c charges from amber values 
 have been observed in conformationally averaged esp charges derived 
 from dipeptide fragments .
the largest 
 discrepancy between the two charge sets is in the terminal nitrogen 
 atom , which is positively charged in amber ( 0.16 e ) but negatively charged in ddec ( 0.52 e ) .
several esp studies have indicated that our ddec 
 value is more appropriate for the charge of the n - terminus .
( a ) correlation 
 between ddec partial atomic charges and the amber 
 ff99sb force field for ubiquitin .
( b ) electrostatic contribution to 
 the backbone hydrogen bond energies in the x - ray crystal structure 
 of ubiquitin .
backbone 
 hydrogen bonds calculated by ddec / onetep are more polarized than in 
 the amber force field .
following ji et al . , we define the electrostatic contribution to the hydrogen bond between 
 backbone nh and co groups as the coulombic interaction between atoms 
 carrying point charges q:12figure 5b reveals that 
 the resulting dipole 
it is possible that the enhancement of the interaction 
 is a systematic bias in the potential of the ddec charges , though 
 it does support previous hypotheses that proteins are stabilized in 
 their native environment via electronic polarization .
we investigate this effect further in section 4 by comparing the dynamics of backbone 
 hydrogen bonds with experimental nmr observables .
the correct treatment of electrostatics 
 is vital in the accurate determination of molecular interactions in 
 biological systems .
in particular , if ddec charges are to be useful 
 as polarized protein - specific charges , they should reproduce ab initio electrostatics for large molecules better than 
 standard force fields that are based on esp charges .
we have calculated 
 the electrostatic potential and the dipole moments of the three proteins 
 studied in the previous section with full dft , the derived ddec charges , 
 and with a standard amber force field ( supporting 
 information figure s4 gives the same information for the amber 
 ff03 force field ) . for all charge distributions ,
the electrostatic 
 potential was computed in vacuum , in order to better observe the discrepancy 
 between ab initio and point charge values that would 
 otherwise be dampened by dielectric screening .
we note , however , that 
 the ab initio , ddec and amber charge distributions 
 have already incorporated solvation effects in their calculation , 
 and the electrostatic comparisons are , therefore , equivalent .
table 1 shows that ddec / onetep aims reproduce the total ab initio potential and dipole moments dft for all three proteins , vastly outperforming the standard 
 mm charges .
figure 6 shows the difference in 
 electrostatic potential between the two charge sets and the dft calculation . 

although the standard mm charges are fit to closely reproduce the 
 esp of small molecules , they neglect large - scale electronic polarization 
 that is present in protein - specific charges , and hence , the error 
 in the computed esp is large .
the computed errors in the mm electrostatic 
 potential are consistent with literature studies of errors in interaction 
 energies between charged biotin ligands and the avidin protein in 
 vacuum ( 11 kcal / mol on average ) . for 
 more realistic simulations , the same study found that the error was 
 reduced to approximately 2 kcal / mol for neutral ligands and to 0.7 
 kcal / mol if solvent screening is also included .
the esp difference is calculated at 
 the solvent - accessible surface of the protein .
also shown as gray , 
 blue , and red arrows are the calculated dft , ddec , and amber dipole 
 moments , respectively .
also shown are 
 errors in the 
 protein dipole moment vectors , which have been calculated relative 
 to the center of mass of each protein to remove ambiguity for charged 
 systems .
the ddec charges reproduce electrostatic properties of large 
 molecules much better than do the mean field classical point charges .
we noted earlier that isa charges 
 are particularly good at reproducing 
 the esp , since the electron density is divided into spherically symmetric 
 partitions .
this point is reinforced by figure 7 , which reveals that the deviation from the ab initio dipole moment of ubiquitin decreases with decreasing  ( increasing 
 isa fraction ) .
however , one of the motivations for including ih reference 
 densities in the ddec information entropy functional ( eq 8) was to alleviate problematic convergence for aim charges 
 of embedded atoms due to the shallow optimization landscapes experienced 
 in such regions .
we find that the ddec 
 method with  = /14 , as recommended for 
 all systems by manz and sholl , is a good 
 compromise between reproducing accurately the esp of the system , reducing 
 the number of self - consistency iterations required for convergence 
 ( figure 7 ) and , as we shall see in the next 
 section , improving the transferability of the charges .
further comparisons 
 between ddec , ih , and isa charges are presented in supporting information figure s5 .
dependence of the ddec 
 dipole moment error on the mixing parameter 
 . increasing the isa content reduces the dipole error but also 
 increases the number of iterations required for charge convergence 
 and , hence , the computational time . in 
 order to be of use in parametrizing flexible force fields , atomic 
 charges , even those of buried atoms
, must not be overly sensitive 
 to small conformational changes . to study conformational dependencies 
 of the ddec charges for large systems
, we have performed ddec analysis 
 of nine experimental nmr conformers of bovine pancreatic trypsin inhibitor 
 ( bpti , pdb : 1pit ) in implicit solvent ( 892 atoms ) .
each residue in figure 8 is colored by the average of the standard deviation 
 in the calculated atomic charges over the nine conformations .
the 
 first point to note is that the ddec charges are much more stable 
 to conformational change than the isa charges .
while some fluctuation 
 due to the changing environment is expected , the protein conformers 
 are all close to the native state and the ddec charges reflect this . 

this is in stark contrast to resp charges , for which a recent review 
 of atomic charge models found strong fluctuations in partial charges 
 for an ensemble of polypeptide chain conformations .
second , there is no discernible difference between charge 
 fluctuations in buried and exposed residues , as would be observed 
 for shallow optimization landscapes .
finally , the width of the protein 
 chains in figure 8 reflects the average positional 
 fluctuation of each residue .
it is noticeable that ddec charges vary 
 more strongly in residues with a high degree of flexibility , especially 
 when they form intramolecular contacts ( between -sheets , within 
 -helices , and in sections linked by disulfide bridges ) , as 
 one might expect .
ddec / onetep charge and conformation fluctuations for nine 
 nmr models 
 of bpti .
average positional fluctuation per residue is represented 
 by its putty width and average charge fluctuation ( standard deviation ) 
 per residue by its color .
isa charge fluctuations have been clipped 
 at 0.04 e for visual clarity ( ten residues in the 
 isa method have values greater than 0.04 e , with 
 a maximum of 0.06 e ) .
isa charges are ddec / onetep 
 charges derived with  = 0 ( eq 7 ) .
we showed in the previous section 
 that ddec charges are able to 
 accurately reproduce the esp around the structure that they are fit 
 to . here
, we investigate whether average charges derived for an ensemble 
 of native state conformations are able to accurately describe the 
 electrostatic properties of each conformation close to that native 
 state .
such a test is vital if the charges are to be used in the construction 
 of force fields for simulating the dynamics of proteins .
figure 9 shows the discrepancy between dipole moments calculated 
 via atomic point charges and dft . to obtain the best agreement with 
 the dft dipole moment of a particular structure , it is preferable 
 to fit the ddec atomic point charges to that particular structure . 

however , encouragingly , there is very little degradation in performance 
 if ddec charges that have been averaged over the full ensemble are 
 used instead .
this reflects the transferability of ddec charges between 
 conformations close to the native state .
the force field charges perform 
 relatively well for this structure but not as well as the polarized 
 protein - specific charges .
this behavior is in agreement with previous 
 work showing that esp charges give the best fit to ab initio electrostatics for a given structure , but that ddec charges are 
 generally more transferable to different , similar structures . deviation between charge model and ab initio dipole 
 moment vectors for nine nmr conformers of bpti .
ddec / onetep charges 
 have been calculated separately for each structure ( dark blue ) and 
 averaged over all structures ( light blue ) .
as we have demonstrated in the previous section , 
 ddec charges derived 
 from the dft electron density of entire proteins have properties that 
 make them desirable for use in flexible force fields .
we now put this 
 into practice by performing molecular dynamics ( md ) simulations of 
 the three proteins described previously ( pdb : 1ubq , 1mhn , 6lyt ) using force fields 
 based on our calculated ddec charges .
we have followed the procedure 
 of tong et al . by taking the bonded and 
 lennard - jones parameters directly from the amber ff99sb force field 
 but replacing the atom - centered point charges by the ddec / onetep charges . 
 it should be emphasized that torsional and lennard - jones parameters 
 are as important as the charges in determining protein dynamics and 
 a large number of alternative force fields exist , many of which outperform 
 the ff99sb force field in comparison with experiment .
nevertheless , here we concentrate on demonstrating the 
 feasibility of using aim approaches to supplement an existing , commonly 
 used force field .
figure 10 shows the 
 backbone rms deviations 
 of the three proteins from their x - ray crystal structures over the 
 course of 10 ns simulations , performed at 300 k in explicit water 
 ( supporting information ) . for both force 
 fields ,
backbone c root - mean - square deviation 
 ( rmsd ) 
 with respect to the initial experimental pdb structure for 1ubq , 1mhn , and 6lyt , with running averages 
 taken over 50 ps windows .
the first and last four residues have been 
 excluded due to their high flexibility . to assess the quality of simulated mm dynamics , it is common 
 to 
 use , as metrics , experimental nmr - derived quantities .
the square of 
 the generalized order parameter s , based 
 on the lipari 
szabo model - free approach to analyzing nuclear 
 spin - relaxation measurements , is the 
 plateau value of the time - correlation function c2(t ) of the second - order legrendre polynomial 
 of backbone n  h bond unit vectors e(t):13where c2(t ) is time - averaged over 
 a trajectory .s represents the 
 spatial restriction of the motion of n  h bond vectors , and 
 a lower s value corresponds to a more 
 flexible residue .
s values were calculated 
 from our 10 ns md simulations using the ired method as implemented in the amber11 ptraj module . 
 in this method ,
h unit bond vector pairs ei , j in a fixed reference frame and 
 averaged over the entire trajectory:14s of residue i is then obtained 
 from solving the eigenvalue equation mm = mm as15where the sum is 
 taken over n  5 eigenvectors excluding five 
 with the largest eigenvalues , 
 and mi is the ith component of eigenvector m. figure 11 compares simulated and experimentals profiles for the standard amber ff99sb 
 force field and the same force field supplemented by ddec / onetep atomic 
 charges .
the only difference between the simulation protocols is in 
 the point charges and so any improvement in the calculated order parameters 
 is due to the inclusion of native state polarization in their calculation
. 
 table 2 compares the agreement between simulated 
 and experimental s data ( for all residues 
 where experimental data is available ) .
it indicates that ddec aim 
 charges perform better than mean field force field charges in providing 
 a suitable electrostatic environment that maintains protein stability 
 throughout the 10 ns trajectory while remaining dynamically consistent 
 with experimental observations .
h 
 order parameter s for 1ubq , 1mhn , and 6lyt using amber ff99sb 
 and ddec atomic charges .
the agreement of the generalized 
 order parameter with experiment is improved by using the ddec / onetep 
 protein - specific charges .
in addition to backbone rigidity , we also compared the nmr scalar 
 coupling jnc of the n 
backbone 
 hydrogen bonds , which provides a measure of hydrogen bond dynamics 
 within a protein .
it was shown by barfield that jnc values , being dependent on the wave function overlap across the 
 hydrogen bond , are strongly geometry dependent and can be parametrized 
 as such by several formulas , with the simplest based on only the ho 
 distance roh and ho = c 
 angle :16equation 16 was used to compute 
 the time - averaged jnc of backbone 
backbone 
 hydrogen bonds using the amber ff99sb and ddec / onetep 
 mm snapshots over the 10 ns trajectories .
the rms deviations between simulation and experiment 
 ( figure 12 , in brackets ) once again indicate 
 that ddec / onetep charges perform at least as well as amber charges , 
 illustrating that backbone n  h and c = o bond polarization 
 is also suitably described by the electron density partitioning approach 
 to charge derivation .
the 
 ddec atoms - in - molecule scheme for the assignment of atomic 
 charges through partitioning of the ab initio electron 
 density has a number of features that makes it desirable for charge 
 derivation for flexible force fields .
however , its 
 use in the analysis of realistic biomolecular systems has , until now , 
 been limited by the computational expense of the underlying qm calculation . 
 in this paper , we have implemented the ddec scheme in the linear - scaling 
 dft code onetep .
strict localization of the ngwf basis in onetep onto 
 atomic regions allows derivation of the electronic density of systems 
 comprising thousands of atoms with affordable computational expense . 
 optimizing the ngwfs in situ allows for a minimal 
 number of atom - centered orbitals to be used while maintaining near - complete 
 basis set accuracy .
systematic improvement in accuracy is achievable 
 via tuning of a very small number of parameters and , along with a 
 minimal parameter solvation model , allows the user to derive charges 
 for any molecular system with no prior knowledge of its chemical characteristics . 

the ddec method is also suitable for implementation in other dft codes , 
 as well as for example fragment molecular orbital approaches to generating 
 the electron density of large systems , for which we hope our implementation 
 will provide an accuracy benchmark . we have validated our implementation 
 of the ddec method in onetep 
 against the standard chargemol package with quantum chemistry calculations 
 in gaussian09 . in agreement with previous observations ,
ddec charges 
 perform less well than esp charges in reproducing the ab initio potential of small molecules .
nevertheless , 
 our interest here is in the calculation of polarized charges for entire 
 proteins .
indeed , we have found that the errors in the electrostatic 
 potential at the surfaces of three proteins are very similar to those 
 for small molecules ( errors of 1.31.5 kcal / mol ) .
furthermore , 
 the system - specific ddec charges outperform standard amber charges 
 for these molecules ( errors of 57 kcal / mol ) .
we would expect 
 similar performances for other force fields that have been fit to 
 small molecules without incorporating large - scale polarization .
although the esp method is unsuited to the direct calculation of 
 atomic charges from a single qm calculation due to the sensitivity 
 of the charges to buried atoms and small conformational changes , a number of studies have determined protein - specific 
 charges based on this method . by separating the system into a large 
 number of dipeptide fragments , fitting esp charges to the resulting 
 structures and averaging them over several conformations of the protein , 

stable and transferable charges have been obtained . by also including the effects of the environment , represented by 
 point charges , in an iterative esp fit
, it should be possible to further 
 derive polarized protein - specific charges that are stable with respect 
 to small changes in conformation . although 
 the esp charges are preferable to ddec for reproducing the ab initio potential of small molecules , it is not obvious 
 that these conformation - averaged esp charges are the optimum choice 
 for large molecules .
in fact , manz and sholl have shown that ddec 
 and conformation - averaged esp charges give essentially the same agreement 
 with the ab initio electrostatic potential of small 
 molecules .
it would be interesting , in 
 future work , to make this comparison for larger molecules .
further 
 advantages of the ddec charges are that , since they are 
 derived from a single dft calculation of entire biomolecules , environmental 
 polarization is naturally included .
moreover , no prior , system - specific 
 chemical knowledge or specialized treatment is required to ensure 
 proper charge derivation , rendering the method extremely versatile . 
 in agreement with previous observations of smaller systems ,
the charges 
 are insensitive to small conformational changes and buried atoms require 
 no special treatment . in this paper , we are interested in the application 
 of ddec charges to standard mm force fields .
therefore , we have restricted 
 the ddec expansion to monopole order in our electrostatic potential 
 comparisons . as discussed by manz and sholl
, the electrostatic potential description of ddec charges can be trivially 
 improved by including higher order multipole terms in the partitioned 
 aim densities , without requiring a new fitting procedure .
we 
 have incorporated the ddec charges of three proteins into a 
 classical force field and run molecular mechanics simulations to compute 
 nmr order parameters and scalar couplings .
our observations indicate 
 that the ddec aim charges perform at least as well , if not better 
 than amber ff99sb in replicating protein dynamics .
backbone hydrogen bonds are more polarized in ddec / onetep 
 than in the standard amber force field .
nevertheless , they are able 
 to replicate dynamics consistent with experiments when used in mm 
 force fields . our results are therefore consistent with previous hypotheses 
 that electronic polarization is important in stabilization of a protein s 
 native state .
it should be noted that 
 these ddec aims were obtained without explicit fitting to replicate 
 electrostatic effects nor parametrization to experimental properties .
the potential applications of linear - scaling ddec analysis span 
 a wide range of problems in which accurate determination of electrostatics via atomic point charges is important .
ddec charges have 
 already been used to study molecular adsorption inside metal  organic 
 frameworks . in the biomolecular sciences , 
 the optimization of protein - inhibitor interactions for drug design , 
 pka calculation , study of ion channel 
 conductivity and elucidation of enzymatic reaction mechanisms all 
 depend critically on an accurate description of the electrostatic 
 potential , for which ddec perform markedly better than standard classical 
 force field charges .
 atomic partial charges for use in 
 traditional force fields for 
 biomolecular simulation are often fit to the electrostatic potentials 
 of small molecules and , hence , neglect large - scale electronic polarization . 
 on the other hand , recent advances in atoms - in - 
 molecule charge derivation 
 schemes show promise for use in flexible force fields but are limited 
 in size by the underlying quantum mechanical calculation of the electron 
 density . here 
 , we implement the density derived electrostatic and 
 chemical charges method in the linear - scaling density functional theory 
 code onetep . 
 our implementation allows the straightforward derivation 
 of partial atomic charges for systems comprising thousands of atoms , 
 including entire proteins . 
 we demonstrate that the derived charges 
 are chemically intuitive , reproduce ab initio electrostatic 
 potentials of proteins and are transferable between closely related 
 systems . 
 simulated nmr data derived from molecular dynamics of three 
 proteins using force fields based on the onetep charges are in good 
 agreement with experiment . 
Introduction
Computational
Methodology
Properties of DDEC Charges
Protein
Dynamics
Conclusions