Divide-and-Conquer-Type Density-Functional ... - ACS Publications

Article pubs.acs.org/JPCB

Divide-and-Conquer-Type Density-Functional Tight-Binding Simulations of Hydroxide Ion Diffusion in Bulk Water Aditya Wibawa Sakti,† Yoshifumi Nishimura,‡ and Hiromi Nakai*,†,‡,§,∥ †

Department of Chemistry and Biochemistry, School of Advanced Science and Engineering, ‡Research Institute for Science and Engineering, Waseda University, Shinjuku, Tokyo 169-8555, Japan § Core Research for Evolutional Science and Technology (CREST), Japan Science and Technology Agency (JST), Chiyoda-ku, Tokyo 102-0075, Japan ∥ Elements Strategy Initiative for Catalysts and Batteries (ESICB), Kyoto University, Kyoto 615-8520, Japan S Supporting Information *

ABSTRACT: The diffusion of the hydroxide ion in bulk water was examined by linear-scaling divide-and-conquer densityfunctional tight-binding molecular dynamics (DC-DFTB-MD) simulations using three different-sized unit cells that contained 522, 1050, and 4999 water molecules as well as one hydroxide ion. The repulsive potential for the oxygen−oxygen pair was improved by iterative Boltzmann inversion, which adjusted the radial distribution function of DFTB-MD simulations to that of the reference density functional theory-MD one. The calculated diffusion coefficients and the Arrhenius diffusion barrier were in good agreement with experimental results. The results of the hydroxide ion coordination number distribution and potential of mean force analyses supported a dynamical hypercoordination diffusion mechanism.

1. INTRODUCTION The hydroxide ion (OH−) in bulk water has been thoroughly studied over this decade because of its anomalous diffusion process and controversial solvation structure. Such anomalous diffusion occurs via proton transfer between OH− and water molecules through the hydrogen bond network, which is characterized by low activation energy.1−4 Experimental studies on the reorientation dynamics,5,6 dielectric relaxation,7 structural properties,8−18 and thermochemical properties19,20 of hydrated OH− are important for explaining its anomalous diffusivity. Furthermore, the dynamical properties of water due to the presence of OH− in solution were also studied to examine the importance of the OH− solvation shell for the diffusion process.21−24 In addition to experimental studies, previous theoretical approaches provide details of hydrogen bond networks in the first solvation structure of OH− via Car− Parinello and ab initio molecular dynamics (AIMD) simulations.25−38 For larger systems, a multistate empirical valence bond MD simulation was used to deal with the hydrated hydroxide ion.39,40 Theoretical aspects of the OH− diffusion mechanism in aqueous environments explain the importance of solvation shells for OH− diffusion.41,42 Static ab initio electronic structure calculations were also performed to describe the energetics and spectral properties of the hydrated OH− ion.43−46 Three theoretical solvation models have been proposed to interpret the above diffusion: the static hypercoordination © XXXX American Chemical Society

mechanism (SHM), mirror image mechanism (MIM), and dynamical hypercoordination mechanism (DHM).25 In the SHM model, OH− is surrounded by four water molecules, three of which act as hydrogen bond donors to the oxygen atom of OH− and the fourth as an acceptor for the hydrogen atom of OH−. The SHM scenario does not allow other water molecules to form hydrogen bonds with OH−. As a result, OH− diffuses as OH−(H2O)4 with a low diffusion coefficient. In the MIM model, OH− binds with three water molecules acting as hydrogen bond donors, whereas the hydrogen atom of OH− is not involved in hydrogen bonding. Thus, this model considers OH−(H2O)3, which is analogous to the eigenform of the hydronium ion, H3O+(H2O)3.47 The diffusion coefficient of deuterated OH− predicted by this model is approximately 6 times larger than the experimental value.25 The DHM model, being an intermediate between SHM and MIM, suggests a dynamical change of the first solvation structure of the OH− ion, namely, OH−(H2O)4 and OH−(H2O)3. Although this model was considered to reflect the observed behavior, the corresponding AIMD simulation underestimated the experimental diffusion coefficient.25 As a result, the most suitable model for interpreting the diffusion of OH − remains controversial. Received: October 22, 2016 Revised: January 22, 2017 Published: January 23, 2017 A

DOI: 10.1021/acs.jpcb.6b10659 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Heaviside step function (η(x)) with Fermi level εF and orbital energy εp is introduced. The Mulliken charge is given by

More accurate simulations are urgently required to resolve this issue, one important factor being the system size. The previous AIMD simulation used less than 128 water molecules in a unit cell, whereas the classical MD simulation suggested that the use of a larger system size is essential for obtaining an accurate diffusion coefficient value.48 In fact, reasonable results for proton transfer in bulk water were reported when at least 500 water molecules were used at the density-functional tightbinding (DFTB) level.49 As the DFTB-MD simulation is more cost-effective than AIMD, a larger unit cell and/or longer simulation time are feasible. The drawback of DFTB-MD is its parameter- or system-dependent accuracy. For example, the DFTB-MD simulation resulted in an overbinding behavior of water molecules both in the presence and absence of OH−.50 The present study examines the diffusion of OH− in bulk water by DFTB-MD simulation. To enable a large unit cell calculation, we combine DFTB-MD with a divide-and-conquer (DC) technique, achieving a linear-scaling computational cost with a low prefactor. To overcome the overbinding problem, the DFTB parameter is improved using the iterative Boltzmann inversion (IBI) protocol to fit the O−O pair radial distribution function (RDF) obtained by the AIMD simulation. The structure of this article is as follows. Section 2 describes the theoretical aspects of the present study, that is, the backgrounds of the DFTB method and the IBI protocol for modifying the O−O pair repulsion parameter. Section 3 describes the computational details of the DC-DFTB-MD simulation for OH− diffusion in bulk water and the related analyses to extract the dynamical, energetic, and structural properties. Section 4 provides a discussion of the obtained results, including the improvement of the O−O pair repulsive potential and the dynamical, energetic, and structural properties of OH− in bulk water. Finally, the concluding remarks for the present study are drawn in Section 5.

AO AO

ΔQ A =

μ∈A

1 rep ∑ UAB + 2 AB

(H − εpS)Cp = 0

+

2

