electrostatic interactions 
 comprise one of the principle interatomic 
 forces , along with exchange - repulsion , dispersion , and polarization 
 or induction .
the importance of electrostatic interactions is paramount 
 at long range and for polar molecules .
much development effort has 
 been focused on computational treatment of long - range electrostatics , 
 e.g. , the development of particle - meshed ewald ( pme ) methods .
electrostatic interactions at short range have received less consideration 
 until recently . at close distances ,
a spherical approximation of atomic 
 charge distributions is insufficiently accurate and use of atomic 
 multipole expansions provides much greater flexibility in modeling 
 complex electrostatic potentials near a molecular surface , an insight 
 which inspired the development of the amoeba force field .
nonetheless , at very close interatomic 
 distances , when electron clouds overlap , a point multipole approximation 
 becomes inadequate .
the electrostatic potential within a spherical 
 electron cloud no longer behaves as a simple 1/r interaction potential at small separation distances .
such 
 deviation from a simple coulomb potential is referred to as a penetration 
 effect . while the charge penetration effect leads to a negative correction 
 to energy at typical molecular interaction distances , where the electron  electron 
 penetration is dominant
, it can be repulsive at very short range . a recent study by lewis and co - workers reported 
 the counterintuitive result that any ring substitutions of the benzene 
 dimer ( parallel ) with electron - withdrawing or electron - donating groups 
 yield more favorable electrostatic contributions than the unsubstituted 
 benzene  benzene dimer itself .
this 
 result is contrary to the conventional thought that such interactions 
 are correlated with the ability to withdraw or donate electrons to 
 the  cloud as described by the hunter  sanders rules .
sherrill and co - workers suggested this is because 
 the electrostatic interactions in such systems at the  
 stacking distance exhibit a significant charge penetration effect . 

the multipole model , which the hunter  sanders rules are based 
 upon , can not correctly account for such effects .
moreover , in a recent study of aromatic crystals , a charge 
 penetration corrected amoeba - like model predicted better crystal properties 
 than the uncorrected model .
it was shown 
 that point atomic multipoles consistently predict positive ( repulsive ) 
 electrostatic interactions between stacked or t - shaped benzene dimers 
 while symmetry - adapted perturbation theory ( sapt ) suggests an opposite trend toward attractive interactions . 

the current amoeba force field seemingly compensates for penetration 
 with a less repulsive van der waals interaction , so the total interaction 
 energy is reasonable at certain dimer configurations .
however , explicit 
 incorporation of the penetration effect provides much better anisotropy 
 in crystal packing and makes the overall force field more transferable . 
 in another study of organochlorine compounds using the amoeba model
, 
 it was found that the transferability of chlorine van der waals parameters 
 was unsatisfactory , likely due to lack of an explicit penetration 
 correction .
there have been previous 
 attempts to incorporate the charge penetration 
 effect into implicit solvent models , multipole - based electrostatic 
 models , charge - density - based ( including gaussian multipole ) models , and combined quantum and molecular mechanics ( qm / mm ) 
 models . generally , the charge penetration correction involves 
 breaking the atom - centered point charge into an effective core and 
 a valence electron density , as suggested by gordon et al . and piquemal et al .
in this way , the electrostatic interaction between two atoms is 
 described as a sum of interactions between core and valence charge 
 densities , which can be modeled with empirical exponential functions . 

alternatively , rigorous integration over the two charge densities 
 can be used to model short - range electrostatic interactions , with a significantly greater 
 expenditure of computational effort .
others have explored incorporating 
 charge penetration effects into the qm / mm framework , using either 
 screened molecular mechanics ( mm ) charges or simple empirical damping corrections . such screened mm charges are typically parametrized 
 for qm / mm applications and may be not directly applicable in full 
 mm calculations , e.g. , to reproduce the attractive sapt electrostatic 
 energy in a stacked benzene  benzene conformation , unless an 
 explicit term for interactions of the valence charge densities is 
 included , as in the model recently proposed by wang and truhlar . in this study , the charge 
is revisited , 
 implemented , and extensively tested in the context of the amoeba force 
 field and using a new parametrization strategy .
