Nuclear Quantum Effects in Sodium Hydroxide Solutions from Neural

Subscriber access provided by UNIV OF LOUISIANA

B: Liquids, Chemical and Dynamical Processes in Solution, Spectroscopy in Solution

Nuclear Quantum Effects in Sodium Hydroxide Solutions from Neural Network Molecular Dynamics Simulations Matti Hellström, Michele Ceriotti, and Jörg Behler J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b06433 • Publication Date (Web): 18 Oct 2018 Downloaded from http://pubs.acs.org on October 24, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Nuclear Quantum Effects in Sodium Hydroxide Solutions from Neural Network Molecular Dynamics Simulations Matti Hellström,†,‡ Michele Ceriotti,¶ and Jörg Behler∗,†,‡ †Lehrstuhl f¨ ur Theoretische Chemie, Ruhr-Universit¨ at Bochum, 44780 Bochum, Germany ‡Universit¨ at G¨ ottingen, Institut f¨ ur Physikalische Chemie, Theoretische Chemie, Tammannstr. 6, 37077 G¨ ottingen, Germany ´ ¶Laboratory of Computational Science and Modeling, Institute of Materials, Ecole Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland E-mail: [email protected] Phone: +49(551) 39-23133

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Nuclear quantum effects (NQEs) cause the nuclei of light elements like hydrogen to delocalize, affecting numerous properties of water and aqueous solutions, such as hydrogen-bonding and proton transfer barriers. Here, we address the prototypical case of aqueous NaOH solutions by investigating the effects of quantum nuclear fluctuations on radial distribution functions, hydrogen bonding geometries, power spectra, proton transfer barriers, proton transfer rates, water self-exchange rates around the Na+ cations, and diffusion coefficients, for the full room-temperature solubility range. These properties were calculated from classical and ring polymer molecular dynamics simulations employing a reactive highdimensional neural network potential based on dispersion-corrected density functional theory reference calculations. We find that NQEs have a very small impact on the solvation structure around Na+ , slightly strengthen the water-water and water-hydroxide hydrogen bonds, and lower the peak positions in the power spectra for the HOH bending and OH stretching modes by about 50 cm−1 and 100 cm−1 , respectively. Moreover, NQEs significantly lower the proton transfer barriers, thus increasing the proton transfer rates, resulting in an increase of the diffusion coefficient in particular of OH – , as well as a decrease of the mean residence time of molecules in the first hydration shell around Na+ at high NaOH concentrations.


Page 2 of 34

Page 3 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

Introduction

Nuclear quantum effects (NQEs), such as tunneling, nuclear delocalization, and zero-point energy, can influence the structure and dynamics of systems containing light elements like hydrogen. Due to the high hydrogen content, water and aqueous systems are among the most important examples for the influence of NQEs. 1 One important concept to rationalize the impact of NQEs on the hydrogen-bond network of water is that on the one hand, the hydrogen bonds become stronger as a result of proton delocalization along the direction of the bond, but on the other hand, the hydrogen bonds become weaker from proton delocalization perpendicular to the bond. 2–4 Which of these two effects is dominant can be understood based on simple models, 5 and depends primarily on the intermolecular distance between the hydrogen-bonding partners. This balance is sometimes referred to as the “competing effects” of NQEs in water. 1 In molecular simulations, it is possible to explicitly model nuclear quantum effects, i.e., to model the nuclei as quantum particles, using path integrals. The simulation techniques have progressed such that the explicit treatment of NQEs, via for example ring polymer molecular dynamics (RPMD), can now routinely be employed, 6 albeit still at an increased computational expense compared to modeling classical nuclei. Several simulations based on path integral molecular dynamics techniques have investigated the influence of NQEs in pure water. 2,4,7–13 Unfortunately, to date much less is known about NQEs in the context of aqueous electrolyte and salt solutions. Wilkins et al. 14 found that NQEs only have a weak effect on the hydration of LiF. In another study, Videla et al. 15 explored hydrogen bonding to halide ions dissolved in water, and found that NQEs favor HOD · · · F – over DOH · · · F – , but the opposite result was found for the other halide ions. Two ions for which NQEs are particularly pronounced are the hydronium (H3 O+ ) and hydroxide (OH – ) ions. 16–18 These species can diffuse in water via the Grotthuss mechanism, 17 in which protons are transferred to or from neighboring water molecules. Such proton transfer reactions can be strongly influenced by NQEs. Indeed, NQEs have also been shown to have a decisive impact on proton diffusion not only in aqueous solution 19 but also in water wires. 20,21 In the present work, we use RPMD and classical molecular dynamics (MD) simulations to identify the influence that NQEs have on the structure and dynamics of NaOH solutions of different concentrations. Dilute and more concentrated NaOH solutions are among the most widely used chemical reagents, and have been characterized in classical molecular simulations using a variety of ab initio, 16,22–28 force-field, 29–32 and machine learning potential 33–35 techniques. Here, we use a high-dimensional neural network potential 33,36 fitted to reference data calculated using the dispersion-corrected generalized-gradient-approximation (GGA) DFT functional RPBE-D3. 37,38 Both the proton transfer mechanism, 26 as well as the influence of NQEs, 10



can be sensitive to the choice of reference method. We have previously shown that our neural network potential works well for classical proton transfer reactions in NaOH solutions; 33 here, we will evaluate the effects of NQEs. Two particularly interesting features of the hydroxide ion are its diffusion mechanism and vibrational fingerprint in vibrational (IR and Raman) spectroscopy. A number of computational works have studied one or both of these properties for hydroxide ions in various environments: in bulk water, 16,19,26,27,39,40 in water nanodroplets, 41–44 in water wires, 20,45 at solid/liquid interfaces, 46–51 in nanoconfinement between two surfaces, 52 and for different concentrations of electrolytes, e.g., NaOH. 23,28,33 Both the proton transfer mechanisms and vibrational frequencies are intimately coupled to the hydrogen bonding environment around the OH – ion. For example, in bulk water at low concentrations, OH – ions are normally found in a “hypercoordinated” state, in which they accept four hydrogen bonds and donate none. 16,40,53 Through thermal fluctuations, the OH – can lose one of its accepted hydrogen bonds and instead donate a hydrogen bond, which activates the OH – for proton transfer (lowers the proton transfer barrier). 16,26,33,40 This mechanism is often referred to as a presolvation mechanism for proton transfer. For OH – ions in other environments, e.g., at a solid/liquid interface 49 or in highly concentrated NaOH solutions, 33 the “default” hydrogen-bonding environments of hydroxide ions and water molecules can be very different, and, consequently, other types of hydrogen bond fluctuations may aid the proton transfer reactivity. For example, we recently showed that for high concentrations of NaOH(aq), many of the water molecules are coordinated to one or more Na+ and therefore mostly accept only one or zero hydrogen bonds. Fluctuations which cause these water molecules to accept more hydrogen bonds are more important for proton transfer events than fluctuations around OH – . 33 A possible experimental probe for these hydrogen-bond fluctuations is spectroscopic. The OH – (aq) stretching frequency is higher than the frequencies for the water stretching modes. 28,54,55 This is, in part, because the hydroxide ion is quite a bad hydrogen bond donor. Simulations have shown that also the number of accepted hydrogen bonds determines the stretching frequency, with the frequency becoming red-shifted (lower) the more hydrogen bonds are accepted. 42,44 Finally, another process in which NQEs may play a role is the water exchange mechanism around the Na+ ions, i.e., events in which one water molecule leaves the first hydration shell and another takes its place. We recently studied this process for different concentrations of NaOH using classical MD simulations, 35 and demonstrated the existence of a “proton transfer pathway” for water exchange, in which water molecules in the first hydration shell transfer a proton to a neighboring OH – ion, thus forming a labile direct-contact ion pair Na+ – OH – , which we found can dissociate very easily.


Page 4 of 34

Page 5 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1:

Aqueous NaOH solutions modeled in this work: mole fraction xNaOH , weight-% (w-%) NaOH, number of NaOH and H2 O molecules N (NaOH) and N (H2 O) in the periodic cubic simulation box with side length L, and the shear viscosities η. xNaOH

w-%

0.016 3.5 0.049 10 0.085 17 0.143 27 0.185 34 0.255 43 0.333 53 a Experimental values at

η (mPa·s)a 8 496 25.35 1.1 24 464 24.84 1.7 40 432 24.39 2.9 64 384 23.76 7.8 80 352 23.37 15 104 304 22.81 34 128 256 22.31 56 T = 298 K interpolated and extrapolated from Ref. 60 N (NaOH)

N (H2 O)

L (˚ A)

The present work addresses the influence of NQEs on the structure and dynamics of aqueous NaOH solutions of different concentrations. NQEs have previously been shown to lower the proton transfer barrier in dilute OH – solutions, 16 but the extent to which this effect depends on the concentration is still unknown. We will address that issue, and, in addition, we will characterize how NQEs influence other fundamental properties of these solutions, such as radial distribution functions, hydrogen-bonding environments around water and hydroxide, power spectra (vibrational density of states, VDOS), hydroxide ion lifetimes (proton transfer rates), water residence times in the first solvation shell of Na+ ions, as well as diffusion coefficients. We will show that NQEs in general favor hydrogen-bonding environments that are more active for proton transfer, lower the OH stretching frequencies, and increase the OH – diffusion coefficient by about 80% for all concentration. At high concentrations, where the proton transfer pathway for water exchange is particularly prominent, NQEs also lower the mean residence time of molecules in the first solvation shell of Na+ .

