weak acid ionization 
 plays a central role in many chemical , biological , 
 and industrial processes . measuring and predicting the free energy 
 of this process ( i.e. , pka values )
have 
 become important experimental and theoretical tools , enabling control 
 over and insight into acid / base chemistry and ph - dependent biomolecular 
 transformations .
although a large number of theoretical methods have 
 been developed to predict pka values ( see 
 recent reviews ( 1 ) and ( 2 ) and references therein ) , 
 less attention has been given to understanding the molecular - scale 
 mechanism of acid dissociation .
this is largely due to the complexity 
 and cost of analyzing the full free energy surface , which involves 
 not only a chemical reaction but also significant solvent and solute 
 reorganization . instead , most pka methods 
 rely on the use of thermodynamic cycles that break the complex process 
 of deprotonation into computationally tractable steps . for small molecules , 

this cycle typically combines the gas phase dissociation energy , calculated 
 with high - level quantum mechanical ( qm ) calculations , with de / solvation 
 free energies of the reactant and product species , calculated with 
 continuum solvent approaches . while the gas phase energies are quite accurate , the simple continuum 
 description of solvation free energies often requires parameter adjustment 
 to reproduce all of the entropic and enthalpic complexities that are 
 inherent in molecular solvation .
they have been improved by hybrid 
 approaches that include some number of explicit water molecules in 
 the solvation free energy analysis . in conjunction 
 with methodological refinements
, these qm - continuum solvent models 
 have enabled pka calculations of small 
 acids to within a few units of the experimentally determined values .
in contrast , the suite of 
 methods that have been developed to predict 
 biological pka values have typically used 
 a thermodynamic cycle in which the solute remains in a solvated environment . 

each has advantages and drawbacks as well , discussed in 
 ref ( 2 ) , however none 
 has yet demonstrated the ability to predict pka values within a few pka units 
 of the experimental value for all test cases .
it has been suggested 
 that challenging moieties ( e.g. , buried residues that have large pka shifts or those with strong coupling to conformational 
 changes or other titration sites ) will require methods that better 
 capture the underlying physics .
these 
 conclusions highlight how interesting and nuanced a biomolecular system 
 can be . while it is clear that protonation states influence everything 
 from stability to binding and reactivity , the details of their influence 
 is still sometimes poorly understood .
few methods to date have 
 simulated the actual process of interest , 
 the exchange of an excess positive charge from an acid to bulk water . 

to do so involves calculating the free energy surface for proton dissociation 
 in the bulk phase environment , from which one learns about the molecular 
 mechanism and rate .
this requires three things : ( 1 ) a method that 
 describes the underlying physics of bond formation and cleavage in 
 the bulk phase environment , ( 2 ) a method that is efficient enough 
 that when combined with enhanced sampling approaches can capture the 
 distribution of phase space that defines a free energy in the condensed 
 phase ( where entropy and energy are often a similar magnitude ) , and 
 ( 3 ) a reaction coordinate that is flexible enough to track the transferring 
 excess charge defect .
the first of these requirements is most 
 obviously filled by qm 
 approaches , which have been used in a handful of helpful studies to 
 probe acid dissociation .
parrinello 
 molecular dynamics ( cpmd ) to calculate the free energy profile of 
 histidine dissociation .
although the chosen reaction 
 coordinate limited their analysis to stopping just beyond the transition 
 state , they were able to calculate the 
 correct relative pka value by subtracting 
 the equivalent profile for water autodissociation .
similarly , park et al . used cpmd in combination with metadynamics 
 to sample the configurational space of acetic acid in bulk water .
while their trajectories were too short to sufficiently 
 converge the free energy surface , they were able to observe and characterize 
 multiple protonation and deprotonation events .
these reactions occurred 
 along a well - characterized pathway where the excess proton is briefly 
 shared between the acid and solvating water .
they also noted the existence 
 of a shallow minimum in the free energy surface corresponding to the 
 existence of a contact ion pair ( cip ) that is separated from the protonated 
 state by a small energy barrier . in an alternative approach that still 
 describes the bond dissociation at the ab initio level , uddin et al . 

recently used qm / mm to calculate the free energy surfaces for several 
 small molecules .
although they obtained 
 quite accurate pka values , their method 
 can not describe the mechanism of dissociation since they include only 
 one water molecule and the acid in the qm region .