the charge penetration 
 corrected amoeba point multipole model ( multipoles + cp ) is developed 
 using a comprehensive set of small molecule complexes , and the parameters 
 are determined for h , c , n , o , p , s , f , cl , and br to cover the elements 
 commonly found in organic and biological molecules . to facilitate 
 model development in this and future studies , a new database of sapt2 + decomposed quantum mechanical energies constructed for 101 small 
 molecule pairs , each at 7 different intermolecular distances ( the 
 s1017 database ) , is presented .
in order to systematically 
 examine the electrostatic and other components of intermolecular forces , 
 the s101 and s1017 databases of homo- and heterodimers of common 
 organic molecules have been constructed .
the first 
 66 pairs , which cover the majority of the typical organic interactions 
 of h , c , n , and o atoms , are taken from the s66 database from hobza 
 et al .
in addition , 15 complexes containing 
 halogen atoms ( f , cl , and br ) , six complexes containing sulfur , and 
 four complexes containing phosphorus have been added .
furthermore , 
 10 monomer  water complexes , which encompass amino acid side 
 chain analogs ( including the charged ones ) missing in the s66 data 
 set , have also been added , yielding a total of 101 pairs . to construct 
 the s1017 database , definitions of the intermolecular distance 
 vectors from the s668 database of hobza et al . were used .
unlike s668 , each of the 101 model complexes 
 were placed at seven separation distances , corresponding to 0.70 , 
 0.80 , 0.90 , 0.95 , 1.00 , 1.05 , and 1.10 times the equilibrium intermolecular 
 distances .
compared against the s668 database , the s1017 
 set includes more dimer configurations at very short separations , 
 which have been rarely investigated but are essential to the study 
 of penetration effects and exchange  repulsion interactions . 

we have selected 0.7 times the equilibrium distance as the lower bound 
 because the sapt calculations below show that at this close distance 
 the electrostatic energy is about 50% of the exchange 
repulsion 
 energy or higher ; i.e. , both electrostatic and van der waals ( vdw ) 
 components are important in the total interaction energy . as these 
 short distances are being sampled in molecular dynamics simulations 
 at room temperature and
even more so at higher temperatures , their 
 contributions to the simulated bulk properties are nonnegligible . 

thus , it is essential to ensure the charge penetration model behaves 
 correctly at these short distances . the newly added structures 
 among the 101 complexes were optimized at the mp2/cc - pvtz level of 
 theory with counterpoise correction using the gaussian09 program . for each of the resulting 707 dimer configurations
, 
 the interaction energy has been decomposed using sapt2 + analysis provided by the psi4 program .
the sapt2 + 
 calculation returns electrostatic , exchange  repulsion , induction , 
 and dispersion energies , all to second order with respect to intramolecular 
 electron correlation .
exact definitions of each component can be found 
 in figure 1 of sherrill et al .
it should 
 be noted that dispersion energy can only be separated from other effects 
 in long range when two molecules do not overlap .
thus , at van der 
 waals distances , it may be more appropriately to refer to this as 
  dispersion - like  or  damped dispersion  
 energy .
this should be kept in mind even though for simplicity the 
 term  dispersion  is used throughout the discussion . 

all sapt calculations were carried out using dunning s correlation 
 consistent basis sets at both aug - cc - pvdz and aug - cc - pvtz 
 levels .
the complete basis set ( cbs ) limits of the sapt2 + energies 
 were also estimated .
( data can be found in the supporting information ) schematic view of monomers and dimers 
 in the s101 data set .
different 
 configurations of the same dimers , e.g. , menh2water , 
 phenol  water , benzene  benzene , and mecl  mecl , 
 are included in the data set to take into account the orientational 
 effect .
a two - point extrapolation 
 strategy has been used to estimate the complete basis set limit of 
 the exchange  repulsion and dispersion energy at the sapt2 + 
 level of theory .
this is similar to helgaker s scheme but with an optimized p value 
 ( eq 1 ) .
such a protocol has previously been 
 applied to extrapolate the dispersion energy of dft - sapt calculation 
 in earlier study.1 different p values , 
 3.0 for exchange  repulsion and 4.3 for dispersion energy , were 
 obtained using the small pairs and subsequently applied for extrapolation 
 over the full s1017 database .
