Partial Molar Solvation Volume of the Hydrated Electron Simulated Via DFT

Different simulation models of the hydrated electron produce different solvation structures, but it has been challenging to determine which simulated solvation structure, if any, is the most comparable to experiment. In a recent work, Neupane et al. [J. Phys. Chem. B2023, 127, 5941–5947] showed using Kirkwood–Buff theory that the partial molar volume of the hydrated electron, which is known experimentally, can be readily computed from an integral over the simulated electron–water radial distribution function. This provides a sensitive way to directly compare the hydration structure of different simulation models of the hydrated electron with experiment. Here, we compute the partial molar volume of an ab-initio-simulated hydrated electron model based on density-functional theory (DFT) with a hybrid functional at different simulated system sizes. We find that the partial molar volume of the DFT-simulated hydrated electron is not converged with respect to the system size for simulations with up to 128 waters. We show that even at the largest simulation sizes, the partial molar volume of DFT-simulated hydrated electrons is underestimated by a factor of 2 with respect to experiment, and at the standard 64-water size commonly used in the literature, DFT-based simulations underestimate the experimental solvation volume by a factor of ∼3.5. An extrapolation to larger box sizes does predict the experimental partial molar volume correctly; however, larger system sizes than those explored here are currently intractable without the use of machine-learned potentials. These results bring into question what aspects of the predicted hydrated electron radial distribution function, as calculated by DFT-based simulations with the PBEh-D3 functional, deviate from the true solvation structure.


Uncertainty Calculations
Uncertainties for the calculated radial distribution functions (RDFs) and values pertaining to the partial molar volume (V M ) were calculated using block averaging to remove any serial correlations from our time-series molecular dynamics data.All block averaging uncertainty analysis was done using pyblock 1 version 0.4 with Python version 3.8.17.Reported uncertainties were those associated with the optimal block given by pyblock's find_optimal_block() functionality.When an optimal block was not found, the block uncertainty with the largest uncertainty value was reported.The convergence of the standard error as a function of block iteration is plotted in Figure S1.Additional analysis of the uncertainty using bootstrapping, shown in Figure S3, is done using the SciPy package, 2 with N = 1000 re-samples.

waters
Figure S1: Standard error as a function of block iteration for the 47-, 64-, and 128-water DFT-simulated electron-water center of mass radial distribution function bins.Each curve represents the convergence of the standard error for a single bin of the radial distribution function at the respective system size.The optimal block was chosen using pyblock's f ind_optimal_block function.Due to the limited amount of statistics available to AIMD simulations, the standard errors do not plateau for certain bins and thus serve only as a best estimate.Further details are discussed in the "Comparison of Uncertainty Measures" section below.The raw standard deviation serves as the most conservative possible limit on the magnitude of the uncertainties and indicates that our estimated standard errors from block averaging are not underestimated by an order of magnitude, so that the differences between the simulated and experimental V M 's for this system are real outside the error; see also Fig. S4.

Comparison of Uncertainty Measures
Due to the limited statistics available to ab initio molecular dynamics (AIMD) simulations, the uncertainties of some bins of our radial distribution function do not plateau in Figure S1, which means these errors could be underestimated.To determine the degree of this possible underestimation, we compare the standard error computed with block averaging, multiplied by a factor of 10, with the raw standard deviation of the counts in each bin (which serves as the upper limit of the possible uncertainty) in Figure S2.This comparison indicates that the computed uncertainties we show in the main text with block averaging are not underestimated by more than an order of magnitude.To explore this further, we provide a comparison of various methods to compute the uncertainty of the radial distribution function in Figure S3.The methods tested include block averaging, bootstrapping, confidence intervals from the Student's t-distribution (the method used by Neupane et al. 3 ), the raw standard error, and the raw standard deviation.The uncertainties presented in our main text use the 95% confidence interval from block averaging (top left panel of Fig. S3), the most conservative estimate with the exception of the raw standard deviation.Given the fact that all radial distribution functions of the hydrated electron simulated with AIMD (and ML potential extensions thereof) presented in the literature to date qualitatively resemble the ones presented here, the uncertainties we use on the main text appear to be a quite conservative estimate.was the method presented in the text, as block averaging is a known method for removing possible serial correlation in time series data, and is most pertinent to use given the limited simulation times currently accessible via ab initio DFT-based (PBEh-D3) methods.We note that even with the most generous possible uncertainty, panel (c), the calculated V M for all system sizes does not match experiment within error.

