studies of complex 
 biomolecular systems are increasingly relying 
 on the information obtained by generating dynamical trajectories of 
 large - scale detailed atomic models . to achieve computational efficiency , 
 these molecular dynamics ( md ) simulations rely on an approximate representation 
 of the born  oppenheimer potential energy surface based on simple 
 mathematical functions , the so - called molecular mechanical force fields . the vast majority of these force fields represent the electrostatic 
 potential energy as pairwise additive coulomb interactions between 
 effective fixed atomic charges that are empirically adjusted to reflect 
 the average polarization of molecular moieties in a realistic environment . examples of nonpolarizable additive potential functions that have 
 been widely used to study biological systems include the charmm , opls , 
 gromos , and amber force fields . despite the success of the classical additive force fields 
 in modeling 
 a wide range of critical chemical and physical events , the limitations 
 associated with a fixed charge model are frequently observed . this 
 became increasingly apparent , as additional efforts to optimize the 
 parameters of the force fields failed to improve the accuracy of the 
 models , especially for the binding energetics of multivalent ions 
 in complex electrostatic environments . by comparing with quantum mechanical ( qm ) calculations , ji et al . 
 pointed out the importance of many - body polarization effects on the 
 pka value of asp26 of thioredoxin , which 
 indicates the lack of nonadditivity in the classical additive force 
 fields . because induction is approximated 
 in a mean - field fashion , the nonpolarizable additive models are incapable 
 of describing situations where electrostatic polarization effects 
 are important . to address this issue , 
 substantial efforts have been made to introduce a more elaborate representation 
 of electrostatic interactions in the molecular mechanical force fields 
 as early as in the 1970s ; although , the 
 first complete polarizable force fields for biomolecular systems have 
 emerged only in the past few years . these 
 include the amoeba force field , in which the static and induced charge 
 distribution within a molecule is explicitly represented with multipoles 
 associated with the nuclei together with induced point dipoles determined 
 through a self - consistent field ( scf ) condition . applications of amoeba 
 include modeling of gas - phase clusters and simulations of condensed - phase 
 systems from a few to tens of nanoseconds . the principal computational 
 cost of this model lies in the scf conditions , limiting the type of 
 systems that can be simulated practically . another polarizable force 
 field that has emerged is the classical drude oscillator model . sometimes called the shell model or the spring - on - charge 
 model , the drude model consists of adding 
 an auxiliary charged particle attached to the nucleus via a harmonic 
 spring that is treated as an additional dynamical degree of freedom 
 ( for more details , see refs ( 13 , 14 , 18 , 22 , and 23 ) ) . despite its deceptive simplicity , 
 the classical drude oscillator model can effectively capture the charge 
 redistribution in response to changes in local electric fields , allowing 
 for a realistic representation of molecular interactions in a variety 
 of condensed - phase systems including solvated proteins , nucleic acids , 
 and membranes . importantly , the drude model allows for computationally efficient 
 md simulations by virtue of its particle - based functional form . in the present study , our efforts are 
 specifically focused on the 
 development of the drude polarizable force field toward the general 
 goal of accurately describing ion  protein interactions . the 
 selective binding of small ions ( na , k , ca , and cl ) often plays a decisive role 
 in the function of many enzymes . ions also play a critical role in the functions of numerous membrane - bound 
 proteins such as transporters and channels , where preferential binding 
 and ion - selective permeation serve to regulate ionic gradients and 
 action potentials across cell membranes . yet , despite their biological importance , it is well understood 
 that ion protein interactions are among the least accurately 
 represented by current additive force fields . to ensure the accuracy 
 and transferability of the final model , the parameters of the drude 
 polarizable force field for ion  protein interactions were optimized 
 on the basis of experimental data as well as a series of ab initio 
 calculations for a set of 30 representative metalloproteins for which 
 high resolution crystal structures are available . the article is organized 
 as follows : we begin by introducing the methodologies employed for 
 developing the drude polarizable force field , then report and conclude 
 on its performance and quality , and finally discuss perspectives and 
 outlooks for its usage and further development . the optimized 
 ion parameters for the drude polarizable force field are summarized 
 in table 1 . a rigorous 
 theoretical framework for modeling condensed matter must reproduce 
 a representative set of target data , and provide accurate descriptions 
 and prospective predictions of molecular mechanisms at both macroscopic 
 and microscopic scales . parameter optimizations for molecular mechanics 
 often suffer from the lack of available experimental data at the molecular 
 level and the large sampling gaps between detailed geometries and 
 bulk properties . to tackle these issues , 
 we implemented a strategy that combines experimental observables with 
 ab initio calculations to yield models that have a consistent performance 
 for various testing systems . at the microscopic level , where geometric 
 properties such as ion  protein coordination are essential , 
 ab initio optimized structures and interaction energies of ion - model 
 compound clusters are used as benchmarks . the variables used to define 
 the geometry of the clusters are shown in figures 1 , 2 , 3 , and 4 . at this stage of gas - phase cluster optimization , 
 potential parameter sets were selected with a multidimensional grid 
 search in an iterative fashion . nma 
 interactions , we started by considering a dimer structure , in which 
 three objective functions xe2 = ( eqm  emodel ) , xr2 = ( rqm  rmodel ) , and x2 = ( qm  model ) are minimized 
 simultaneously , where e is the k nma dimer interaction energy , r represents 
 the distance between the cation and the nma carbonyl oxygen atom , 
  is the intermolecular angle ( c = o  k in table 3 ) formed by the carbonyl group 
 and the cation . in the second , third , and fourth iterations , we searched 
 for the parameter set that minimizes the objective function xr2 = i=1n ( ri , qm ri , model ) , in which n indicates the 
 number of nma molecules in the cluster , and i represents 
 the index of each nma molecule . given the geometric properties , at 
 the macroscopic level , where thermodynamics properties are critical , 
 observables such as numerical solvation free energies of ions in condensed - phase 
 liquids ( sfe ) are strictly matched with experimental values . at this 
 stage , we evaluated the objective function xsfe2 = ( sfeqm  sfemodel ) by using the candidate 
 parameter sets obtained in the stage of gas - phase cluster optimization . the best parameter sets were then examined by applying them to compute 
 the ion the parameter sets with the 
 smallest standard deviations from the qm data were selected as the 
 final models . schematic illustrations 
 of ion - model compound clusters employed 
 in the gas - phase cation the dashed lines indicate the reference 
 properties ( listed in table 3 ) employed in 
 the target functions of the parameter optimization process . schematic illustrations of ion - model compound clusters employed 
 in the gas - phase cl the dashed lines indicate the reference ion - model compound distances 
 ( listed in table 3 ) employed in the target 
 functions of the parameter optimization process . schematic illustrations of ion - model compound clusters employed 
 in the gas - phase calculations of the cation  propanoic acid 
 ( pa ) clusters . the green and black dashed lines indicate the reference 
 ion  pa distances ( listed in table 5 ) employed in the target functions of the parameter optimization 
 process for the monovalent cation - pa clusters and the divalent cation schematic illustrations of ion - model compound 
 clusters employed 
 in the gas - phase calculations of the cation  ethanol ( etoh ) 
 clusters . the dashed lines indicate the reference ion - etoh distances 
 ( listed in table 4 ) employed in the target 
 functions of the parameter optimization process . a cation is represented 
 by a pink sphere . in the drude force field , 
 an auxiliary particle carrying a negative 
 charge qd any addition partial atomic 
 charges needed to represent the static charge distribution of the 
 molecule is simply superimposed to the charge of the nucleus , whose 
 single - point charge becomes qi  qd . the force constant 
 of the harmonic springs is typically kd = 1000 kcal / mol /  . the 
 value of  is determined by fitting the electrostatic field 
 from ab initio calculations in the presence of small perturbing charges . the polarizabilities of the ions ( na , k , ca , and cl ) experimentally determined in 
 aqueous solutions are used in this study . the interactions between 
 the ions and proteins are classified into three categories : ( 1 ) interactions 
 between cations and backbone carbonyl groups of asparagine ( asn ) , 
 glutamine ( gln ) acids , and model - compound n - methylacetamide 
 ( nma ) representing peptide linkage structures ; ( 2 ) interactions between 
 cations and hydroxyl groups of serine ( ser ) and threonine ( thr ) acids , 
 and model - compound ethanol ( etoh ) ; ( 3 ) interactions between cations 
 and carboxylate groups of glutamic ( glu ) and aspartic ( asp ) acids , 
 and model - compound propanoic acid ( pa ) . the underlying drude ion models 
 were published in our previous study . to modulate the interactions between the ions and proteins , we introduced 
 the following terms listed in table 1:(i)pair specific nonbonded 
 lennard - jones 
 ( lj ) interaction parameters emin and rmin ( referred to as nbfix in the charmm force 
 field ) , which override the default values determined by the lorentz we found that two pairs of nbfix parameters are 
 necessary when determining the interactions between cations and carbonyl 
 groups of protein backbones . they are the parameters for the cations 
 and carbonyl oxygen atoms , and for the cations and the lp particles 
 ( the lone pair particles in the force field files ) . since cl mostly interacts with the amide groups of protein backbones , the 
 force field contains only nbfix parameters for cl and hydrogen atoms of amide groups , nbfix and nbthole ( see below ) 
 parameters for cl and nitrogen atoms of amide groups.(ii)pair - specific thole charge 
 shielding 
 factors ( referred to as nbthole in the charmm force field ) , designed 
 to control the strength of electrostatic induction for a pair of charges . pair specific nonbonded 
 lennard - jones 
 ( lj ) interaction parameters emin and rmin ( referred to as nbfix in the charmm force 
 field ) , which override the default values determined by the lorentz we found that two pairs of nbfix parameters are 
 necessary when determining the interactions between cations and carbonyl 
 groups of protein backbones . they are the parameters for the cations 
 and carbonyl oxygen atoms , and for the cations and the lp particles 
 ( the lone pair particles in the force field files ) . since cl mostly interacts with the amide groups of protein backbones , the 
 force field contains only nbfix parameters for cl and hydrogen atoms of amide groups , nbfix and nbthole ( see below ) 
 parameters for cl and nitrogen atoms of amide groups . pair - specific thole charge 
 shielding 
 factors ( referred to as nbthole in the charmm force field ) , designed 
 to control the strength of electrostatic induction for a pair of charges . the absolute solvation 
 free energies of k , na , ca , and 
 cl in model - compound liquid nma and liquid etoh 
 were calculated by the following equation:1 the term 
 gelec arises from the electrostatic 
 interactions . the terms 
 grep and gdisp represent the repulsive and dispersive components 
 of van der waals ( vdw ) interactions , respectively . each of the three components was sampled in independent 
 md simulations in which a single ion is placed at the center of a 
 cubic box containing 144 nma or 125 etoh molecules . the contributions 
 from the electrostatic and vdw interactions were computed using the 
 thermodynamic integration ( ti ) approach . we implemented a soft - core 
 scheme for the calculation of grep to avoid an end - point singularity , together 
 with the weighted histogram analysis method ( wham ) . the term kbt ln(1/vm ) is an entropic contribution 
 accounting for compression from the ideal gas at 1 atm to an idealized 
 bulk solution at 1.0 molar concentration . the last term qf represents the contribution from the interface potential of a liquid 
 relative to that of vacuum , where  is the interface potential 
 and q is an ionic charge , and f is 
 the faraday constant = 23.06 kcal / mol / v . to estimate this term , we constructed slab simulations of a cubic 
 box containing nma or etoh liquid , which is continuous in the x and y dimensions . the z dimension is equal to three times the x  y dimension , so there are two neat - liquid / vacuum interfaces forming 
 along the z - axis . from the slab simulations , the 
 charge density of a liquid system is calculated as a function of the z - coordinate , and the interface potential  is determined 
 by numerically integrating the corresponding poisson equation . for na , the solvation 
 free energies were calculated by the charmm program , 
 version c38b1 . the drude force field 
 describing the model - compounds is based on recent work on polarizable 
 amino acids by lopes et al . periodic 
 boundary conditions ( pbc ) were applied and the system pressure was 
 maintained at 1 atm by using a modified anderson to active maximum computational efficiency , the system was propagated 
 with a dual - thermostat nos hoover ( nh ) extended lagrangian 
 algorithm . to closely reproduce the self - consistent 
 field ( scf ) induced polarization condition , a first nh thermostat 
 at 1 k is coupled to the relative distance of all pairs of reduced - mass 
 oscillators comprising the ( heavy ) nucleus and its ( light ) drude particle . 
 to maintain the global temperature of the system , a second nh thermostat 
 at 303.15 k is coupled to the centers - of - mass of all nucleus a time step of 1 fs was used , and all bonds 
 involving hydrogen atoms were constrained using the shake algorithm . long - range electrostatic interactions were evaluated 
 using the particle - mesh ewald method , with a smooth real - space cutoff 
 applied between 10 and 12  . to optimize the force field 
 parameters , 30 high - resolution crystal structures of enzymatic proteins 
 from the protein data bank ( pdb ) were chosen in which the binding 
 of na , k , and ca is functionally 
 important . the surrounding 
 bulk solution , with the tip3p water model , with 0.15 m of the corresponding 
 salt , was built with the charmm - gui web - portal . the information on the simulation systems , including the pdb identification 
 codes , system dimensions and numbers of atoms is given in table 2 . each system 
 is solvated with 
 tip3p water and has 0.15 m of salt . the program namd and 
 the charmm c36 force field were used to carry out md simulations . energy minimization by means of the conjugate 
 gradient algorithm was performed to relax the initial structures and 
 remove any high - energy contacts . each of the solvated proteins was 
 equilibrated for 4 ns using npt ensemble coupling with p = 1 atm . the system temperature was 
 kept at 300 k using langevin dynamics with a damping constant of 5 
 ps and a time step of 1 fs . the purpose of this 
 md simulation step is to sample different conformations of the nearest 
 protein ligands around the bound ions , so the numerical interaction 
 energies between the ions and the ligands via ab initio methods , the 
 drude polarizable force field , and the charmm c36 additive force field 
 possess thermal effects ( see supporting information 
 figure s1 ) . we found that most of the systems 
 have the ions stably bound in the binding sites , with a few exceptions , 
 which have the binding ions migrating out of the binding site during 
 the 4 ns md simulations . in those few cases , we applied a soft harmonic 
 potential with spring constant equal to 1.0 kcal / mol /  to the bound ions . this restraint allows a similar sampling of the 
 ligand conformations around the ions ( see figure 
 s2 ) . from the 4 ns conformations , we extracted 20 snapshots 
 to build a truncated model of the ion - binding sites . each ion - binding 
 site in the enzymatic proteins was truncated within a sphere of 5.5 
  from the ions , which is larger than the distance from the first 
 hydration shell to the ions (  3.0  ) . the truncated 
 ion binding pockets formed by amino - acid residues have c- and n - termini 
 defined in the charmm force field : ( 1 ) acetylated n - terminus ( ace ) 
 and n - methylamide c - terminus ( ct3 ) ; ( 2 ) n - terminus 
 ( nter ) and c - terminus ( cter ) . the two sets of termini allow testing 
 if there is a dependence of binding energies on the terminated groups . since they yield almost the same results , we discuss the results for 
 the first set of termini in the main text , and the other results are 
 provided in the si . when the truncated 
 procedure introduced several small and sparse segments of the original 
 undisrupted protein sequences , we added one or two original amino 
 acids to the segments in order to avoid disruptions , thus reducing 
 the numbers of end - amino - acid termini , even though the added amino 
 acids of the proteins were not found in the truncation . however , if 
 the disruptions required adding more than two original amino acids 
 between the sparse segments , we treated the sparse segments as separate 
 ligands . we also kept about two to five water molecules nearest to 
 the pockets of the binding sites . as a result , the truncated models 
 of the binding sites contain about 200  50 atoms , including 
 the binding - ions , two or five water molecules and the number of ligands 
 varies from 1 to 6 in different systems . all ab initio calculations using 
 basis set superposition error ( bsse ) corrections for the gas - phase ion - model compound clusters and ion - binding 
 proteins were carried out using the gaussian 09 suite of programs . another ab initio software , called general atomic 
 and molecular electronic structure system ( gamess ) , was also used to ensure the consistency between the two 
 commonly used softwares ( we tested on small complexes in si ) . for 
 qm calculations , it is essential to choose basis sets and types of 
 approximations , whose detailed analysis 
 are provided in the si from tables s17 
 and figures s36 , to compute ion - binding energies of metalloprotein 
 sites , which are defined as interaction energies between the cations 
 and surroundings using the counterpoise option in gaussian . we found 
 that the ion  protein binding energies computed by the b3lyp 
 functional with basis set cep-121 are consistent and reliable in terms 
 of accuracy and computational cost . nma 
 clusters were characterized at various coordination numbers from 1 
 to 4 as illustrated in figure 1 . nma 
 clusters , because the system of ca and one nma leads 
 to severely distorted electron distributions . the geometries of small 
 gas - phase clusters were optimized at the b3pw91 level with the lanl2dz 
 basis set . table 1 lists pair - specific 
 interactions between the ions 
 and protein atoms / particles : ( 1 ) interactions between k , na , ca , cl and protein 
 backbones tested in the model compound nma ; ( 2 ) interactions between 
 k , na , ca , cl and protein side - chains having hydroxyl groups tested in the model 
 compound etoh ; and ( 3 ) interactions between k , na , ca and protein side - chains having carboxylate 
 groups tested in the model compound pa . in the case of the ion - nma 
 interactions , we introduced extra pair - specific interactions between 
 the ions and the lone - pair particles attached to the carbonyl oxygen 
 atoms ( see methodologies ) . these interactions 
 effectively restrain the c  o - ion angles in the ion nma 
 dimers ( table 3 ) 
 to a reasonable range of about 180 , thus preventing the angular 
 geometry from being suddenly distorted , where c and o are for carbon 
 and oxygen atoms of carbonyl groups , respectively . reproducing the 
 stability of this angle is critical for modeling at molecular levels ; 
 otherwise any instability indicates imbalanced interactions between 
 the ions and the amide group , thus leading to incorrect coordination 
 of the ions . in the case of the cl amide 
 interactions , in addition to the interactions between cl and its immediate neighbors of amide hdp1a atoms , we found it to 
 be necessary to include a thole screening function between the charges 
 on the cl and nd2a2 atoms on the amide groups to 
 reproduce the qm energies more accurately ( see below ) . schematic representations 
 for 
 the geometric variables used in this table are shown in figures 1 , 2 , 3 , and 4 . the ab initio calculations are performed 
 at the b3pw91/lanl2dz level with the bsse corrections , using the gaussian 
 09 program suite . experimental values 
 of the binding enthalpies are taken from the study by klassen et al . o1 to o4 indicated the 
 oxygen atoms in the nma molecule , arranged by their distances to the 
 cation . the interaction energies in brackets are not employed as benchmarks 
 in the parameterization process . all other properties in the table 
 are employed as target values for parameterizing the drude model . schematic representations for 
 the geometric variables used in this table are shown in figures 1 , 2 , 3 , and 4 . the ab initio calculations are performed 
 at the b3pw91/lanl2dz level with the bsse corrections , using the gaussian 
 09 program suite . experimental values 
 of the binding enthalpies are taken from the study by klassen et al . the interaction energies in brackets are not 
 employed as benchmarks in the parameterization process . all other 
 properties in the table are employed as target values for parameterizing 
 the drude model . schematic 
 representations for 
 the geometric variables used in this table are shown in figures 1 , 2 , 3 and 3 . the ab initio calculations are performed 
 at the b3pw91/lanl2dz level with the bsse corrections , using the gaussian 
 09 program suite . experimental values 
 of the binding enthalpies are taken from the study by klassen et al . the interaction energies in brackets are not 
 employed as benchmarks in the parameterization process . all other 
 properties in the table are employed as target values for parameterizing 
 the drude model . to optimize the interactions 
 of the ion - nma clusters with different sizes ( figures 14 ) , we examined multiple properties 
 from the interaction energies to the coordination geometries that 
 are key features for predicting and reproducing realistic features 
 in biological modeling . table 3 lists the properties with both the ab initio calculations 
 of the monovalent ion  nma dimers and their experimentally measured 
 binding enthalpies . using a gaussian-3 ( g3 ) protocol with the geometry - corrected 
 counterpoise method , siu et al . and tsang et al . reported the binding enthalpies of 30.3 kcal / mol and 38.1 
 kcal / mol for k and na with single nma molecule , 
 respectively . using lanl2dz basis set and b3pw91 functional , we obtained 
 the binding enthalpies of 30.8 kcal / mol for the k nma dimer and 38.4 kcal / mol for the na - dimer , both within 0.5 kcal / mol of the g3 calculations . nma dimers differ by less than 
 3.0 kcal / mol from the first - principles calculations . although the 
 binding enthalpy for the na  nma dimer obtained 
 from the drude model is 2.8 kcal / mol larger than that from the ab 
 initio calculations , it deviates by only 0.1 kcal / mol from the value 
 measured in the collision - induced - disassociation experiments by klassen 
 et al . in terms of fundamental geometries , 
 the drude model successfully reproduces the following characteristics 
 of the qm calculations : using the drude 
 force field , both k and na have the c  o  ion 
 angles within 10 of the qm angles . a difference of about 
 5 is found for the n  h  cl angles . 
 in the case of ca , we used a reference structure having 
 a trimer that contains one cation and 2 nma molecules to mimic the 
 condensed phase , where ca is usually coordinated with 
 at least four water molecules or other electron - rich protein groups . the optimized distances , 2.48  between k and nma , and 2.15  between na and nma , 
 are 0.05  longer than the values by siu et al . , which converge to 2.43 and 2.10  , respectively . the 
 differences between the drude and ab initio optimized ion - nma distances 
 are less than 0.1  for all of the monovalent ions . using the 
 drude force field , ca has the optimized distances from 
 nma of about 0.12  lower than that from the ab initio calculations . 
 in this case , compensation on the nma  ion distance ( 0.12  
 in the ca cluster vs 0.10  in monovalent cation 
 nma cluster ) is necessary to achieve a balanced description of the 
 geometry and interaction energies . for this geometry , the charmm27 
 model , which was developed to reproduce the solvation free energies 
 of alkali cations in the condensed phase but , in the gas - phase , leads 
 to significant underestimations for the absolute binding enthalpies , 
 and thus overestimates the optimized distances between the ions and 
 amide groups . for example , charmm27 estimated an nma  k distance to be 2.68  in the dimer cluster , in comparison 
 with the qm predicted value of 2.48  . we found 
 that , as the number of nma molecules increases from 1 to 4 , the ion nma 
 distances in the drude models increase by 0.09  for k and 0.06  for na , reproducing well the values 
 found in the ab initio calculations . for ca , the distance 
 increases by 0.08  in the ab initio calculations and by 0.06 
  in the drude models . this phenomenon is of particular importance 
 because it reflects the behavioral correlations between ions and protein 
 structures , which dynamically change local electrostatic environments , 
 for example , in the processes of ion permeation through channel structures 
 such as the ion - channel kcsa protein . from the bulk water environment to the channel interior , the coordination 
 number of water molecules and the oxygen atoms of the negatively charged 
 carbonyl groups in the immediate solvation shell of a cation vary 
 noticeably as the cation enters deeper into the selectivity filter 
 of kcsa . the interactions and geometries between the ions and 
 other typical amino - acid side - chains , which are modeled as pa and 
 etoh molecules ( figures 3 and 4 ) , are listed in tables 4 and 5 . since the first coordination shell of the solvent molecules 
 contributes approximately half of the ion - electrochemical solvation 
 free energies , the size effects of the 
 coordination shell were examined by varying the number of etoh molecules 
 from 1 to 4 , and the number of pa molecules from 1 to 2 for the monovalent 
 cations , and 2 to 3 for ca . the numbers of pa molecules 
 are chosen based on the survey of the carboxylate groups within 5.5 
  of the cations in the 30 ion - binding sites ( see the sampling of ion - protein binding section ) . pa is 
 a natural carboxylic acid having one negative charge ( |e| ) , 
 whose electrostatic attractive forces with the cations can have large 
 contributions due to the polarization effects , since a 0.1  
 deviation of the ion pa distance can result in roughly 610 
 kcal / mol changes in the interaction energies for the dimer structures , 
 and as large as 30 kcal / mol for the ca  pa trimer . 
 to account for these contributions , we parametrized the ion pa 
 interactions in an iterative fashion by using the qm results as the 
 targets . table 4 summarizes an assessment 
 of the ion  pa interactions , which shows the drude force field 
 yielding an absolute difference of < 0.1  from the qm calculations 
 for all of the ion - oxygen distances , and the energies are different 
 from the qm ones by < 5 kcal / mol . etoh interaction energies , 
 which differ from the qm results by less than 2 kcal / mol ( < 7% ) 
 ( table 5 ) . the optimized ion - oxygen distances 
 show deviations of 0.02  from the qm optimized values , 
 while the optimized angles have such deviations of < 5. in 
 fact , these deviations are less than those among various ab initio 
 approximations ( table s1 ) . to examine the quality of 
 the drude force field in terms of thermodynamics , we computed the 
 solvation free energies and solvent  solvent transfer free energies 
 of natural salts in several liquids , to directly compare with the 
 experiments by yu et al . and klassen et al . in tables 6 and 7 , 
 we summarize the numerical and experimental free energies of the ions 
 in the liquids of nma and ethanol transferring from liquid water simulation 
 systems . the accurate prediction of the drude model for these thermodynamic 
 properties is remarkable because the model can be used to describe 
 changes of not only local electrostatic environments ( see above ) , 
 but also bulk - solvent environments having different electrostatic 
 environments and substantially different polarizable effects . the absolute 
 solvation free energies 
 of single ions include the contributions of the interfacial potentials 
 for ion translocation across the vacuum liquid  nma interface . the 
 solvation free energy 
 for nacl and kcl salts in all of the liquids are different from the 
 experimental values by only 0.3 kcal / mol . in contrast , the charmm 
 param27 nonpolarizable force field ( using the nbfix parameters ) yields 
 the free energies of 140.7 kcal / mol for kcl and 157.0 
 kcal / mol for nacl in nma , which are at least 10 kcal / mol larger than 
 the experimental values , thus indicating inaccurate interactions between 
 the ions and the condensed phase liquids . using the drude model , we 
 also estimated the solvation free - energy differences to be 17.7 , 
 18.6 , and 15.2 kcal / mol for the substitution of the 
 ions ( k  na ) in liquid water , nma 
 and etoh , respectively . these values are almost identical to the experimental 
 values reported in many experiments , whereas the free energy difference 
 for the substitution in liquid nma computed by the charmm c36 force 
 field is 2.3 kcal / mol larger than the experimental value . for cacl2 solvation in liquid nma , 
 the solvation free energy was experimentally estimated to be 526.7 
 kcal / mol in liquid water environment at 293 k , while the drude force 
 field accurately predicted this value to be 526.6 kcal / mol 
 in liquid nma . for the ion solvation 
 free energies in this liquid , the drude force field reproduces the 
 values of 145.5 kcal / mol for kcl and 160.5 kcal / mol 
 for nacl , which deviate by 0.5 kcal / mol from the experiments . for cacl2 in this liquid , we computed 
 the solvation free energy by considering the transfer free energy 
 of a salt between two saturated solutions using2where a and b refer to two types of solvents 
 with molarity a and b , a and b are the activity coefficients of 
 the salt in two solvents , and t is the temperature . based on experimental data , this equation can estimate the solvation 
 free energy of cacl2 in ethanol to be 524.9 kcal / mol . however , the drude model predicts this value to be 473.2 kcal / mol , 
 which is about 50 kcal / mol ( 10% ) different from the experimental 
 estimation . this difference , even though it is large , is commonly 
 encountered in qm calculations for the interactions between calcium 
 and proteins ( see figure s5(a ) ) . as an 
 approximation of electron densities , the drude model perhaps suffers 
 the same difficulty in describing the divalent ion interactions with 
 surroundings . as we will show later , the drude model yields an acceptable 
 accuracy ( uncertainty < 10% ) for a number of calcium - binding sites . table 7 shows the free energy differences 
 for transferring salts from one solvent to another . a free energy 
 of transferring between two salts can indicate the capability of a 
 force field to describe changes of electrostatic environment , which 
 can be quantified by dielectric permittivity . for example , at 293.15 
 k the dielectric permittivity is 80.1 for water , 37.8 for liquid nma , 
 only 24.5 for liquid ethanol , and in proteins the dielectric 
 permittivity ranges from 80 to as low as 4 . indeed , the drude polarizable force field can accurately 
 reflect the changes between the different electrostatic environments 
 with the solvent  solvent free energies less than 0.5 kcal / mol 
 away from the experimental values . for the final and critical 
 test , we examined how the drude force field performs for ion - binding 
 sites in 30 biologically important enzymes ( see methodologies ) . for protein atoms , table 8 shows the types 
 and numbers of ion - coordinating oxygen atoms in the various binding 
 sites found in these enzymes . those oxygen atoms include backbone 
 amide linkages , hydroxyl groups of ser and thr , carboxylate groups 
 of asp and glu , and carboxamide groups of asn and gln . table 8 also shows that k mainly interacts 
 with backbone carbonyl groups , na tends to have a preference 
 for both carbonyl groups and the carboxylate oxygen atoms , and ca displays favorable interactions with the carboxylate groups . snapshots of enzymes 1ni4 ( k ) , 1l5b ( na ) , and 3li3 ( ca ) are illustrated in 
 figure 5 . figures 6 and s79 show the binding energies computed 
 by using qm calculations , the drude force field , and the c36 force 
 field ( the numerical data are provided in tables 
 s89 ) . in order to illustrate how different force fields 
 influence the ion - protein binding , figure 7 shows the structure alignments of simulation snapshots obtained 
 from the two force fields for the binding sites having k , na , and ca . snapshots of toy - model 
 ( truncated ) ion - binding sites employed in 
 qm calculations . the backbone 
 structures of the protein segments are highlighted in orange with 
 the ribbon model . the red , 
 green , and blue bars with statistic errors represent calculations 
 using qm , c36 , and the drude model , respectively . ( b ) na binding proteins ( c ) ca binding proteins . structure alignments for proteins 1ni4 ( a ) , 193l ( b ) and 1a4v ( c ) with k , na and ca cations , respectively . the structures are the 
 x - ray structures ( colored by atom name ) and snapshots from md full 
 system simulations using the charmm c36 force field ( yellow licorice 
 stick models ) and the drude force field ( light purple licorice stick 
 models ) . the red vdw 
 sphere represents position of water molecules in the crystal structures . 
 in panels a and b , the rmsd values indicate the deviations between 
 the backbone carbonyl groups of the crystal structure and of the equilibrated 
 structures from 4 ns md calculations using the c36 and the drude force 
 field . in panel c , the rmsd indicates the deviation between the crystal 
 structure and the md equilibrated structure , by aligning the carbonyl 
 group oxygen as well as the side chain of residue asp87 and asp88 . each element ( types of oxygen 
 atoms , h2o , ion and amino acids components ) is counted 
 within a 3.5  sphere centered at the bound ion . consistent with the experimental solvation 
 energies of k in the neat liquids , the qm calculations 
 yield the binding affinities 
 for 10 binding sites in a range of 60 to 200 kcal / mol 
 with standard deviations of about 10 kcal / mol ( figure 6a ) . overall , the drude force field reproduces all 10 of the 
 affinities with the same uncertainties as in the qm calculations . the largest deviations between the qm and drude calculations are found 
 in enzymes 1v3z and 4ls7 , 
 in which the binding affinities have been underestimated and overestimated 
 by 12 kcal / mol , respectively . for comparison , the charmm36 force field 
 with or without the nbfix parameters for k systemically 
 overestimates the absolute binding affinities by 437 kcal / mol 
 ( see figures s79 in the si ) . this 
 overestimation and data in table 6 , where we 
 showed that the nonpolarizable force field underestimates the solvation 
 free energy of neutral salt kcl in nma by 12.6 kcal / mol , suggest that 
 the c36 force field is unable to reflect the highly polarizable nature 
 of the cl anion and the interactions between k and nma . in particular , to quantify the consistency between 
 the developed models and the qm benchmarks , we computed a concordance 
 correlation coefficient c between the qm calculated 
 binding energies ( eiqm ) and the energies evaluated using 
 the c36 and drude models ( eimodel ) , for k - binding protein 1jf8 ( see figure 8a)3where  and  are 
 the variance and the average of the calculated energies , respectively . the parameter r , which is the correlation coefficient 
 between two sets of energies , is evaluated by4we found that the drude model for 
 k not only yields binding energies close to those of the 
 qm calculations , 
 but also improves the consistency with qm , from c = 0.012 with the c36 force field to 0.317 with the drude force field . the improved accuracy is also reflected in the ability of the force 
 field to maintain the crystallographic structures . for k binding in protein 1ni4 ( figure 7a ) , both the drude and the c36 model 
 have carbonyl group oxygen root mean square deviation ( rmsd ) from 
 the crystal structure of 0.41  relative to the crystal structure 
 ( obtained by aligning the backbone carbonyl group oxygen atom of the 
 ion - binding site ) . a closer look to the binding site reveals that 
 the deviations are primarily associated with the backbone carbonyl 
 groups and a few water molecules . in the c36 and the drude simulations , the k ion is in direct contact with one water molecule 
 and the backbones of ile161 , arg163 , asn165 , and ile112 . a subtle 
 difference is observed in the orientation of the pro110 backbone . 
 in the md simulations based on the drude force field , the pro110 carbonyl 
 group points toward the interior of the binding site in agreement 
 with the crystal structure , while in the simulation based on c36 this 
 carbonyl group flipped and points toward the exterior of the binding 
 site . the structural difference decreases the binding energies from 67 
  4 kcal / mol in the drude simulation to 83  11 
 kcal / mol in the c36 simulation . correlation of binding energies in ( a ) 
 k - binding protein 1jf8 , ( b ) na - binding protein 1s36 , and ( c ) ca - binding protein 1rwy . in each system , we used 20 snapshots 
 from 4 ns md simulations for computing the energies ( see methodologies ) via qm ( horizontal axis ) , and c36 
 and the drude force field ( vertical axis ) . the green dots and blue 
 squares represent the c36 and drude energies versus qm energies , respectively . the solid red line represents the agreement between molecular mechanics 
 and qm binding energies . in panel c , the black line in c presents 
 the qm binding energies ; the red line presents the reference values 
 employed in the parametrization process , which has been increased 
 by 10 kcal / mol with respect to the black line . na displays stronger binding energetics 
 than k , as reflected by the higher absolute na - binding affinities of 60 to 300 kcal / mol in the 
 qm calculations ( figure 6b ) . the largest affinity 
 appears in the case of 1qus , which has one hydroxyl oxygen atom , one 
 carbonyl oxygen atom and nearly 4 carboxylate oxygen atoms on average 
 within a cutoff of 3.0  . overall , the drude model shows absolute 
 differences of the affinities from 1 to 27 kcal / mol with respect to 
 the qm calculations . in contrast , the c36 force field systematically overestimates 
 the binding affinities , especially for 1qnj and 1qus , where the binding affinities are overestimated 
 by > 50 kcal / mol , with two exceptions , where the same affinity as 
 computed 
 from qm calculations was found for 193l and an underestimate by 6 kcal / mol with 
 respect to the qm energy was found for 3n0u . in particular , we found that the concordance 
 correlation coefficient between the c36/drude binding energies and 
 the qm energies for na - binding protein 1s36 ( see figure 8b ) differs significantly with c = 0.365 of the drude model from c = 0.027 of the 
 c36 force field . although both force fields yield the same binding 
 energies of 64 kcal / mol for na binding in protein 
 193l , the packing of the ser72 side chain differs in the drude and 
 the c36 simulations . as shown in figure 7b , 
 the c36 simulation predicts a more compact structure , which deviates 
 further from the crystal structure than that from the drude simulation . the ser72 hydroxyl group is 0.2  closer to the cation than in 
 the drude simulation . in terms of the packing of the backbone carbonyl 
 group , the drude simulation yields an rmsd of 0.49  relative 
 to the crystal structure . however , the deviation between the c36 structure 
 and the crystal structure is as large as 0.86  . among the three cations examined in the present 
 study ( na = 0.157  ; k = 0.830  ; ca = 0.490  ) , ca has the largest net charge and therefore the strongest polarizing 
 power . the ab initio estimates ( figure 6c ) 
 indicate the binding affinities with an average of approximately 700 
 kcal / mol and a broad span of nearly 350 kcal / mol from the most negative 
 value ( 2aaa , 
 973 kcal / mol ) to the least negative one ( 1rwy , 629 kcal / mol ) . 
 overall , although the drude force field seems to underestimate the 
 binding affinities by about 1040 kcal / mol for all 10 of the 
 binding sites , it yields the binding affinities more accurately than 
 the c36 additive force field , which overestimates the binding affinities 
 by more than 80 kcal / mol compared to the qm calculations ( figure 9 ) . in particular , for ca - binding protein 1rwy , the c36 model yields 
 c = 0.003 , which is much lower 
 than c = 0.578 computed from the drude force field . 
 in the three cases ( 1bli , 1exr , and 1rwy ) , the c36 force 
 field misrepresents the binding energies over 170200 
 kcal / mol . interestingly , we found that in the three proteins , ca interacts with more than five oxygen atoms of the carboxylate 
 groups ( see figure 10 ) . the overestimated binding 
 energies with the c36 force field are also found to increase almost 
 linearly as a function of the number of carboxylate oxygen atoms in 
 the crystal structures ( see table 8) . consequently , 
 this indicates significant overbinding of ca as well 
 as potentially distorted structures and inaccurate coordination numbers 
 by the c36 force field due to too strong attraction between ca and electron - rich groups . however , the drude force field 
 yields the binding affinities for ca in agreement with 
 the ab initio calculations , especially for 1a4v and 4kts with the differences being only 8 and 
 9 note that b3lyp approximation used for the 
 qm calculations can overestimate the absolute binding free energies 
 for ca by about 10 kcal / mol compared with the mp2 method ; 
 and the cep-121 basis set employed in the benchmark calculation can 
 yield an overestimate of 1020 kcal / mol compared with the full - electron 
 calculations using the 6 - 31++g * * basis set ( see table s7 and figures s36 ) . these large uncertainties 
 even in the qm calculations imply that such a simpler approximation 
 as the drude model for describing approximate electron densities around 
 the divalent ion can not be more accurate than 1020 kcal / mol . therefore , an uncertainty of 1040 kcal / mol ( < 6% of the 
 averaged binding affinity ) for all 10 of the ca - binding 
 sites can be considered as acceptable . consistent with the clear differences 
 between the c36 and drude / qm binding energies , we observed several 
 distinguishable features of ca binding in 1a4v during the 4 ns 
 full - system simulation using the drude force field and the c36 additive 
 force field . in figure 7c , we aligned the backbone 
 carbonyl group oxygen atoms as well as the side chain of asp87 and 
 asp88 of the equilibrated structure with the crystal structure . in 
 the drude simulation , the side chain of residue asp87 and asp88 both 
 have one oxygen atom in direct contact with the ca cation , 
 consistent with the crystal structure . however , we found that the equilibrated c36 structure 
 deviates from the crystal structure by an rmsd of more than 1.0  . 
 the backbone of residue asp82 , and the side chain of residue asp87 
 and asp88 pack tightly around the divalent cation . interestingly , 
 the latter two residues are bound to ca through both 
 of the carboxylate oxygen atoms , showing a completely different coordination 
 from the crystal structure , and thus explaining why the c36 model 
 significantly overestimates the ca binding affinity in 1a4v by 164 kcal / mol . difference 
 between qm ion - binding energy and classical ion - binding 
 energy versus number of carboxylate ( co ) oxygen atoms in 30 ion - binding 
 sites . illustration of ligands for the analysis of the nonadditive effects 
 in ca - binding protein 1bli . one of the key 
 features of a polarizable 
 force field is its ability to account for qm nonadditive effects , 
 which follow from the observation that the sum over individual interaction 
 energies of each ion - ligand pairs ( esum ) does not match the total interaction energy of the ion for the 
 complete binding site ( eoverall ) . hence , nonadditivity is simply quantified as the difference ( e ) between esum and eoverall . by default , nonadditivity is absent 
 with a nonpolarizable additive force field . to evaluate this feature 
 of the models , we compared the k - nma cluster interaction 
 energies computed by qm , the charmm c36 nonpolarizable force field , 
 the drude force field , and the amoeba polarizable force field . table 9 shows that as 
 the number of nma molecules increases from 1 to 4 , the magnitude of 
 e = esum  eoverall accordingly 
 increases from 0 to 14.6 kcal / mol , where esum is a sum of interaction energies for all k ion - nma pairs , and eoverall is 
 the interaction energy of the overall cluster . the deviations 
 between the drude and qm results are + 0.1 , 1.8 , and 3.3 
 kcal / mol for the systems of 2 , 3 , and 4 nma molecules , respectively . these deviations are smaller than those in the amoeba force field , 
 which overestimates nonadditive effects by + 2.7 , + 3.9 , and + 5.5 kcal / mol 
 for the systems of 2 , 3 , and 4 nma molecules , respectively . esum is a sum of interaction energies for all 
 individual k - ion  nma pairs . the 0.1 kcal difference is due to round off error . to further assess nonadditive effects 
 produced by the new drude 
 force field , we applied the same analysis to some ion - binding sites , 
 which have 46 individual ligands . each ion - binding site was 
 then dissected into smaller fragments for computing the binding interaction 
 energies between a cation and all individual fragments . one fragment 
 can be a peptide , or a single water molecule , or another binding ion 
 as illustrated in figure 10 ( specifically , 
 the binding site of 1e43 has two ca ions , the binding site of 1bli has another 
 na , and the binding sites of 1gen and 1s36 have a chloride anion ( see table 8) ) . tables 10 , 11 , and 12 report the nonadditive effects 
 in 1jf8 , 1ljl , 1e43 , 1l5b , 1bli , and 3l3i . for example , in 
 the k - binding enzyme 1jf8 , the qm binding energy eoverall is 196 kcal / mol and esum is 239 kcal / mol , thus rendering 
 a value of the nonadditivity , e = 43 
 kcal / mol , which accounts for 18% of esum . remarkably , nonadditive effects with the 
 drude force field for this binding site differs by only 6 kcal / mol 
 from the qm nonadditivity value . similar consistency of the nonadditive 
 effects is found for other monovalent ion - binding sites . nonadditive 
 effects become more pronounced in the ca binding proteins , 
 where they contribute approximately 2030% of the accumulated 
 binding energies in both the qm and drude calculations . the magnitude 
 of nonadditive effects explains the difficulties in approximating 
 electron densities around a divalent ion on the basis of effective 
 fixed - charge models . the ability of the drude model to successfully 
 reproduce the many - body electrostatic nonadditive effects present 
 in the qm is a real breakthrough in the efforts to accurately represent 
 ion  protein interactions . esum is the sum 
 of interaction energies for all individual ion  ligand , ion eoverall is the overall qm interaction energy by treating all the ligands , 
 water molecules and other ions as a whole system . all values in kcal / mol are averaged for 20 frames extracted from 
 4 ns simulations . all conclusions derived 
 from md simulations of cation interactions 
 with proteins are dependent on the accuracy of the underlying potential 
 functions . the conformational sampling , selection of low energy states 
 and therefore the time - evolution of solvated metalloproteins are essentially 
 limited by the accuracy of the potential energy surfaces and by statistical 
 availability . new simulation techniques , software and hardware solutions 
 are pushing time - scales into the micro- and even milli - second range , 
 but one must seek a balance between more accurate but more expensive 
 representations of the energy surfaces , and the obligation to adequately 
 sample the relevant configurations . we believe that the present study 
 clarifies the limitations of existing additive force fields in studies 
 of metalloproteins and exposes possible sources of systematic error . according to the energy mapping of 30 protein sites , there is an extremely 
 sharp variation in accuracy between an additive force - field and ab 
 initio data with an increase in ( a ) cation and/or site net charge , 
 ( b ) complexity of the electrostatic environment and chemical types 
 of coordinating ligands . interestingly , the same can be concluded 
 from comparative studies of simple liquid mimetic systems such as 
 ethanol and n - methylacetamide as reported in this 
 and previous studies . this can be attributed 
 in part to the nonadditive effects , but another source of error is 
 an assumption on the transferability of ion parameters developed for 
 reproducing the solvation in water to simulations of metalloprotein 
 sites . the findings reported in this paper warrant the revisitation 
 and validation of some previous conclusions about the transport mechanisms 
 or selective binding to metalloprotein sites . we have refined the drude polarizable force field to accurately 
 represent ion  protein interactions . pairwise - specific nonbonded 
 parameters and charge - screening functions were systematically optimized 
 to reproduce a representative set of target data comprising ab initio 
 calculations and experimental thermodynamic properties . the resulting 
 force field accurately models the interactions of k , na , ca , and cl with the different 
 polar groups present in all proteins , including the backbone amide 
 linkage , and the amino acid side chains that contain hydroxyl , carboxylate 
 or carboxamide groups . by carefully calibrating the force field parameters 
 for a set of model systems on the basis of geometrical and thermodynamic 
 properties and ab initio calculations , the resulting force field is 
 likely to remain accurate for a wide range of situations . protein models successfully 
 reproduce consistent energetics in 30 different ion - bound enzymatic 
 proteins in comparison with ab initio calculations , 
 whereas the additive models such as the c36 additive force field can 
 fail . the accuracy is clearly improved in the case where the many - body 
 polarization effect is pronounced , for example , in the ca binding proteins . remarkably , we showed that the many - body nonadditive 
 effects found in qm are successfully accounted for in the drude model 
 for more accurate descriptions of the ion protein interactions 
 than the nonadditive force fields , and the drude model is substantially 
 less computationally expensive than qm calculations for simulating 
 a variety of biological systems . since the ion  protein 
 interactions are significantly improved , 
 the simulation approaches using the drude force field will be viable , 
 inexpensive , and robust . the force field maintains the format of the 
 charmm force field , so it can be used in many simulation softwares 
 such as namd , amber , etc . , while effectively modulating the pair - specific 
 interactions and describing well the dynamical changes of electron 
 densities . our study leaves little doubt that even the most carefully 
 parametrized additive force - field will meet with increasing difficulties 
 in reproducing correctly the ion  site interaction energetics 
 for diverse electrostatic environments . the nonadditive effects in 
 the description of divalent ion / protein site interactions lead to 
 up to 15% adjustment in interaction energetics compared to that of 
 well - parametrized force fields such as charmm c36 . in the situation 
 where the electrostatic environment has to dynamically change in response 
 to ion binding , these effects become crucial to the understanding 
 of metalloprotein function . systematically assessing and enhancing 
 the accuracy , we believe 
 that the present study will facilitate more advanced improvements 
 and innovations in the fields of computational biophysics and biochemistry 
 for studies of ion channels , transporters , enzymes , and so on . this 
 study also outlines a novel protocol for comprehensively constructing 
 and evaluating a force field for modeling macromolecules such as a 
 force field for qm / mm scc - dftb embedded in charmm . the protocol incorporates 
 training of parameter sets on a large number of biological systems , 
 whose properties are diverse . in that protocol , additional parameters 
 can be introduced based on specific physical properties , either by 
 utilizing customized parameter optimization tools or with readily 
 available gadgets , for example , the general automated atomic model 
 parameterization ( gaamp ) web portal ( url : http://gaamp.lcrc.anl.gov/ ) .
small metal ions play critical roles 
 in numerous biological processes . 
 of particular interest 
 is how metalloenzymes are allosterically regulated 
 by the binding of specific ions . 
 understanding how ion binding affects 
 these biological processes requires atomic models that accurately 
 treat the microscopic interactions with the protein ligands . 
 theoretical 
 approaches at different levels of sophistication can contribute to 
 a deeper understanding of these systems , although computational models 
 must strike a balance between accuracy and efficiency in order to 
 enable long molecular dynamics simulations . in this study , we present 
 a systematic effort to optimize the parameters of a polarizable force 
 field based on classical drude oscillators to accurately represent 
 the interactions between ions ( k+ , na+ , ca2 + , and cl ) and coordinating amino - acid 
 residues for a set of 30 biologically important proteins . by combining 
 ab initio calculations and experimental thermodynamic data , 
 we derive 
 a polarizable force field that is consistent with a wide range of 
 properties , including the geometries and interaction energies of gas - phase 
 ion / protein - like model compound clusters , and the experimental solvation 
 free - energies of the cations in liquids . 
 the resulting models display 
 significant improvements relative to the fixed - atomic - charge additive 
 charmm c36 force field , particularly in their ability to reproduce 
 the many - body electrostatic nonadditivity effects estimated from ab 
 initio calculations . 
 the analysis clarifies the fundamental limitations 
 of the pairwise additivity assumption inherent in classical fixed - charge 
 force fields , and shows its dramatic failures in the case of ca2 + binding sites . 
 these optimized polarizable models , amenable 
 to computationally efficient large - scale md simulations , set a firm 
 foundation and offer a powerful avenue to study the roles of the ions 
 in soluble and membrane transport proteins .