⎡ 1 atom 0 + Sμν⎢ { ∑ ΔQ C(γAC + γBC) + ΔQ A2 Γ diag Hμν = Hμν AA ⎢⎣ 2 C ⎤ ⎥ + ΔQ B2 Γ diag BB } ⎥⎦ 0 =Hμν + Sμν ΩAB

(5)

subsystem DC Dμν ≈ Dμν =

α α Pμν Dμν

(6) α

The subsystem density matrix, D , is expanded with a local basis set, {χμ, μ∈L(α)}, which consists of mutually exclusive S(α) and environmental B(α) sets. The spatial local regions spanned by S(α) and B(α) are referred to as the central and buffer regions of subsystem α, respectively. Pαμν in eq 6 is the partition matrix containing the following elements α Pμν

(1)

⎧1 (μ ∈ S(α) ∧ ν ∈ S(α)) ⎪ = ⎨1/2 (μ ∈ S(α) ∧ ν ∈ B(α), or vice versa) ⎪ ⎪0 (otherwise) ⎩

Cαp ,

(7)

εαp ,

The subsystem coefficient vector, and orbital energy, are obtained by solving the eigenproblem for the local DFTB Hamiltonian (H α − εpαS α)C pα = 0 α

(8)

α

where H and S are the submatrices of H and S, respectively, in basis set L(α). Subsequently, the subsystem density matrix, Dα, is constructed as follows MO(α) α Dμν =2

∑ p

fβ (εF − εpα)CμαpCνα* p (9)

MO

In eq 9, the step function of eq 2 is modified by continuous Fermi function fβ(x) = [1 + exp(−βx)]−1, where β denotes the inverse temperature parameter. Note that the universal Fermi level, εF, is determined from the constraint on the total number of electrons in the DC-DFTB method.

Dμν = 2∑ CμiCν*i = 2∑ η(εF − εp)CμpCν*p i

∑ α

where Urep AB is a pair repulsive potential between atoms A and B; Dμν and H0μν are the density and zero-order Hamiltonian matrix elements, respectively, for atomic orbitals χμ and χν; γAB is the distance-dependent function describing the interaction between induced Mulliken charges ΔQA and ΔQB; and Γdiag AA is a diagonal term of the charge derivative of γAB. The second-order level of DFTB energy, denoted DFTB2, is described up to the third term of eq 1. The simplified density matrix for a closed-shell system is given by occ

(μ ∈ A , ν ∈ B )

For the DFTB2 formalism, the terms including ΔQ2AΓdiag AA and ΔQ2BΓdiag are deleted in eq 5. BB 2.2. DC-DFTB Method. In the DC-DFTB scheme,54 the density matrix in eq 2 is constructed as the sum of subsystem contributions

AB

1 3 ∑ Γ diag AA ΔQ A 6 A

(4)

where H is the DFTB Hamiltonian. In the case of DFTB3-diag, the explicit formula is given by

∑ DμνHμν0 1 ∑ γABΔQ AΔQ B μν

(3)

ν

where Sμν is the overlap integral, Sμν = ⟨χμ|χν⟩, and ZA is the core charge of atom A. The coefficient vector, Cp, and orbital energy, εp, are obtained by solving the eigenproblem

2. THEORETICAL ASPECTS 2.1. DFTB Method. This subsection briefly explains the fundamentals of the DFTB method.51−53 The formalism of DFTB energy up to the diagonal (or on-site) third-order term, denoted DFTB3-diag, is as follows EDFTB3 ‐ diag =

∑ ∑ DμνSνμ − ZA

p

(2)

where Cμi is a molecular orbital (MO) coefficient and i runs over occupied MOs. When arbitrary MOs are used (φp), the B


Article


O pair repulsive potential was first improved using the IBI protocol, which required the reference RDF obtained by accurate simulations. The AIMD simulation was performed using the Becke exchange58 and Lee−Yang−Parr correlation59 functional (BLYP), together with the valence double-zeta plus polarization (DZVP) basis set and Goedecker−Teter−Hutter (GTH) pseudopotentials60,61 in the CP2K code.62 The BLYP functional was chosen for the reference due to its capability to describe the DHM model for the hydroxide ion diffusion process25 compromising its limitation to describe an experimental value of water density.63,64 The used system featured 31 water molecules with one OH− ion in a unit cell of 10.794 × 10.682 × 10.606 Å3, which was determined by pre-equilibration using the classical force field. The reference O−O pair RDF was estimated by AIMD simulation for 100 ps with a time step of Δt = 0.5 fs under the Nosé−Hoover NVT ensemble. All AIMD simulations were performed using the CP2K program.62 The present study utilized the mio parameter set at the DFTB3-diag level. Similarly to several previous studies,50,55,65,66 the modified form of γAB was used in the present study to improve hydrogen-bonding interactions (denoted γh). Only the O−O pair repulsive potential was improved by the IBI protocol. DFTB-MD simulations were performed for 31 water molecules and one OH− ion in a unit cell of 10.794 × 10.682 ×10.606 Å3, which was the same as the one used in the reference AIMD simulation for updating the O−O pair repulsive potential. Each IBI protocol utilized the O−O pair RDF obtained by 100 ps DFTB-MD simulation with a time step of Δt = 0.5 fs under the Andersen NVT ensemble. The O− O pair repulsive potential, UOO(r), was updated according to eq 13. The cutoff distance of UOO(r) was extended above the value of 2.22 Å utilized for the original mio parameter51 to ensure sufficient flexibility. The shape of UOO(r) up to the original cutoff distance was generally not modified. All DFTB-MD simulations were carried out using the in-house DC-DFTB-K package,54 in connection with the IBI protocol. 3.2. DC-DFTB-MD Simulations for OH− Diffusion in Bulk Water. DC-DFTB-MD simulations using the IBIimproved O−O pair repulsive potential were performed to investigate the OH− diffusion in bulk water, which was characterized by three-dimensional periodic boundary conditions. Three different unit cells were examined, containing 522, 1050, and 4999 water molecules and one OH− ion. The box sizes of 25.044 × 25.042 × 25.046, 31.542 × 31.486 × 31.535, and 53.028 × 53.016 × 53.004 Å3 determined by the classical force field were used for 522, 1050, and 4999 water systems, which corresponded to the densities of 0.995, 1.003, and 1.003 g/cm3. The O−O pair RDFs for the 522 and 1050 water molecule systems were extracted from 10 ps trajectories with Δt = 1.0 fs under the NVT ensemble at 300 K, where the temperature was controlled using the Nosé−Hoover thermostat, whereas those for the 4999 system were obtained from a 5 ps trajectory. The temperature dependence of the 522 system was investigated by additional 10 ps simulations at 320, 360, and 400 K. Here, the box sizes of 25.109 × 25.107 × 25.111, 25.345 × 25.343 × 25.347, and 25.587 × 25.585 × 25.589 Å3 were adopted for 320, 360, and 400 K, which corresponded to the densities of 0.988, 0.960, and 0.933 g/cm3. For comparison, DC-DFTB-MD simulations for the original mio parameter set at DFTB3-diag level were additionally performed for the 522 system, where the simulation setup was