2 2.1

Methods Molecular dynamics simulations

We performed classical and quantum molecular dynamics simulations of NaOH(aq) solutions of different concentrations at T = 300 K. The potential energy surface was described by a high-dimensional neural network potential previously developed for NaOH solutions for the entire room-temperature solubility range. 33 Neural networks constitute one of the most versatile machine learning techniques, and can be fitted to reproduce ab-initio-calculated potential energy surfaces very accurately. 36,56 In our previous work, 33 the neural network was parameterized to reproduce the dispersion-corrected density functional RPBE-D3; 37,38 this DFT functional, and the closely related revPBE-D3 functional, have previously been shown to yield a good description of liquid water. 10,57–59



Page 6 of 34

We express the concentration of NaOH(aq) in terms of the mole fraction xNaOH ,

xNaOH =

N (NaOH) N (NaOH) + N (H2 O)

,

(1)

where N (NaOH) is the number of NaOH formula units, and N (H2 O) is the number of H2 O molecules in the periodic simulation box. Table 1 gives N (NaOH) and N (H2 O) at each of the different mole fractions considered in this work up to the solubility limit of about 20 moles per liter (xNaOH ≈ 0.333), together with the cubic periodic box length L, which determines the density. The value of L was taken from Ref. 33, where the density had been equilibrated from classical MD simulations in the N pT ensemble at p = 1 bar, T = 300 K. For the classical simulations, we reused the trajectories from Ref. 35, where for each concentration, at least 20 independent 1 ns long N V E simulations were initialized from structures and velocities obtained from an N V T simulation (at T = 300 K). Here, we have performed additional analyses of those trajectories. In order to explicitly model nuclear quantum effects, we employed the methodology of thermostatted ring polymer molecular dynamics (TRPMD) 61 with 16 beads (see section 3.1 for a convergence test made with respect to the number of beads), together with the PILE-G thermostat. 62 The TRPMD simulations were at least 200 ps long, with a timestep of 0.25 fs. Figure 1 shows a snapshot of a part of the simulation box for xNaOH = 0.185, with all 16 beads depicted simultaneously. The simulations were performed using a custom module 58 for neural network potentials implemented in LAMMPS, 63 used as a driver for the i-PI universal force engine. 64,65

Figure 1: Snapshot of a ring polymer molecular dynamics simulation for a NaOH(aq) solution with a mole fraction xNaOH = 0.185. The positions of each of the 16 beads for each atom are shown. Colors: Na black, hydroxide O orange, water O red, H white. The image was rendered using VMD. 66,67


Page 7 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2.2

Trajectory analysis

In the TRPMD simulations, we used the positions of the individual beads to calculate radial distribution functions and structural observables, but the positions of the centroids for all other types of analysis. Each H was assigned to be covalently bound to its nearest O atom; O atoms that coordinated only one H were considered to be hydroxide oxygen atoms, and O atoms that coordinated two H atoms water oxygen atoms. When this analysis was applied to the positions of the individual beads, some of the O atoms were also considered to be hydronium (H3 O+ ) oxygens resulting from transient spontaneous water autolysis events (2 H2 O −−→ OH – + H3 O+ ), similar to the observations made in Ref. 4. These events were quite rare; for example, at the lowest NaOH concentration xNaOH = 0.016, at least one bead for about 0.09% of O atoms would be classified as hydronium. Consequently, OH – ions formed through such autolysis events were also rare, and no distinction in the analysis of, for example, radial distribution functions, was made between the “explicit” OH – ions and those formed through autolysis events. We applied a geometric definition of hydrogen bonds, considering a hydrogen bond as a donor-H-acceptor group Od Hd · · · Oa for which the distance Od Oa < 3.5 ˚ A and the angle 6 Oa Od Hd < 30◦ . 68 The proton transfer free-energy barrier, ∆F ‡ , was evaluated from the proton transfer free-energy landscape as a function of the one-dimensional proton transfer coordinate δmin . For each hydrogen bond that an OH – accepts, the proton sharing coordinate δ is calculated as the difference between the Hd · · · Oa and Od Hd distances. The smallest value of the proton sharing coordinate is the proton transfer coordinate δmin for that particular OH – ion. The distribution of δmin -values is then binned into a histogram W (δmin ), and the relative Helmholtz free energy ∆F calculated as

∆F = −kB T ln W (δmin ).

(2)

Four different types of time correlation functions were calculated, in order to estimate the power spectra, hydroxide ion lifetimes (the inverse of the proton transfer rate constants), mean residence times for molecules in the first solvation shell around Na+ ions, and diffusion coefficients. We calculated both “total” and “partial” power spectra. The total power spectra were calculated as the Fourier transform of the velocity autocorrelation function, which was calculated over 5 ps intervals as an average over all time origins. In order to compare the vibrational frequencies of molecules in different hydrogen bonding environments, we also calculated partial power spectra, by calculating the velocity autocorrelation function only for H atoms fulfilling particular criteria, such as belonging to a hydroxide ion accepting 3 hydrogen bonds, at the time origin. The condition was put only at the time origin, meaning that the hydrogen-bonding environment



Page 8 of 34

could change during the correlation time. In order to minimize the effect of hydrogen bonding fluctuations during the correlation time on the partial power spectra, the correlation time was chosen to be rather short, namely 1 ps. The hydroxide ion lifetimes were calculated using the “stable states picture” (SSP), 69 where the time correlation function C(∆t) gives the probability that a “stable” hydroxide ion at time t = t1 has not turned into a “stable” water molecule at any time in the interval t1 ≤ t ≤ t1 + ∆t. This treatment minimizes the influence of proton “rattling”, in which the proton repeatedly jumps back and forth between the same two O atoms before settling on one of them. In this work, we define a stable water molecule to be an O that coordinates two H within a sphere of radius 1.1 ˚ A, and a stable hydroxide ion to be an O that coordinates only one H within a sphere of radius 1.4 ˚ A. The characteristic lifetime τPT is extracted from C(∆t) by fitting the function to a double exponential of the form C(∆t) = A exp(−∆t/τ1 ) + (1 − A) exp(−∆t/τ2 ), and using the weighted average τPT = Aτ1 + (1 − A)τ2 . For the calculation of the mean residence time (MRT) of a molecule in the first solvation shell around Na+ , we similarly employed an SSP-type time correlation function, which gives the probability that a molecule which is “definitely” within the first solvation shell of a particular Na+ ion (Na+ – O distance < 2.7 ˚ A) at time t = t1 , has not “definitely” entered the second solvation shell of the same Na+ (Na+ – O distance > 3.6 ˚ A) at any time in the interval t1 ≤ t ≤ t1 + ∆t. These are the same criteria that we used in Ref. 35. The mean residence time, τexchange , is calculated from the time correlation function in the same way as τPT . The mean squared displacement (MSD) was calculated for the element I according to

MSDI (∆t) = h|ri (t1 + ∆t) − ri (t1 )|2 i

(3)

where the average h·i is taken over all time origins t1 and over all atoms i of element I. We also calculated the mean squared displacement for OH – , taking both mass transport and proton transfer events into account. The position of a hydroxide ion was followed from one timestep to the next, by minimizing the sum of squared distances moved by all hydroxide ions, i.e., the positions of the hydroxide oxygens at timestep t were assigned so that the sum N (NaOH)

J=

X

|ri (t) − ri (t − 1)|2

(4)

i=1

was minimized. The minimal sum was found using the Hungarian algorithm. 70 By fitting a straight line to the MSD in the interval 4 ps ≤ ∆t ≤ 20 ps, the diffusion coefficient D(L),


Page 9 of 34

where L refers to the side length of the cubic periodic simulation box (Table 1) was calculated as

D(L) =

slope(MSD) 6

(5)

The estimated diffusion coefficients that would be obtained for an infinitely large periodic box, D(∞), was calculated by the following finite-size correction: 71,72

D(∞) = D(L) + ξ

kB T 6πηL

(6)

where ξ ≈ 2.837 for a cubic simulation box, 71,72 and where we used the experimentally measured shear viscosities η given in Table 1.

3 3.1

Results Convergence with respect to the number of beads (a) gOH(r)

2

1.5

5 0

1

0.8

1

16 beads 64 beads

7 Intensity (arb. units)

10

(c) Power spectra

8

16 beads 64 beads

1.2

0.5

6 5 4 3 2 1 0

0

1

2

3

4

0

500

1000 1500 2000 2500 3000 3500 4000 ν (cm−1)

r (Å) (b) gHH(r)

−6

16 beads 64 beads

1.6 1.4

ΔF/kBT (arbitrarily shifted)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1.2 1 0.8 0.6 0.4 0.2 0

2

3

4

5

6

16 beads 64 beads

−7 −8 −9 −10 −11

1

(d) Proton transfer free energy landscape

−0.6

r (Å)

−0.4

−0.2

0 0.2 δmin (Å)

0.4

0.6

