Nuclear Quantum Effects on Autoionization of Water Isotopologues Studied by Ab Initio Path Integral Molecular Dynamics
NNuclear Quantum Effects on Autoionization of Water Isotopologues Studied by AbInitio Path Integral Molecular Dynamics
Bo Thomsen a) and Motoyuki Shiga b) CCSE, Japan Atomic Energy Agency, 178-4-4, Wakashiba, Kashiwa, Chiba,277-0871, Japan (Dated: 1 February 2021)
In this study we investigate the nuclear quantum effects (NQEs) on the acidity con-stant (p K A ) of liquid water isotopologues at the ambient condition by path inte-gral molecular dynamics (PIMD) simulations. We compared simulations using afully explicit solvent model with a classical polarizable force field, density functionaltight binding, and ab initio density functional theory, which correspond to empirical,semiempirical, and ab initio PIMD simulations, respectively. The centroid variablewith respect to the proton coordination number of a water molecule was restrainedto compute the gradient of the free energy, which measures the reversible work ofthe proton abstraction for the quantum mechanical system. The free energy curveobtained by thermodynamic integration was used to compute the p K A value basedon probabilistic determination. This technique not only reproduces the p K A value ofliquid D O experimentally measured (14.86) but also allows for a theoretical predic-tion of the p K A values of liquid T O, aqueous HDO and HTO which are unknowndue to its scarcity. It is also shown that the NQEs on the free energy curve can resultin a downshift of 4 . ± . K A units in the case of liquid water, which indicates thatthe NQEs plays an indispensable role in the absolute determination of p K A . Theresults of this study can help to inform further extensions into the calculation of theacidity constants of isotope substituted species with high accuracy. a) Electronic mail: [email protected] b) Electronic mail: [email protected] a r X i v : . [ phy s i c s . c h e m - ph ] J a n . INTRODUCTION The acidity constant, p K A , plays a fundamental role in acid-base chemistry. There is nodoubt about the importance of quantitatively estimating p K A for a given functional groupin a molecular species. Computational chemistry rooted in molecular theory has the meritthat it is possible to evaluate p K A independent of experiments. However, computationalevaluation of p K A has yet to reach the quantitative level of accuracy, even for the most basiccase of liquid water, i.e. , the water autoionization constant, p K W .Previous computational evaluations of p K A have been based on either an implicit solventmodel, in which a solute molecule is embedded in a solvent described by polarizable con-tinuum medium, or an explicit solvent model, in which both the solute and the solventare described explicitly as molecules. In general, the implicit solvent model is able toprovide p K A values with small computational effort but limited accuracy. The intrinsic errorof the implicit solvent model arises from the lack of complexity to describe solute-solventinteractions.The inclusion of explicit solvent molecules and extended conformational sampling have recently been shown to partially solve these issues. Both methods, however, requirecareful (re)-parametrization of the implicit solvent model and selection of the conformers tobe considered in the calculation. It is therefore preferable to deal with the solvent moleculesexplicitly for quantitatively estimating p K A , where solute-solvent interactions are takeninto account properly. In this case p K A is estimated directly from the free energy changeupon the protonation of a solute molecule by molecular dynamics (MD) techniques, such ascoordination-constrained MD (often to referred to as the blue-moon ensemble method). p K A is known to have a strong dependence on thermodynamic conditions such as thetemperature and pressure. For liquid water under a pressure of 15 MPa, the p K W value of14 at the ambient temperature decreases to about 11 in subcritical conditions at 300 ◦ C,and then increases to about 20 in supercritical conditions at 400 ◦ C. It is also known thatp K A clearly differs between the hydrogen isotopologues, e.g. , p K W of D O under ambientconditions is 14.86, which is larger than that of H O, 14.00. Constrained MD simulationsallow the estimation of p K A under different thermodynamic conditions. However, they willproduce identical results for the hydrogen isotopologues, because the free energy changecalculated by classical MD does not depend on the nuclear masses. The isotope effect on2he free energy can be traced back to the quantum nature whereby the kinetic and potentialenergy operators in the Boltzmann density do not commute. Therefore, it follows that thenuclear quantum effects (NQEs) should not be ignored for the quantitative estimation ofp K A . In fact, hydrogen is generally known to exhibit quantum behaviors such as zero-pointvibration and tunneling because of its light mass.The p K W of water has been a long-standing interest for theoretical studies, since the longtime-scale on which the autodissociation takes place makes it difficult to sample efficiently.Furthermore the strong solvation effects of the produced OH − and H + ions make it importantto consider the solvation structure explicitly.The computation of p K W of liquid water has been done with explicit solvent in severalstudies. The reaction mechanism and kinetics of the autoionization of liquid water havebeen the subject in extensive studies.
However, none of these studies have explored theNQEs on the p K W of water in ab initio simulations based on first principles.In this paper we introduce a coordination-restrained path integral molecular dynamics(PIMD) method based on an explicit solvent model to study the p K W of liquid water and itsisotopologues. PIMD is a rigorous approach that takes account of both nuclear quantum andtemperature effects based on the imaginary-time path integral theory for quantum statisticalmechanics . PIMD has been used to study the NQEs on several properties of water assummarized in a recent review. In this review the NQEs on the p K W of water is estimated by extrapolating the effect onthe water isotopologues to the case of infinite mass nuclei, which corresponds to the resultof a classical MD simulation. The change from classical to quantum nuclei was found tocause a downshift of around 3 p K A units. This empirical shift was subsequently used in thestudy by Wang et al. to correct the p K W calculated from classical MD. Two important references in the context of this paper are the study on the solvated protonand hydroxide ion by Marx et al. and the recent study on the NQEs in proton transportin water under an electrical field by Cassone. Here we will extend the PIMD methodfor p K W /p K A estimations, allowing the computations of quantum free energies upon theprotonation of the solute, by restraining the centroid variable of the coordination number(CN).In Section II we will outline the theory of coordination-restrained PIMD and the cal-culation of p K W in terms of probabilistic and absolute methods. Section III contains the3omputational details for the simulations conducted in this study. The results of MD andPIMD will then be discussed in Section IV. Furthermore we will investigate the isotopeeffects in pure D O and T O, and the solvated HDO and HTO isotope substituted waterlike molecules. Finally in Section V we will summarize our findings.
II. THEORYA. Coordination Number Restrained Path Integral Molecular Dynamics
Consider the quantum Hamiltonian for an N atom system within the Born Oppenheimerapproximation ˆ H = N (cid:88) I =1 ˆ P I M I + V ( ˆ R , . . . , ˆ R N ) , (1)where M I is the mass of the I th particle, ˆ P I is the momentum operator (cid:16) ˆ P I,x , ˆ P I,y , ˆ P I,z (cid:17) ,ˆ R I is the position operator (cid:16) ˆ R I,x , ˆ R I,y , ˆ R I,z (cid:17) and V is the potential. The quantum partitionfunction for this system is given as Z = (cid:90) d R · · · (cid:90) d R N (cid:68) R · · · R N (cid:12)(cid:12)(cid:12) exp (cid:16) − β ˆ H (cid:17)(cid:12)(cid:12)(cid:12) R · · · R N (cid:69) = Tr exp (cid:16) − β ˆ H (cid:17) . (2)Here β = k b T , where k b is the Boltzmann constant and T is the temperature. It is assumedthat the thermal de Broglie wavelength is smaller than the separation of two identical atoms,thereby allowing us to ignore the possibility for exchange of fermion/boson positions.By dividing the Boltzmann operator in Equation (2) into P terms, inserting closurerelations in coordinate spaces between each term of the Boltzmann operator, and applyingthe second-order Suzuki-Trotter expansion one can derive the following expression for thepartition function, Z = lim P →∞ Z P (3)where Z P = N (cid:89) I =1 (cid:34)(cid:18) M I P πβ (cid:126) (cid:19) P (cid:90) d R (1) I (cid:90) d R (2) I · · · (cid:90) d R ( P ) I (cid:35) × exp (cid:34) − β (cid:40) P (cid:88) s =1 N (cid:88) I =1 M I ω P (cid:16) R ( s ) I − R ( s − I (cid:17) + P (cid:88) s =1 P V (cid:16) R ( s )1 , . . . , R ( s ) N (cid:17)(cid:41)(cid:35) . (4)4he constant ω P is given as √ Pβ (cid:126) . The factor (cid:81) NI =1 (cid:16) M I P πβ (cid:126) (cid:17) P is a constant and will be omittedin the following, as it does not alter the relative free energy differences. Rearranging theexponent in the above equation we arrive at the following, Z P ∝ N (cid:89) I =1 (cid:20)(cid:90) d R (1) I (cid:90) d R (2) I · · · (cid:90) d R ( P ) I (cid:21) exp ( − βV eff ( { R } )) , (5)where V eff ( { R } ) = P (cid:88) s =1 (cid:40) N (cid:88) I =1 M I ω P (cid:16) R ( s ) I − R ( s − I (cid:17) + 1 P V (cid:16) R ( s )1 , . . . , R ( s ) N (cid:17)(cid:41) . (6)These two equations show that the quantum behaviour of an N particle system can bemimicked by considering an N P particle system. The individual terms in P are oftenreferred to as beads on a chain, where each bead corresponds to a single copy of the classicalsystem. This system is evolved according to its classical forces and its coupling to a polymerchain. The cost of this method thus scales steeply with the number of beads, and byextension accuracy desired. Convergence with respect to number of beads does howeverscale inversely with temperature, i.e. a few tens of beads can be sufficient to accuratelydescribe quantum effects for systems under standard conditions.The blue-moon sampling approach is here used to compute the free energy curve ofwater autoionization. Following a study by Sprik , the CN of an oxygen atom (labeled as“O ∗ ”) is used for studying this reaction. We introduce a rational function for the CN as n O ∗ ( { R } ) = (cid:88) j ∈ H − (cid:16) r O ∗ j d OH (cid:17) − (cid:16) r O ∗ j d OH (cid:17) , (7)where { R } is the set of atomic coordinates of the system, r ij are the distances to all hydrogenatoms in the system, and d OH is a constant set to 1.35 ˚A.The CN is in this study is restrained to vary the coordination of hydrogen in a singleOH − moiety from one to zero. These conditions correspond to the moiety forming a H Omolecule and a solvated OH − ion respectively. In practice, for a random oxygen atom, oneof its attached hydrogens is chosen to remain bound, while the other attached hydrogen andall other hydrogens are subject to the CN restraint. The same method is used in the caseof D O and T O for the deuterium and tritium atoms in the simulation, respectively. ForHDO and HTO in solution all hydrogens are restrained, while the core OD − or OT − is keptunrestrained. 5o derive the free energy difference in the blue-moon ensemble one has to consider thefree energy of a system restrained to a fixed coordination number ¯ n O ∗ , A (¯ n O ∗ ) = − β − log ρ (¯ n O ∗ ) + A . (8) A is a constant term taking care of the constants dropped in Equation (5). This term isonly necessary to consider when calculating the absolute value of the free energy, as it willcancel for relative free energy differences. The distribution ρ (¯ n O ∗ ) is the scaled probabilityfor finding the system with the given coordination number ¯ n O ∗ , which is expressed as ρ (¯ n O ∗ ) = (cid:42) δ (cid:32) ¯ n O ∗ − P P (cid:88) s =1 n O ∗ (cid:0) R ( s ) (cid:1)(cid:33)(cid:43) = lim P →∞ Z − P N (cid:89) I =1 (cid:20)(cid:90) d R (1) I · · · (cid:90) d R ( P ) I (cid:21) δ (cid:32) ¯ n O ∗ − P P (cid:88) s =1 n O ∗ (cid:0) R ( s ) (cid:1)(cid:33) exp ( − βV eff ( { R } ))(9)in the PIMD formalism. The brackets indicate the ensemble average of the contained func-tion, which is equal to its time average assuming ergodicity. To further simplify the ex-pression a narrow peaked Gaussian function, with an inverse variance given as κ , is used toapproximate the delta function, giving ρ (¯ n O ∗ ) ≈ lim P →∞ Z − P (cid:114) βκ π N (cid:89) I =1 (cid:20)(cid:90) d R (1) I · · · (cid:90) d R ( P ) I (cid:21) exp ( − βV conseff ( { R } , ¯ n O ∗ )) (10)where V conseff ( { R } , ¯ n O ∗ ) = V eff ( { R } ) + κ (cid:32) ¯ n O ∗ − P P (cid:88) s =1 n O ∗ (cid:0) R ( s ) (cid:1)(cid:33) . (11)It is exceedingly difficult to calculate the free energy value itself through the above equations.Calculating the derivative of the free energy with respect to the restrained value is, however,less of a challenge. This derivative is given as f (¯ n O ∗ ) = ∂A (¯ n O ∗ ) ∂ ¯ n O ∗ = (cid:42) κ (cid:32) ¯ n O ∗ − P P (cid:88) s =1 n O ∗ (cid:0) R ( s ) (cid:1)(cid:33)(cid:43) eff , (12)where the subscript “eff” stands for the sampling by a PIMD simulation with the restraint.This derivative corresponds to the time-averaged force which will be used in the followingto calculate the free energy surface as a function of the CN.6 . Methodology for Calculating the Autoionization Constant of Water Using the time-averaged forces, f (¯ n O ∗ i ), from Equation (12) and the time-averaged CNs, n O ∗ i , obtained from a number of CN restrained simulations it is possible to estimate the freeenergy difference through thermodynamic integration,∆ A (¯ n O ∗ i ) = (cid:90) ¯ n O ∗ i ¯ n O ∗ f (¯ n ) d ¯ n. (13)The numerical integral is evaluated by spline interpolation between each of the simulatedrestrained CNs (¯ n O ∗ , . . . , ¯ n O ∗ M ). The numerical integral is here calculated using a thirdorder B-spline which passes through all the calculated points.Using the free energy surface one can employ a probabilistic (PROB) method to calculatethe p K W of water and p K A of an acid, as suggested by Davies et al. , based on the work ofChandler. This method relies on the relative probability of finding the system in a boundstate. For this purpose we define a cutoff bond distance, R c , at which the O-H bond breaksand the OH − and H O + ions are formed. The probability ratio between the bound anddissociated states is given by γ ( R c ) = (cid:82) R c exp( − β ∆ A (¯ n O ∗ ( r )))4 πr dr (cid:82) R max exp( − β ∆ A (¯ n O ∗ ( r )))4 πr dr , (14)where the factor 4 πr arises from the Jacobian of polar coodinates. The mapping ¯ n O ∗ ( r ) iscarried out by assigning the time averaged closest distance between the central oxygen andthe restrained hydrogens, r i , for the given restraint ¯ n O ∗ i . The assigned free energy differencecalculated in Equation (13), ∆ A (¯ n O ∗ ( r )), can then be linearly interpolated to numericallyevaluate the probability. R max is the time averaged distance when the CN restraint is set tothe lowest coordination number, ¯ n O ∗ M .The value of R c is determined so that it gives p K W = 14 .
00 for liquid H O. For waterautoionization H O(l) (cid:10) OH − (aq) + H + (aq) , (15)the autoionization constant isp K W = − log (cid:0) [OH − (aq)][H + (aq)] (cid:1) . (16)Rewriting this using the probabilities of water dissociation ( γ W ( R c )) from Equation (14)leads to p K W ( R c ) = − (cid:18) (1 − γ W ( R c )) N W c V (cid:19) , (17)7here N W and V are the number of water molecules and the volume of the simulation box,respectively, and c is the standard concentration (1 M). To calculate R c , the reference valueat standard conditions is used, resulting in γ W ( R c ) = 1 − − c W (18)where c W = N W c V (which is about 55.6 in ambient conditions). This can then be used to find R c for Equation (14). It is then commonly assumed that this distance is also applicable tobreaking the A-H bond to form A − and H O + for a common acid A.The p K A of an acid, or acid group of a molecule, (A) in water can be found by consideringAH(aq) (cid:10) A − (aq) + H + (aq) , (19)and the following expression for the acidity constant (p K A ),p K A = − log (cid:18) [A − ][H + ][AH] (cid:19) . (20)Here we have assumed a dilute aqueous solution, where the activities of all species are deter-mined by their concentration in the solution. The solvated protons, H + , can in principle stemfrom both water autoionization and the dissociation from A leaving us with the followingexpression for their concentration[H + ] = [H + ] acid + [H + ] water = (1 − γ A ( R c )) c A + (1 − γ W ( R c )) c (cid:48) W , (21)where c (cid:48) W = N (cid:48) W c V and c A = N A c V . N (cid:48) W = N W − N A , and N A is the number of acid moleculesreplacing water molecules in the box. It is implicitly assumed that the p K W of water isunchanged in the solution of AH, a reasonable assumption given that the solution is dilute.Using the probabilities of dissociation for the acid and water we can rewrite Equation (20)as p K A (PROB) = − log (cid:18) (1 − γ A ( R c )) c A ((1 − γ A ( R c )) c A + (1 − γ W ( R c )) c (cid:48) W ) c A γ A (cid:19) = − log (cid:18) (1 − γ A ) γ A (cid:18) (1 − γ A ) N A c V + 10 − N W − N A N W (cid:19)(cid:19) . (22)Generally one would find that (1 − γ A ( R c )) (cid:29) (1 − γ W ( R c )), so the above equation can thenbe reasonably approximated asp K A (PROB) ≈ − log (cid:18) (1 − γ ( R c )) ( γ ( R c )) N A c V (cid:19) . (23)8quation (23) was used in previous studies to predict the p K A of acidic substances wherethe proton concentration are mainly from the solute dissociation. In this study we willresort to using Equation (22), since the p K A of HDO and HTO in aqueous solution isexpected to be close to that of the solvent H O itself. This approach is however general toall very weak acids, with p K A values around their solvent. These types of acids have notpreviously been studied in this context, which is why most studies have relied on Equation(23) for calculating p K A values.Another way of calculating the probabilities used in the equations above is to employa basic two state model. That is, we assume that the minimum of the free energy surfacecorresponds to the free energy of the bound state of water and the free energy of the lowestCN restraint corresponds to the dissociated state. From this model we can formulate theprobability of finding a water molecule in a dissociated state as1 − γ = exp ( − β ∆ A )1 + exp ( − β ∆ A ) , (24)where ∆ A is taken as the difference between the maximum (dissociated state) value andthe minimum (bound state) value of ∆ A (¯ n O ∗ i ) in Equation (13). For the autoionization ofwater and its isotopologues we can approximate this probability as exp ( − β ∆ A ), as the freeenergy difference is very large. Inserting this into Equation (17) results inp K W (ABS) = − (cid:18) exp ( − β ∆ A ) N W c V (cid:19) = 2 β ∆ A ln(10) − (cid:18) N W c V (cid:19) . (25)This method allows us to compare the calculated p K W of H O to that of D O and T O, asopposed to the probabilistic method, in which the H O simulation is used as a reference. Inthe following this method will be referred to as the absolute (ABS) method, as it dependson the absolute free energy difference between two states.
III. COMPUTATIONAL DETAILS
MD and PIMD simulations of liquid H O and its isotopologues were undertaken in thecanonical ensemble at temperature 300 K. The interatomic potential and the associatedforce were computed on the fly by ab initio DFT, semiempirical DFTB, or empirical OSS2methods. The simulation conditions are listed in Table I. Systems of 32 −
64 water moleculeswere contained in cubic boxes with periodic boundary conditions. The box size was set9uch that the number density was 29.86 molecules/˚A , which amounts to 1.00 g/cm forH O. An example structure is depicted in Figure 1 along with snapshots of the constrainedwater molecule and its surroundings under different CN restraints recorded from the DFTtrajectories.The numerical integration schemes for MD and PIMD were based on the reversible refer-ence system propagation algorithm (RESPA) as implemented in the
PIMD software package. The temperature was strongly controlled by attaching massive Nos´e-Hoover chain (MNHC)thermostats to each degree of freedom in the MD simulations. The MNHC thermostatswere also attached to each normal mode representing the ring polymer in the PIMD simu-lations. The fifth-order Suzuki-Yoshida factorization was used for the numerical integrationof the MNHC thermostats. The simulations were run for 12 . − . . − .
43 fs depending on the system and the interatomic potential. The restraints witha force constant κ = 4 −
10 hartree were applied at 15 points in the range 0 . ≤ n i ≤ . VASP ). We employed the Perdew-Burke-Ernzerhof (PBE) exchangecorrelation functional and Grimme’s D3 van der Waals correction. The core electrons weretaken into account using the projector augmented wave (PAW) method. The valence elec-trons are expressed in terms of a linear combination of plane wave basis functions with acutoff at 400 eV. Only the Γ-point of the Brillouin zone was computed.The semiempirical DFTB energy calculations were carried out using
DFTB+ . We em-ployed the third-order self-consistent-charge density-functional tight-binding (SCC-DFTB) method with the 3ob Slater-Koster parameter set. The empirical energy calculations based on the OSS2 potential were implemented inan in-house version of the PIMD software package. The OSS2 potential can be categorizedas a polarizable force field composed of short-range intramolecular bonds and long-rangeinteractions between point charges and induced point dipoles with damping. The parametersin the OSS2 potential are fitted so as to reproduce the ab initio energies of neutral andprotonated water clusters at the level of second-order Møller-Plesset perturbation theory(MP2). The OSS2 potential was originally developed for small water clusters in the free10oundary condition, but it can be applied to bulk water by adaption for periodic boundaryconditions. We used the Ewald sum of point charges and induced point dipoles for this,where the point dipoles in the electrostatic field were determined by the matrix inversionmethod. The resulting restraints, forces and trajectories were analysed using a locally developedpython script, where the
MDAnalysis library was used for calculating the O-H dis-tances. The errors reported here were calculated by the method outlined by Flyvbjerg andPetersen. The errors of the cut off distance R c and the probabilistic p K A or p K W are notconsidered to be correlated, i.e. , the errors in p K W or p K A are calculated using a fixed valueof R c . All figures depicting the molecular systems were visualized using the VMD software. The numerical integrals needed for the calculations outlined in the theory section were allcarried out using linear interpolation and a step size of 1 · − in both CN and O-H distancespace. IV. RESULTS
The free energy curves, A ( n ( r )), obtained from ab initio MD and PIMD simulations aredisplayed in Figure 2(a). These represent in our opinion the most reliable results in thepresent study. The following features become clear when comparing the results presentedin Figure 2(a). The ab initio PIMD results in a much lower free energy for the dissociationprocess when compared to that of ab initio MD. This effect can be explained by the NQEs,which is expected to lead to lower free energies of dissociation due to the delocalization ofthe proton. Comparing the results for the isotopologues H O and T O we find that the freeenergy curves are different in the ab initio PIMD. The calculation of these two species usingab initio MD would result in the same curves due to the absence of NQEs in classical MD .An important consequence of this difference is that the cut off distance R c determined by abinitio MD (1 . ± .
03 ˚A) must be adapted to that by ab initio PIMD (1 . ± .
02 ˚A(H O)1 . ± .
04 ˚A(D O), see Table II or SI, respectively) in order to reproduce the referencep K W value of liquid water. We note in passing that the result obtained from ab initio MD isconsistent with values used in previous studies of 1 . − . One way of interpretingthe increase in cut off distance in the PIMD simulations is that the proton stays bonded tothe restrained oxygen longer due to NQEs. 11n Figures 2(b) and 2(c) we display the results obtained from semiempirical MD/PIMDand empirical MD/PIMD, respectively. Comparing Figures 2(b) and 2(c) with Figure 2(a),the free energy curves of semiempirical MD and PIMD look very different from that of abinitio MD and PIMD in the absolute values. Accordingly, the cutoff distances R c are quitedifferent from one another. On the other hand, the reduction of the free energy curvesbehaves similarly with respect to the NQEs. Therefore, we speculate that the isotope effectspredicted by empirical and semiempirical PIMD can be as reliable as those from ab initioPIMD.The free energy curves of the classical simulation show an anharmonic behavior sincethe mean force upon proton dissociation is a nonlinear function of r . In addition, thedifference between the free energy curves of the classical and quantum simulations vary along r , which implies the influence of anharmonicity on the NQEs. In general, the magnitude ofthe NQEs tends to be large where the potential curvature ω is larger than 1 /β (cid:126) , which canbe understood from the formula of the quantum harmonic correction to the free energy, A qhc = − β − log (cid:110) β (cid:126) ω/ β (cid:126) ω/ (cid:111) . The computational effort needed to obtain the free energy curve increases proportionallywith respect to the number of restraints. We therefore setup the ab initio MD and PIMDsimulations with a smaller system and number of beads ( N = 32 and P = 12) compared toour previous work on the same system without restraints ( N = 64 and P = 16). To verifythis setup, we checked the size- and bead-dependence of the MD and PIMD simulationsusing the semiempirical DFTB potential. Figure 3(a) shows the free energy curves obtainedfrom semiempirical PIMD simulation for a larger number of beads with P = 32, while Figure3(b) shows the free energy curves obtained from semiempirical MD and PIMD simulationfor a larger system size with 64 water molecules. It can be clearly seen that Figures 3(a)and (b) follow the same trend as Figures 2(b), in the sense that the free energy is reducedby the NQEs within the semiempirical simulations. We can therefore expect that the resultsobtained from ab initio MD and PIMD simulations shown in Figure 2(a) are reasonable withrespect to the nuclear quantum and isotope effects on the free energy curves.In Table II, we display the autoionization constants, p K W , of liquid D O and T O cal-culated using Equation (17) and the free energy curves obtained from PIMD simulations.The cutoff parameter R c was determined for a particular setup of PIMD simulations suchthat the p K W of liquid H O is 14.00. We see the trend that the p K W value of liquid T O12s higher than 14.00, which is consistent with experimental expectations. However, thedifference between the p K W values of liquid H O and D O in the ab initio PIMD simula-tions lies within the error bars. All the PIMD simulations were able to predict that thep K W value of liquid T O is larger than that of liquid H O with statistical significance. Wealso tried calculating p K W using the D O p K W as a reference, results given in Table SI ofthe supplementary information (SI), and we found that the results agree well with the onespresented in Table II.The autoionization constants calculated using the absolute method outlined around Equa-tion (25) are given in Table III. This method makes it possible to compare the p K W cal-culated using MD and PIMD simulations directly, as no reference value is required for thismethod. Comparing the ab initio MD and PIMD results for H O reveals that the NQEs playan important role in determining the correct p K W value. We note that the values predictedfrom the semi empirical MD and PIMD simulations are far from the correct value of p K W ,see Table SII in the SI, as can be expected from their free energy profiles. The results usingthe empirical OSS force field are, however, comparable with those of the ab initio method,with a slightly smaller difference between MD and PIMD results. As shown above, theabsolute method of ab initio PIMD simulation can correct for the large overestimation ofthe unbinding energy of a proton from water of ab initio MD simulation. Additionally, theabsolute method produces a reasonable result for the isotope effect for all pure isotopologuesstudied here.In Table IV, we display the acidity constants, p K A , associated with the reaction,HXO(aq) (cid:10) OX − (aq) + H + (aq) , (26)where X = D or T. The values were calculated using Equation (22) and the free energycurves obtained from PIMD simulations.Experimental values do not exist for the p K A of HDO and HTO. However, the “ratio-nal” p K A of H O, which is 15 .
74 = p K W + log ( c W ), would be a reference. The rationalp K A is obtained by using an activity for the H O molecules when the dissociating watermolecule is assumed misleadingly to be distinguishable from the rest of the water moleculesin solution.
It can however serve well as a reference in this case where the HDO andHTO molecules are in fact distinguishable from the solvent molecules.It is expected that the p K A of HDO and HTO in water are either similar to or larger13han the “rational” p K A of H O, assuming that the isotope effect follows the same trend asthe case of the p K W of H O, D O and T O.Here it turned out that thep K A value of HDO predicted from empirical PIMD is close to the “rational” p K A of H O,while the predicted p K A value of HTO is larger than that of H O with statistical significance.For the semiempirical calculations we find that both the p K A of HDO and HTO are largerthan the “rational” p K A of H O. We note that these p K A values are difficult to measureexperimentally, but they are important in determining the concentration of OD − and OT − ions in liquid water. V. CONCLUSION
In this study we established a first-principles approach to compute the autoionizationand acidity constants of water taking account of NQEs by PIMD with CN restraints. Thesimulations were carried out using different potential energy surfaces, i.e. , ab initio DFT,semiempirical DFTB and empirical OSS2 methods.The findings presented here are in line with previous empirical results. The currentstudy does differentiate itself from previous studies, by targeting the autoionization processdirectly, without any empirical factors to estimate its free energy curves.It was found that the free energy curve in the proton dissociation obtained from thequantum PIMD simulation is downshifted significantly compared with that obtained fromthe classical MD simulation, thus showing the importance of the NQEs on the autoionizationand acidity constants. The p K W values of water isotopologues, liquid D O and T O, wereestimated based on a probabilistic method using shifts in the free energy curves of D O andT O with respect to that of H O. The results agree well with experimental values, accountingfor the statistical uncertainties of our simulations.We went on to compute the p K A values of aqueous HDO and HTO molecules which aredifficult to measure experimentally. The results predict that p K A of aqueous HTO (HDO)is larger than (close to) that of H O.The work presented here opens the possibility for accurate calculations of p K A for morecomplex systems, such as small organic molecules in solution. It furthermore makes itpossible to to predict the isotope effect on these system by direct calculation. These goals14nd calculating the temperature and pressure dependence of the autoionization constant ofwater will be a subject of future studies.We finally note that the probabilistic method requires the reference p K W value of H O(14.00) while the absolute method does not. For the latter, however, an accurate estimationof absolute p K A values remains a difficult challenge. In fact the p K W estimated using theabsolute method with the present simulations have very different values to those from theprobablistic method, see Tables II and III. This is because the standard free energy curveis strongly dependent on the potential models and the system sizes. These issues should bestudied more carefully in the future. It is however clear that the inclusion of the NQEs isimportant for determining autoionization and acidity constants.In fact, we do find a difference of 4 . ± . K A units between the ab initio MD and PIMDresults, where the PIMD result is clearly closer to the true value of the p K W of water. Thisdifference does seem to be in line, to some extent, with the experimental extrapolation of 3p K A units which was suggested earlier. VI. SUPPLEMENTARY MATERIAL
See the supplementary material for Tables showing the autoionization constants of waterisotopologues calculated by the probabilistic method using the experimental p K W of D Oto calculate R c , and the autoionization constants of water isotopologues calculated by theabsolute method using the semiempirical DFTB method. VII. ACKNOWLEDGEMENTS
This work was completed under the project “Hydrogenomics” in Grant-in-Aid for Scien-tific Research on Innovative Areas, MEXT, Japan. The computations were mostly run on thesupercomputer facilities at Japan Atomic Energy Agency and the Institute for Solid StatePhysics, The University of Tokyo. M.S. thanks JSPS KAKENHI (18H05519, 18H01693,18K05208) and MEXT Program for Promoting Research on the Supercomputer Fugaku(Fugaku Battery & Fuel Cell Project) for financial support. We thank Prof. Nikos Doltsinisin Universit¨at M¨unster for his advice on the coordination number constraints. We thankDr. Alex Malins in JAEA for proofreading the text.15
III. DATA AVAILABILITY
The data that support the findings of this study are available within the article and itssupplementary material.
REFERENCES K. S. Alongi and G. C. Shields, “Chapter 8 - Theoretical Calculations of Acid DissociationConstants: A Review Article,” in
Annual Reports in Computational Chemistry (R. A.Wheeler, ed.), vol. 6, pp. 113–138, Elsevier, Jan. 2010. E. Alexov, E. L. Mehler, N. Baker, A. M. Baptista, Y. Huang, F. Milletti, J. E. Nielsen,D. Farrell, T. Carstensen, M. H. M. Olsson, J. K. Shen, J. Warwicker, S. Williams, andJ. M. Word, “Progress in the prediction of pKa values in proteins,”
Proteins , vol. 79,pp. 3260–3275, Dec. 2011. J. Ho and M. Coote, “A universal approach for continuum solvent pK a calculations: arewe there yet?,”
Theor. Chem. Acc. , 2010. J. Tomasi, B. Mennucci, and R. Cammi, “Quantum Mechanical Continuum SolvationModels,”
Chem. Rev. , vol. 105, pp. 2999–3094, Aug. 2005. S. Miertuˇs, E. Scrocco, and J. Tomasi, “Electrostatic interaction of a solute with a con-tinuum. A direct utilizaion of AB initio molecular potentials for the prevision of solventeffects,”
Chem. Phys. , vol. 55, pp. 117–129, Feb. 1981. J. B. Foresman, T. A. Keith, K. B. Wiberg, J. Snoonian, and M. J. Frisch, “Solvent Effects.5. Influence of Cavity Shape, Truncation of Electrostatics, and Electron Correlation on abInitio Reaction Field Calculations,”
J. Phys. Chem. , vol. 100, pp. 16098–16104, Jan. 1996. J. L. Pascual-ahuir, E. Silla, and I. Tu˜non, “GEPOL: An improved description of molecularsurfaces. III. A new algorithm for the computation of a solvent-excluding surface,”
J.Comput. Chem. , vol. 15, no. 10, pp. 1127–1138, 1994. S. Miertuˇs and J. Tomasi, “Approximate evaluations of the electrostatic free energy andinternal energy changes in solution processes,”
Chem. Phys. , vol. 65, pp. 239–245, Mar.1982. M. Cossi, N. Rega, G. Scalmani, and V. Barone, “Energies, structures, and electronicproperties of molecules in solution with the C-PCM solvation model,”
J. Comput. Chem. ,16ol. 24, no. 6, pp. 669–681, 2003. A. Klamt and G. Sch¨u¨urmann, “COSMO: a new approach to dielectric screening in solventswith explicit expressions for the screening energy and its gradient,”
J. Chem. Soc., PerkinTrans. 2 , pp. 799–805, Jan. 1993. A. Klamt, “Conductor-like Screening Model for Real Solvents: A New Approach to theQuantitative Calculation of Solvation Phenomena,”
J. Phys. Chem. , vol. 99, pp. 2224–2235, Feb. 1995. A. V. Marenich, C. J. Cramer, and D. G. Truhlar, “Universal Solvation Model Based onSolute Electron Density and on a Continuum Model of the Solvent Defined by the BulkDielectric Constant and Atomic Surface Tensions,”
J. Phys. Chem. B , vol. 113, pp. 6378–6396, May 2009. I. Soteras, C. Curutchet, A. Bidon-Chanal, M. Orozco, and F. J. Luque, “Extension ofthe MST model to the IEF formalism: HF and B3lyp parametrizations,”
J. Mol. Struct.(Theochem) , vol. 727, pp. 29–40, Aug. 2005. N. L. Doltsinis and M. Sprik, “Theoretical pKa estimates for solvated P(OH)5 from co-ordination constrained Car–Parrinello molecular dynamics,”
Phys. Chem. Chem. Phys. ,vol. 5, pp. 2612–2618, June 2003. M. Schilling and S. Luber, “Determination of pKa Values via ab initio Molecular Dynamicsand its Application to Transition Metal-Based Water Oxidation Catalysts,”
Inorganics ,vol. 7, p. 73, June 2019. Y.-L. Chen, N. L. Doltsinis, R. C. Hider, and D. J. Barlow, “Prediction of AbsoluteHydroxyl pKa Values for 3-Hydroxypyridin-4-ones,”
J. Phys. Chem. Lett. , vol. 3, pp. 2980–2985, Oct. 2012. J. E. Davies, N. L. Doltsinis, A. J. Kirby, C. D. Roussev, and M. Sprik, “Estimating pKaValues for Pentaoxyphosphoranes,”
J. Am. Chem. Soc. , vol. 124, pp. 6594–6599, June2002. N. Sandmann, J. Bachmann, A. Hepp, N. L. Doltsinis, and J. M¨uller, “Copper(II)-mediated base pairing involving the artificial nucleobase 3h-imidazo[4,5-f]quinolin-5-ol,”
Dalton Trans. , vol. 48, pp. 10505–10515, July 2019. A. K. Tummanapelli and S. Vasudevan, “Dissociation Constants of Weak Acids from abInitio Molecular Dynamics Using Metadynamics: Influence of the Inductive Effect andHydrogen Bonding on pKa Values,”
J. Phys. Chem. B , vol. 118, pp. 13651–13657, Nov.17014. A. K. Tummanapelli and S. Vasudevan, “Ab Initio MD Simulations of the Brønsted Acidityof Glutathione in Aqueous Solutions: Predicting pKa Shifts of the Cysteine Residue,”
J.Phys. Chem. B , vol. 119, pp. 15353–15358, Dec. 2015. A. K. Tummanapelli and S. Vasudevan, “Ab Initio Molecular Dynamics Simulations ofAmino Acids in Aqueous Solutions: Estimating pKa Values from Metadynamics Sam-pling,”
J. Phys. Chem. B , vol. 119, pp. 12249–12255, Sept. 2015. C. D. Daub and L. Halonen, “Ab Initio Molecular Dynamics Simulations of the Influenceof Lithium Bromide Salt on the Deprotonation of Formic Acid in Aqueous Solution,”
J.Phys. Chem. B , vol. 123, pp. 6823–6829, Aug. 2019. L. Xu and M. L. Coote, “Methods To Improve the Calculations of Solvation Model Den-sity Solvation Free Energies and Associated Aqueous pKa Values: Comparison betweenChoosing an Optimal Theoretical Level, Solute Cavity Scaling, and Using Explicit SolventMolecules,”
J. Phys. Chem. A , vol. 123, pp. 7430–7438, Aug. 2019. C. C. R. Sutton, G. V. Franks, and G. da Silva, “First Principles pKa Calculations onCarboxylic Acids Using the SMD Solvation Model: Effect of Thermodynamic Cycle, ModelChemistry, and Explicit Solvent Molecules,”
J. Phys. Chem. B , vol. 116, pp. 11999–12006,Oct. 2012. S. Zhang, “A reliable and efficient first principles-based method for predicting pKa values.III. Adding explicit water molecules: Can the theoretical slope be reproduced and pKavalues predicted more accurately?,”
J. Comput. Chem. , vol. 33, no. 5, pp. 517–526, 2012. B. Thapa and H. B. Schlegel, “Calculations of pKa’s and Redox Potentials of Nucleobaseswith Explicit Waters and Polarizable Continuum Solvation,”
J. Phys. Chem. A , vol. 119,pp. 5134–5144, May 2015. B. Thapa and H. B. Schlegel, “Density Functional Theory Calculation of pKa’s of Thiols inAqueous Solution Using Explicit Water Molecules and the Polarizable Continuum Model,”
J. Phys. Chem. A , vol. 120, pp. 5726–5735, July 2016. N. L. Haworth, Q. Wang, and M. L. Coote, “Modeling Flexible Molecules in Solution: ApKa Case Study,”
J. Phys. Chem. A , vol. 121, pp. 5217–5225, July 2017. E. A. Carter, G. Ciccotti, J. T. Hynes, and R. Kapral, “Constrained reaction coordinatedynamics for the simulation of rare events,”
Chem. Phys. Lett. , vol. 156, pp. 472–477, Apr.1989. 18 M. Sprik and G. Ciccotti, “Free energy from constrained molecular dynamics,”
J. Chem.Phys. , vol. 109, pp. 7737–7744, Nov. 1998. A. V. Bandura and S. N. Lvov, “The Ionization Constant of Water over Wide Ranges ofTemperature and Density,”
J. Phys. Chem. Ref. Data , vol. 35, pp. 15–30, Dec. 2005. N. Mora-Diez, Y. Egorova, H. Plommer, and P. R. Tremaine, “Theoretical study of deu-terium isotope effects on acid–base equilibria under ambient and hydrothermal conditions,”
RSC Adv. , vol. 5, pp. 9097–9109, Jan. 2015. D. W. Shoesmith and W. Lee, “The ionization constant of heavy water (D2o) in thetemperature range 298 to 523 K,”
Can. J. Chem. , vol. 54, pp. 3553–3558, Nov. 1976. M. Sprik, “Computation of the pK of liquid water using coordination constraints,”
Chem-ical Physics , vol. 258, pp. 139–150, Aug. 2000. E. Perlt, M. v. Domaros, B. Kirchner, R. Ludwig, and F. Weinhold, “Predicting the IonicProduct of Water,”
Sci Rep , vol. 7, pp. 1–10, Aug. 2017. M. ˇStrajbl, G. Hong, and A. Warshel, “Ab Initio QM/MM Simulation with Proper Sam-pling: “First Principle” Calculations of the Free Energy of the Autodissociation of Waterin Aqueous Solution,”
J. Phys. Chem. B , vol. 106, pp. 13333–13343, Dec. 2002. R. Wang, V. Carnevale, M. L. Klein, and E. Borguet, “First-Principles Calculation ofWater pKa Using the Newly Developed SCAN Functional,”
J. Phys. Chem. Lett. , vol. 11,pp. 54–59, Jan. 2020. M. Moqadam, A. Lervik, E. Riccardi, V. Venkatraman, B. K. Alsberg, and T. S. v. Erp,“Local initiation conditions for water autoionization,”
PNAS , vol. 115, pp. E4569–E4576,May 2018. P. L. Geissler, C. Dellago, D. Chandler, J. Hutter, and M. Parrinello, “Autoionization inLiquid Water,”
Science , vol. 291, pp. 2121–2124, Mar. 2001. B. L. Trout and M. Parrinello, “Analysis of the Dissociation of H2o in Water Using First-Principles Molecular Dynamics,”
J. Phys. Chem. B , vol. 103, pp. 7340–7345, Aug. 1999. M. Shiga, “Path Integral Simulations,” in
Reference Module in Chemistry, Molecular Sci-ences and Chemical Engineering , Elsevier, Jan. 2018. M. Parrinello and A. Rahman, “Study of an F center in molten KCl,”
J. Chem. Phys. ,vol. 80, pp. 860–867, Jan. 1984. R. W. Hall and B. J. Berne, “Nonergodicity in path integral molecular dynamics,”
J.Chem. Phys. , vol. 81, pp. 3641–3643, Oct. 1984.19 M. Ceriotti, M. Parrinello, T. E. Markland, and D. E. Manolopoulos, “Efficient stochasticthermostatting of path integral molecular dynamics,”
J. Chem. Phys. , vol. 133, p. 124104,Sept. 2010. M. E. Tuckerman, B. J. Berne, G. J. Martyna, and M. L. Klein, “Efficient moleculardynamics and hybrid Monte Carlo algorithms for path integrals,”
J. Chem. Phys. , vol. 99,pp. 2796–2808, Aug. 1993. M. Ceriotti, W. Fang, P. G. Kusalik, R. H. McKenzie, A. Michaelides, M. A. Morales, andT. E. Markland, “Nuclear Quantum Effects in Water and Aqueous Systems: Experiment,Theory, and Current Challenges,”
Chem. Rev. , vol. 116, pp. 7529–7550, July 2016. D. Marx, M. E. Tuckerman, J. Hutter, and M. Parrinello, “The nature of the hydratedexcess proton in water,”
Nature , vol. 397, pp. 601–604, Feb. 1999. G. Cassone, “Nuclear Quantum Effects Largely Influence Molecular Dissociation and Pro-ton Transfer in Liquid Water under an Electric Field,”
J. Phys. Chem. Lett. , vol. 11,pp. 8983–8988, Nov. 2020. D. Chandler,
Introduction to Modern Statistical Mechanics . Oxford: Oxford UniversityPress, 1987. M. Shiga, “PIMD,” 2020. S. Nos´e, “A unified formulation of the constant temperature molecular dynamics methods,”
J. Chem. Phys. , vol. 81, pp. 511–519, July 1984. W. G. Hoover, “Canonical dynamics: Equilibrium phase-space distributions,”
Phys. Rev.A , vol. 31, pp. 1695–1697, Mar. 1985. G. J. Martyna, M. L. Klein, and M. Tuckerman, “Nos´e–Hoover chains: The canonicalensemble via continuous dynamics,”
J. Chem. Phys. , vol. 97, pp. 2635–2643, Aug. 1992. G. Kresse and J. Hafner, “Ab initio molecular dynamics for liquid metals,”
Phys. Rev. B ,vol. 47, pp. 558–561, Jan. 1993. G. Kresse and J. Hafner, “Ab initio molecular-dynamics simulation of the liquid-metal–amorphous-semiconductor transition in germanium,”
Phys. Rev. B , vol. 49, pp. 14251–14269, May 1994. G. Kresse and J. Furthm¨uller, “Efficiency of ab-initio total energy calculations for metalsand semiconductors using a plane-wave basis set,”
Comput. Mater. Sci. , vol. 6, pp. 15–50,July 1996. 20 J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized Gradient Approximation MadeSimple,”
Phys. Rev. Lett. , vol. 77, pp. 3865–3868, Oct. 1996. S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, “A consistent and accurate ab initioparametrization of density functional dispersion correction (DFT-D) for the 94 elementsH-Pu,”
J. Chem. Phys. , vol. 132, p. 154104, Apr. 2010. B. Hourahine, B. Aradi, V. Blum, F. Bonaf´e, A. Buccheri, C. Camacho, C. Cevallos,M. Y. Deshaye, T. Dumitric˘a, A. Dominguez, S. Ehlert, M. Elstner, T. van der Heide,J. Hermann, S. Irle, J. J. Kranz, C. K¨ohler, T. Kowalczyk, T. Kubaˇr, I. S. Lee, V. Lutsker,R. J. Maurer, S. K. Min, I. Mitchell, C. Negre, T. A. Niehaus, A. M. N. Niklasson, A. J.Page, A. Pecchia, G. Penazzi, M. P. Persson, J. ˇRez´aˇc, C. G. S´anchez, M. Sternberg,M. St¨ohr, F. Stuckenberg, A. Tkatchenko, V. W.-z. Yu, and T. Frauenheim, “DFTB+,a software package for efficient approximate density functional theory based atomisticsimulations,”
J. Chem. Phys. , vol. 152, p. 124101, Mar. 2020. M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai, andG. Seifert, “Self-consistent-charge density-functional tight-binding method for simulationsof complex materials properties,”
Phys. Rev. B , vol. 58, pp. 7260–7268, Sept. 1998. M. Gaus, A. Goez, and M. Elstner, “Parametrization and Benchmark of DFTB3 for Or-ganic Molecules,”
J. Chem. Theory Comput. , vol. 9, pp. 338–354, Jan. 2013. L. Ojam¨ae, I. Shavitt, and S. J. Singer, “Potential models for simulations of the solvatedproton in water,”
J. Chem. Phys. , vol. 109, pp. 5547–5564, Sept. 1998. J. Sala, E. Gu`ardia, and M. Masia, “The polarizable point dipoles method with electro-static damping: Implementation on a model system,”
J. Chem. Phys. , vol. 133, p. 234101,Dec. 2010. N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein, “MDAnalysis: Atoolkit for the analysis of molecular dynamics simulations,”
J. Comput. Chem. , vol. 32,no. 10, pp. 2319–2327, 2011. R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, J. Doma´nski,D. L. Dotson, S. Buchoux, I. M. Kenney, and O. Beckstein, “MDAnalysis: A PythonPackage for the Rapid Analysis of Molecular Dynamics Simulations,”
Proceedings of the15th Python in Science Conference , pp. 98–105, 2016. H. Flyvbjerg and H. G. Petersen, “Error estimates on averages of correlated data,”
J.Chem. Phys. , vol. 91, pp. 461–466, July 1989.21 W. Humphrey, A. Dalke, and K. Schulten, “VMD: Visual molecular dynamics,”
J. Mol.Graph. , vol. 14, pp. 33–38, Feb. 1996. M. Shiga and H. Fujisaki, “A quantum generalization of intrinsic reaction coordinate usingpath integral centroid coordinates,”
J. Chem. Phys. , vol. 136, p. 184103, May 2012. M. Machida, K. Kato, and M. Shiga, “Nuclear quantum effects of light and heavy waterstudied by all-electron first principles path integral simulations,”
J. Chem. Phys. , vol. 148,p. 102324, Dec. 2017. T. P. Silverstein and S. T. Heller, “pKa Values in the Undergraduate Curriculum: WhatIs the Real pKa of Water?,”
J. Chem. Educ. , vol. 94, pp. 690–695, June 2017. E. C. Meister, M. Willeke, W. Angst, A. Togni, and P. Walde, “Confusing QuantitativeDescriptions of Brønsted-Lowry Acid-Base Equilibria in Chemistry Textbooks – A CriticalReview and Clarifications for Chemical Educators,”
Helv. Chim. Acta , vol. 97, no. 1, pp. 1–31, 2014. 22 igure Captions
Figure 1: (A) Initial structure containing 32 water molcules. The four inserts on theright show the water molecules within 4 ˚A of the central oxygen molecule during the simu-lation with a coordination number restraint of (B) 0.98, (C) 0.80, (D) 0.60, (E) 0.31. Thedistance from the central oxygen to the nearest constrained hydrogen is (B) 1.00 ˚A, (C) 1.16˚A, (D) 1.29 ˚A, (E) 1.47 ˚A. The O-H bonds are for all figures drawn if the O-H distance isless than 1.3 ˚A, with the exception of the O ∗ -H bond in which case the bonds are drawn fordistances up to 1.5 ˚A. The central oxygen is in (B-E) marked with orange, and the teal hy-drogen atom is the only hydrogen atom not subject to any constraints during the simulation.Figure 2: The free energy curves, A ( n ( r )), for H O using MD (black), H O using PIMD(orange) and T O using PIMD (green). (a) was calcultated using the ab initio potential,while (b) uses the semiempirical potential and (c) uses the empirical potential. For the MDsimulations using OSS2, DFTB and DFT we found the values of R c to be 1 . ± .
02 ˚A,1 . ± .
00 ˚A, 1 . ± .
01 ˚A, respectively. See Table II for the corresponding values for thePIMD simulations. The R c values are shown in this figure as vertical lines for H O usingMD (black) and PIMD (orange).Figure 3: The free energy curves of semiempirical simulations for H O and T O in abox containing (a) 32 and (b) 64 molecules. In both figures the black curve representsclassical MD on H O, and the corresponding black vertical line is the calculated R c distancefor this simulation, (a) 1 . ± .
00 ˚A and (b) 1 . ± .
00 ˚A, respectively. The orange andgreen curves represent H O and T O respectively in a PIMD simulation with (a) P = 32or (b) P = 12. The orange vertical lines represent the calculated R c of the two H O PIMDsimulations, (a) 1 . ± .
00 ˚A and (b) 1 . ± .
00 ˚A, respectively.23
A) (B)(D) (C)(E)
FIG. 1. F r ee E n e r g y ( k c a l / m o l ) (b) Semiemperical(b) Semiemperical(b) Semiemperical1.0 1.1 1.2 1.3 1.4 1.5R (Aangstrom)0510 (c) Emperical(c) Emperical(c) Emperical FIG. 2. F r ee E n e r g y ( k c a l / m o l ) FIG. 3. able I. Simulation setup of water isotopologues.Method System Molecules P ∆ t [ps] Length [ps] RestraintsDFT a liq. H O 32 1 0.25 24.0 15DFT b liq. H O 32 12 0.25 25.0 15DFT b liq. D O 32 12 0.25 25.0 15DFT b liq. T O 32 12 0.25 25.0 15DFTB c liq. H O 32 1 0.25 25.0 15DFTB d liq. H O 32 12 0.25 25.0 15DFTB d liq. H O 32 32 0.25 25.0 15DFTB d liq. D O 32 12 0.25 25.0 15DFTB d liq. T O 32 12 0.25 25.0 15DFTB d liq. T O 32 32 0.25 25.0 15DFTB d aq. HDO 32 12 0.25 25.0 15DFTB d aq. HTO 32 12 0.25 25.0 15DFTB c liq. H O 64 1 0.25 25.0 15DFTB d liq. H O 64 12 0.25 25.0 15DFTB d liq. T O 64 12 0.25 25.0 15OSS2 e liq. H O 64 1 0.25 12.5 15OSS2 f liq. H O 64 12 0.25 12.5 15OSS2 f liq. D O 64 12 0.35 17.5 15OSS2 f liq. T O 64 12 0.43 21.5 15OSS2 f aq. HDO 64 12 0.25 12.5 15OSS2 f aq. HTO 64 12 0.25 12.5 15 a Ab initio MD. b Ab initio PIMD. c Semiempirical MD. d Semiempirical PIMD. e Empirical MD. f Empirical PIMD. able II. Autoionization constants of water isotopologues calculated bythe probabilistic method.Method System Molecules P p K W (PROB) R c [˚A]DFT a liq. H O 32 12 14.0 e ± a liq. D O 32 12 13.5 ± ± a liq. T O 32 12 15.4 ± ± b liq. H O 32 12 14.0 e ± b liq. D O 32 12 15.0 ± ± b liq. T O 32 12 15.0 ± ± b liq. H O 32 32 14.0 e ± b liq. T O 32 32 15.9 ± ± b liq. H O 64 12 14.0 e ± b liq. T O 64 12 14.4 ± ± c liq. H O 64 12 14.0 e ± c liq. D O 64 12 15.3 ± ± c liq. T O 64 12 15.6 ± ± d liq. H O 14.00Exptl d liq. D O 14.86 Exptl d liq. T O 15.2 a Ab initio PIMD. b Semiempirical PIMD. c Empirical PIMD. d Experimental values. e Reference value to determine R c . able III. Autoionization constants of water isotopologues calculated bythe absolute method.Method System Molecules P p K W (ABS)DFT a liq. H O 32 1 19.7 ± b liq. H O 32 12 15.2 ± b liq. D O 32 12 14.8 ± b liq. T O 32 12 16.0 ± c liq. H O 64 1 15.0 ± d liq. H O 64 12 13.5 ± d liq. D O 64 12 14.9 ± d liq. T O 64 12 14.6 ± a Ab initio MD. b Ab initio PIMD. c Empirical MD. d Empirical PIMD. able IV. Acidity constants of water isotopologues, calculated using Equation (23).Method System p K A (PROP) R c [˚A]DFTB a aq. HDO 16.1 ± ± a aq. HTO 16.3 ± ± b aq. HDO 15.6 ± ± b aq. HTO 17.1 ± ± a Semiempirical PIMD of 32 water molecules with P = 12. b Empirical PIMD of 64 water molecules with P = 12.= 12.