since the 
 truncated terms in the sapt2 + dispersion energy make a considerable 
 contribution to the total interaction energy , dispersion energies 
 obtained at the sapt2+/cbs level are scaled by a factor f in order to match the sapt total interaction energy to those obtained 
 at the ccsd(t)/cbs level of theory ( eq 2).2 by 
 minimizing eq 2 using the 66 pairs in the s66 
 data set , a scale factor of f = 0.89 has been determined 
 and used to construct the s1017 
 database . in order 
 to model the charge penetration effect , the method of piquemal et 
 al . is revisited .
their original model 
 corrects the charge  charge and charge  dipole interactions . 
 here , we propose to retain the charge 
, each atomic point charge is divided into an effective 
 core and a damped valence electron distribution .
thus , the electrostatic 
 energy between two atomic charges can be written as3where r is the interatomic 
 distance ; z is the positive effective core charge , 
 which is set to be equal to the number of valence electrons of each 
 atom ; q is the net charge of the atom , thus ( z  q ) can be considered as the 
 magnitude of the ( negatively charged ) electron cloud ; and  
 and  are two parameters controlling the magnitude of the damping 
 of the electron cloud when the atom is interacting with the core and 
 with electrons from other atoms , respectively .
thus , the total electrostatic 
 energy between two atoms now involves three components , the core 