hence , this method 
 changes the physical process from the dissociation of a delocalized 
 charge defect to proton transfer followed by separation of a localized 
 hydronium cation and acid anion in a solvating mm environment .
just 
 as for thermodynamic cycles , their results show that sampling an alternative 
 pathway between reactant and product states can yield accurate pka values for small molecules .
another 
 approach that describes the bond dissociation process is 
 reactive md ( rmd ) .
the efficiency of these models allows them to sample reactive dynamics 
 on much longer length and time scales than ab initio md .
since they 
 are parametrized to work with mm force fields , they also avoid the 
 boundary issues suffered by qm / mm calculations .
both of these factors 
 make rmd models ideal for simulating reactions in large , complex systems 
 with slow dynamics .
proteins and proton exchange membranes are two 
 important examples of such systems . without sufficient sampling in these types of systems , the final 
 result can be strongly influenced by the initial conditions .
a common 
 criticism of rmd models is that they are only approximations to the 
 true quantum mechanical behavior of the system .
while this is certainly 
 true , if the functional form is both flexible enough and parametrized 
 correctly , then the reactive model will retain the ability to accurately 
 reproduce the potential and free energy surfaces from which it was 
 parametrized . the multistate empirical 
 valence bond ( ms - evb ) method is a rmd 
 force field that has been used to characterize proton transport in 
 many condensed phase and biological environments .
it is especially well suited to the challenge of acid ionization 
 since it inherently defines a center of excess charge that can be 
 used to track reactions involving proton transport , thereby avoiding 
 issues with geometric reaction coordinates .
although 
 ms - evb has been used to study weak acids and amino acid dissociation , these efforts involved 
 a larger degree of empiricism since the pka value was needed to fit the msevb potential .
it has been a long 
 sought goal to find a procedure that would minimize the empiricism 
 in rmd models , such that the bulk property of interest was derived 
 entirely from quantum data in combination with extensive sampling .
this work reports the achievement of that goal , a multiscale approach 
 to calculate the free energy profile for acid ionization that is based 
 entirely on qm data and yields an accurate absolute pka value .
an iterative procedure , in the spirit of adaptive 
 force matching , is used to sample the reactive 
 phase space , obtain reference ab initio data ( qm / mm forces ) , and fit 
 a rmd potential .
care has been taken to modify the functional form 
 such that it is sufficiently flexible to reproduce the qm / mm reference 
 data and , in principle , capture the essential physics .
our approach 
 is applied to aspartic acid ( asp ) and to glutamic acid ( glu ) .
both 
 of these amino acids play crucial roles in proton and ion transport 
 in biological channels and pumps and are also commonly involved in 
 salt bridges that influence protein dynamics , structure , and binding . 

we show that accurate estimates of the pka values of both asp and glu are obtained from the resulting reactive 
 models , which is the first time a reactive model has been able to 
 accurately predict the correct pka without 
 being explicitly parametrized to do so .
finally , we discuss the mechanistic 
 description our models provide of acid deprotonation in water . to emphasize that these models were fit entirely to qm / mm data 
 and without adjustments to match the experimental pka
, we will hereafter refer to them as multiscale reactive 
 molecular dynamics ( ms - rmd ) models .
it 
 should be noted , however , that the ms - rmd asp and glu models are parametrized 
 to work with a refined version of the ms - evb3 hydronium model that 
 accurately describes a proton solvation and transport in water .
ms - rmd is similar to the ms - evb 
 framework that has been successfully used to model proton solvation 
 and transport in many aqueous systems .
the power in this approach lies 
 in its ability to describe electronic delocalization as a linear combination 
 of topologically distinct states . for a given configuration ,
for example , state |i 
 could be a protonated amino acid in a box of water ( figure 1a ) .
state |j would have 
 the same coordinates but have a topology describing a cip between 
 a solvated hydronium ion and a deprotonated amino acid ( figure 1b ) . at each step in the simulation
, a state search 
 algorithm first determines all the possible states based on geometric 
 criteria and then calculates the energy of each state .
the lowest 
 energy state is termed the pivot state and is the start of the state 
 search algorithm at the next step . since the number of waters , and 
 therefore the number of possible states , grows exponentially with 
 distance from the pivot state , only states in the first three solvation 
 shells
several different bonding 
 topologies are shown for a given set 
 of nuclear coordinates , r. dashed lines surround the 
 protonated molecule , as defined by the bonding topology .