Radial Distribution Function Bin Size Analysis
Figure S5 shows the impact of bin count on the calculated partial molar volume (V M ).We note that the final value of V M changed less than 1% for a wide range of bin size choices.For the analysis presented in the main text, we chose to use bin counts of 35, 25, and 60 for the 47-, 64-and 128-water DFT-simulated radial distribution functions (RDFs), respectively, as these values were converged with the predicted V M and allowed the most statistics in each bin.

Dependence of Simulation Time on Radial Distribution Function Convergence
Since the partial molar volume is quite sensitive to changes in the RDF, we calculated RDF's using different fractions of our simulation data to better judge the convergence.Figure S6 shows RDF's for the DFT-simulated hydrated electron for each system size calculated using 25%, 50%, 75%, and 100 % of our data.The results indicate that the RDFs are converged using only half the simulation data for the 47-water simulation and 75% of the data for the 64-and 128-water simulations.The figures presented in the main text used 100% of the data.

Compressibility of Water (κ T )
Since the K-B integral used to compute V M from the RDF depends on the isothermal compressibility, it is possible to obtain different values of V M depending on whether one uses an experimental or simulated value of κ T .Figure S7  The data show that our chosen bin size is adequate for all simulation sizes and that changing the RDF bin size has a negligible impact on the calculated partial molar volume (on the order of 0.01-0.04cm 3 /mol).

Testing the Convergence of V M Calculated via the K-B Method with a Classical Aqueous Chloride Ion
To validate our work for a test case where both the experimental and simulated V M is known, we simulated a classical chloride ion with 47, 64, and 128 waters and calculated the partial molar volume using the same method we used for the hydrated electron in the main text.
This test case is especially useful as the chloride ion is an anion with a very similar RDF to that of the DFT-simulated hydrated electron (i.e., with highly structured solvent shells with similar peak heights and widths).Moreover, the V M of chloride is experimentally known to be 23.7 cm 3 /mol, while the converged simulation V M for the particular model we chose is 21.2 cm 3 /mol. 6These values are of the same sign and magnitude as the (experimental) hydrated electron V M .This provides an excellent opportunity to validate that our analysis (using the same system sizes and trajectory lengths) correctly predicts V M for a similar anion to the DFT-simulated hydrated electron.
The classical Lennard-Jones parameters for Cl − that we chose are from the work of Imai et al., 6 which predict a V M of chloride that is close to its experimental value.The water model we used in these simulation is SPC/fw with a 0.5-fs time step.We ran these simulations in the N, V, T ensemble at 298 K with box sizes, numbers of water and trajectory lengths chosen to match the AIMD simulation parameters for the hydrated electron.What follows is a reproduction of the analysis done in the manuscript for the hydrated electron but on these additional aqueous chloride simulations.
Figure S8 shows the K-B-calculated V M values for our simulated aqueous Cl − as a function of simulation size, with the dashed horizontal line indicating the fully converged simulation value from Imai et al. 6 These results indicate that the systems sizes and trajectory lengths used in our simulations are indeed sufficient to give a reasonable estimate of the true partial molar volume.
Figure S9 shows the chloride-water RDFs for all system sizes.There is an increase in first solvent shell peak intensity going from 47 to 64 waters, along with a shift in peak location for the first and second solvent shells.From 64 to 128 waters there is little change in the first solvent shell peak, but there is a shift in second solvent shell peak location.
Figure S10 shows the different solvent shell contributions to the partial molar volume of aqueous chloride, defined in the same way as in the main text for the hydrated eletron.The cavity contribution shows a monotonic increase with system size that is outside the error, while the first shell contribution shows a monotonic decrease of comparable magnitude.The  (blue), and 128 (red) water molecules using the same trajectory lengths as used for our DFT simulations of the hydrated electron.The dashed green line indicates the converged simulation value for V M . 6The partial molar volume values for each system size are the same within error, indicating a convergence of this quantity with respect to the number of waters.All of the values are slightly smaller than the fully converged simulation value, although within 10%, and the 128-water value is nearly within error of the converged value.The small slope, along with a y-intercept value that is close that of the converged simulation V M , indicates that the partial molar volume is relatively converged at these system sizes.As with the hydrated electron, the calculated V M 's are not monotonic with system size.volume should be at infinite system size.The relatively small slope of the fit line, as well as the close correspondence of our simulated V M values to the converged simulated value, indicates that the K-B-calculated V M is fairly converged at these system sizes.This analysis