2.3. IBI Protocol. The accuracy of the DFTB method remarkably depends on the quality of the repulsive potential, 0 Urep AB , the zero-order Hamiltonian matrix, Hμν, and the distancedependent functions, γAB and Γdiag , that include semiempirical AA parameters. This method has become increasingly accurate because of the continuous improvement of parameters. Nonetheless, the previous DFTB-MD simulation failed to describe the dynamical properties of hydrated OH− in bulk water, mainly due to overbinding behavior.50 To overcome the overbinding problem, we have adopted the IBI protocol, inspired by two previous studies on the improvement of bulk structure and diffusion properties of aqueous systems.55,56 The main idea of the IBI technique is the improvement of the pair repulsive potential using the related structural properties of the bulk system. The RDF, gAB(r), is obtained from the MD trajectory as follows V p (r ) 4πr 2 AB

gAB(r ) =

(10)

In eq 10, pAB(r) is the probability density at distance r between particles A and B in volume V. In the present case, A and B were oxygen atoms. The pair repulsive potential, UAB(r), is derived with the consideration that gAB(r) is strongly correlated with the target two-body interaction. With the Boltzmann inversion scheme, the relationship between free energy FAB(r) and gAB(r) is given as57 FAB(r ) = −kBT ln gAB(r ) + c

(11)

where c is an arbitrary constant, kB is the Boltzmann constant, and T is the temperature. In the real implementation of IBI, the difference in free energy is initially computed as (0) (r ) ΔF AB

⎡ g(0)(r ) ⎤ = kBT ln⎢ AB ⎥ ⎢⎣ g (r ) ⎥⎦ AB

(12)

g(0) AB (r)

where is the A−B pair RDF obtained from the DFTBMD simulation trajectory using the original parameter set, whereas gAB(r) is the reference RDF. The new repulsive potential, U(i+1) AB (r), is iteratively generated from the previous U(i) AB(r) (i + 1) (i) (i) U AB (r ) = U AB (r ) + ΔF AB (r )

where

ΔF(i) AB(r)

(13)

is defined as

⎡ g(i) (r ) ⎤ (i) ΔF AB (r ) = kBT ln⎢ AB ⎥ ⎢⎣ g (r ) ⎥⎦ AB

(14)

In eq 14, g(i) AB(r) is the A−B pair RDF, which is calculated from the ith iteration of the IBI-generated pair repulsive potential in the DFTB-MD simulation. The process is repeated until the reference RDF is well reproduced. The present study adopts the following merit function59 to evaluate the fitting quality f(i) =

(i) (r )|2 dr ∫ w(r)|gOO(r) − gOO

(15)

Here, we adopt w(r) = exp(−r) as the weighting function to impose larger weights in the small-distance region.

3. COMPUTATIONAL DETAILS 3.1. AIMD and DFTB-MD Simulations for Improving O−O Pair Repulsive Potential. In the present study, the O− C


Article


Figure 1. (a) O−O pair RDF obtained at the DFTB3-diag/mio (red line) and BLYP/DZVP-GTH (black line) levels, and (b) O−O pair RDF at the DFTB3-diag level during updating by the IBI protocol. The reference RDF (black line) was calculated using BLYP/DZVP-GTH.

Note that δh(t) values of 1 and −1 correspond to the forward and backward shuttling processes. The Mulliken charge of oxygen atoms was used for identifying OH− in the DC-DFTB-MD trajectories. The threshold charge for the oxygen atom of the hydroxide ion was determined to be −1.000. Finally, the overall diffusion coefficient of OH−, DOH, was defined as the sum of Dv and DG. Second, we estimated the diffusion barrier of OH− as the energetic property. The energy barrier of OH− diffusion, ΔE, was estimated from the logarithmic form of the Arrhenius equation

the same as that in the case of the IBI-improved O−O pair repulsive potential. The DC calculation setup was similar to that in our previous study on proton diffusion.49 The entire system was automatically divided into several subsystems with grid boxes of 3 × 3 × 3 Å3. The buffer radius, which could include the effect of adjacent atoms into the subsystem density matrix, was set to 6 Å. Note that an automatic grid box with a reasonable buffer radius enabled the description of the dynamical change of chemical bond formation and cleavage during the DC-DFTBMD simulation. All DC-DFTB-MD simulations were carried out using the DC-DFTB-K package,54 in which the AUTO option was turned on to define the subsystems automatically. 3.3. Evaluation of Dynamical, Energetic, and Structural Properties of OH− in Bulk Water. The DC-DFTB-MD simulations described in the previous subsection were analyzed to extract the dynamical, energetic, and structural properties of OH− in bulk water. First, we evaluated the diffusion coefficient of OH− as a dynamical property. The vehicular diffusion coefficient, Dv, was estimated using the Einstein relation67 Dv = lim

t →∞

⟨|rO(t ) − rO(0)|2 ⟩ 6t

ln D = ln D0 −

l2 rp 6

(16)

intra δ = |rOinter ‐ H − rO ‐ H|

(17)

PMF(δ , n) = −kBT ln P(δ , n)

(22)

where P(δ, n) is the probability to find a value of δ and coordination number n, which is counted from the whole trajectories. When summing up P(δ, n) with respect to the coordination number, n, the probability becomes a function of δ

(18)