the reactive system is described by the following hamiltonian:1where r is the vector describing 
 the nuclear coordinates , hij is the sum of mm potential terms for state i , and hij is the coupling 
 between states |i and |j. the details of these terms will be provided in the force field section . with
the hamiltonian defined 
 this way , the relative weight of each state is described by the components 
 of the eigenvector , ci:2finally , forces are calculated using 
 the hellmann  feynman 
 theorem .
these forces are then passed to a standard md integrator 
 such as velocity verlet .
free energy calculations require a 
 continuously varying reaction 
 coordinate ( rc ) to describe the process of interest .
the evb formalism 
 provides a convenient reference for the rc ; the center of the excess 
 charge ( cec):3 in this equation , ri is the geometric center - of - charge 
 for the protonated molecule in state i , and ci is the relative weight of 
 that state .
we then use the distance between the cec and the center 
 of mass of the asp or glu carboxyl group to describe the progress 
 of the deprotonation reaction .
the 
 diagonal elements of the ms - rmd hamiltonian 
 matrix are derived from the standard charmm mm force field .
the protonated form uses all of the standard 
 parameters with the exception of the acidic o  h bond . this 
 bond
is replaced with a morse potential to allow dissociation:4where r is the bond length 
 and d ,  , and r0 are parameters fit to reproduce the bond length and frequency of 
 the harmonic potential it replaces while still allowing for a reasonable 
 bond dissociation energy . the deprotonated form of amino acid uses 
 all of the standard charmm bonded and nonbonded terms and parameters . 
 since the charmm force field was not originally intended for use in 
 reactive simulations , additional repulsive terms are used between 
 the carboxyl oxygens of the deprotonated model and the hydronium oxygen 
 and hydrogens .
these terms help to correct for the overattraction 
 at close distances due to the use of point charges .
the functional forms for the repulsions are the same as 
 those used in the ms - evb3 hydronium model:56where b , b , b , c , and c are fitted parameters .
following 
 the logic in previous 
 work , and to avoid parameter redundancy , doo and doh were 
 fixed the same values used in ms - evb3 .
the sum of gaussians term in 
 eq 5 is included to help the fitted potential 
 to reproduce the correct asymmetric solvation structure around the 
 hydronium , and the same switching function is used as given in eq 
 9 of ref ( 28 ) .
the asp and glu models are designed to work with a refined version 
 of the ms - evb3 force field , herein referred to as msevb 3.1 .
this 
 refinement uses the same functional forms for the hamiltonian but 
 revised parameters . using a genetic algorithm ,
the fitting parameters 
 were optimized using the original ab initio data sets as well as potential 
 energy surface scans for the eigen ion at the same level of the original 
 ab initio calculations .
furthermore , the oxygen partial charge of 
 the hydronium is an adjustable parameter ( the charge for the hydrogens 
 is automatically determined to maintain the hydronium + 1 total charge ) . 

parameters for this 
 force field are presented in the supporting information . in this study ,
the limiting species are the de / protonated 
 asp or 
 glu models and the ms - evb 3.1 solvated proton in spc - fw water . the ms - rmd force field is parametrized to smoothly switch 
 from one protonated species to another .
the asymmetry of the reaction 
 necessitates several modifications to the original ms - evb equations , 
 which were designed with water in mind .
the off diagonal coupling 
 term , hij , is designed 
 to reproduce the correct transition state geometry .
previous ms - evb 
 models used a complex off - diagonal term to describe a symmetric reaction 
 pathway between two waters .
however , the asymmetry of the asp - water 
 system was better fit by a simple gaussian potential:7where roh is the 
 distance between the deprotonated asp oxygen and the excess proton . 

the other parameters determine the shape and size of the gaussian 
 and are fit to the qm / mm data . in practice ,
protonation reactions 
 occurring between different species require a constant energy term 
 ( vii ) to be added to 
 the deprotonated state .
this term accounts for the difference in energy 
 resulting from different mm force fields . without it ,
the protonated 
 and deprotonated asp potential energy surfaces would be offset by 
 hundreds of kilocalories per mole , and no reaction would ever take 
 place .
this difference is due primarily to electrostatics and can 
 be estimated by subtracting the coulomb energy of the most favorable 
 hydronium state ( containing the favorable electrostatic interactions 
 between the hydronium cation and amino acid anion ) from the coulomb 
 energy of the protonated ( and neutral ) amino acid state .