Figure S2 :
Figure S2: Comparison of possible measures of the uncertainty for the 128-water DFTsimulated electron-water center of mass radial distribution function.The blue shading corresponds to the computed standard errors from block averaging multiplied by a factor of 10 (i.e., five times larger than what is shown in the main text), while the error bars represent the raw standard deviation of each bin.The raw standard deviation serves as the most conservative possible limit on the magnitude of the uncertainties and indicates that our estimated standard errors from block averaging are not underestimated by an order of magnitude, so that the differences between the simulated and experimental V M 's for this system are real outside the error; see also Fig.S4.

FigureFigure S3 :
Figure S4 panels (a-d) show the partial molar volume running integrals plotted using 4 different measures of uncertainty.Even with the largest possible choice for the error bars, the predicted V M 's from the DFT-based simulations do not agree with experiment within error.

Figure S5 :
FigureS5: Variation of the predicted partial molar volume of the hydrated electron with the bin count for the 47-, 64-, and 128-water DFT-based simulations.The number of bins used in our RDF and the calculated V M are 35, 25, and 60 for the 47, 64, and 128 water simulations, respectively.The data show that our chosen bin size is adequate for all simulation sizes and that changing the RDF bin size has a negligible impact on the calculated partial molar volume (on the order of 0.01-0.04cm 3 /mol).8

Figure S6 :Figure S7 :
FigureS6: RDFs of the DFT-simulated hydrated electron for all 3 simulation sizes as a function of the fraction of the available trajectory data used to compute the RDF.The 47-water simulation is converged using only 50% of the available data, while the 64-and 128-water simulations are clearly converged using only 75% of the data; 100% of the data was used in the main text.

Figure
FigureS8: K-B-calculated partial molar volumes for chloride simulated with 47 (black), 64 (blue), and 128 (red) water molecules using the same trajectory lengths as used for our DFT simulations of the hydrated electron.The dashed green line indicates the converged simulation value for V M .6The partial molar volume values for each system size are the same within error, indicating a convergence of this quantity with respect to the number of waters.All of the values are slightly smaller than the fully converged simulation value, although within 10%, and the 128-water value is nearly within error of the converged value.

Figure S9 :Figure S10 :Figure S11 :
FigureS9: Chloride-water RDFs for the 47-, 64-, and 128-water simulations.Going from 47 to 64 waters causes a change in first solvent shell peak height as well as slight shifts in the first and second solvent shell peak locations.From 64 to 126 waters there is little change in first solvent shell peak height, however, there is a shift in second solvent shell peak location.These changes are qualitatively comparable to those we observed in the DFT simulations of the hydrated electron, indicating comparable convergence.

Figure
Figure S11 shows an extrapolation of the calculated partial molar volume with the inverse of the simulation cell length for aqueous Cl − , giving an estimate of what the partial molar