Figure 2: Results from the convergence test with respect to the number of beads in the ring polymer molecular dynamics simulations modeling quantum nuclei at xNaOH = 0.085. (a) gOH (r), with the inset showing the region between 0.8–1.2 ˚ A, (b) gHH (r), (c) power spectrum, and (d) one-dimensional proton transfer free energy landscape (calculated from the position of the centroid).



We performed convergence tests with respect to the number of beads for one of the intermediate concentrations, xNaOH = 0.085, with up to 64 beads. Those simulations were performed in a box one eighth the size of the corresponding production simulation in Table 1 (the simulation box size L was halved in each of the three spatial dimensions). The first 10 ps were used for the equilibration period and the next 65 ps were used to calculate the properties shown in Figure 2, namely, the O-H radial distribution function (Figure 2a), the H-H radial distribution function (Figure 2b), the total power spectrum (Figure 2c), and the proton-transfer free energy landscape calculated as a function of the one-dimensional proton transfer coordinate δmin (Figure 2d). Results are shown for ring polymers containing 16 beads and 64 beads. The results obtained with the different number of beads are only marginally different; with 16 beads, the first peaks in gOH (r) and gHH (r) are slightly higher than with 64 beads, the HOH bending frequency is blue-shifted by about 10 cm−1 and the OH stretching frequency by about 5 cm−1 as compared to 64 beads, and the proton transfer free energy landscapes are practically identical. Because of the very similar results obtained with 16 and 64 beads, the production simulations were carried out using 16 beads. We expect the errors associated with this choice of the number of beads not to affect the qualitative results discussed in this work.

3.2

Radial distribution functions

The radial distribution functions gNaO (r), gNaOOH− (r), gOO (r), gOH (r), and gHH (r), as calculated for all systems in Table 1 using both classical and quantum nuclei, are given in tabular format in the Supporting Information. Here, we just highlight the most interesting features, using the three mole fractions xNaOH = 0.049, 0.185, and 0.333 as examples covering the full concentration range. Figure 3a–c reveal that NQEs only have a marginal effect on the gNaO (r) radial distribution function. The first peak is only slightly smaller in the quantum simulations compared to the classical simulations (similar to the Li+ – O radial distribution functions reported by Wilkins et al. 14 for dilute solutions of LiF). The average number of molecules in the first solvation shell, calculated by integrating gNaO (r) up to the first minimum, is around 5.5 at the lowest concentration, in agreement with available experimental data. 73,74 At higher concentrations, this number increases slightly to become 5.8. 34 The degree of Na+ – OH – ion pairing is barely influenced by NQEs; only at the highest concentration (xNaOH = 0.333) is the ion pairing slightly decreased (Figure 3d–f). There is a bigger, although still small, effect for the amount of OH – in the second solvation shell at the lower concentrations xNaOH = 0.049 and 0.185, where there are more OH – ions in the second solvation shell for classical nuclei compared to quantum nuclei. In general, the quantum nuclei give a less “structured” gNaOOH− (r) than the classical nuclei in the first and second solvation shells. The


Page 10 of 34

Page 11 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


5

gNaO(r)

5 (a) xNaOH = 0.049

4

4

classical nuclei quantum nuclei (beads)

3

5

(b) xNaOH = 0.185

3

3

2

2

2

1

1

1

0

2

4

6 r (Å)

8

10

0

2

4

6 r (Å)

gNaO (d) xNaOH = 0.049

2

1

0

2

4

6 r (Å)

8

10

0

2

1

1

2

4

4

6 r (Å)

8

10

6 r (Å)

8

(f) xNaOH = 0.333

3

2

0

2

(r)

(e) xNaOH = 0.185

3


10

−

OH

3

8

(c) xNaOH = 0.333

4

10

0

2

4

6 r (Å)

8

2

(i) xNaOH = 0.333

10

gOO(r) 2

2

(g) xNaOH = 0.049

1.8


1.6

(h) xNaOH = 0.185

1.8

1.8

1.6

1.6

1.4

1.4

1.2

1.2

1.2

1

1

1

0.8

0.8

0.8

0.6

0.6

1.4

2

4

6 r (Å)

8

10

2

4

6 r (Å)

8

10

0.6

2

4

6 r (Å)

8

10

Figure 3: The radial distribution function gNaO (r), gNaOOH− (r), and gOO (r) calculated for different mole fractions of NaOH(aq) solutions, using classical and quantum nuclei (calculated from the positions of the individual beads). The radial distribution functions for all concentrations are given in the Supporting Information.

function gOO (r) (Figure 3g-i) is also only marginally affected by quantum nuclei; at low concentration the first minimum becomes somewhat deeper, but the second maximum larger, in agreement with the pure water MD simulations of Marsalek and Markland. 10 We also observe that the first peak becomes slightly smaller with quantum nuclei; this feature was not observed in the revPBE-D3 simulations of Marsalek and Markland, but instead in the recent force-field simulations of pure water by Wilkins et al. 11 The radial distribution functions gOH (r) (Figure 4a–c) and gHH (r) (Figure 4d–f) are more affected by the NQEs, because of the small mass of H. With quantum nuclei, the first minimum of gOH (r) is “lifted” (Figure 4a–c). The effect becomes more pronounced as the concentration of NaOH increases, because of the



2.5

gOH(r)

2.5 (a) xNaOH = 0.049

2 1.5

2.5

(b) xNaOH = 0.185

2


1.5

1

1

1

0.5

0.5

0.5

0

2

4 6 r (Å)

8

3

0

10

2


2

1

2

4

4 6 r (Å)

8

gHH(r)

3 (d) xNaOH = 0.049

0

0

0

10

2

1

1

4

r (Å)

2

4 6 r (Å)

8

10

(f) xNaOH = 0.333

2

2

0

3

(e) xNaOH = 0.185

0

(c) xNaOH = 0.333

2

1.5

0

Page 12 of 34

0

r (Å)

2

4 r (Å)

Figure 4: The radial distribution function gOH (r) and gHH (r) calculated for different mole fractions of NaOH(aq) solutions, using classical and quantum nuclei (calculated from the positions of the individual beads). The radial distribution functions for all concentrations are given in the Supporting Information. greater concentration of OH – in the solution. The same is true for gHH (r) (Figure 4d–f).

3.3

Hydrogen bonding

Figure 5 shows the average hydrogen bond angles

6

Oa Od Hd and noncovalent R(Hd · · · Oa ) distances, for

both OH – and OH2 as hydrogen bond acceptors. The averages were calculated over all hydrogen bonds, using the hydrogen bond criterion 6 Oa Od Hd ≤ 30◦ and R(Oa · · · Od ) ≤ 3.5 ˚ A. The hydrogen bonds accepted by OH – are on average stronger than those accepted by OH2 , as shown by the smaller average HB angles and intermolecular distances. NQEs give, for all concentrations, slightly stronger hydrogen bonds (smaller angles and intermolecular distances) as compared to the classical nuclei, for both OH – and OH2 as hydrogen bond acceptors. Figure 6 gives the number of accepted (A) and donated (D) hydrogen bonds for both H2 O and OH – as a function of mole fraction xNaOH for both classical (solid lines) and quantum (dashed lines) nuclei. The data for classical nuclei were originally reported in Ref. 33. Figure 6a shows that OH – ions are most likely to accept 4 hydrogen bonds for all concentrations except at the saturation limit (xNaOH = 0.333), where the populations of A(OH− ) = 4 and A(OH) = 3 are approximately equal. NQEs decrease the population


Page 13 of 34

(a)

(b) donor

2.1 acceptor class. nuclei H 2O 2 quant. nuclei avg. R(Hd...Oa) (Å)

class. nuclei

14 13

acceptor

acceptor donor

15 avg. hydrogen bond angle (deg.)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


quant. nuclei

12 11 10

class. nuclei

9 quant. nuclei

8 0

0.1 0.2 xNaOH

1.9 class. nuclei

1.8

1.7 acceptor OH–

0.3

quant. nuclei 0

0.1 0.2 xNaOH

0.3

Figure 5: Average (a) hydrogen bond angles, and (b) intermolecular hydrogen-bonding distances, calculated as a function of the mole fraction xNaOH , for both H2 O and OH – acceptors, for both classical and quantum nuclei. The schematic illustrations show two of the hydrogens as semi-transparent, indicating that both the hydrogen bond donors and acceptors can be either H2 O or OH – . of A(OH− ) = 5 and instead increase the populations of the lower-coordinated A(OH− ) = 3 and, at high concentration, A(OH− ) = 2. At low concentration, NQEs also decrease the population of A(OH− ) = 4. The effect of NQEs on A(H2 O) is generally smaller (Figure 6b). There is a slight increase in the fraction of A(H2 O) = 3, for which the proton transfer barrier is the smallest, 33 and a slight decrease in the fraction of A(H2 O) = 0, for which the proton transfer barrier is the largest. 33 NQEs increase the fraction of OH – donating a hydrogen bond (Figure 6c), such that with the present structural hydrogen bond definition, even at low concentration the population of D(OH− ) = 1 is slightly larger than the population of D(OH− ) = 0, which is not the case for the classical nuclei. However, for both classical and quantum nuclei, the fraction of OH – donating a hydrogen bond is around 50%; the exact number is likely to depend on how a hydrogen bond is defined. Figure 6d shows that NQEs barely influence the number of hydrogen bonds donated by H2 O molecules. It should be noted that – given the near-perfect balance between competing quantum effects that characterizes water – the overall impact of NQEs depends strongly on the choice of density functional. 75 Given that charged water species tend to form stronger hydrogen bonds, and that NQEs are known to accentuate the effects of strong bonds, 3 we expect that the changes due to the flavor of DFT would be quantitative but not qualitative in these highly-concentrated models. This issue is discussed further in section 4.2.