where δt is the time step of the MD simulation and δh(t) is defined in eq 19 ⎧ 0 (if proton does not hop) ⎪ ⎪ δh(t ) = ⎨1 (if proton hops to a new donor) ⎪ ⎪−1 (if proton hops to the last donor) ⎩

(21)

− where rinter O‑H is the distance between the OH oxygen atom and the hydrogen atom of the closest water molecule, whereas rintra O‑H is the distance between the oxygen and hydrogen atoms of the closest water molecule. The numerical value of PMF was estimated as

where l is the proton transfer length and rp is the proton transfer rate, determined from the slope of the time-course accumulation function h(t) presented in eq 18 and having an initial value of zero49,50,65 h(t ) = h(t − δt ) + δh(t )

(20)

where D corresponds to either the vehicular or overall OH− diffusion coefficient, T denotes temperature, D0 is the preexponential factor, and R is the universal gas constant (8.314 J mol−1 K−1). Third, we evaluated the potential of mean force (PMF) with respect to the coordination number of OH− as a structural property. Note that the PMF is closely related to the OH− diffusion mechanism. We defined the reaction coordinate for the proton transfer between OH− and the closest water molecule as follows

where rO(t) is the position vector of the oxygen atom at time t and ⟨···⟩ is an ensemble average obtained from NVE trajectories of DC-DFTB-MD simulations. The Grotthuss-type diffusion coefficient, DG, was estimated as41,49,68 DG =

ΔE RT

P(δ) =

∑ P(δ , n) n

(23)

On the other hand, the integral of P(δ, n) with respect to δ provides P(n) P(n) =

(19) D

∫ P(δ , n) dδ

(24) DOI: 10.1021/acs.jpcb.6b10659 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article


Figure 2. Bulk water structure obtained from the last snapshot of the 10 ps trajectory under the NVT ensemble at the DFTB3-diag/mio level (a) and the DFTB3-diag/mio-IBI-BLYP level (b).

Figure 3. (a) System-size dependence and (b) temperature dependence of the O−O pair RDF calculated at the DFTB3-diag/mio-IBI-BLYP level. The integral values of gOO(r) (NOO(r)) are shown by dashed lines.

Utilizing P(δ) and P(n), the corresponding PMF becomes a function of the distance difference, δ, and coordination number, n.

4. RESULTS AND DISCUSSION 4.1. Improvement of the O−O Pair Repulsive Potential. This subsection describes the improvement of the O−O pair repulsive potential by the IBI protocol. Before discussing the IBI results, the performance of the original mio parameter is demonstrated. Figure 1a compares the O−O pair RDF (gOO(r)), obtained by the DFTB-MD simulation at the DFTB3-diag/mio level, with that obtained by the AIMD simulation at the BLYP/DZVP-GTH level. Note that the O−O pair RDF was plotted for O−O distances of less than 5.0 Å because both DFTB-MD and AIMD simulations adopted a small system with 31 water molecules and one OH− ion. The RDF determined by the AIMD simulation has two peaks corresponding to the first and second solvation, with a minimum in between. On the contrary, the RDF of the DFTB-MD simulation has only one peak in this region, with a very shallow minimum around 3.7 Å. In particular, this peak is larger than that of AIMD and corresponds to a shorter distance, that is, 2.74 and 2.78 Å for DFTB-MD and AIMD values,

Figure 4. Time-course accumulation function h(t) obtained for DFTB3-diag/mio (black) and DFTB3-diag/mio-IBI-BLYP (blue).

respectively. This indicates that the original mio parameter shows some overbinding character. E


Article

The Journal of Physical Chemistry B Table 1. Proton Transfer Rates and Diffusion Coefficients Estimated from the DC-DFTB-MD Simulations of 522 Water Molecules with One OH− Ion at 300, 360, and 400 K at the DFTB3-diag/mio-IBI-BLYP Level exptl.a

calc.

a

T [K]

rp [ps−1]

DG [Å2 ps−1]

Dv [Å2 ps−1]

DG/Dv

DOH [Å2 ps−1]

DOH [Å2 ps−1]

300 360 400

0.38 0.48 0.79

0.40 0.50 0.82

0.27 0.82 1.27

1.48 0.61 0.65

0.67 1.32 2.09

0.53

Ref 71.

Figure 6. Probability of finding hydrogen bond donors around the OH− ion (P(n)), estimated using DFTB3-diag/mio (black) and DFTB3-diag/mio-IBI-BLYP (blue).

underestimates the water density of the liquid−vapor interface system. Using the converged O−O pair repulsive potential, we performed DFTB-MD simulations for larger unit cell sizes (522, 1050, and 4999 water molecules and one OH− ion) with the help of the DC technique. The O−O pair RDFs were obtained for distances up to 12.0 Å. Figure 3a shows the dependence of the O−O pair RDF on unit cell size, obtained using the improved mio-IBI-BLYP potentials. The height of gOO(r) at the first maximum decreases for larger unit cell sizes, being equal to 3.26, 3.24, and 3.10 for 522, 1050, and 4999 systems, respectively. On the other hand, the height of gOO(r) at the first minimum increases with increasing unit cell size, being equal to 0.34, 0.35, and 0.45 for 522, 1050, and 4999 systems, respectively. Next, the temperature dependence of the O−O pair RDF was investigated for the 522 system. Figure 3b shows the O−O pair RDFs obtained using the improved mio-IBI-BLYP potentials, demonstrating the tendency of the first gOO(r) solvation peak to decrease at higher temperatures. This behavior is interpreted by the easier dissociation of hydrogen bond networks at higher temperatures, which decreases the probability of finding the O−O pair at short distances. To compensate such events, the adjacent minimum is expected to become deeper with increasing temperature. In addition, a less pronounced second solvation peak is observed at elevated temperatures, showing the same temperature dependence as the first one. Note that a similar temperature dependence was found in a previous simulation of water using the AMOEBA water force field.69 4.2. Dynamical Properties of OH− in Bulk Water. The time-course accumulation functions (h(t) in eq 18) obtained from the first 40 ps trajectory at the DFTB3-diag/mio and DFTB3-diag/mio-IBI-BLYP levels are shown in Figure 4. Grotthuss shuttling was suppressed for DFTB3-diag/mio, where the first proton transfer occurred at 8.72 ps, with the proton immediately hopping to the last donor at 8.73 ps. The maximum resting time of the proton transfer event was 17.27 ps, with OH− diffusivity during this time being pure vehicular diffusion. The rare proton transfer event was a clear indication