figure 2 shows a plot of the average value of this term 
 as a function of rc . since vii controls the relative stability of the protonated and deprotonated 
 states , it was used in past work to tune the behavior of the model 
 to reproduce experimental quantities such as pka . in the current work ,
an iterative 
 fitting approach is used to rigorously determine the value of vii along with all other ms - rmd 
 parameters .
plot of the difference in coulomb energy between the protonated 
 asp state and the lowest energy non - asp state . at short distances , 

it is similar in some 
 ways to that used in previous work in which ms - rmd models were force 
 matched from aimd data .
there are , however , 
 a few important differences , including the use of qm / mm instead of 
 aimd , the use of empirical functional forms instead of tabulated potentials , 
 and the use of an iterative scheme for sampling configuration space 
 starting from a guessed hamiltonian .
these choices were partially 
 motivated by the success that has been demonstrated with the adaptive 
 force matching method by wang and co - workers .
reactive 
 hamiltonian in order to generate a set of configurations along the 
 reactive pathway .
next , high - level qm / mm calculations are performed 
 to collect the atomic forces for those configurations .
the new reactive 
 hamiltonian is then used to generate new configurations via umbrella 
 sampling , and the process is repeated until the parameters reach the 
 desired convergence .
details of each step are discussed in the text . due to the nonlinearity of several 
 of the reactive force field 
 parameters ,
the reactive model is fit to the qm / mm reference data 
 using a genetic algorithm to minimize 
 the residual:8here , nc and na are the number of configurations and the number 
 of atoms in each configuration , respectively .
w(rij ) is the weight for each atom 
 that , unless specified otherwise , is set to 1 .
fij is the atomic force from the current ms - rmd 
 parameter set and fijqm / mm is the reference force from 
 the qm / mm calculation .
the fitting calculations were run in 
 parallel and used 4000 genomes 
 ( parameter sets ) per generation . to lessen the danger of overfitting
the first step followed the 
 uniformly distributed mutation scheme outlined in section 3.1 of ref ( 40 ) .
the range of parameters 
 was chosen to be as broad as possible without allowing unreasonable 
 values .
for example , c3 in the off - diagonal 
 term would be unreasonable if it were greater than the aspo 
once a converged parameter set was reached with 
 the uniform distribution , the mutation scheme was switched to the 
 normally distributed scheme described in section 3.2 of ref ( 40 ) .
this scheme is very good 
 at quickly finding a local minimum and was therefore used to find 
 the best genome in the vicinity of the result from the first step . 

since the second scheme is easily trapped in local minima , it was 
 not ideal for use on its own .
the combination described here was found 
 to quickly and reliably converge on the optimum solution .
the discrete 
 recombination scheme described in section 2.2 of ref ( 40 ) is used with both mutation 
 schemes .
the nature of the reactive model introduces additional 
 complexity 
 into this otherwise simple scheme .
vii , which directly controls the depth of the protonated well 
 relative to the deprotonated cip , is a constant that is only included 
 for states in which the amino acid is protonated .
the off - diagonal term acts 
 over a broader range spanning the protonated amino acid and cip configurations 
 but also has a negligible effect once the cip dissociates .
in contrast , 
 the repulsive terms , which are only applied to the deprotonated amino acid , are significant at longer distances and negligible at 
 distances below the transition state .
furthermore , since both vii and the repulsive terms 
 are added to the diagonal matrix elements , they tend to have correlated 
 behavior if fit simultaneously . thus , the value of vii for the first iteration after generating new qm / mm data 
 is chosen to be the average value described in figure 2 .
this value has been found to be within a few kilocalories 
 per mole of the final value of vii and 
 is therefore a good initial approximation .
holding the value of vii fixed , the off - diagonal and repulsive terms 
 can be fit to the full set of qm / mm reference forces .
once the off - diagonal 
 and repulsive terms reach their optimum value for the current value 
 of vii , umbrella sampling is used to generate 
 a pmf with the model . since vii s 
 effect is most important at and around the transition state , a range 
 of rc values immediately surrounding the transition state
is chosen 
 from the pmf , and configurations from that range are then used to 
 fit a new value of vii . using this new 
 value of vii , the off - diagonal and repulsive 
 terms are refit while vii is again held 
 constant .