(a) A (OH−) 0.7

0.3 0.2

0.4 0.3

A=1

0.1

A=2 0

0.05

0.1

0.15 0.2 xNaOH

0.25

0.3

0

0.35

A=3 0

0.05

0.1

−

0.15 0.2 xNaOH

0.25

0.3

0.35

0.3

0.35

(d) D (H2O)

(c) D (OH ) 0.7

1 D=1

0.6

Fraction

0.4

D=0

0.3 0.2

0

0.05

0.1

0.15 0.2 xNaOH

0.25

0.6 0.4

1 2

D=1

0 1

0.1

D=2

0.8

0.5

0

A=0

0.2

0.1 0

A=2

0.5

A=5

A=3

0 1 2 3

0.6

Fraction

dashed: quantum

0.4

0.7

2 3 4 5

A=4

0.5 Fraction

(b) A (H2O)

solid: classical

0.6

Fraction

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 34

0.2

0.3

0.35

0

0

0.05

0.1

0.15 0.2 xNaOH

0.25

Figure 6: Distributions of the number of accepted (A) and donated (D) hydrogen bonds for (a,c) OH – and (b,d) H2 O, as a function of the mole fraction xNaOH , for both classical (solid dark lines) and quantum (dashed light lines) nuclei. For clarity, species that form rarely [with a fraction < 5% for every mole fraction, e.g., A(OH – ) = 1 and D(H2 O) = 0] have been omitted from the figure. The data for classical nuclei were originally reported in Ref. 33.


Page 15 of 34

3.4

Power spectra (a) power spectra (classical nuclei)

Intensity (arb. units)

10 xNaOH = 0.016 0.049 0.085 0.143 0.185 0.255 0.333

8 6 4 2 0

0

500 1000 1500 2000 2500 3000 3500 4000 ν (cm−1) (b) power spectra (quantum nuclei)

10 Intensity (arb. units)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


xNaOH = 0.016 0.049 0.085 0.143 0.185 0.255 0.333

8 6 4 2 0

0

500 1000 1500 2000 2500 3000 3500 4000 ν (cm−1)

Figure 7: Power spectra calculated from the Fourier transform of the velocity autocorrelation function for (a) classical nuclei, and (b) quantum nuclei, for different NaOH mole fractions xNaOH . Figures 7a–b show the calculated total power spectra, using both classical and quantum nuclei, for different mole fractions of NaOH. The power spectra contain contributions from both infrared and Raman active vibrations. In general, NQEs lower (red-shift) the vibrational frequencies by about 50 cm−1 for the HOH bend, and by about 100 cm−1 for the OH stretch, regardless of the NaOH concentration. This redshift is less pronounced than that observed by Marsalek and Markland 10 in their simulations for pure water, about 100 cm−1 for the HOH bend and 200-300 cm−1 for the OH stretch when using the revPBE-D3 GGA functional. Moreover, NQEs give a more pronounced tail towards low frequencies for the OH stretch around 2500–3300 cm−1 . This tail is associated with vibration along the proton transfer coordinate between H2 O and OH – , 28,54 and with increasing concentration, the tail towards low frequencies becomes more pronounced, in agreement with Raman experiments. 28,54 In infrared absorption experiments 55 of NaOH solutions, the HOH bending peak was found at ∼1640 cm−1 mostly independent of the concentration. In the present simulations, the peak position of the HOH bending vibration decreases as xNaOH increases, from ∼1660 cm−1 to ∼1610 cm−1 with classical nuclei, and from ∼1600 cm−1 to about ∼1550 cm−1 with quantum nuclei, as xNaOH increases from 0.016 to 0.333. Such a concentration dependence was not found in the IR experiments of Mandal et al. 55 This discrepancy 15 ACS Paragon Plus Environment


is possibly related to the fact that the intensity in the power spectra is not weighted by the change in dipole moment during the vibration. For the HOH bend, the frequency obtained using classical nuclei better matches the experimental data. In Raman experiments 28,54 on NaOH solutions, a sharp feature at ∼3620 cm−1 , assigned to the OH – stretch, was found to become more prominent with increasing concentration. Such a sharp feature is absent from the total power spectra in Figure 7, since the intensity in the power spectra is not the same as in Raman spectra. In order to evaluate the stretching frequencies for H2 O and OH – separately, and, additionally, evaluate the frequency dependence on the hydrogen-bonding environment, we calculated partial power spectra for H2 O and OH – molecules accepting different number of hydrogen bonds (an example is shown in Figure 8a for xNaOH = 0.085, extracted peak positions for different mole fractions are given in Figure 8b-c), as well as for hydrogens that either donate a hydrogen bond or that are left dangling (example for xNaOH = 0.085 in Figure 8d, extracted peak positions for different mole fractions in Figure 8e–f). Figure 8b shows that at low concentrations, a higher value of A(OH− ) gives a lower vibrational frequency, in agreement with the study of microsolvated OH – by Crespo and Hassanali. 44 However, at higher concentrations, the difference in stretching frequency for different values of A(OH− ) decreases; in fact, for the classical nuclei, OH – ions with A = 3 vibrate at somewhat smaller wavenumbers than OH – ions with A = 4 or A = 5. The reason for this is that the undercoordinated OH – ions are more likely to also donate hydrogen bonds, which causes a red-shift of the vibrational frequency (see below). Figure 8c reveals that, especially at low concentrations, the stretching frequencies for water molecules is more sensitive to the number of accepted hydrogen bonds than what was the case for hydroxide ions; at low concentration, the stretching frequency for A(H2 O) = 3 is about 200 cm−1 smaller than for A(H2 O) = 0. This can easily be rationalized, since H2 O molecules with high values of A are more likely to participate in proton transfer reactions, 33 which contribute to the long low-frequency tail of the stretching region around 2500-3500 cm−1 . 28 Figures 8e,f show that hydrogens donating a hydrogen bond vibrate at significantly lower wavenumbers (about 100 cm−1 for OH – and 200 cm−1 for H2 O) than those that do not donate a hydrogen bond.


Page 16 of 34

Page 17 of 34

(a) Power spectra at xNaOH = 0.085

classical nuclei

8


4

A=3

A=5 decr. − A(OH )

6

A=0

0 −2

OH− (dangl.) OH (don. HB) HOH (dangl.)

classical nuclei

8 6

decr. A(H2O)

2 A=3

(d) Power spectra at xNaOH = 0.085

10


10

−

HOH (don. HB)

4 2 0 −2

−4

−4 quantum nuclei

−6

3000

3200

3400

3600

3800

quantum nuclei

−6

4000

3000

3200

3400

−1

A=4

A=5

3700

A=5

3600

−1

3500

3400

3300

3300 total (class.) total (quant.) 0

0.1

total (class.) total (quant.)

0.2 xNaOH

3200 0.3 0

4000

dan

3800

do

n.

A=3

A=0

A=1 A=2 A=3 0.1

0.2

total (class.) total (quant.)

3800 HB

dangl. (class.)

3700 A=0 A=1 A=2

(f) H2O stretching

3900

gl.

3500

3400

3200

3900

3800

A=3 A=4

3700

(e) OH– stretching

(c) H2O stretching

0.3

ν (cm−1)

A=3

3800

3600

3900

ν (cm )

3900

3800

ν (cm )

da

3700 ng

l.

3600 do

n.

3500

ν (cm−1)

(b) OH– stretching

3600 −1

ν (cm )

ν (cm−1)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


HB

3500

3400

3400

3300

3300

3200

total (class.) total (quant.) 0

0.1

xNaOH

0.2 xNaOH

dangl. (quant.)

3600

don. HB (class.)

3200 0.3 0

don

. HB

0.1

( qu

ant. )

0.2 xNaOH

0.3

Figure 8: (a) Partial power spectra (VDOS) for OH – (brown/orange) and H2 O (purple/pink) accepting different numbers of hydrogen bonds at xNaOH = 0.085. The power spectra are plotted as negative values for the quantum nuclei. Each curve has been scaled slightly in order to enhance the clarity of the figure. (b-c) Extracted peak positions for the OH stretch for different numbers of accepted hydrogen bonds (A) as a function of xNaOH . For clarity, the points are slightly shifted horizontally. (d-f) The same as a, b, c, but for hydrogen atoms donating a hydrogen bond or not (“dangling”).



3.5

Proton transfer barriers and hydroxide ion lifetimes