Figure 5. Arrhenius plot of diffusion coefficients for the vehicular (red) and overall OH− ion (blue) diffusion at 300−400 K.

Table 2. Activation Energies of Vehicular and Hydroxide Ion Diffusions [kJ mol−1] Estimated from the Arrhenius Plota ΔEv ΔEOH a

calc.

exptl.b

15.62 11.23 (−0.90)

12.13

Deviation from the experimental value is given in parentheses. bRef 1.

Next, we improved the O−O pair repulsive potential using the IBI method, which adopted the RDF of the AIMD simulation as a reference. Figure 1b shows the iterative update of the O−O pair RDF. The overbinding character was gradually reduced, as shown by the decrease of the first peak. The convergence of the IBI update, which satisfied the condition of f(i) ≤ 0.001, was achieved at the 14th iteration. The converged O−O pair repulsive potential was denoted mio-IBI-BLYP. The shape of the converged repulsive potential is given in the Supporting Information. Upon adopting the NVT ensemble, the overbinding character of the original mio parameter led to unphysical voids (Figure 2a), which is a drawback that was pointed out in a previous study.50 On the other hand, no voids were predicted when the mio-IBI-BLYP parameter was used, as shown in Figure 2b. The removal of voids in the bulk water system is likely to destabilize the hydrogen bond networks in the overall system. In the Supporting Information, the performance of the DFTB3-diag/mio and DFTB3-diag/mio-IBI-BLYP parameters to describe a liquid−vapor system72,73 is discussed. As expected, the present DFTB3-diag/mio-IBI-BLYP method, however, still F


Article


Figure 7. (a) Calculation-method dependence and (b) coordination dependence of PMF. The definition of reaction coordinate δ is described in the text.

(18 kJ mol−1),3,4 indicating that vehicular diffusion and proton transfer jumps are mutually distinguishable in the hydrated hydroxide diffusion process. 4.4. Structural Properties of OH− in Bulk Water. As described in Section 1, three different hydroxide diffusion mechanisms were proposed: SHM, MIM, and DHM. This subsection discusses these mechanisms from the structural point of view. Figure 6 shows the probabilities of an OH− ion accepting n hydrogen bonds, that is, P(n), which were evaluated at the DFTB3-diag/mio and DFTB3-diag/mio-IBI-BLYP levels. In agreement with the overbinding behavior of its O−O pair RDF, the DFTB3-diag/mio results tend to form OH− ions fivefold coordinated by water molecules. On the other hand, as the repulsive potential was improved to reduce the O−O pair RDF overbinding, the DFTB3-diag/mio-IBI-BLYP results tended to form fourfold coordinated OH− ions. The average coordination numbers obtained using DFTB3-diag/mio and DFTB3-diag/ mio-IBI-BLYP are 4.76 and 4.03, respectively. The fourfold coordination of the hydroxide ion at the DFTB3-diag/mio-IBIBLYP level is in good agreement with the previous AIMD simulation,25 in which both DHM and SHM models showed fourfold coordinated OH− ions. Figure 7a shows plots of PMF(δ) versus the distance difference, δ, at the DFTB3-diag/mio and DFTB3-diag/mioIBI-BLYP levels. Here, the coordination number, n, was summed up. PMF(δ) values exhibit maxima at δ = 0 and minima at δ = 0.61 and 0.58 for DFTB3-diag/mio and DFTB3diag/mio-IBI-BLYP, respectively. The energy difference between the maximum and minimum corresponds to the activation barrier of proton transfer in the hydrated hydroxide system. The estimated proton transfer barrier for DFTB3-diag/ mio-IBI-BLYP is ∼1.3 kcal mol−1, whereas the corresponding value for DFTB3-diag/mio is ∼2.8 kcal mol−1. Hence, the improved O−O pair repulsive potential decreases the proton transfer barrier by up to ∼1.5 kcal mol−1. As another comparison, the proton transfer barrier estimated by Car− Parrinello MD is 1.15 kcal mol−1,25,26 which is close to the present DFTB3-diag/mio-IBI-BLYP result.

of the SHM model. Conversely, the maximum resting time for DFTB3-diag/mio-IBI-BLYP calculations was 5.37 ps. The vehicular (Dv) and Grotthuss (DG) diffusion coefficients were obtained from the DC-DFTB-MD simulations of 522 water molecules and one OH− ion at the DFTB3-diag/mio-IBIBLYP level under the NVE ensemble at 300, 360, and 400 K. The estimated diffusion coefficients at individual temperatures are shown in Table 1, along with the experimental result at 300 K. The calculated value of Dv at 300 K is reasonably close to the experimental value of neutral water diffusion (0.23 Å2 ps−1).70 The OH− diffusion (DOH) was estimated to be 0.67 Å2 ps−1 at 300 K, which is also in good agreement with the experimental value of 0.53 Å2 ps−1.71 Overall, the diffusion coefficients increased with increasing temperature. The DG/Dv ratio, representing the contribution of Grotthuss diffusion, decreased as the temperature is increased from 300 to 360 K. A further increase from 360 to 400 K does not considerably affect the ratio due to the destabilization of hydrogen-bonding interactions around the first solvation structure, manifested by the decrease of the first RDF solvation peak in Figure 3b. Conversely, our previous simulation yielded similar DG/Dv ratios for the hydrated hydronium system.49 Hence, the above ratio implicitly shows that the configuration of water molecules in the first solvation structure is crucial to force the proton transfer event for the hydrated OH− ion system. A more detailed analysis of the number of hydrogen bond donors around OH− is discussed in Section 4.4. 4.3. Energetic Properties of OH− in Bulk Water. Figure 5 shows the Arrhenius plot of diffusion coefficients for vehicular and overall OH− diffusion. As both plots were reasonably fitted by straight lines, the activation barriers of OH− diffusion were estimated using eq 20. Table 2 summarizes the estimated activation barriers and compares them with the experimental value for the overall diffusion.1 The obtained results could reproduce the experimental value within an error of 1 kJ mol−1. The barrier for vehicular diffusion is about 4.4 kJ mol−1 larger than that for overall diffusion, indicating that the proton transfer event significantly contributes to the latter. Moreover, the vehicular diffusion barrier is close to that of water molecules G