this procedure fits an empirical model from solely qm / mm 
 data and results in models that correctly predict the bulk phase pka values for both asp and glu .
solution simulations of asp and 
 glu include one solute molecule solvated by 486 and 467 waters , respectively . 
 to reduce interactions with the charged backbone protonatable sites , 
 the n - termini were capped with an acetyl group and the c - termini were 
 capped with methylamine .
the systems were equilibrated in the npt 
 ensemble at 1 atm and 310 k for 200 ps before changing to the nvt 
 ensemble with a cubic box of 24.55486  on a side for asp and 
 24.24938  for glu .
all reactive simulations were run with a 
 modified version of the lammps simulation package .
the rc was chosen to be the distance between the cec and the center - of - mass 
 of the amino acid s side chain carboxylic acid .
umbrella sampling 
 was used to evenly sample the rc . for the asp system ,
25 windows were 
 spaced every 0.36  starting at a rc value of 1.0  .
the glu system used 36 
 windows spaced every 0.25  with 20.0 kcal / mol restraints . each 
 production window
was further equilibrated for 100 ps , and then statistics 
 were collected for 1.7 ns .
the pka is calculated from the resulting pmf using the following 
 equation:9where w(r ) is the free energy from the pmf and  
 is the product of the 
 simulation temperature and the boltzmann constant .
the integral is 
 calculated from zero to the transition state , as denoted by .
c is the standard state concentration whose 
 value is 1660  / molecule and results from the entropic 
 freedom that is gained by the proton when it dissociates from the 
 acid .
atomistic configurations for the qm / mm 
 calculations were selected from the umbrella sampling trajectories . 
 in order to ensure a uniform distribution along the rc ,
individual 
 configurations were sorted by their rc value into 60 0.067  
 wide bins .
five configurations were then selected at random from each 
 bin , giving a total of 300 configurations for each iteration .
qm / mm 
 was performed in cp2k using b3lyp with the double - zeta basis set .
unrestricted dft was used to allow 
 the most accurate calculation of forces near the dissociation barrier . 

b3lyp was chosen over mp2 because recent studies have shown comparable 
 accuracy with a reduced computational cost for water and the solvated 
 excess proton .
calculations were run on 256 intel sandy bridge cores and on average 
 completed within 10 min .
the qm region included the entire amino acid 
 residue and all of the waters found by the ms - rmd state search . because 
 proximity to classical point charges reduces the accuracy of forces 
 calculated on qm atoms , an additional 3  buffer layer of qm 
 water was added around the qm atoms . all together , this results in 
 a qm region containing an average of 208 atoms , or one asp / glu and 
 61 waters . in total ,
four iterations were needed for asp parameters 
 to converge , while glu parameters converged in three iterations .
parameters 
 for the new asp 
 and glu models are presented in table 1 .
these 
 parameters were fit to qm / mm reference forces using the previously 
 described iterative scheme to ensure convergence .
it must be noted 
 that the asp and glu models were parametrized using the same method , 
 but entirely independent data sets .
the fact that the final parameters 
 for each model are so similar is not surprising , however , given the 
 molecular similarity and small difference in experimental pka values .
as the next section discusses , the 
 pka predicted by each model is a very 
 good estimate of the experimental value .
the fact that the models 
 are able to capture such subtle differences is evidence that the methodology 
 works well .
the potential of 
 mean force ( pmf ) of deprotonation describes the 
 free energy of the process as a function of the rc . to obtain a free 
 energy change between any two points on the pmf , one simply integrates 
 the pmf .
figure 4a shows the pmf for deprotonation 
 of asp , and figure 4b shows the pmf for glu 
 ( black lines ) .
as depicted in figure 5 , the 
 rc gives the distance between the center - of - mass of the carboxylic 
 moiety and the cec .
rc values around 1.3  correspond to the 
 protonated amino acid , while values between 2.1 and 5  correspond 
 to the stable cip .
values above 5.5  represent an excess proton 
 separated from the anionic acid by bulk solvent .
pmfs and cmax distributions 
 from the final ( a ) asp and ( b ) glu models .
the colormap in the background gives 
 the cmax probability distribution 
 normalized for each rc value .
this distribution describes the extent 
 of charge delocalization and , by extension , the solvation structure .
snapshots from umbrella sampling trajectories 
 illustrating the 
 types of structures suggested by the cmax distributions and the 2d - rdfs .