Figure 9: Proton transfer free energy landscapes (arbitrarily shifted) calculated at xNaOH = 0.185 as a function of (a) the proton transfer coordinate δmin , and (b–c) both δmin and the intermolecular donoracceptor Od Oa distance. Figure 9a shows the calculated proton transfer free energy landscapes (PTFELs) using both classical and quantum nuclei at xNaOH = 0.185. Two effects are visible: the proton-transfer barrier decreases, and the free-energy minimum is shifted towards a smaller value of δmin . Overall, NQEs give a flatter PTFEL. Figure 9b–c show the same free-energy landscapes using the intermolecular donor-acceptor OO distance as an additional coordinate; a shortening of the OO distance gives a smaller PT barrier, as has also been observed for many other PT reactions. 17,46


Page 18 of 34

Page 19 of 34

2.98

2.85

3 classical nuclei

2.5

2.65 2.41

2.08 1.85

2 ΔF‡/kBT

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1.60

1.50

1.5 1.06

1 0.5 0

0.54 0

1.23

0.63

0.05

1.61

1.76

1.69

1.37

0.71

0.91

1.04

1.22

1.16

quantum nuclei (centroid)

0.1 0.15 0.2 0.25 mole fraction xNaOH

0.3

0.35

Figure 10: Calculated proton transfer free energy barriers in units of kB T for NaOH solutions of different mole fractions. The barriers are obtained from one-dimensional proton transfer free energy landscapes like the one in Figure 9a. The black arrows and black text indicate the difference in proton transfer barriers calculated using classical and quantum nuclei.

The aforementioned trends are valid for all concentrations of NaOH. Figure 10 gives the PT barriers (calculated from the one-dimensional PTFELs, as in Figure 9a) for all of the studied concentrations of NaOH. At low concentration (xNaOH = 0.016), quantum nuclei give a PT barrier ∆F ‡ that is 1.06 kB T smaller than what is obtained in the classical simulation. However, this difference increases with concentration to 1.76 kB T in the saturated solution (xNaOH = 0.333). Thus, NQEs play a bigger role for the proton transfer barriers at high concentrations. Figure 11a shows two example time correlation functions C(∆t) for a proton transfer reaction, calculated for classical and quantum nuclei at xNaOH = 0.185. The correlation function gives the probability that a “stable” OH – at time t = t1 has not turned into a “stable” H2 O molecule within the interval t1 ≤ t ≤ t1 +∆t (see also Methods). As can be seen, the function decays faster for quantum nuclei, in accordance with the smaller PT barrier calculated in the quantum case (Figure 10). The fitted double exponential functions that are used to extract the characteristic OH – lifetimes, τPT , are shown as lines in the example in Figure 11a. The values of τPT for all mole fractions xNaOH are given in Figure 11b, as calculated using both classical and quantum nuclei. Indeed, when NQEs are taken into account the lifetime of OH – is calculated to be about half of that obtained with classical nuclei, throughout the NaOH room-temperature solubility range. This can also be stated in terms of the proton transfer rate, that is roughly doubled when NQEs are taken into account, for all concentrations of NaOH.



(a) OH− → H2O

1

xNaOH = 0.185

0.9 classical nuclei τPT = 5.50 ps

0.8 0.7 C(Δt)

0.6 0.5 0.4

quantum nuclei (centroid) τPT = 2.57 ps

0.3 0.2 0.1 0

0

2

4

6

(b)

8

10 12 Δt (ps)

14

16

18

20

OH− → H2O

16

13.77

14 12 classical nuclei

10 τPT (ps)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 34

8.92

8

4 2

7.29

5.60

6 3.79 1.33 2.33 1.74 0.65 0.84 1.00

0 0

0.05

4.36

2.57 1.76 quantum nuclei (centroid)


0.3

0.35

Figure 11: (a) Calculated time correlation functions C(∆t) at xNaOH = 0.185 for both classical and quantum nuclei, that give the probability that a stable OH – ion at time t = t1 has not turned into a stable H2 O molecule within the interval t1 < t ≤ t1 + ∆t. The lines indicate fits to double exponential functions. (b) Calculated average hydroxide ion lifetimes τPT . The error bars for the classical nuclei simulations are ±2.09×(standard error of the mean), calculated over 20 independent NVE trajectories.


Page 21 of 34

3.6

Mean residence times and diffusion coefficients Na−[OH−/OH2] → Na ⋅ ⋅ ⋅ [OH−/OH2] 51.3

50

48.3 τexchange (ps)

40 30

19.4

20

0

29.2

classical nuclei

14.9

9.8 10.7 12.0

10

11.4 9.4 10.7 0

0.05

13.3

24.2 17.1

quantum nuclei (centroid)


0.3

0.35

Figure 12: Calculated mean residence times τexchange for molecules in the first solvation shell around Na+ . The error bars for the classical nuclei simulations are ±2.09×(standard error of the mean), calculated over 20 independent NVE trajectories. Figure 12 shows the calculated mean residence times (τexchange ) for a ligand around Na+ to stay in the first solvation shell before it escapes to the second solvation shell. At low NaOH concentration (xNaOH < 0.1), NQEs have only marginal impact on the calculated exchange rates. However, at higher concentrations, NQEs increase the rate of exchange (decrease the mean residence time of molecules in the first solvation shell). The mean residence time around Na+ is quite short and therefore difficult to determine experimentally. 76 Our results agree well with previous molecular dynamics simulations, which have given estimates between 10 and 100 ps for dilute aqueous Na+ solutions. 77–82 (a) Na

(b) O

(c) H

(d) Hydroxide 14

class. quant.

9

1

1

2

D × 109 m2 s−1

2

D × 109 m2 s−1

2 −1

2

12 D × 10 m s

D × 109 m2 s−1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

10 8 6 4 2

0

0

0.1 0.2 0.3 xNaOH

0

0

0.1 0.2 0.3 xNaOH

0

0

0.1 0.2 0.3 xNaOH

0

0

0.1 0.2 0.3 xNaOH

Figure 13: Calculated diffusion coefficients for (a) Na, (b) O, (c) H, and (d) OH – as a function of mole fraction xNaOH for both classical and quantum nuclei. The arrows indicate the experimental room temperature ionic diffusion coefficients at infinite dilution (from Ref. 83) and the water self-diffusion coefficient (from Ref. 84). Figure 13 gives the obtained diffusion coefficients for Na, O, H, and OH – , calculated from the slope



Page 22 of 34

of least-squares straight-line fit to the MSD in the interval 4–20 ps, and corrected for finite-size effects according to equation 6, from our simulations using both classical and quantum nuclei. In general, NQEs result in marginally slower diffusion of Na, O, and H at low concentration, but faster at high concentration. The calculated value of D(OH− ) takes both mass transport and proton transfer events into account. In −

OH accordance with the smaller PT barriers when NQEs are considered, Dquant. is considerably larger, by about

a factor of 1.8 for all concentrations, than the corresponding classical value. +

In the dilute limit, the experimental room-temperature diffusion coefficients of Na+ and OH – are DNa = −

1.33 × 10−9 m2 ·s−1 and DOH = 5.27 × 10−9 m2 ·s−1 , 83 and the water self-diffusion coefficient is 2.299 × 10−9 m2 ·s−1 . 84 Thus, the present neural network potential slightly overestimates the diffusion coefficient of Na+ and self-diffusion coefficient of water, but significantly overestimates the diffusion coefficient of OH – . In Ref. 40, it was shown that ab initio simulations employing the dispersion-corrected density functional PBE-TS also overestimated the diffusion coefficient of OH – . When NQEs are taken into account, the problem is exacerbated. Nevertheless, the present simulations illustrate that NQEs increase the diffusion coefficient of OH – by about 80% for all concentrations.

4 4.1

Discussion Structure and dynamics of NaOH solutions

The solvation structure around OH – has been debated for a long time, with studies proposing that the OH – ion either accepts three hydrogen bonds, 39 or that it is “hypercoordinated” and accepts four hydrogen bonds, 16,26,40,53 with proton transfer primarily occurring after a hydrogen bond fluctuation in which one of these hydrogen bonds is lost. Our present simulations reaffirm the notion that OH – ions at low concentrations are primarily “hypercoordinated”, i.e., that they accept four hydrogen bonds. This is the case for both the classical and quantum simulations (Figure 6). However, NQEs increase the fraction of threecoordinated OH – , which is consistent with the lower proton transfer barriers in the quantum simulations (Figure 10). Whether the NQEs primarily cause a change in the hydrogen-bonding environment, thus indirectly decreasing the proton transfer barriers, or whether the NQEs primarily decrease the proton transfer barriers, thus indirectly changing the hydrogen-bonding environment, is still an open question. Most likely, it is a combination of the two effects. In general, we find that NQEs in aqueous NaOH solutions have only a small impact on the calculated radial distribution functions involving the heavier elements Na and O (Figure 3). For the lighter element H, a greater degree of delocalization is clearly observed (Figure 4), which cause the HOH bending and OH


Page 23 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