the first method involves fitting the damped potential to 
 the qm electrostatic potential at short range , near or within the 
 molecular surface . by considering a probe charge of + 1e as a particle with an effective core charge of + 1e and having no valence electrons , ( z2 
thus , the 
 electrostatic potential can be written as4 once z and q are determined , 
  can be obtained easily by fitting eq 4 to the qm electrostatic potential . in the second method ,  

is intuitively set to be the same 
 as the number of valence electrons ( except the hydrogen atom):5 when z and 
  are fixed in eq 3 , the electrostatic 
 energy is more attractive when  
 is greater .
this is in accordance with the intuition that atoms having 
 a larger electron cloud may exhibit a stronger penetration effect . 
 although the final parameters for h , c , n , o , p , s , f , cl , and br 
 were derived based on the second method , the performance of both methods 
 is examined for h , c , n , and o containing molecules in later sections . as the distance between two atoms increases ,
thus , the electrostatic interaction 
 at medium and long distances can still be accurately modeled via a 
 multipole expansion , as the penetration correction diminishes rapidly 
 with distance . as the data will show , the penetration correction is 
 only significant when atomic separation is shorter than the sum of 
 atomic van der waals radii and thus does not affect the reciprocal 
 space portion of an ewald summation approach such as particle mesh 
 ewald ( pme ) .
in addition , to ensure the continuity between the real 
 and reciprocal space , a switching function is used near the real space 
 ewald cutoff distance ( typically 7  for atomic multipole pme ) 
 to ensure the penetration correction completely disappears:6where r is the interatomic 
 distance and rl and ru are the lower and upper bounds of the switching function .
quadrupole 
 interactions . following previously detailed procedures , an initial set of atomic multipole 
 moments for each molecule was obtained from distributed multipole 
 analysis ( dma ) at the mp2/6 - 311 g * * level 
 of theory .
then the dipole and quadrupole moments were further optimized 
 by fitting to the electrostatic potential calculated at the mp2/aug - cc - pvtz 
 level .
the same strategy has been used 
 in developing the amoeba force field for small molecules , proteins , and organochloroine 
 compounds previously . with the addition 
 of the penetration correction described here , monopole 
monopole 
 ( charge  charge ) interactions are calculated using eq 3 while all other terms in the amoeba retain their 
 original form .
the sapt and amoeba multipole based intermolecular 
 interaction energies were compared on exactly the same dimer structures . since h , c , 
 n , and o are the most common elements in organics and biomolecules , 
 their parameters were determined first within the new charge penetration 
 formalism .
a training set of 357 molecule pairs was used , consisting 
 of 157 pairs of hydrogen - bonded complexes , 127 pairs of 
 dispersion - dominant complexes , and 87 pairs with mixed features 
 of both .
for example , an initial set of parameters for sp or sp carbon , nonpolar hydrogen , and hydrogen 
 attached to sp c were obtained by selecting a smaller 
 number of pairs ( e.g. , 87 ) from the dispersion - dominant and 
 mixed complexes .
similarly , other parameters such as those for sp oxygen and polar hydrogen were obtained initially from water 
 dimers .
then , one or more sets of initial parameters 
 for all h , c , n , and o atom types were optimized together using the 
 entire training set . once the parameters for these four elements were 
 finalized , further parametrization for p , s , f , cl , and br was carried 
 out using the subsets in s1017 . in each of these subsets , an 
 80/20 ratio for the training and testing complexes
was maintained 
 to ensure a sufficient amount of data points for each atom type .
an 
 optimization program written in python , using the quasi - newton and 
 nelder  mead simplex methods from the scipy library , was applied 
 to all of the parametrization work .
the first derivative of the sum 
 of unsigned errors with respect to each parameter was calculated numerically .
in order to examine the convergence of the sapt2 + energy toward the 
 basis set limit , five small molecule pairs were selected and calculated 
 using the aug - cc - pvxz ( x = d , t , q , or 5 , abbreviated as axz in the 
 following paragraphs ) basis sets .
ethyne 
 ( t - shaped ) , and ethyne  water ( cho ) .
the energy 
 difference between different axz ( x = d , t , q , and 5 ) basis sets for 
 the total interaction energy and each energy component , including 
 electrostatic , induction , exchange  repulsion , and dispersion 
 energy , are compared . in general , a steady decrease in the energy 
 gaps between adz 
the difference between 
 aqz and a5z basis sets of all of the energy components as well as 
 the total interaction energy are already well below 0.05 kcal / mol . 

in particular , the differences in electrostatic and induction energies 
 are even smaller , at 0.003 and 0.005 kcal / mol , respectively .
this 
 implies that the difference between a5z and a bigger basis set , e.g. , 
 a6z , should be even smaller and negligible .
therefore , the results 
 obtained using the a5z basis set were used to approximate the complete 
 basis set limit .
since sapt2 + calculations 
 are computationally expensive , the practical 
 size of the basis set has been limited to aug - cc - pvtz for most molecule 
 pairs in the s101 database . to obtain an estimate of the sapt2 + energy 
 at the cbs limit , extrapolation
as shown in figure 2 , the electrostatic and induction energy components 
 converge quickly to the cbs limit ( approximated by a5z results ) .
the 
 mean unsigned errors between atz and a5z of five pairs are 0.018 and 
 0.010 kcal / mol for the two components , respectively .
therefore , for 
 electrostatic and induction energies , the results obtained with the 
 atz basis set are considered a reasonable approximation of the cbs 
 limit . for exchange  repulsion and dispersion energy ,
a two - point 
 scheme was applied to extrapolate the energy calculated at adz and 
 atz to the cbs limit .
difference between the sapt2 + energy components calculated 
 using 
 aug - cc - pvxz ( x = d , t , and q ) basis set with the value obtained using 
 aug - cc - pv5z . as mentioned in 
 the previous section , different energy components converge at different 
 rates with respect to the basis set size .
electrostatic and induction 
 energies calculated using the atz basis set are sufficiently converged , 
 while the exchange  repulsion and dispersion energy terms are 
 not .
thus , a two - point extrapolation scheme is used to extrapolate 
 the exchange  repulsion and dispersion energy at atz basis set 
 to the complete basis set limit . the extrapolated dispersion energy 

is further scaled by a factor of 0.89 to compensate for higher order 
 dispersion terms missing in the sapt2 + approach ( refer to methods for details of the extrapolation and scale 
 factor determination ) .
then the total sapt2 + interaction energy is 
 obtained by summing up the individual energy components .
finally , 
 the quality of the sapt2 + interaction energy at different basis set 
 levels was examined by comparing with the ccsd(t)/cbs interaction 
 energy of the s66 data set ( figure 3 ) . in general , the sapt2+/cbs estimates with scaled 
 dispersion energy ( which will be referred to as sapt2+/cbs / scaled 
 to distinguish from the sapt2+/cbs values without scaling ) have the 
 smallest mean unsigned error ( mue ) among tested combinations , given 
 mue values of 0.16 , 0.68 , 0.47 , and 0.17 kcal / mol ( and rmse of 0.25 , 
 0.80 , 0.56 , and 0.21 kcal / mol ) for sapt2+/cbs / scaled , sapt2+/cbs , 
 sapt2+/atz , and sapt2+/adz , respectively .
it is not surprising that 
 the sapt2+/adz combination yields very reasonable total interaction 
 energy as this is consistent with a previous report .
however , the small error in the total interaction energy 
 of sapt2+/adz is due to the error cancellation of individual energy 
 components .
to best estimate individual energy components , sapt2+/cbs / scaled 
 results remain the most accurate choice in this study .
however , for 
 the larger systems where aug - cc - pvtz calculations are not practical , 
 sapt2+/adz may be considered as an alternative .
for all 707 pairs 
 in the s1017 database reported in this study , the same strategy 
 ( sapt2+/cbs / scaled ) is applied to calculate the individual energy 
 components .
errors of sapt2 + interaction 
 energy compared to ccsd(t)/cbs estimation for dimers in the s66 data set .
sapt2+/aug - cc - pvdz 
 energy is shown in solid circles , sapt2+/cbs is shown in triangles , 
 and sapt2+/cbs / scaled is shown in hollow circles .
the electrostatic 
 interaction energy due to amoeba point multipoles as well as the charge 
 penetration correction ( multipoles + cp ) are calculated and compared 
 with sapt2 + data for the s1017 data set ( excluding 77 
 complexes containing an ethyne molecule due to a lack of amoeba parameters ) 
 ( supporting information table s3 ) . in our 
 current model , the parameters z and  for each 
 atom are uniquely determined by the element type .
 is set equal to z , or if z is less than 2 , then  
 is set to 2 .
the only parameter to be determined for the penetration 
 correction is  in eq 3 . for each of the 
 h , c , n , and o elements ,
carbon , nitrogen , and oxygen all have three  
 values representing sp , sp , and aromatic cases . 
 for sulfur , distinct  values are used for sulfide and sulfur 
 iv , while p , f , cl , and br have only a single  value per element 
 in current parametrization . in general , after fitting of  
 parameters , the new electrostatic model with charge penetration correction 
 shows excellent agreement with the sapt2 + results ( figure 4 ) . taking the valence- parameter set as an 
 example , for dimers near the equilibrium distances ( rmin ) , i.e. , 0.90 , 0.95 , 1.00 , 1.05 , and 1.10 of rmin , the mean unsigned error ( mue ) of the original 
 point multipoles is 3.16 kcal / mol , which is reduced about 5-fold to 
 0.57 kcal / mol after inclusion of the charge penetration correction 
 ( table 3 ) . for the dimers at very short separation , 

i.e. , 0.70 and 0.80 of rmin , the mues 
 for the corrected and uncorrected electrostatic energy are 3.28 and 
 19.16 kcal / mol , respectively .
as shown in figure 4 , it is striking that point - multipole - based electrostatic 
 energy alone yields very large errors for dimers in close contact , 
 and the simple charge penetration correction applied here is able 
 to systematically improve agreement with sapt - derived electrostatics . 
 based upon the mean unsigned errors , the charge penetration corrected 
 model results in a percentage error of 13.6% and 13.4% at near - equilibrium 
 and very short separations , respectively .
in contrast , the uncorrected 
 model has errors of 53% and 69% for these same two distance ranges . 

it is clear the charge penetration corrected model not only reduces 
 the magnitude of absolute and relative errors compared to sapt but 
 also provides consistent performance over a range of distances . in 
 the uncorrected model ,
the percentage of error at very short distances 
 is larger than at near - equilibrium distances , due to the increased 
 effect of short - ranged charge penetration .
plots of multipole 
 electrostatic energy ( kcal / mol ) against the 
 reference sapt2+/aug - cc - pvtz calculation for ( a ) near - equilibrium 
 ( 0.90 , 0.95 , 1.00 , 1.05 , and 1.10 ) complexes taken from the s1017 
 data set , ( b ) expanded plot of the boxed region in a , and ( c ) short - range 
 ( 0.70 and 0.80 ) complexes in the s1017 data set . the uncorrected 
 amoeba point multipole energy ( multipoles only ) is shown in red circles , 
 and the charge penetration corrected point multipole energies using 
 the valence- parameter set ( multipoles + cp ) are denoted by 
 blue crosses .
charge penetration corrected model 
 using the fitted- parameter set . for s66 dimers at near - equilibrium separations and 
 using uncorrected 
 amoeba multipoles , the hydrogen - bonded complexes exhibit the largest 
 mean unsigned error of 4.41 kcal / mol , compared to mues of 3.08 and 
 2.08 kcal / mol for dispersion - dominant and mixed complexes ( table 4 ) .
this is not surprising since the hydrogen - bonded 
 complexes generally have the strongest electrostatic interactions . 

however , in terms of relative errors , the dispersion - dominant complexes 
 carry the largest error at 105% , while the hydrogen - bonded and mixed 
 complexes have the mean percentage of errors of 30% and 58% , respectively . 

it is somewhat surprising the dispersion - dominant complexes have such 
 absolute and relative errors , as they are normally considered to have 
 the weakest electrostatic interaction among the three types .
charge penetration corrected model 
 using the fitted- parameter set . to help understand why dispersion - dominant complexes 
 have such 
 large relative errors , the electrostatic energies of benzene dimers 
 and  stacked and t - shaped complexes , as well 
 as hydrogen - bonded water dimers , are shown in figure 5 . for the uncorrected amoeba model ,
the calculated electrostatic 
 energy is positive for the  benzene pairs yet 
 qm calculations suggest the interaction is attractive with a negative 
 electrostatic energy . taking the electrostatic energy for this pair 
 at the equilibrium distance as an example , the sapt2+/cbs / scaled calculation 

yields a value of 2.6 kcal / mol , while the uncorrected amoeba 
 multipoles give + 1.0 kcal / mol , an error of 3.6 kcal / mol or 138% . for 
 the t - shaped benzene dimer ,
the sapt2+/cbs / scaled and the uncorrected 
 amoeba multipoles have values of 2.2 and 0.4 kcal / mol , 
 respectively .
the unsigned error is 1.8 kcal / mol or 82% of the sapt 
 values , both somewhat less than for the  complex . 

in contrast , although the 
 hydrogen - bonded water dimer has the larger electrostatic energy of 
 7.2 kcal / mol , the uncorrected model has an unsigned error 
 of 1.6 kcal / mol and a relative error of only 22% .
this trend is in 
 accordance with the averaged errors reported in table 4 and suggests the electrostatic interaction in dispersion - dominant 
 complexes is the most charge penetration dependent .
first , in the nonpolar molecules , the electron 
 distribution is more  balanced  ; i.e. , there is more 
 electron density on the hydrogen atoms , hence a stronger penetration 
 effect for hydrogens .
second , in the stacked benzene dimer , interactions 
 between heavier atoms , carbon 
carbon for example , suffer stronger 
 charge penetration effect , thus weight more in electrostatic energy . 
 for hydrogen - bonded pairs
, although the percentage of error is relatively 
 low for the uncorrected atomic multipoles , the absolute error remains 
 significant .
therefore , a correction is still necessary in order to 
 achieve better accuracy in the force field .
it is notable that , after 
 the charge penetration correction , the mean unsigned errors of all 
 three types of complexes are reduced to 0.50.6 kcal / mol near 
 the equilibrium distances , which is approaching the possible error 
 of the qm calculation itself .
plots of the electrostatic energy profiles of 
 water  water 
 and benzene  benzene dimer complexes .
the valence- parameter 
 set was used in calculations of the charge penetration corrected model 
 ( multipoles + cp ) .
the vertical line indicates the equilibrium distance . as 
 mentioned in methods , an alternative way to 
 derive
the  parameter is to fit the penetration - damped electrostatic 
 potential ( eq 4 ) to the target qm values .
an 
 attempt to use this fitting strategy has been also made , and the resulting 
 parameters have been compared .
the parametrization of  is restricted 
 to a single unique value for each element type , as before .
a brute 
 force scanning of the parameter using a grid size of 0.1  was used to search for the global minimum since the 
  parameter is less sensitive than . all 13 monomers 
 ( excluding ethyne ) in the s66 data set were used in fitting of the 
  for h , c , n , and o elements .
then  parameters were 
 determined as before with  values fixed to their potential - fitted 
 values .
the penetration parameter set obtained this way will be referred 
 to as the fitted- set , while the parameter set with  
 based on eq 5 will be referred as the valence- 
 set . with the fitted- parameters ,
the rmse of the electrostatic 
 potential of the 13 monomers calculated using eq 4 is greatly reduced to 0.07 kcal / mol , compared to an rmse of 0.95 
 kcal / mol for the valence- ( see supporting 
 information table s2 ) .
the performance of the two sets 
 of parameters has been compared using the s101 data set .
the overall 
 performance of the two parameter sets is very similar to mean unsigned 
 errors of 1.35 and 1.33 kcal / mol for the valence- and fitted- 
 sets , respectively ( table 3 ) . for near - equilibrium 
 pairs ,
the valence- set has a marginally better mue of 0.57 
 kcal / mol against 0.72 kcal / mol for the fitted- set . for short - ranged 
 pairs ,
the fitted- set with a mue of 2.84 kcal / mol yet is slightly 
 better than a mue of 3.28 kcal / mol of the valence- set .
similar 
 trends in rmses of the two sets of parameters are also observed . however
, 
 the valence- parameter set tends to have more balanced performances 
 for hydrogen - bonded , dispersion - dominant , and mixed complexes , giving 
 the mues of 0.50 , 0.53 , and 0.61 kcal / mol for the three groups , respectively 
 ( table 4 ) .
in contrast , the fitted- 
 set , with a mue of 0.47 kcal / mol for the hydrogen - bonded complexes , 
 shows subtly better agreement with sapt results yet has slightly worse 
 performances for the aromatic compounds .
the mues of the dispersion - dominant 
 and mixed complexes are 0.69 and 0.74 kcal / mol , respectively ( electrostatic 
 energy of individual pairs can be found in table s3 in the supporting information ) .
nonetheless , the two 
 sets of parameters all have excellent agreement with the sapt results 
 for the whole s1017 database , while the fitted- set yields 
 better electrostatic potential than the valence- set .
the charge penetration model also exhibited good transferability 
 during the fitting of  parameters .
although three atom types 
 are used for h , c , n , and o in the current parametrization , restriction 
 to a single  for each element also results in reasonable accuracy . 
 simply applying the arithmetic mean of the three  parameters 
 in valence- parameter
set for each element ( supporting information table s1 ) increases the mue by only 
 0.1 to 0.65 kcal / mol for the near - equilibrium pairs in the s66 set 
 ( table 5 ) . for pairs with shorter distances , 
 the mue increases by 0.8 to 3.52 kcal / mol in the same set .
only marginal 
 improvements in mues were found after optimizing the  parameters 
 for each element starting from the averaged value .
we believe this 
 demonstrates the robustness and transferability of the charge penetration 
 correction and the parametrization strategy . for the purpose of retaining 
 flexibility , we recommend the use of three atom types for each of 
 the h , c , n , and o elements in our final model .
an additional 
 set of parameters 
 which has a unique  for each element is also presented .
the 
 charge penetration effect is usually overlooked in molecular 
 mechanical models and traditional force fields .
our results show that 
 dma - derived point multipoles systematically underestimate the sapt 
 electrostatic interaction energy at typical molecular interaction 
 distances based on the 1017 dimers studied here ( see the supporting information ) .
an exponential damping 
 function providing a simple charge  charge penetration model 
 suitable for force field incorporation has been revisited , along with 
 a new parametrization strategy . the s1017 sapt - decomposed quantum 
 mechanical energy database is developed as a reference for parameter 
 training and for use in future force field comparison .
the database 
 is an extension of the s66 and s668 data set previously developed 
 by hobza and co - workers , with additional 
 prototype molecular complexes .
the decomposed energies are calculated 
 at the sapt2+/aug - cc - pvtz level of theory , with exchange 
the dispersion energy is further scaled to compensate for missing 
 higher order terms in the sapt2 + method . the total sapt interaction 
 energy is in excellent agreement with ccsd(t)/cbs results , which are 
 currently considered to be the  gold standard  for estimation 
 of intermolecular interactions , with a mean unsigned error of 0.16 
 kcal / mol for the s66 data set .
thus , the sapt results should provide 
 a reliable reference for force field development . by replacing 
 the idealized charge 
charge ( coulomb ) interaction 
 with the charge penetration corrected model ( eq 3 ) in the amoeba framework , the accuracy of calculated electrostatic 
 energies for the s1017 database is improved by 5-fold . for the 
 five distance pairs near the equilibrium distances ( i.e. , 0.901.10 
 times the equilibrium distance ) , the mean unsigned error of the charge 
 penetration corrected and uncorrected point multipole models are 0.57 
 and 3.16 kcal / mol , respectively ; for the extremely close distance 
 separations ( i.e. , 0.70 and 0.80 times the equilibrium distance ) , 
 the mean unsigned errors of the two models are 3.28 and 19.16 kcal / mol , 
 respectively .
the improvement for the corrected model is significant 
 and shows a consistent agreement with the quantum mechanics data at 
 both long and short distances . the robustness and transferability 
 of this model is also reflected in the use of very limited ( element - based ) 
 parameters .
the charge penetration correction is short - ranged and 
 rapidly converges to the classical coulomb interaction beyond 67 
  .
thus , it can be completely incorporated into the real space 
 of ewald summation without any additional computational cost in reciprocal 
 space . because simulations including penetration correction are clearly 
 feasible
, there is ongoing work dedicated to the optimization of parallel 
 scaling the coupled penetration / smooth particle mesh ewald approach . 

quadrupole penetration ) are also possible and have 
 been implemented in models such as sibfa .