Article


Computational Science through the HPCI System Research project (Project ID: hp150267).

The proton transfer barrier estimated from PMF(δ) can be subdivided into coordination-number-dependent components that provide information on each contribution of coordination number dynamical change to the barrier height. Figure 7b shows plots of PMF(δ, n) versus the distance difference, δ, at the DFTB3-diag/mio-IBI-BLYP level. The barrier height decreases with decreasing coordination number, being equal to ∼1.8, ∼1.6, and ∼1.0 kcal mol−1 for n = 5, 4, and 3, respectively. This indicates that the decrease of coordination number promotes the proton transfer event. Consequently, the results of DFTB3-diag/mio-IBI-BLYP calculations support the DHM mechanism of OH− diffusion.

■

(1) Luz, Z.; Meiboom, S. The Activation Energies of Proton Transfer Reactions in Water. J. Am. Chem. Soc. 1964, 86, 4768−4769. (2) Loewenstein, A.; Szoke, A. The Activation Energies of Proton Transfer Reactions in Water. J. Am. Chem. Soc. 1962, 84, 1151−1154. (3) Easteal, A. J.; Price, W. E.; Woolf, L. A. Diaphragm Cell for HighTemperature Diffusion Measurements. Tracer Diffusion Coefficients for Water to 363 K. J. Chem. Soc., Faraday Trans. 1 1989, 85, 1091− 1097. (4) Halle, B.; Karlström, G. Prototropic Charge Migration in Water. J. Chem. Soc., Faraday Trans. 2 1983, 79, 1031−1046. (5) Petersen, C.; Tho̷ gersen, J.; Jensen, S. K.; Keiding, S. R. Electron Detachment and Relaxation of OH− (aq). J. Phys. Chem. A 2007, 111, 11410−11420. (6) Tho̷ gersen, J.; Jensen, S. K.; Petersen, C.; Keiding, S. R. Reorientation of Hydroxide Ions in Water. Chem. Phys. Lett. 2008, 466, 1−5. (7) Buchner, R.; Hefter, G.; May, P. M.; Sipos, P. Dielectric Relaxation of Dilute Aqueous NaOH, NaAl(OH)4, and NaB(OH)4. J. Phys. Chem. B 1999, 103, 11186−11190. (8) Megyes, T.; Bálint, S.; Grósz, T.; Radnai, T.; Bakó, I.; Sipos, P. The Structure of Aqueous Sodium Hydroxide Solutions: A Combined Solution X-ray Diffraction and Simulation Study. J. Chem. Phys. 2008, 128, No. 044501. (9) De Paz, M.; Giardini, A. G.; Friedman, L. Tandem-MassSpectrometer Study of Solvated Derivatives of OD−. Total Hydration Energy of the Proton. J. Chem. Phys. 1970, 52, 687−692. (10) Yang, X.; Castleman, A. W. Production and Magic Numbers of Large Hydrated Anion Clusters X−(H2O)n=0‑59 (X = OH, O, O2, and O3) under Thermal Conditions. J. Phys. Chem. 1990, 94, 8500−8502. (11) 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. (12) Chaudhuri, C.; Wang, Y.-S.; Jiang, J. C.; Lee, Y. T.; Chang, H.C.; Niedner-Schatteburg, G. Infrared Spectra and Isomeric Structures of Hydroxide Ion-Water Clusters OH− (H2O)1−5: A Comparison with H3O+ (H2O)1−5. Mol. Phys. 2001, 99, 1161−1173. (13) 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. (14) Imberti, S.; Botti, A.; Bruni, F.; Cappa, G.; Ricci, M. A.; Soper, A. K. Ions in Water: The Microscopic Structure of Concentrated Hydroxide Solutions. J. Chem. Phys. 2005, 122, No. 194509. (15) Robertson, W. H.; Diken, E. G.; Price, E. A.; Shin, J.-W.; Johnson, M. A. Spectroscopic Determination of the OH− Solvation Shell in the OH−·(H2O)n Clusters. Science 2003, 299, 1367−1372. (16) McLain, S. E.; Imberti, S.; Soper, A. K.; Botti, A.; Bruni, F.; Ricci, M. A. Structure of 2 Molar NaOH in Aqueous Solution from Neutron Diffraction and Empirical Potential Structure Refinement. Phys. Rev. B 2006, 74, No. 094201. ́ (17) Smiechowski, M.; Stangret, J. Hydroxide Ion Hydration in Aqueous Solutions. J. Phys. Chem. A 2007, 111, 2889−2897. (18) Cappa, C. D.; Smith, J. D.; Messer, B. M.; Cohen, R. C.; Saykally, R. J. Nature of the Aqueous Hydroxide Ion Probed by X-ray Absorption Spectroscopy. J. Phys. Chem. A 2007, 111, 4776−4785. (19) Meot-Ner, M.; Speller, C. V. Filling of Solvent Shells about Ions. 1. Thermochemical Criteria and the Effects of Isomeric Clusters. J. Phys. Chem. 1986, 90, 6616−6624. (20) Wang, Y.-S.; Tsai, C.-H.; Lee, Y. T.; Chang, H.-C.; Jiang, J. C.; Asvany, O.; Schlemmer, S.; Gerlich, D. Investigations of Protonated and Deprotonated Water Clusters Using a Low-Temperature 22-Pole Ion Trap. J. Phys. Chem. A 2003, 107, 4217−4225. (21) Nienhuys, H.-K.; Lock, A. J.; van Santen, R. A.; Bakker, H. J. Dynamics of Water Molecules in an Alkaline Environment. J. Chem. Phys. 2002, 117, 8021−8029.

5. CONCLUSIONS We modified the DFTB repulsive potential for the O−O pair using the IBI method to improve the description of the proton transfer event in the DFTB-MD simulation of the hydrated hydroxide ion system. Parameterization was performed to mitigate the overbinding issue of the first solvation peak in the O−O pair RDF by improving the two-body interaction tail to reproduce the BLYP reference. The modified parameter performance was inspected via a large-scale DC-DFTB-MD simulation of a hydrated hydroxide ion system containing up to 4999 water molecules. Fitting to structural properties caused a dramatic improvement of dynamical and energetic properties, such as the overall hydroxide ion diffusion and Arrhenius diffusion barriers. The present results support the DHM of hydroxide ion diffusion, as previously described in an AIMD simulation.