stretching frequencies to become smaller (Figure 7). At every concentration of NaOH, NQEs yield hydrogen bonds that, from a geometric viewpoint, are slightly more “ideal” (shorter intermolecular distances and smaller Oa Od Hd angles) and consequently stronger, regardless of whether the hydrogen bond acceptor is a hydroxide ion or a water molecule (Figure 5). These stronger hydrogen bonds slow down the diffusion of Na, O, and H at low NaOH concentration somewhat (Figure 13). However, at high concentrations (xNaOH ≥ 0.143), NQEs increase the diffusion coefficients of Na, O, and H (Figure 13), and also shorten the mean residence time of molecules in the first solvation shell around Na+ (Figure 12). Because the NQE-influenced changes to the average hydrogen-bonding structure are roughly independent of the concentration (Figure 5), we attribute the increase in D to the significant increase in proton transfer rates that NQEs give rise to. After successful proton transfer events, the local environment around the newly-formed H2 O and OH – relax, resulting in the breaking and forming of new hydrogen bonds. 40 Indeed, in a recent work, we also showed that proton transfer events in the first hydration shell around Na+ were coupled to an increased rate of exchange of molecules between the first and second hydration shells, 35 which explains the decreased mean residence time around Na+ when NQEs are taken into account at the higher concentrations (Figure 12). Figure 11b shows that NQEs approximately double the proton transfer rate (halve the average OH – lifetime) for all concentration of NaOH. This is because the delocalization of the H atoms result in a significantly decreased proton transfer barrier ∆F ‡ (Figure 10). Because OH – (aq) ions mainly diffuse via a Grotthuss-like mechanism, in which proton transfer events cause the diffusion of the ions, the simulations including NQEs give an OH – diffusion coefficient that is approximately 80% greater than those obtained in classical simulations. Figure 10 shows that the PT barriers are more affected by NQEs at high concentrations than at low concentrations. Despite this, the hydroxide ion lifetimes and diffusion coefficients are about equally affected by NQEs regardless of the concentration. We attribute this to the fact that a greater decrease of the PT barrier at high concentration, is partially counteracted by the resulting increased “flatness” of the free-energy potential energy surface for proton transfer (as illustrated in Figure 9a at xNaOH = 0.185). A flatter potential energy surface means that the protons will have an increased propensity for “rattling” multiple times between the same pair of donor/acceptor O atoms before settling on one of them. Such rattling events clearly do not contribute to the long-scale diffusion of OH – , and the lifetimes in Figure 11b were calculated using the stable states formalism in order to explicitly exclude them.



4.2

Limitations of the DFT reference

This work explored the role of nuclear quantum effects on aqueous NaOH solutions, by means of molecular dynamics simulations employing a high-dimensional neural network potential. The potential had been trained to reproduce a potential energy surface calculated at the RPBE-D3 level of theory. RPBE is a GGA-type functional, and the results presented in this work therefore cannot be expected to be more accurate than what can be obtained at the GGA level of theory. Many authors have concluded that nuclear quantum effects can have very different impact in aqueous systems depending on the method used to calculate the potential energy surface (e.g. force-fields, GGA DFT, hybrid DFT, etc.) 1,2,7,10 Marsalek and Markland 10 recently calculated a number of properties of liquid water using both classical and quantum nuclei, for both dispersion-corrected GGA (revPBE-D3) and dispersion-corrected hybrid (revPBE0-D3) functionals. Using a GGA DFT led to a softening of the vibrational modes, and so their classical spectra were in better agreement with experiments than those with quantum nuclei, due to error cancellation. Using a hybrid functional with exact exchange corrected the softening, and led to better agreement between quantum spectra and experimental measurements. Similarly, they found that the diffusion coefficients were better reproduced using classical nuclei at the GGA level, but using quantum nuclei at the hybrid DFT level. Our simulation results in this work, based on the dispersion-corrected RPBE-D3 density functional, give mostly better agreement to experiment when using classical nuclei. This is particularly true for the frequency of the HOH bending vibration (Figure 7), and the diffusion coefficient for OH – (Figure 13). Because the OH – diffusion coefficient is better reproduced (although overestimated) using classical nuclei, it is reasonable to believe that also the proton transfer barriers (Figure 10) are more accurate from the classical simulations as opposed to the quantum simulations. However, in some cases, we find that the difference between classical and quantum simulations is quite small; for example, in the Na-O radial distribution functions (Figure 3), the average hydrogen-bonding angles and intermolecular distances (Figure 5), and the diffusion coefficients for Na, O, and H (Figure 13). Finally, in the case of the OH stretching vibrations, we find that for the lowconcentration limit (approaching pure water, which in experiments has a peak around 3404 cm−1 ), the peak in the H2 O VDOS is ∼3500 cm−1 for classical nuclei, but ∼3400 for quantum nuclei (lines in Figure 6f). Similarly, at low concentration, the OH – stretching peak is ∼3850 cm−1 for classical nuclei, but ∼3680 cm−1 for quantum nuclei (which can be compared with ∼3610 cm−1 in Raman experiments 28,54 ). Thus, using classical nuclei, the OH stretching frequency, for both H2 O and OH – , is somewhat blue-shifted with respect to the experimental reference. We obtain better agreement in the quantum simulations. Marsalek and Markland 10 found near-perfect agreement for pure water between experiment and the GGA revPBE-D3 functional using classical, not quantum, nuclei. There are two possible reasons for the discrepancy between


Page 24 of 34

Page 25 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


our results and theirs: firstly, the neural network potential is fitted to reproduce a slightly different density functional, namely RPBE-D3; secondly, our simulations were performed at the neural-network-equilibrated densities (at T = 300 K, p = 1 bar), which are about 10% smaller than the experimentally observed densities. Consequently, the OH stretches are less “hindered” by neighboring hydrogen bond acceptors, potentially resulting in a blue-shift of the vibrational frequency. In summary, our work mostly supports the notion that fairly accurate simulations can be performed using RPBE-D3 with classical nuclei (with biggest exception being the overestimation of the OH – diffusion coefficient; as shown in Ref. 40, a better value for this quantity can be obtained using a hybrid density functional). This is useful, since the underlying electronic structure theory is relatively computationally inexpensive, and using classical nuclei is still much less computationally demanding than modeling quantum nuclei. However, the real-life nuclei are quantum particles, and the employed RPBE-D3 functional therefore gives (approximately) the right answer for the wrong reason in the classical case. A more accurate description of the physics of both the electrons and light nuclei necessitates a more advanced description of the electronic structure theory, coupled with an explicit treatment of nuclear quantum effects.

5

Conclusions

This work investigated the role of nuclear quantum effects (NQEs) in NaOH solutions at T = 300 K by comparing thermostatted ring polymer and classical molecular dynamics simulations employing a highdimensional neural network potential parameterized to reproduce a dispersion-corrected density-functionaltheory-calculated potential energy surface. 33 TRPMD simulations are computationally very demanding, which typically limits their applicability to small systems and short time scales. Here, by using a reactive neural network potential, that is computationally inexpensive but nevertheless faithfully reproduces the underlying reference calculations, the influence of NQEs on NaOH solutions could be studied using long TRPMD simulations for large systems (∼1000-1500 atoms, 7 different concentrations, ring polymers using 16 beads with each simulation lasting several hundred picoseconds). The simulations show that NQEs, which are most prominently manifested as a delocalization of the H atoms, have a very small effect on the solvation shell structure around Na+ , slightly strengthen the hydrogen bonds accepted by both OH – and OH2 , lower the vibrational frequencies by about 50 cm−1 for the HOH bend, and by about 100 cm−1 for the OH stretch, decrease the proton transfer barriers, and increase the proton transfer rates. The greater proton transfer rates in turn lead to an increased amount of protontransfer-driven water exchange events in the first hydration shell around Na+ ions, thus decreasing the mean residence time of molecules around Na+ . The impact of NQEs on the PT barrier increases with increasing 25 ACS Paragon Plus Environment


concentration of NaOH, but the smaller PT barriers also give a flatter free-energy landscape, thus increasing the amount of proton “rattling” that occurs between two O atoms before the proton settles on one of them. Moreover, the increased PT rate from NQEs leads to a roughly 80 % increase of the diffusion coefficient of OH – (aq) for all NaOH concentrations, exacerbating the overestimation of DOH− that is due to the limitations of the underlying GGA DFT calculations. This study provides a qualitative assessment of the importance of NQEs in concentrated basic solutions. It appears that, as observed in neat water, gradient-corrected density functional theory leads to an excessive softening of the covalent OH bond, which leads to a better agreement with experiments in classical simulations than when fully accounting for nuclear quantum effects. This underscores the importance of performing simulations that fully account for all physical effects, as otherwise error cancellation can lead to deceptively accurate results that do not fully capture the nature of the phenomenon. Refining the reference data set using a more accurate electronic structure method, such as hybrid DFT, shall make it possible to achieve closer agreement with experimental observables, and to reach the goal of a predictive model that treats at the highest level of theory the quantum nature of both electrons and light nuclei. Finally, we comment that NaOH solutions contain water molecules and hydroxide ions in many different local solvation environments, with respect to hydrogen-bonding and cation coordination. This makes it an excellent system for evaluating how, for example, vibrational frequencies and proton transfer barriers depend on the solvation environment. As an example, in this work we explored how the vibrational frequencies depend on the local hydrogen-bonding environment around H2 O and OH – . However, there are still many possibilities to explore in this regard, for example, evaluating the vibrational fingerprint for water molecules as a function of the “tetrahedrality” of the solvation environment, distinguishing between hydrogen bonds donated to H2 O molecules and those donated to OH – ions, and evaluating the IR and Raman spectra including the correct intensities. Such analyses would help to provide a complete picture of the structure and dynamics of NaOH solutions, and would serve as excellent avenues to explore in future work.