arrows indicate the region 
 of the pmf where the windows were run , while the distances in the 
 box are the center of the harmonic restraint .
the position of the cec is shown as a yellow sphere . 
 in eigen structures ,
the cec is nearly aligned with the central water s 
 oxygen . in zundel structures ,
the protonated forms of asp and glu are much more stable 
 than the 
 deprotonated forms , which is consistent with their status as weak 
 acids . evaluating the integral in eq 9 estimates 
 the pka values for asp and glu to be 3.82 
  0.17 and 4.56  0.18 , respectively .
these are both in 
 very good agreement with the experimental ranges of 3.71 to 3.90 for 
 asp and 4.15 to 4.31 for glu .
the barrier that prevents protonation is due to the cost of 
 forming a zundel - like arrangement between the acid and hydronium . 

the depth and extent of the cip varies with the system ; however its 
 existence has been predicted in a variety of weak acid systems simulated 
 with different methods . when preparing a simulation to 
 calculate a pmf of deprotonation and a pka , care must be taken to ensure that the rc is sampled far enough 
 beyond the cip to ensure bulk - like properties .
if this is not done , 
 then the pmf from which the pka is calculated 
 will be shifted and the resulting  pka  will neither be accurate nor reflect the quality 
 of the model .
fortunately , sampling well beyond the cip is not a problem 
 with ms - rmd models due to their efficiency .
however , care needs to 
 be taken when running qm / mm simulations to ensure that a large enough 
 qm box and enough solvating waters are used to not only sample beyond 
 the cip but also to avoid boundary issues . with the pmf in hand , one
the 
 ms - rmd methodology conveniently provides a measure of the structure 
 of the solvated complex .
it has been shown that the square of the 
 largest component of the ground state eigenvector ( cmax , from c in eq 3 ) describes the solvated structure .
a value of 1 indicates the limiting case of an entirely 
 localized excess charge , while values less than that quantify the 
 degree of charge delocalization .
eigen structures , where the excess 
 charge is shared among a central hydronium and three h - bound water 
 molecules , have a cmax value 
 around 0.65 .
zundel structures , which share the excess charge more 
 equally between two water molecules , have a cmax value closer to 0.5 .
it should be emphasized 
 that there is always additional delocalization to other surrounding 
 water molecules , even for these limiting eigen- and zundel - like cations . 
 by plotting the cmax probability 
 density as a function of the rc , we gain insight into the structural 
 changes that accompany deprotonation .
this probability density is 
 plotted as a colormap in figure 4 . at low rc 
 values ,
both asp and glu have an average cmax greater than 0.9 , which indicates that the structure 
 and dynamics mostly resemble the limiting case of a protonated acid 
 in bulk water with little charge transfer to the surrounding water 
 molecules .
this is notably different than the protonated water molecule , 
 which always involves a significant amount of charge delocalization .
as the proton begins to dissociate , the cmax decreases as the proton goes 
 through a clear zundel - like transition between the acid and solvating 
 water . beyond the transition state ,
the system begins to resemble 
 an eigen cation where one of the solvating waters is replaced with 
 the acid anion .
however , the eigen structure seen in the cip is distinctly 
 different from that seen in bulk water , as evidenced by a cmax of 0.73 .
the close proximity 
 of the anionic acid helps to stabilize the hydronium , leading to a 
 more localized charge .
further evidence of the uniqueness of the cip 
 is that there is a well - defined zundel transition around 4  
 linking the cip to the bulk - like structure at longer distances .
eigen 
 transitions in bulk , the water around the anionic amino acid is more 
 structured , as evident in the distinct zundel - like transition at a 
 rc value around 4  .
figure 5 shows snapshots 
 of the typical structures found at different rc values and characteristic 
 of significant regions in the pmf .
the bonding topology reflects o  h distances 
 less than or equal to 1.6  , which is commonly considered to 
 be the length of a hydrogen bond .
it is also interesting to 
 look at the structural changes that occur directly around the cec . 

to do so , the two - dimensional radial distribution functions ( 2d - rdfs ) 
 between the cec and surrounding water oxygens ( ow ) and 
 water hydrogens ( hw ) are shown in figure 6 .
a slice along the first dimension of the 2d - rdfs 
 gives a traditional 1d - rdf , which describes the probability of finding 
 the indicated pair a given distance apart .