■

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.6b10659. The shape of O−O pair repulsive potential obtained from the IBI protocol and the discussion of the liquid− vapor interfacial system (PDF)

■

REFERENCES

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel: +81-3-5286-3452. ORCID

Hiromi Nakai: 0000-0001-5646-2931 Author Contributions

The manuscript was written through contributions of all authors. All authors have given their approval to the final version of the manuscript. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This study was supported in part by the FLAGSHIP2020 program of the Ministry of Education, Culture, Sports, Science and Technology (MEXT) within the priority study 5 (Development of new fundamental technologies for highefficiency energy creation, conversion/storage, and use). One of the authors (A.W.S.) acknowledges financial support from the Yoshida Scholarship Foundation (YSF). Some simulations were partially conducted using the computer resources of the Research Center for Computational Science (RCCS) and of the K computer provided by the RIKEN Advanced Institute for H


Article

The Journal of Physical Chemistry B (22) Corridoni, T.; Sodo, A.; Bruni, F.; Ricci, M. A.; Nardone, M. Probing Water Dynamics with OH−. Chem. Phys. 2007, 336, 183−187. (23) Aziz, E. F.; Ottoson, N.; Faubel, M.; Hertel, I. V.; Winter, B. Interaction between Liquid Water and Hydroxide Revealed by CoreHole De-Excitation. Nature 2008, 455, 89−91. (24) Roberts, S. T.; Petersen, P. B.; Ramasesha, K.; Tokmakoff, A.; Ufimtsev, I. S.; Martinez, T. J. Observation of a Zundel-Like Transition State During Proton Transfer in Aqueous Hydroxide Solutions. Proc. Natl. Acad. Sci. U.S.A. 2009, 106, 15154−15159. (25) 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. (26) Tuckerman, M. E.; Chandra, A.; Marx, D. Structure and Dynamics of OH− (aq). Acc. Chem. Res. 2006, 39, 151−158. (27) Ludwig, R. New Insight into the Transport Mechanism of Hydrated Hydroxide Ions in Water. Angew. Chem., Int. Ed. 2003, 42, 258−260. (28) Renault, J. P.; Vuilleumier, R.; Pommeret, S. Hydrated Electron Production by Reaction of Hydrogen Atoms with Hydroxide Ions: A First-Principles Molecular Dynamics Study. J. Phys. Chem. A 2008, 112, 7027−7034. (29) Tuckerman, M. E.; Laasonen, K.; Sprik, M.; Parrinello, M. Ab Initio Simulations of Water and Water Ions. J. Phys.: Condens. Matter 1994, 6, A93−A100. (30) Tuckerman, M.; Laasonen, K.; Sprik, M.; Parrinello, M. Ab Initio Molecular Dynamics Simulation of the Solvation and Transport of H3O+ and OH− Ions in Water. J. Phys. Chem. 1995, 99, 5749−5752. (31) Tuckerman, M.; Laasonen, K.; Sprik, M.; Parrinello, M. Ab initio Molecular Dynamics Simulation of the Solvation and Transport of Hydronium and Hydroxyl Ions in Water. J. Chem. Phys. 1995, 103, 150−161. (32) Chandra, A.; Tuckerman, M. E.; Marx, D. Connecting Solvation Shell Structure to Proton Transport Kinetics in Hydrogen-Bonded Networks via Population Correlation Functions. Phys. Rev. Lett. 2007, 99, No. 145901. (33) Tuckerman, M. E.; Marx, D.; Parrinello, M. The Nature and Transport Mechanism of Hydrated Hydroxide Ions in Aqueous Solution. Nature 2002, 417, 925−929. (34) Asthagiri, D.; Pratt, L. R.; Kress, J. D.; Gomez, M. A. The Hydration State of HO− (aq). Chem. Phys. Lett. 2003, 380, 530−535. (35) 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. (36) Chen, B.; Ivanov, I.; Park, J. M.; Parrinello, M.; Klein, M. L. Solvation Structure and Mobility Mechanism of OH−: A CarParrinello Molecular Dynamics Investigation of Alkaline Solutions. J. Phys. Chem. B 2002, 106, 12006−12016. (37) 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. (38) Zhu, Z.; Tuckerman, M. E. Ab Initio Molecular Dynamics Investigation of the Concentration Dependence of Charged Defect Transport in Basic Solutions via Calculation of the Infrared Spectrum. J. Phys. Chem. B 2002, 106, 8009−8018. (39) Wick, C. D.; Dang, L. X. Investigating Hydroxide Anion Interfacial Activity by Classical and Multistate Empirical Valence Bond Molecular Dynamics Simulations. J. Phys. Chem. A 2009, 113, 6356− 6364. (40) Ufimtsev, I. S.; Kalinichev, A. G.; Martinez, T. J.; Kirkpatrick, R. J. A Charged Ring Model for Classical OH− (aq) Simulations. Chem. Phys. Lett. 2007, 442, 128−133. (41) Agmon, N. Mechanism of Hydroxide Mobility. Chem. Phys. Lett. 2000, 319, 247−252. (42) Wannier, G. Die Beweglichkeit des Wasserstoff- und Hydroxylions in wäßriger Lösung. I. Ann. Phys. 1935, 24, 545−568. (43) Tunon, I.; Rinaldi, D.; Ruiz-Lopez, M. F.; Rivail, J. L. Hydroxide Ion in Liquid Water: Structure, Energetics, and Proton Transfer Using a Mixed Discrete-Continuum Ab Initio Model. J. Phys. Chem. 1995, 99, 3798−3805.