Acknowledgement This work was supported by the Cluster of Excellence RESOLV (EXC 1069) funded by the Deutsche Forschungsgemeinschaft, and the DFG Heisenberg professorship Be3264/11-2. M.C. acknowledges financial support by the Swiss National Science Foundation (project ID 200021-159896).


Page 26 of 34

Page 27 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Supporting Information Available Calculated radial distribution functions for all concentrations of NaOH(aq) in this work, using both quantum and classical nuclei.

References (1) Ceriotti, M.; Fang, W.; Kusalik, P. G.; McKenzie, R. H.; Michaelides, A.; Morales, M. A.; Markland, T. E. Nuclear quantum effects in water and aqueous systems: experiment, theory, and current challenges. Chem. Rev. 2016, 116, 7529–7550. (2) Habershon, S.; Markland, T. E.; Manolopoulos, D. E. Competing quantum effects in the dynamics of a flexible water model. J. Chem. Phys. 2009, 131, 024501. (3) Li, X.-Z.; Walker, B.; Michaelides, A. Quantum nature of the hydrogen bond. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 6369–6373. (4) Ceriotti, M.; Cuny, J.; Parrinello, M.; Manolopoulos, D. E. Nuclear quantum effects and hydrogen bond fluctuations in water. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 15591–15596. (5) McKenzie, R. H. A diabatic state model for donor-hydrogen vibrational frequency shifts in hydrogen bonded complexes. Chem. Phys. Lett. 2012, 535, 196–200. (6) Markland, T. E.; Ceriotti, M. Nuclear quantum effects enter the mainstream. Nat. Rev. Chem. 2018, 2, 0109. (7) Miller III, T. F.; Manolopoulos, D. E. Quantum diffusion in liquid water from ring polymer molecular dynamics. J. Chem. Phys. 2005, 123, 154504. (8) Cheng, B.; Behler, J.; Ceriotti, M. Nuclear quantum effects in water at the triple point: using theory as a link between experiments. J. Phys. Chem. Lett. 2016, 7, 2210–2215. (9) Kapil, V.; Behler, J.; Ceriotti, M. High order path integrals made easy. J. Chem. Phys. 2016, 145, 234103. (10) Marsalek, O.; Markland, T. E. Quantum dynamics and spectroscopy of ab initio liquid water: the interplay of nuclear and electronic quantum effects. J. Phys. Chem. Lett. 2017, 8, 1545–1551.



Page 28 of 34

(11) Wilkins, D. M.; Manolopoulos, D. E.; Pipolo, S.; Laage, D.; Hynes, J. T. Nuclear quantum effects in water reorientation and hydrogen-bond dynamics. J. Phys. Chem. Lett. 2017, 8, 2602–2607. (12) Machida, M.; Kato, K.; Shiga, M. Nuclear quantum effects of light and heavy water studied by allelectron first principles path integral simulations. J. Chem. Phys. 2018, 148, 102324. (13) Ojha, D.; Henao, A.; K¨ uhne, T. D. Nuclear quantum effects on the vibrational dynamics of liquid water. J. Chem. Phys. 2018, 148, 102328. (14) Wilkins, D. M.; Manolopoulos, D. E.; Dang, L. X. Nuclear quantum effects in water exchange around lithium and fluoride ions. J. Chem. Phys. 2015, 142, 064509. (15) Videla, P. E.; Rossky, P. J.; Laria, D. Isotope effects in aqueous solvation of simple halides. J. Chem. Phys. 2018, 148, 102306. (16) Tuckerman, M. E.; Marx, D.; Parrinello, M. The nature and transport mechanism of hydrated hydroxide ions in aqueous solution. Nature 2002, 417, 925–929. (17) Marx, D. Proton transfer 200 years after von Grotthuss:

insights from ab initio simulations.

ChemPhysChem 2006, 7, 1848–1870. (18) Agmon, N.; Bakker, H. J.; Campen, R. K.; Henchman, R. H.; Pohl, P.; Roke, S.; Thämer, M.; Hassanali, A. Protons and hydroxide ions in aqueous systems. Chem. Rev. 2016, 116, 7642–7672. (19) Hassanali, A.; Giberti, F.; Cuny, J.; K¨ uhne, T. D.; Parrinello, M. Proton transfer through the water gossamer. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 13723–13728. (20) Rossi, M.; Ceriotti, M.; Manolopoulos, D. E. Nuclear quantum effects in H+ and OH− diffusion along confined water wires. J. Phys. Chem. Lett. 2016, 7, 3001–3007. (21) Litman, Y.; Donadio, D.; Ceriotti, M.; Rossi, M. Decisive role of nuclear quantum effects on surface mediated water dissociation at finite temperature. J. Chem. Phys. 2018, 148, 102320. (22) Chen, B.; Park, J. M.; Ivanov, I.; Tabacchi, G.; Klein, M. L.; Parrinello, M. First-principles study of aqueous hydroxide solutions. J. Am. Chem. Soc. 2002, 124, 8534–8535. (23) Chen, B.; Ivanov, I.; Park, J. M.; Parrinello, M.; Klein, M. L. Solvation structure and mobility mechanism of OH− : a Car-Parrinello molecular dynamics investigation of alkaline solutions. J. Phys. Chem. B 2002, 106, 12006–12016. 28 ACS Paragon Plus Environment

Page 29 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(24) Tuckerman, M. E.; Chandra, A.; Marx, D. Structure and dynamics of OH− (aq). Acc. Chem. Res. 2006, 39, 151–158. (25) Megyes, T.; Bálint, S.; Grósz, T.; Radnai, T.; Bakó, I.; Sipos, P. The structure of aqueous sodium hydroxide solutions: a combined solution x-ray diffraction and simulation study. J. Chem. Phys. 2008, 128, 044501. (26) Marx, D.; Chandra, A.; Tuckerman, M. E. Aqueous basic solutions: hydroxide solvation, structural diffusion, and comparison to the hydrated proton. Chem. Rev. 2010, 110, 2174–2216. (27) Hermansson, K.; Bopp, P. A.; Sp˚ angberg, D.; Pejov, L.; Bak, I.; Mitev, P. D. The vibrating hydroxide ion in water. Chem. Phys. Lett. 2011, 514, 1 – 15. (28) Stefanski, J.; Schmidt, C.; Jahn, S. Aqueous sodium hydroxide (NaOH) solutions at high pressure and temperature: insights from in situ Raman spectroscopy and ab initio molecular dynamics simulations. Phys. Chem. Chem. Phys. 2018, 20, 21629–21639. (29) Bruni, F.; Ricci, M. A.; Soper, A. K. Structural characterization of NaOH aqueous solution in the glass and liquid states. J. Chem. Phys. 2001, 114, 8056–8063. (30) Bonthuis, D. J.; Mamatkulov, S. I.; Netz, R. R. Optimization of classical nonpolarizable force fields for OH− and H3 O+ . J. Chem. Phys. 2016, 144, 104503. (31) Zhang, W.; van Duin, A. C. T. Second-generation ReaxFF water force field: improvements in the description of water density and OH− anion diffusion. J. Phys. Chem. B 2017, 121, 6021–6032. (32) Rimsza, J.; Jones, R.; Criscenti, L. Interaction of NaOH solutions with silica surfaces. J. Coll. Interf. Sci. 2018, 516, 128 – 137. (33) Hellström, M.; Behler, J. Concentration-dependent proton transfer mechanisms in aqueous NaOH solutions: from acceptor-driven to donor-driven and back. J. Phys. Chem. Lett. 2016, 7, 3302–3306. (34) Hellström, M.; Behler, J. Structure of aqueous NaOH solutions: insights from neural-network-based molecular dynamics simulations. Phys. Chem. Chem. Phys. 2017, 19, 82–96. (35) Hellström, M.; Behler, J. Proton-transfer-driven water exchange mechanism in the Na+ solvation shell. J. Phys. Chem. B 2017, 121, 4184–4190.