nonetheless , the simple empirical charge penetration model presented 
 in this work provides us with an efficient approach to achieve accurate 
 electrostatic energy that is systematically modeled after sapt quantum 
 mechanical energy decomposition .
the change in the electrostatic component 
 requires re - examination of the van der waals interaction to arrive 
 at a balanced representation of the total energy .
overall we expect 
 this improvement in the electrostatic component will alleviate the 
 need for error compensation via other components and lead to more 
 balanced and transferable potential energy functions in general
. the 
 comprehensive sapt database we developed in this work will also be 
 useful for many others who are interested in understanding intermolecular 
 forces or evaluating different empirical models .
 classical 
 molecular mechanics force fields typically model interatomic 
 electrostatic interactions with point charges or multipole expansions , 
 which can fail for atoms in close contact due to the lack of a description 
 of penetration effects between their electron clouds . 
 these short - range 
 penetration effects can be significant and are essential for accurate 
 modeling of intermolecular interactions . in this work we report 
 parametrization 
 of an empirical charge  
 charge function previously reported 
 ( piquemalj .- p . ; 
 j. phys . 
 
 chem . 
 a2003 , 107 , 1035326313624 ) to correct for the missing penetration term in standard 
 molecular mechanics force fields . for this purpose , 
 we have developed 
 a database ( s1017 ) of 101 unique molecular dimers , each at 7 
 different intermolecular distances . 
 electrostatic , induction / polarization , 
 repulsion , and dispersion energies , as well as the total interaction 
 energy for each complex in the database are calculated using the sapt2 + 
 method ( parkert . 
 m. ; j. chem . 
 
 phys.2014 , 140 , 09410624606352 ) . 
 this empirical penetration model 
 significantly improves agreement between point multipole and quantum 
 mechanical electrostatic energies across the set of dimers and distances , 
 while using only a limited set of parameters for each chemical element . 
 
 given the simplicity and effectiveness of the model , we expect the 
 electrostatic penetration correction will become a standard component 
 of future molecular mechanics force fields . 
Introduction
Methods
Results
Conclusion