Zifan
Ye†
a,
Cunzhi
Zhang†
a and
Giulia
Galli
*abc
aPritzker School of Molecular Engineering, The University of Chicago, Chicago, IL 60637, USA. E-mail: gagalli@uchicago.edu
bDepartment of Chemistry, The University of Chicago, Chicago, IL 60637, USA
cMaterials Science Division, Argonne National Laboratory, Lemont, IL 60439, USA
First published on 27th January 2022
Determining the electronic structure of aqueous solutions at extreme conditions is an important step towards understanding chemical bonding and reactions in water under pressure (P) and at high temperature (T). We present calculations of the photoelectron spectra of water and a simple solution of NaCl under pressure at conditions relevant to the Earth’s interior (11 GPa and 1000 K). We combine first-principles and deep-potential molecular dynamics with electronic structure calculations with dielectric-dependent hybrid functionals. These functionals are defined with a fraction of exact exchange determined from the dielectric constant of the liquid computed in extreme conditions. We find a broadening of the spectra relative to ambient conditions, particularly prominent in the merging of the two main peaks below the onset of the spectra. Furthermore we find an overall red shift at high pressure and temperature, which is however not constant over the whole energy range and varies between 1.1 and 2.4 eV. Our results also show that the anion energy levels are closer to the valence band maximum of the liquid than at ambient conditions, indicating that as P and T are increased, the defect levels of Cl− and OH− in water may eventually lie below the valence band maximum of water. Finally, we characterize the ionization potential of hydrated species deriving from rapid water dissociation, e.g. hydrated hydroxide and hydronium, and we elucidate the electronic states associated with proton transfer events at high pressure. Our results represent a first, important step in predicting the electronic properties of solutions in super-critical conditions.
While most of the current knowledge of deep-earth aqueous fluids comes from geophysical models14–23 and from vibrational spectroscopic measurements, progress has been made in recent years in understanding the properties of these fluids at the microscopic level, using first principles simulations. The latter have been used to investigate fundamental properties at high P and T (HPT) including the dielectric constant of water,24,25 ionic speciation,5,6 ion solvation26,27 and ionic conductivity of simple solutions.28–31 However, relatively little is known about the electronic structure32 of aqueous fluids under extreme conditions, which is important to understand, for instance, electron and proton transfer processes in redox reactions occurring33,34 in water and solutions in contact with rocks in the Earth’s interior.
Many electronic structure probes used at ambient conditions are difficult to employ at high temperature and pressure. For example, photoelectron (PE) spectroscopy has been widely used to probe the electronic properties of hydrogen bonded liquids, providing information on their occupied energy levels,35–37 including the ionization potentials (IP) of solvated ions,38–41e.g. hydrated OH−, H3O+, Cl−, Na+, thus elucidating their structural and chemical environment. However, obtaining PE spectra for HPT water and solutions is a challenging task not only because of the reactivity of water at high P and T, but also for the difficulty in preparing super-critical liquids in contact with vacuum and detecting PE signals with high time-resolution in the presence of rapid thermal expansion.42–44
First principles theoretical methods based on quantum mechanical calculations have been applied to predict the PE spectra of several aqueous solutions at ambient P and T (APT),40,41,45–47 providing results in excellent agreement with experiments. In particular, recent investigations have shown that a combination of first-principles molecular dynamics (FPMD) and density functional theory calculations with dielectric-dependent hybrid (DDH) functionals48 is a robust, predictive approach to study photoelectron spectra of aqueous fluids.40,41 It is hence interesting to explore how to predict PE under extreme conditions and gain insight into the electronic properties of aqueous fluids at high P and T.
Here, we generalize the approach recently used to compute PE spectra at ambient conditions to the study of aqueous fluids under pressure, at high temperature. We carry out first principles simulations with semilocal and hybrid functionals, and we consider two specific systems at conditions (11 GPa and 1000 K) relevant to the Earth upper mantle, where pressure can reach ∼13 GPa and temperature ∼1700 K;49,50 in particular we investigate pure water and a 0.68 M NaCl solution, whose structural and bonding properties have been previously studied.27,31 We report PE spectra on an absolute scale, by using a technique to refer energy levels to vacuum which is appropriate for hot compressed fluids; we discuss results obtained with dielectric-dependent hybrid functionals, and with a fraction of exact exchange determined for the specific conditions studied here.32 We also analyze the contribution of dissociated species to the PE spectra, in particular hydroxide (OH−), hydronium (H3O+) ions and H4O2, which play an important role in proton transfer processes and acid–base chemistry,51–53 and contribute to the increased ionic conductivity31 of water under extreme conditions.
The rest of the paper is organized as follows: the methods used here are described in the next section, followed by a discussion of our results in Section 3. Our conclusions are presented in Section 4.
To predict PE spectra comparable to experiments, the calculation of absolute energies referenced to vacuum is required, and the strategy chosen to do so in our work is described in detail in the next section.
![]() | (1) |
Fig. 1 illustrates three possible strategies to obtain ΔV, which are summarized below.
• Vacuum: a water sample at a chosen, initial density is in contact with vacuum within a periodically repeated slab that is equilibrated using MD at a given temperature. This scheme was used for calculations of PE spectra at APT and is applicable if the density of water is such that immediate evaporation is not observed, over the time scale of the simulations. Hence it is not applicable to simulate water at high P and T, where the sample would tend to expand into vacuum over short time scales.
• Wall: similar to the previous approach, a water sample is in contact with vacuum within a periodically repeated slab that is equilibrated using MD at given temperature. However in this case a repulsive potential is added to the Hamiltonian so as to maintain the water sample at a density corresponding to the desired pressure. Such a repulsive potential prevents water from evaporating into the vacuum portion of the slab. This method can in principle be used at any P–T conditions, however the results may be affected by the choice of the specific form of the repulsive potential and the approach requires the use of relatively large systems to obtain a converged electrostatic potential.
• Slice: a relatively large bulk water model is equilibrated in a periodically repeated supercell at desired P and T conditions. A slice of the sample is then extracted after equilibration and placed in contact with vacuum, with no further MD simulations of the whole slab. To avoid dangling bonds at the interface, when taking the slice, we retain the H–O covalent bonds. This scheme is applicable to any P–T conditions, provided a slice of sufficient thickness is chosen.
To reduce the computational cost in preparing the samples required to obtain ΔV, we carried out simulations with empirical potentials first, for benchmarking purposes, and then we used the deep-MD potential61,62 implemented in the DeePMD-kit package.63 The deep-MD potential (DP)64 is trained on first principle data obtained with the SCAN functional65 at HPT and allows for the proper description of dissociated species at extreme conditions, unlike empirical potentials. In addition, using the DP potential, simulations with cells larger than those affordable with FPMD can be carried out and longer time scales can be sampled with first-principle accuracy, as described below.
In order to assess the robustness of our methodology, we first carried out simulations with the TIP3P66 potential and we compared results obtained with the Vacuum and Slice procedures at ambient conditions; we also compared our results with those previously reported in the literature. We then further compared the results obtained with the Wall and Slice methods at HPT. MD simulations were carried out in the NVT ensemble with the LAMMPS package.67,68 A timestep of 2 (1) fs was used for TIP3P simulations at APT (11 GPa & 1000 K). When using the Vacuum method, our samples consisted of 108 water molecules interfaced with vacuum in a cell of dimensions 12.77 × 12.77 × 80 Å. When using the Slice method, we equilibrated water samples with 435 and 683 molecules at APT and 11 GPa & 1000 K, respectively, corresponding to the densities of 1.0 and 1.57 g cm−3 in a cell of dimensions 12.77 × 12.77 × 80 Å. The thickness of a slice was chosen to be 20 (18) Å at APT (11 GPa & 1000 K). When adopting the Wall method, we considered 300 water molecules at 11 GPa & 1000 K, and the wall/lj126 scheme implemented in the LAMMPS code. The distance between two virtual walls simulating a repulsive potential in the Hamiltonian was set to be 39 Å. The repulsive part of Lennard-Jones interactions between the wall and water was described by using carbon–oxygen Lennard-Jones parameters69 to mimic the interaction with water oxygen atoms only. For all the procedures adopted here, the MD simulations with the TIP3P potential were performed for 4 ns following a 2 ns equilibration, with 200 equally-spaced snapshots extracted every 20 ps.
The results for the average potential obtained with the PBE functional on the configurations extracted from our MD simulations are presented in Fig. 2. Our results with the Vacuum method are nearly identical with those reported using the same scheme in the literature: 3.53 eV versus 3.54 eV.40 When using the Slice method we obtain a value of 3.36 eV and we consider the difference relative to the result with the Vacuum method (0.17 eV) to be negligible, given that the position of the valence band maximum (VBM) of water is at approximately 10 eV (ref. 39–41 and 46) relative to vacuum. The difference of 0.17 eV may be ascribed to the specific surface of the chosen slice (the statistical error in our evaluation of ΔV is smaller than 0.1 eV).
![]() | ||
Fig. 2 Plane-average electrostatic potential ![]() |
At 11 GPa and 1000 K, as reported in Fig. 2(b), ΔV is computed to be 5.43 (Slice) and 5.58 (Wall) eV. Again, we consider the difference of 0.15 eV to be small for the purpose of our study.
Next we compared the results obtained with the Slice method and the TIP3P potential with those obtained with the DP potential. When using the DP potential, we performed 210 independent MD simulations (with a time step of 0.25 fs) each of 13 ps long, at the end of which we extracted configurations for our electronic structure calculations. We obtained ΔV = 5.53 eV, in good agreement with the value computed with TIP3P, indicating that the value of ΔV depends rather weakly on subtle structural differences in the liquid samples. These comparisons between different approaches gave us confidence that the Slice method is a reliable technique to use both at APT and HPT.
Therefore, in the following we compute ΔV with the Slice method and the DP potential to obtain the spectra at HPT and to discuss the physical properties of water and the NaCl solution. We note that when using the DDH functional, we could only carry out calculations of ΔV for a small number of configurations (5), much smaller than with the PBE functional. Hence we estimated the difference between ΔV computed with PBE and DDH for those limited number of configurations, and we applied the same difference between the DDH and PBE result (5.31 vs. 5.53 eV) to all eigenvalues obtained at the DDH level of theory for the bulk samples.
![]() | ||
Fig. 3 Photoelectron spectra of a 1 M NaCl solution at ambient conditions (APT)40 and a 0.68 M NaCl solution at high pressure conditions (HPT) (11 GPa & 1000 K) obtained with the PBE (a) and DDH (b) functionals. The intensities are rescaled to those of the water 1b1 peak; the shaded areas show the distribution of ionization potentials of solvated Cl− and Na+ ions. |
Method | 2p (Na+) | 2a1 | 3s (Cl−) | 1b2 | 3a1 | 1b1 | 3p (Cl−) |
---|---|---|---|---|---|---|---|
a Results at APT of a 1 M NaCl solution from ref. 40. b Results at 11 GPa & 1000 K of a 0.68 M NaCl solution. | |||||||
PBE-APTa | 26.90 | 24.71 | 18.06 | 12.63 | 9.03 | 6.99 | 6.20 |
DDH-APTa | 34.04 | 32.57 | 24.15 | 17.14 | 13.56 | 11.34 | 9.46 |
PBE-HPTb | 27.02 | 24.81 | 18.62 | 12.84 | 8.98 | 7.12 | 6.28 |
DDH-HPTb | 31.92 | 30.18 | 22.93 | 15.98 | 12.01 | 10.02 | 8.41 |
The four major peaks (labeled, in order of increasing BE as: 2a1, 1b2, 3a1, and 1b1) are notably broadened at HPT, with the 3a1 and 1b1 features strongly overlapping with each other, consistent with larger structural fluctuations in the HPT samples. The DDH results show that the relative spacing of the water 2a1 band with respect to the 1b1 peak is decreased by ∼1.1 eV at high P, while that of 1b2 and 3a1 bands are less affected at extreme conditions. Interestingly, calculations with the DDH functional predict a reduction of the BE by roughly 1.1–2.4 eV at high pressure, depending on the character of electronic states; however calculation with the PBE or standard hybrid functionals such as PBE0 (ref. 70 and 71) yield negligible differences between the band positions in the PE spectra at high and low pressure. This finding is important and highlights the need to correctly take into account screening effects at HPT and correspondingly adjust the amount of exact exchange used to define hybrid functionals (∼60% in DDH at APT;40,41 ∼42% in DDH at the HPT conditions of this study, corresponding to the dielectric constants ε∞ = 1.78 and 2.37, respectively).
The ionization potentials (IPs) of hydrated Cl− and Na+ at APT and HPT are indicated in Fig. 3 as shaded areas. Similar to the BE of water electronic states, the IPs of ions are underestimated by PBE, relative to DDH, and DDH results show a red shift with increasing P and T which is not reproduced by the PBE functional. In addition, when using the DDH functional, the relative positions of IPs referenced to water is slightly different at low and high P. In particular, the separation between the 3p-level of Cl− and the 1b1 peak of water is decreased by 0.27 eV, suggesting a modified electronic structure near the band edges of water at high P, which in turn may affect the chemical reactivity of solutions under extreme conditions. Our results hint at a possible tendency of anion levels to move below the valence band of water as a function of increasing pressure, at high T. We note that the BE of water and IP of ions are not sensitive to the Cl−–Na+ distance, hence the results reported here should be representative of weakly concentrated solutions.
In our FPMD simulations, we directly observe frequent dissociation events at HPT, indicating that water itself becomes an electrolyte. The short-lived dissociated products27,31 are characterized using a cutoff distance for O–H bonds of 1.25 Å. A detailed analysis reveals that about 2% of water molecules are dissociated, into species including solvated HO−, H3O+ and H4O2. These species play a critical role in proton transfer events and acid–base reactions in water.
As shown in Fig. 5(a), the signal coming from ionic species resulting from water dissociation spans the entire spectrum range. The signal is weak and expected to become more prominent at higher pressure and temperature, where the concentration of dissociated species will increase. The contributions to the overall PE signal of hydrated HO− and H3O+ are presented in Fig. 5(b) and (c), and it has been identified by carrying out projected EDOS calculations. Similar to ambient conditions, the solvated HO− level is just above the VBM of water, and it originates from two degenerate lone-pair orbitals (1e1) localized on the oxygen atom. The solvated H3O+ levels lie instead close to the high-energy edge of the 1b2 feature of bulk water and originates from the ionization from two degenerate covalent orbitals with mainly p-character (1e). However the IPs of hydrated HO− and H3O+ at HPT, peaking at 8.9 and 17.75 eV, are red-shifted relative to ambient conditions, where they are centered at ∼9.2 and ∼20 eV (∼9.99 and ∼19.01 eV) in the experiments of ref. 39 (calculated by many-body perturbation approach46). In addition, the separation between the 1e1 level of HO− and 1b1 peak of water is ∼1.91 eV at APT based on experiment,39 while we find that at HPT it is 1.12 eV, when using DDH. Hence our simulations show that similar to the case of the Cl−, the level of HO− gets closer to the VBM of water as pressure is increased.
To understand the electronic states involved in proton transfer events at high pressure, we plot the IP distributions for three hydrated intermediate complexes H3O2−, H5O2+ and H4O2 in Fig. 5(b–d). With respect to HO−, the IP of H3O2− is blueshifted with a slightly enhanced signal near the 1b2 region of bulk water, similar to what is found in proton transfer events at APT.46 Compared to H3O+, we observe a red shift in the IP of H5O2+ with increased spectral weights near the 1b1 region of bulk water, again similar to what is observed in proton transfer events at APT.46 For solvated H4O2 cluster referenced to bulk water, the major change at high pressure is the double-peak near the 2a1 region of bulk water, giving rise to the 2a1 peaks with lower or higher IP in the resulting species HO− and H3O+. The results reported here for ionization potentials of dissociated species in water may help, at least qualitatively, to understand proton transfer reactions at high pressure.
Footnote |
† These two authors contributed equally, listed in alphabetical order. |
This journal is © The Royal Society of Chemistry 2022 |