(36) Behler, J.; Parrinello, M. Generalized neural-network representation of high-dimensional potentialenergy surfaces. Phys. Rev. Lett. 2007, 98, 146401. (37) Hammer, B.; Hansen, L. B.; Nørskov, J. K. Improved adsorption energetics within density-functional theory using revised Perdew-Burke-Ernzerhof functionals. Phys. Rev. B 1999, 59, 7413–7421. (38) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (39) Asthagiri, D.; Pratt, L. R.; Kress, J. D.; Gomez, M. A. Hydration and mobility of HO− (aq). Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 7229–7233. (40) Chen, M.; Zheng, L.; Santra, B.; Ko, H.-Y.; Jr, R. A. D.; Klein, M. L.; Car, R.; Wu, X. Hydroxide diffuses slower than hydronium in water because its solvated structure inhibits correlated proton transfer. Nat. Chem. 2018, 10, 413–419. (41) Li, X.; Teige, V. E.; Iyengar, S. S. Can the four-coordinated, penta-valent oxygen in hydroxide water clusters be detected through experimental vibrational spectroscopy? J. Phys. Chem. A 2007, 111, 4815–4820. (42) Bankura, A.; Chandra, A. Hydration structure and dynamics of a hydroxide ion in water clusters of varying size and temperature: quantum chemical and ab initio molecular dynamics studies. Chem. Phys. 2012, 400, 154 – 164. (43) Crespo, Y.; Hassanali, A. Unveiling the Janus-like properties of OH− . J. Phys. Chem. Lett. 2015, 6, 272–278. (44) Crespo, Y.; Hassanali, A. Characterizing the local solvation environment of OH− in water clusters with AIMD. J. Chem. Phys. 2016, 144, 074304. (45) Bankura, A.; Chandra, A. Hydroxide ion can move faster than an excess proton through onedimensional water chains in hydrophobic narrow pores. J. Phys. Chem. B 2012, 116, 9744–9757. (46) Tocci, G.; Michaelides, A. Solvent-induced proton hopping at a water-oxide interface. J. Phys. Chem. Lett. 2014, 5, 474–480. (47) Chen, C.; Tse, Y.-L. S.; Lindberg, G. E.; Knight, C.; Voth, G. A. Hydroxide solvation and transport in anion exchange membranes. J. Am. Chem. Soc. 2016, 138, 991–1000. 30 ACS Paragon Plus Environment

Page 30 of 34

Page 31 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(48) Dong, D.; Zhang, W.; van Duin, A. C. T.; Bedrov, D. Grotthuss versus vehicular transport of hydroxide in anion-exchange membranes: insight from combined reactive and nonreactive molecular simulations. J. Phys. Chem. Lett. 2018, 9, 825–829. (49) Quaranta, V.; Hellström, M.; Behler, J. Proton-transfer mechanisms at the water-ZnO interface: the role of presolvation. J. Phys. Chem. Lett. 2017, 8, 1476–1483. (50) Quaranta, V.; Hellström, M.; Behler, J.; Kullgren, J.; Mitev, P. D.; Hermansson, K. Maximally resolved anharmonic OH vibrational spectrum of the water/ZnO(10¯10) interface from a high-dimensional neural network potential. J. Chem. Phys. 2018, 148, 241720. (51) Kebede, G. G.; Mitev, P. D.; Broqvist, P.; Kullgren, J.; Hermansson, K. Hydrogen-bond relations for surface OH species. J. Phys. Chem. C 2018, 122, 4849–4858. (52) Mu˜ noz Santiburcio, D.; Marx, D. On the complex structural diffusion of proton holes in nanoconfined alkaline solutions within slit pores. Nat. Commun. 2016, 7, 12625. (53) Botti, A.; Bruni, F.; Imberti, S.; Ricci, M. A.; Soper, A. K. Solvation of hydroxyl ions in water. J. Chem. Phys. 2003, 119, 5001–5004. (54) Corridoni, T.; Sodo, A.; Bruni, F.; Ricci, M.; Nardone, M. Probing water dynamics with OH− . Chem. Phys. 2007, 336, 183 – 187. (55) Mandal, A.; Tokmakoff, A. Vibrational dynamics of aqueous hydroxide solutions probed using broadband 2DIR spectroscopy. J. Chem. Phys. 2015, 143, 194501. (56) Behler, J. First principles neural network potentials for reactive simulations of large molecular and condensed systems. Angew. Chem. Int. Edit. 2017, 56, 12828–12840. (57) Forster-Tonigold, K.; Groß, A. Dispersion corrected RPBE studies of liquid water. J. Chem. Phys. 2014, 141, 064501. (58) Morawietz, T.; Singraber, A.; Dellago, C.; Behler, J. How van der Waals interactions determine the unique properties of water. Proc. Natl. Acad. Sci. U. S. A. 2016, 113, 8368–8373. (59) Morawietz, T.; Marsalek, O.; Pattenaude, S. R.; Streacker, L. M.; Ben-Amotz, D.; Markland, T. E. The interplay of structure and dynamics in the Raman spectrum of liquid water over the full frequency and temperature range. J. Phys. Chem. Lett. 2018, 9, 851–857.



(60) Sipos, P. M.; Hefter, G.; May, P. M. Viscosities and densities of highly concentrated aqueous MOH Solutions (M+ = Na+ , K+ , Li+ , Cs+ , (CH3 )4 N+ ) at 25.0 ◦ C. J. Chem. Eng. Data 2000, 45, 613–617. (61) Rossi, M.; Ceriotti, M.; Manolopoulos, D. E. How to remove the spurious resonances from ring polymer molecular dynamics. J. Chem. Phys. 2014, 140, 234116. (62) Ceriotti, M.; Parrinello, M.; Markland, T. E.; Manolopoulos, D. E. Efficient stochastic thermostatting of path integral molecular dynamics. J. Chem. Phys. 2010, 133, 124104. (63) Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. J. Comp. Phys. 1995, 117, 1–19. (64) Ceriotti, M.; More, J.; Manolopoulos, D. E. i-PI: a Python interface for ab initio path integral molecular dynamics simulations. Comp. Phys. Comm. 2014, 185, 1019–1026. (65) Kapil, V.; Rossi, M.; Marsalek, O.; Petraglia, R.; Litman, Y.; Spura, T.; Cheng, B.; Cuzzocrea, A.; Meiner, R. H.; Wilkins, D. M. et al. i-PI 2.0: a universal force engine for advanced molecular simulations. arXiv:1808.03824. (66) Humphrey, W.; Dalke, A.; Schulten, K. VMD – Visual Molecular Dynamics. J. Mol. Graph. 1996, 14, 33–38. (67) Stone, J. An efficient library for parallel ray tracing and animation. M.Sc. thesis, Computer Science Department, University of Missouri-Rolla, 1998. (68) Luzar, A.; Chandler, D. Effect of environment on hydrogen bond dynamics in liquid water. Phys. Rev. Lett. 1996, 76, 928–931. (69) Laage, D.; Hynes, J. T. On the residence time for water in a solute hydration shell: application to aqueous halide solutions. J. Phys. Chem. B 2008, 112, 7697–7701. (70) Kuhn, H. W. The Hungarian method for the assignment problem. Nav. Res. Logist. Q. 1955, 2, 83–97. (71) D¨ unweg, B.; Kremer, K. Molecular dynamics simulation of a polymer chain in solution. J. Chem. Phys. 1993, 99, 6983–6997. (72) Yeh, I.-C.; Hummer, G. System-size dependence of diffusion coefficients and viscosities from molecular dynamics simulations with periodic boundary conditions. J. Phys. Chem. B 2004, 108, 15873–15879.


Page 32 of 34

Page 33 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(73) Mähler, J.; Persson, I. A study of the hydration of the alkali metal ions in aqueous solution. Inorg. Chem. 2012, 51, 425–438. (74) Galib, M.; Baer, M. D.; Skinner, L. B.; Mundy, C. J.; Huthwelker, T.; Schenter, G. K.; Benmore, C. J.; Govind, N.; Fulton, J. L. Revisiting the hydration structure of aqueous Na+ . J. Chem. Phys. 2017, 146, 084504. (75) Wang, L. L.; Ceriotti, M.; Markland, T. E. T. Quantum fluctuations and isotope effects in ab initio descriptions of water. J. Chem. Phys. 2014, 141, 104502. (76) Helm, L.; Merbach, A. E. Inorganic and bioinorganic solvent exchange mechanisms. Chem. Rev. 2005, 105, 1923–1960. (77) Rey, R.; Hynes, J. T. Hydration shell exchange kinetics: an MD study for Na+ (aq). J. Phys. Chem. 1996, 100, 5611–5615. (78) Koneshan, S.; Rasaiah, J. C.; Lynden-Bell, R. M.; Lee, S. H. Solvent structure, dynamics, and ion mobility in aqueous solutions at 25 ◦ C. J. Phys. Chem. B 1998, 102, 4193–4204. (79) Kerisit, S.; Rosso, K. M. Transition path sampling of water exchange rates and mechanisms around aqueous ions. J. Chem. Phys. 2009, 131, 114512. (80) Bankura, A.; Carnevale, V.; Klein, M. L. Hydration structure of salt solutions from ab initio molecular dynamics. J. Chem. Phys. 2013, 138, 014501. (81) Kelley, M.; Donley, A.; Clark, S.; Clark, A. Structure and dynamics of NaCl ion pairing in solutions of water and methanol. J. Phys. Chem. B 2015, 119, 15652–15661. (82) Lee, Y.; Thirumalai, D.; Hyeon, C. Ultrasensitivity of water exchange kinetics to the size of metal ion. J. Am. Chem. Soc. 2017, 139, 12334–12337. (83) Haynes, W. M., Ed. CRC Handbook of chemistry and physics, 97th ed.; CRC Press: Boca Raton, 2017. (84) Holz, M.; Heil, S. R.; Sacco, A. Temperature-dependent self-diffusion coefficients of water and six selected molecular liquids for calibration in accurate 1 H NMR PFG measurements. Phys. Chem. Chem. Phys. 2000, 2, 4740–4742.



TOC graphic 10 mol/L NaOH(aq)

OH vibrations OH–

Na+ H 2O

proton transfer

hydrogen bonding


Page 34 of 34

Nuclear Quantum Effects in Sodium Hydroxide Solutions from Neural

Recommend Documents