(44) Hermida-Ramón, J. M.; Karlström, G. Study of the Hydroxyl Ion in Water. A Combined Quantum Chemical and Statistical Mechanical Treatment. J. Phys. Chem. A 2003, 107, 5217−5222. (45) Vassilev, P.; Louwerse, M. J.; Baerends, E. J. Hydroxyl Radical and Hydroxide Ion in Liquid Water: A Comparative Electron Density Functional Theory Study. J. Phys. Chem. B 2005, 109, 23605−23610. (46) Pliego, J. R.; Riveros, J. M. On the Calculation of the Absolute Solvation Free Energy of Ionic Species: Application of the Extrapolation Method to the Hydroxide Ion in Aqueous Solution. J. Phys. Chem. B 2000, 104, 5155−5160. (47) Kirchner, B. Eigen or Zundel Ion: News from Calculated and Experimental Photoelectron Spectroscopy. ChemPhysChem 2007, 8, 41−43. (48) 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. (49) Nakai, H.; Sakti, A. W.; Nishimura, Y. Divide-and-ConquerType Density-Functional Tight-Binding Molecular Dynamics Simulations of Proton Diffusion in a Bulk Water System. J. Phys. Chem. B 2016, 120, 217−221. (50) Choi, T. H.; Liang, R.; Maupin, C. M.; Voth, G. A. Application of the SCC-DFTB Method to Hydroxide Water Clusters and Aqueous Hydroxide Solutions. J. Phys. Chem. B 2013, 117, 5165−5179. (51) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-Consistent-Charge DensityFunctional Tight-Binding Method for Simulations of Complex Materials Properties. Phys. Rev. B 1998, 58, 7260−7268. (52) Seifert, G.; Joswig, J.-O. Density-Functional Tight Binding − An Approximate Density-Functional Theory Method. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2012, 2, 456−465. (53) Yang, Y.; Yu, H.; York, D.; Cui, Q.; Elstner, M. Extension of the Self-Consistent-Charge Density-Functional Tight-Binding Method: Third-Order Expansion of the Density Functional Theory Total Energy and Introduction of a Modified Effective Coulomb Interaction. J. Phys. Chem. A 2007, 111, 10861−10873. (54) Nishizawa, H.; Nishimura, Y.; Kobayashi, M.; Irle, S.; Nakai, H. Three Pillars for Achieving Quantum Mechanical Molecular Dynamics Simulations of Huge Systems: Divide-and-Conquer, Density-Functional Tight-Binding, and Massively Parallel Computation. J. Comput. Chem. 2016, 37, 1983−1992. (55) Goyal, P.; Qian, H.-J.; Irle, S.; Lu, X.; Roston, D.; Mori, T.; Elstner, M.; Cui, Q. Molecular Simulation of Water and Hydration Effects in Different Environments: Challenges and Developments for DFTB Based Models. J. Phys. Chem. B 2014, 118, 11007−11027. (56) Doemer, M.; Liberatore, E.; Knaup, J. M.; Tavernelli, I.; Rothlisberger, U. In Situ Parameterisation of SCC-DFTB Repulsive Potentials by Iterative Boltzmann Inversion. Mol. Phys. 2013, 111, 3595−3607. (57) Reith, D.; Pütz, M.; Müller-Plathe, F. Deriving Effective Mesoscale Potentials from Atomistic Simulations. J. Comput. Chem. 2003, 24, 1624−1636. (58) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. J. Phys. A: Gen. Phys. 1988, 38, 3098−3100. (59) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37, 785−789. (60) Goedecker, S.; Teter, M.; Hutter, J. Separable Dual-Space Gaussian Pseudopotentials. Phys. Rev. B 1996, 54, 1703−1710. (61) Hartwigsen, C.; Goedecker, S.; Hutter, J. Relativistic Separable Dual-Space Gaussian Pseudopotentials from H to Rn. Phys. Rev. B 1998, 58, 3641−3662. (62) Hutter, J.; Iannuzzi, M.; Schiffmann, F.; VandeVondele, J. CP2K: Atomistic Simulations of Condensed Matter Systems. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2014, 4, 15−25. (63) Wang, J.; Román-Pérez, G.; Soler, J. M.; Artacho, E.; FernándezSerra, M.-V. Density, Structure, and Dynamics of Water: The Effect of Van der Waals Interactions. J. Chem. Phys. 2011, 134, No. 024516. I


Article

The Journal of Physical Chemistry B (64) Ho, M.-H.; Klein, M. L.; Kuo, I.-F. W. Bulk and Interfacial Aqueous Fluoride: An Investigation via First Principles Molecular Dynamics. J. Phys. Chem. A 2009, 113, 2070−2074. (65) Maupin, C. M.; Aradi, B.; Voth, G. A. The Self-Consistent Charge Density Functional Tight Binding Method Applied to Liquid Water and the Hydrated Excess Proton: Benchmark Simulations. J. Phys. Chem. B 2010, 114, 6922−6931. (66) Goyal, P.; Elstner, M.; Cui, Q. Application of the SCC-DFTB Method to Neutral and Protonated Water Clusters and Bulk Water. J. Phys. Chem. B 2011, 115, 6790−6805. (67) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Oxford Science: Oxford, UK, 1987. (68) Meiboom, S. Nuclear Magnetic Resonance Study of the Proton Transfer in Water. J. Chem. Phys. 1961, 34, 375−388. (69) Ren, P.; Ponder, J. W. Temperature and Pressure Dependence of the AMOEBA Water Model. J. Phys. Chem. B 2004, 108, 13427− 13437. (70) Mills, R. Self-Diffusion in Normal and Heavy Water in the Range 1−45.deg. J. Phys. Chem. 1973, 77, 685−688. (71) Atkins, P.; de Paula, J. Physical Chemistry, 7th ed.; Freeman: New York, 2002; p 1104. (72) Baer, M. D.; Mundy, C. J.; McGarth, M. J.; Kuo, I.-F. W.; Siepmann, J. I.; Tobias, D. J. Re-examining the Properties of the Aqueous Vapor-Liquid Interface using Dispersion Corrected Density Functional Theory. J. Chem. Phys. 2011, 135, No. 124712. (73) Baer, M. D.; Kuo, I.-F. W.; Tobias, D. J.; Mundy, C. J. Toward a Unified Picture of the Water Self-Ions at the Air-Water Interface: A Density Functional Theory Perspective. J. Phys. Chem. B 2014, 118, 8364−8372.

J


Divide-and-Conquer-Type Density-Functional ... - ACS Publications

Recommend Documents