scanning the second dimension 
 shows how the solvation structure changes as the reaction progresses 
 according to the rc . in combination with the cmax analysis ,
the 2d - rdfs show that the solvent 
 has well - defined rearrangements as the proton dissociates from asp 
 and travels to bulk .
figure 6a shows the 2d - rdf 
 for asp s cec  ow distance .
once asp is deprotonated , 
 the first peak is near 2.5  , which is characteristic of the 
 eigen structure .
when the rc reaches 4  , the peak distance drops 
 to around 1.5  , which is characteristic of a slightly asymmetric 
 zundel cation where the cec resides close to the midpoint of the o  o 
 vector . the 2d - rdf showing asp s cec 
2d - rdfs for ( a ) cec  ow and ( b ) cec  hw with asp .
these plots show 
 how the water structure changes 
 around the cec as deprotonation occurs .
2d - rdfs for ( a ) cec  ow and ( b ) cec  hw with glu .
we 
 have presented a procedure for fitting the parameters of a ms - rmd 
 model to reproduce the forces from qm / mm calculations . due to the 
 efficiency of ms - rmd
, we are able to collect and analyze data from 
 large time and length scales that can reveal intricate behavior and 
 is essential for converged bulk phase thermodynamic properties . this 
 new methodology was used to parametrize new models of both asp and 
 glu that correctly predict the correct experimental pka values .
this is the first time that a ms - rmd model has 
 been able to reproduce the experimental pka values without any fitting to experimental data .
multi - nanosecond 
 trajectories were analyzed to extract important information about 
 the mechanism underlying the deprotonation process .
it is shown that 
 the deprotonation process follows a well - defined series of zundel 
 to eigen transitions as the excess proton travels from asp to bulk 
 solution .
this work is considered continued progress toward 
 the goal of having 
 flexible and systematic algorithms to reliably make a multiscale connection 
 between accurate electronic structure calculations and efficient empirical 
 models to describe complex reactive processes . as more accurate ab 
 initio methodologies ( e.g. ,
those explicitly accounting for electronic 
 correlation ) become computationally tractable for qm / mm of condensed 
 phase systems , the procedure presented herein and related force matching 
 algorithms will be invaluable tools for mapping high level data into 
 efficient rmd models that can significantly extend the time and length 
 scales of reactive simulations to explore complex phenomena .
while the results presented herein are important , the true test 
 of this procedure will be to apply it in biological systems where 
 the pka of individual residues can be 
 shifted by several pka units .
proton transport 
 in proteins is very sensitive to the underlying protonatable residue 
 models , and it is therefore crucial to have a procedure to ensure 
 that those models are physically accurate . furthermore , biological 
 systems are orders of magnitude larger and have much slower dynamics 
 than the systems presented here .
while it is conceivable that qm / mm 
 md simulations could reach the same sort of time scales as the solution 
 simulations presented here , they would be prohibitively expensive 
 to run the length - and time - scales needed to account for the slow dynamics 
 commonly found in biological processes .
the ms - rmd methodology has 
 already been shown to scale well to hundreds of thousands of atoms 
 and can realistically reach tens of nanosecond simulations with large 
 explicit biomembrane systems .
work is already underway to implement 
 these reactive models in protein environments and to use the presented 
 methodology to fine - tune these and other individual amino acid models 
 to behave properly in their specific protein environments .
 accurately 
 calculating a weak acid s pka from 
 simulations remains a challenging task . 
 we report 
 a multiscale theoretical approach to calculate the free energy profile 
 for acid ionization , resulting in accurate absolute pka values in addition to insights into the underlying mechanism . 
 
 importantly , our approach minimizes empiricism by mapping electronic 
 structure data ( qm / mm forces ) into a reactive molecular dynamics model 
 capable of extensive sampling . 
 consequently , the bulk property of 
 interest ( the absolute pka ) is the natural 
 consequence of the model , not a parameter used to fit it . 
 this approach 
 is applied to create reactive models of aspartic and glutamic acids . 
 
 we show that these models predict the correct pka values and provide ample statistics to probe the molecular 
 mechanism of dissociation . 
 this analysis shows changes in the solvation 
 structure and zundel - dominated transitions between the protonated 
 acid , contact ion pair , and bulk solvated excess proton . 
Introduction
Methodology
Results and Discussion
Conclusion