Development of the ReaxFF Methodology for ElectrolyteâWater

Subscriber access provided by WEBSTER UNIV

A: New Tools and Methods in Experiment and Theory

Development of the ReaxFF Methodology for Electrolyte-Water Systems Mark Fedkin, Yun Kyung Shin, Nabankur Dasgupta, Jejoon Yeon, Weiwei Zhang, Diana van Duin, Adri C.T. van Duin, Kento Mori, Atsushi Fujiwara, Masahiko Machida, Hiroki Nakamura, and Masahiko Okumura J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.8b10453 • Publication Date (Web): 18 Feb 2019 Downloaded from http://pubs.acs.org on February 19, 2019

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 49 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

Development of the ReaxFF Methodology for Electrolyte-Water Systems Mark V. Fedkin,† Yun Kyung Shin,† Nabankur Dasgupta,¶ Jejoon Yeon,†,‡ Weiwei Zhang,† Diana van Duin† and Adri C. T. van Duin*,† Kento Mori and Atsushi Fujiwara§ Masahiko Machida, Hiroki Nakamura and Masahiko Okumura∥ †

Department of Mechanical and Nuclear Engineering, Pennsylvania State University, University

Park, Pennsylvania 16802, United States ¶

Department of Engineering Science and Mechanics, Pennsylvania State University, University

Park, Pennsylvania 16802, United States ‡

Center for Composite Materials, University of Delaware, Newark, Delaware 19716, United

States §

Materials Science Department, MOLSIS Inc., 3-19-9, Hatchobori, Chuo-ku, Tokyo 104-0032,

Japan ∥

Center for Computational Science & e-Systems, Japan Atomic Energy Agency, 178-4-

4 Wakashiba, Kashiwa, Chiba, 277-0871, Japan * Corresponding author: [email protected]

ABSTRACT: A new ReaxFF reactive force field has been developed for water-electrolyte systems including cations Li+, Na+, K+ and Cs+, and anions F-, Cl- and I-. The reactive force field parameters have been trained against quantum mechanical (QM) calculations related to water binding energies, hydration energies and energies of proton transfer. The new force field has been validated by applying it to molecular dynamics (MD) simulations of the ionization of different electrolytes in water and comparison of the results with experimental observations and thermodynamics. Radial distribution functions (RDF) determined for most of the atom pairs 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

(cation or anion with oxygen and hydrogen of water) show a good agreement with the RDF values obtained from DFT calculations. Based on the applied force field, the ReaxFF simulations have described the diffusion constants for water and electrolyte ions in alkali metal hydroxide and chloride salt solutions as a function of composition and electrolyte concentration. The obtained results open opportunities to advance ReaxFF methodology to a wide range of applications involving electrolyte ions and solutions. 1. INTRODUCTION Electrolyte aqueous solutions are universal medium that enables myriads of naturally occurring chemical reactions both in the liquid phase and on the interfaces. Understanding those reactions at the atomistic level is a key scientific task that would enable restoring mechanisms of many important processes, such as dissolution and precipitation of solids, ion adsorption and desorption on surfaces, degradation and synthesis of materials, ion transport, and ion exchange. Development of molecular dynamics (MD) approaches for simulating the ion interactions in the aqueous liquid phase and at the water-solid interfaces is essential for understanding mechanisms of many environmental processes. Modeling of such processes as ion solvation, adsorption, or ion interaction with solid interfaces, which typically involve chemical bond forming and bond breaking, is possible by applying quantum mechanical (QM) methods, such as density functional theory (DFT). A number of such studies have addressed the ion behavior in water and structure of ion hydration shells.1-8 However, quantum mechanical methods are usually limited to the systems of a few hundred atoms and times of a few picoseconds and become too computationally expensive when a scale-up to larger bulk systems is needed. Development of the molecular dynamics approaches based on reactive force fields has been proved to be an efficient solution to the scale-up problem. ReaxFF method9, 10 involves an empirical force field that is parameterized versus DFT data and is capable of accurately describing the bond breaking and bond formation in intermediate and large-size systems (≥ 106 atoms). The reliability of the ReaxFF simulation is dependent on the quality of the reactive force field and must be validated through comparison with DFT or experimental results. A number of non-reactive force fields have been previously developed and successfully applied to water11-15 and aqueous electrolyte systems.16-19 Tazi and co-authors17 performed 2 ACS Paragon Plus Environment

Page 2 of 49

Page 3 of 49 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


parameterization of polarizable force field for Na+, K+, Rb+, Cs+, Mg2+, Ca2+, Sr2+ and Cl- by applying DFT calculations to configurations with a single ion surrounded by 32 water molecules and to the ionic crystals containing from 16 to 108 MCl or MCl2 units. This force field was validated via MD simulations over a wide range of concentrations, from infinite dilution (ionwater ratio 1:215) to concentrated solutions (ion-water ratio 1:17) and ionic crystals, and provided good structural evidence for a number of the salts (e.g., NaCl, KCl and CsCl).17 The DFT data upon which the force field development was based appear to be in agreement with ab initio MD results for potassium and sodium ions.1 The reactive force fields for the ReaxFF method have been developed for a number of aqueous systems including metal ions or solids in contact with water.20-24 The present paper continues on that path and describes the development of the reactive force field (which would also work as non-reactive) for key monovalent ionic species, including alkali metal cations Li+, Na+, K+ and Cs+, and halide anions F-, Cl- and I-. For the force field validation, three different criteria are used: degree of aqueous dissociation, radial distribution functions (RDF), and ionic diffusion and water self-diffusion constants. The goal of this work is to provide basis for molecular dynamics modeling of aqueous and interfacial reactions with the intermediate and large size systems (more than 1000 atoms) and to extend the time scale for reaction simulation. 2. COMPUTATIONAL METHODS 2.1. ReaxFF reactive force field method. ReaxFF9, 10 is an empirical force field that employs a bond-length/bond-order/bond-energy relationship to obtain a smooth transition between bonded and non-bonded systems, thus allowing reactions to happen during molecular dynamics simulation. During an MD run, the bond orders are updated every iteration and the calculation of the system potential energy takes the bond-order into account. ReaxFF also considers nonbonded interactions for each pair of atoms, including van der Waals and Coulomb terms. ReaxFF also includes a geometry dependent, polarizable charge calculation. These charges are updated every iteration. We used the general H/O parameters from the ReaxFF water-description23 - these water parameters have been used for a wide range of systems, including proteins,25, oxides23,

27, 28

26

metal

and organic molecules29 that are all fully transferable with the electrolyte

parameters described in this work.

3 ACS Paragon Plus Environment


In this study we have optimized the reactive force field parameters for water-cation (H2O-Li+, H2O-Na+, H2O-K+, H2O-Cs+) and water-anion (H2O-F-, H2O-Cl-, H2O-I-) systems to develop a transferable ReaxFF potential for M (or X)/O/H interactions (M = Li, Na, K, Cs and X = F, Cl, I). The parameter sets were developed by training against QM data or available literature data. These training sets contained the following information: 1. Bond lengths and bond angles of the various conformations of the M+ (or X-)/water clusters or MOH (or HX)/water clusters and MX/water clusters. 2. Water binding energies of the above-described clusters. Since ReaxFF and QM use a different energy scale, ReaxFF was trained against the difference in energy between the water cluster with n water molecules and the cluster with n-1 water molecules. If more than one structure was possible, then the most stable structure was used for the cluster with n-1 water molecules. All clusters were energy-minimized using QM methods. 3. Equations of state for the condensed metal oxide and metal hydroxide (data shown in the Supporting Information), and metal halide phases. ReaxFF was trained against the difference in energy between the phase at a particular volume and the energy of the most stable volume. 4. Heats of formation for metal oxide and metal hydroxide phases. ReaxFF was trained against literature values. We used the general ReaxFF parameter optimization strategy, as first introduced for nonreactive polarizable force fields in van Duin et al.6 This involves a single-parameter linear search method. Parameter correlations – which are quite extensive in ReaxFF – are captured by performing multiple loops over the optimizable force field parameters until the force field error converges. The weights depend on (a) the relevance of a particular training set data point – typically, data points closer to the equilibrium are given higher weights and (b) the expected error in the DFT data; for example for ionic species solvated by incomplete water shells we anticipate significant uncertainty in the DFT data – as such, these points are given lower weights in the training set. We added our ReaxFF training set to the Supporting Information. This contains detailed information regarding the weights of individual training set items. 2.2. Quantum Mechanical Methods. The quantum mechanical DFT calculations of ion-water clusters have been performed with Jaguar30 and ADF2012.31-33


Page 4 of 49

Page 5 of 49 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


In ADF calculations, the PBE-D3 dispersion-corrected GGA functional34,

35

was used with

all-electron Slater type orbital basis sets.36 For alkali cation-water clusters [M+(H2O)n] (M = K, Cs), a triple-zeta + 2 polarization (TZ2P) basis set was used for geometry optimization and a larger basis set of quadruple-zeta + 4 polarization (QZ4P) for single point energies. An augmented TZ2P basis set (ATZ2P)37 was used for halide anion-water clusters [X-(H2O)n; X = F, Cl and I]. Relativistic effects were included with the scalar relativistic ZORA approach.38, 39 The following technical settings were employed: te-Velde numerical integration grid40, 41 of a high quality (keyword in ADF2012: Integration 6), ExactDensity, QZ4P fit for coulomb potential, SCF convergence of 1E-6 Hartree, and gradients convergence of 1E-4 Hartree/Å for geometry optimization. The molecular geometries were relaxed using adapted delocalized coordinates as implemented in QUILD,42 which functions as a wrapper around ADF2012. The optimized geometries were verified as local minima by frequency calculations. Reaction profiles for the selected cation-water and anion-water clusters were simulated with the QUILD ability of adding constraints, which were used to scan a potential energy surface as function of the relevant bond distances (M-O, X-H and M-X). The DFT calculations of molecular fragments, LiCl(H2O)n and KCl(H2O)n (n=1-6) clusters were performed with Jaguar software package using M06-2X functional with the 6-311++G(d,p) basis set. The initial structures of LiCl(H2O)n and KCl(H2O)n clusters were generated by taking geometries from NaCl work by Ghosh et al.43 Full geometry optimizations were performed on these clusters without using symmetry or structural constraints. To obtain the potential energy profiles along the Li-Cl and K-Cl bond dissociation, the constrained geometry optimization was applied in the M-Cl bond distance range from 1.3 Å to 5.0 Å. In addition, the dissociation energies of Li-O bond in an isolated LiOH molecule and in a LiOH molecule solvated by six water molecules were calculated. Periodic DFT calculations were performed with BAND201231, 44, 45 employing very similar settings as in the cluster calculations with ADF2012: the same dispersion-corrected GGA functional (PBE-D3),34,

35

the same scalar relativistic ZORA Hamiltonian,38,

39

the same

numerical integration and the same TZ2P basis set although in this case with a small frozen core. A dense k-space sampling was adopted (keyword in BAND2012: Kspace 5). All atomic positions in a given unit cell were optimized in a regular geometry optimization as implemented 5 ACS Paragon Plus Environment


in BAND2012, where the gradient convergence was chosen as 10-4 Hartree/Bohr. The volume values were varied by a uniform scaling and the atomic coordinates were optimized for each volume cell. Periodic DFT calculations for LiCl, KCl, and CsCl condensed phases were performed using the Vienna ab initio simulation package (VASP).46 The DFT calculations use projector augmented wave (PAW)47 pseudo-potentials and Perdew-Burke-Ernzerhof (PBE)34 exchangecorrelation functional with an energy cutoff of 520 eV. For LiCl and CsCl, the Brillouin zone is sampled with an (4 × 4 × 4) Monkhorst-Pack k-point mesh for rock-salt, wurtzite and sphalerite phases, and (6 × 6 × 6) Monkhorst-Pack k-point mesh for CsCl-type. For KCl, the Brillouin zone is sampled with an (6 × 6 × 6) Monkhorst-Pack k-point mesh for rock-salt, wurtzite and sphalerite phases, and (8 × 8 × 8) Monkhorst-Pack k-point mesh for CsCl-type. Calculations were performed with spin polarization.

3. RESULTS AND DISCUSSION 3.1. Force field development for monovalent electrolyte solutions. 3.1.1. M/O/H force field (M = Li, Na, K, Cs). Figures 13 compare the ReaxFF and the QM data in the training set for Na+, K+ and Cs+ cations in water clusters. While ReaxFF has a tendency to slightly underbind the water molecules except for [Na(H2O)2]+ and [K(H2O)2]+ clusters, the comparison to QM data is reasonably close and the difference is comparable to the error of DFT calculations. The average unassigned deviation of ReaxFF is 1.7 kcal/mol for the [K(H2O)n]+ clusters and 5.5 kcal/mol for the [KOH(H2O)n] clusters. The average unassigned deviation of ReaxFF for the [Cs(H2O)n]+ and [CsOH(H2O)n] clusters is 1.8 kcal/mol and 3.6 kcal/mol, respectively. Figures 1a and 1b describe the ReaxFF and QM water binding energies for Na+ cation and NaOH in (H2O)n clusters (n=2-6). In both systems, the water binding energy decreases as the number of water coordination increases. We note that the data for n=2 was given a relatively low weight in the training since we believe these data points to be of relatively low relevance for fully solvated Na/OH cases. Furthermore, we believe that the DFT data for such partially ionic species may carry a significant higher uncertainty than the fully solvated cases. Similarly, the strength of water binding to K+/KOH was fitted against the DFT data, placing low priority on its partially solvated state. In [KOH(H2O)5-6], it is observed that the water neighbors in the second shell help stabilizing the clusters. In the [Cs(H2O)n]+/[CsOH(H2O)n] systems 6 ACS Paragon Plus Environment

Page 6 of 49

Page 7 of 49 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


(Figure 3), the trend is well reproduced – the strength of water binding to Cs+/CsOH gradually decreases as the water neighbors in the first solvation shell increases. The relative order of water binding strength at n=2 is Cs+ < K+ < Na+.

Figure 1 QM and ReaxFF water binding energies to [Na(H2O)n]+ (a) and [NaOH(H2O)n] (b) clusters. In the cluster notation Na_xw_yw_Z/NaOH_xw_yw_Z, x describes the number of water molecules in the first shell around Na/NaOH, y describes the number of water molecules in the second shell (if available) and Z describes the group symmetry of the molecule (if available). Binding Energies were defined as E([Na(H2O)x'+y']+) + E(H2O) − E([Na(H2O)x+y]+) and E([NaOH(H2O)x'+y']) + E(H2O)− E([NaOH(H2O)x+y]), respectively.



Figure 2 QM and ReaxFF water binding energies to [K(H2O)n]+ (a) and [KOH(H2O)n] (b) clusters. In the cluster notation K_xw_yw_Z/KOH_xw_yw_Z, x describes the number of water molecules in the first shell around K/KOH, y describes the number of water molecules in the second shell (if available) and Z describes the group symmetry of the molecule (if available). Binding Energies were defined as E([K(H2O)x'+y']+) + E(H2O) − E([K(H2O)x+y]+) and E([KOH(H2O)x'+y']) + E(H2O)- E([KOH(H2O)x+y]), respectively.


Page 8 of 49

Page 9 of 49 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


Figure 3 QM and ReaxFF water binding energies to [Cs(H2O)n]+ (a) and [CsOH(H2O)n] (b) clusters. In the cluster notation Cs_xw_yw_Z/CsOH_xw_yw, x describes the number of water molecules in the first shell around Cs/CsOH, y describes the number of water molecules in the second shell (if available) and Z describes the group symmetry of the molecule (if available). Binding Energies were defined as E([Cs(H2O)x'+y']+) + E(H2O) - E([Cs(H2O)x+y]+) and E([CsOH(H2O)x'+y']) + E(H2O) - E([CsOH(H2O)x+y]), respectively.

In order to get a fair view of the energetics and structural stability of [Li(H2O)n]+ clusters, the water association energies of [Li(H2O)n]+ clusters (n ≤ 20) were also included into the force field training set. The structures and energetics of the clusters were collected from the work of Gonzalez et al.48 In their work, the global potential energy minima were searched using a basinhopping global optimization scheme. It has been known that the number of water molecules that comprises the first solvation shell of Li+ cation is not larger than four. The [Li(H2O)4]+ cluster has a tetrahedral structure as shown in Figure 4. In the clusters where the number of water molecule in the solvation shell exceeds four, the extra water molecules are not directly bound to Li+ cation and form hydrogen bonds with the first solvation shell. The association energy is 9 ACS Paragon Plus Environment


Page 10 of 49

calculated to investigate the energetics of the Li-water clusters from the chemical equation as follows.

[Li(H2O)n] + →n(H2O) + Li +

(1)

In Figure 4, it is observed that the association energy gradually decreases as the number of water molecules increases. The major deviations are observed for the n=1 and n=2 cases – these cases, in our view, are not very relevant for solvated Li-cations since these cations will always have at least four water neighbors. If we consider the first solvation shell to extend from n=1 to n=6, then ReaxFF captures the n=3-6 cases and their trends quite well. It also shows the binding energy differential between the first and second solvation shells. We also believe that the DFT data for the n=1 and n=2 cases carries a significantly higher uncertainty than for the more highly solvated cases as seen in the published binding energies obtained from ab initio calculations (Figure 4).49-51 Thus we gave these solvated cases more weight in the training.

Figure 4 The association energy of [Li(H2O)n]+ clusters (n=1-20). The ReaxFF values are compared to the work by Gonzalez et al.48 The magenta bars indicate the mean value and the range of the binding energies obtained from quantum chemical calculations using different levels of theory and basis sets in the literature - the binding energies at n > 1 were computed by us using the literature data.49-51 The structures of Li(H2O), Li(H2O)4, Li(H2O)8, and Li(H2O)12 are displayed. The colors of the atoms in the molecules represent element types – yellow (Li), red (oxygen) and white (hydrogen).


Page 11 of 49 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 bond dissociation of LiOH was investigated in gas phase and in the presence of solvation shell. As shown in Figure 5, an isolated LiOH molecule in gas phase requires a high energy to break the bond. However, a LiOH molecule solvated by six water molecules can be easily dissociated along the almost flat-bottomed potential energy curve with the very low bond dissociation energy less than 2 kcal/mol. The bond distance of Li-OH at the bottom of the potential energy curve is longer than 3.7 Å in both ReaxFF and DFT. Similar to the LiOH bond dissociation curves, when NaOH is solvated by a large number of waters, the Na-O bond of NaOH becomes longer, from 2.1 Å to 3.9 Å and weaker. We note that even when NaOH is partially solvated with only one water molecule, the dissociation energy is already very low, only 12 kcal/mol. With an increasing number of water in the solvation shell, the dissociation of NaOH becomes almost spontaneous. On the other hand, the dissociation energy of Li-OH in a highly solvated state, LiOH(H2O)5 is found to be higher than 10 kcal/mol (data not shown), inducing the difference from other alkali metals in the degree of dissociation in an aqueous solution. Figure 6 shows the potential energy profile for the Cs-OH bond dissociation in the presence of water neighbors. CsOH dissociates along the almost flat-bottomed potential energy curve – the equilibrium bond length of Cs-O in CsOH(H2O)4 cluster is longer than that in gas phase (2.51 Å) and the dissociation energy is much lower than that in gas phase (89.9 kcal/mol).52



Figure 5 Potential energy curves for the bond dissociation of (top) an isolated LiOH molecule in gas phase and a LiOH molecule solvated by six water molecules and (bottom) a NaOH molecule solvated by one and six water molecules in ReaxFF (solid line) and DFT (dashed line). The configurations of LiOH(H2O)0,6 and NaOH(H2O)1,6 clusters are shown. The colors of the atoms in the molecules represent element types – yellow (Li or Na), red (oxygen) and white (hydrogen).

Figure 6 Potential energy profiles for the CsOH bond dissociation in CsOH(H2O)4 cluster and the CsCl bond dissociation in CsCl(H2O)6 cluster in ReaxFF (solid line) and DFT (dashed line). Inset shows the zoomed-in profile of E near 3 Å. The colors of the atoms in the molecules represent element types – yellow (Cs), green (Cl), red (oxygen) and white (hydrogen). 12 ACS Paragon Plus Environment

Page 12 of 49

Page 13 of 49 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


For the condensed crystalline phases of the alkali metal oxides (M2O) and hydroxides (MOH), the equations of state are presented in the Supporting Information (Figures S1-S3). The heats of formation obtained from ReaxFF are in good agreement with literature values (Table 1). Table 1 ReaxFF and literature values for heats of formation for the crystalline phases of metal oxides and hydroxides. System

ReaxFF (kcal/mol)

Literature (kcal/mol)

Li2O

-146.1

-143.0a

LiOH

-123.9

-115.9a

Na2O

-89.4

-99.9a

NaOH

-109.5

-101.8a

K2O

-77.1

-86.8a

KOH

-106.7

-101.5a

Cs2O

-67.2

-82.6b

CsOH

-99.3

-99.6a

a Reference 53, b Reference 54

3.1.2. M/Cl force field (M = Li, Na, K, Cs). The energy-volume relationship of MCl (M = Li, Na, K, Cs) crystal phases is essential for predicting phase stability and structural phase transitions under various pressure and temperature conditions. Accordingly, the equations of state of the crystal phases, e.g., rock-salt, wurtzite, sphalerite and CsCl-type were trained against DFT data in the volume range from -30% to +40% of the ground state volume as plotted in Figure 7. The energy differences between the crystal phases are calculated by using a rock-salt phase as a reference. As one can see, the equations of state for MCl crystal phases are well reproduced in ReaxFF. The cohesive energies of the crystal phases were calculated using atomic metals Li, Na, K and Cs, and atomic Cl as references at 0 K. The cohesive energy for the rock-salt phase of LiCl (156.5 kcal/mol) agrees well with the DFT data (157.5 kcal/mol). We note that the low density phases such as wurtzite and sphalerite with a coordination number of four are also stable in LiCl unlike the larger size of alkali metal chlorides, NaCl, KCl and CsCl. In ReaxFF, the cohesive energies for these phases are 158.7 kcal/mol (wurtzite) and 159.8 kcal/mol (sphalerite), which are also close to the DFT values (158.7 kcal/mol and 158.5 kcal/mol, respectively). Thus, 13 ACS Paragon Plus Environment


the phase transition between rock-salt to wurtzite or sphalerite is much easier in LiCl than in other alkali metal chlorides. The CsCl-type is least stable with the cohesive energy of 144.1 kcal/mol in ReaxFF and 144.5 kcal/mol in DFT.

Figure 7 Equations of state of MCl (M = Li, Na, K, Cs) crystal phases - rock-salt (black), CsCltype (red), wurtzite (green) and sphalerite (blue) in ReaxFF and DFT. The volume of the crystals changes in the range from -30% to +40% of the ground state volume. The energy difference between phases was calculated using rock-salt phase as reference at 0 K. For NaCl and KCl crystal phases, the rock-salt phase is the most stable phase with the cohesive energy of 152.3 kcal/mol and 146.7 kcal/mol, respectively, which are in excellent agreement with the DFT values (152.2 kcal/mol and 146.6 kcal/mol, respectively). From the comparison of the cohesive energy of alkali metal chlorides, it is observed that the cohesive energy decreases as moving down the alkali metal group from Li to K due to the lower charge density and a weak bonding with Cl. For CsCl crystal phases, the energy-volume relationship is reproduced to an acceptable degree in ReaxFF. However, the lattice parameters are slightly underestimated and the cohesive energies are overestimated in ReaxFF – the cohesive energy of the rock-salt phase is found to be 152.1 kcal/mol in ReaxFF and 143.6 kcal/mol in DFT. We note that further improvement in CsCl condensed phase can be required for reliable modeling and 14 ACS Paragon Plus Environment

Page 14 of 49

Page 15 of 49 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


study of crystallization and dissolution. The cohesive energy difference between rock-salt and denser CsCl-type phase becomes smaller from Li to Cs, which means that the phase transition from rock-salt to CsCl-type phase is easier in CsCl crystal. To describe the structures and the stability of MCl-water clusters formed in the process of hydration, we examined the heats of hydration of MCl(H2O)n (n=1-6) clusters at different coordination numbers. The hydration energy is calculated as follows. MCl(H2O)n→MCl(H2O)n - 1 + H2O

(2)

As shown in Figures 8a-10a, the hydration energy of MCl(H2O)n clusters (M = Li, Na, K) is underestimated at low coordination numbers in ReaxFF while it is reproduced properly in the large clusters of a fully solvated state. Note that force field development requires compromises – in particular, the solvation cases with low amounts of water molecules were given a relatively low weight in the force field training. Virtually all of the ReaxFF results are within the DFTerror bar for these systems (in the 4 kcal/mol range). Note, also, that the data in Figures 5, 6 and Figure 11 are of the highest importance as these describe the transition of ion pairs into separated ions, which is critical for the reactive concepts unique for the ReaxFF capability. As such, these data were emphasized in the fitting and individual solvation energies, as depicted in Figures 1, 2, 8 and 10 were given lower weights. We believe that the discrepancies and trends obtained by ReaxFF are quite acceptable. Figures 8b-10b show the hydration energy of the LiCl(H2O)6, NaCl(H2O)6 and KCl(H2O)6 clusters at different water configurations. In LiCl(H2O)6, the configuration at the lowest hydration energy is slightly distorted cubic-like with the corners occupied by Li+ cation, Cl- anion and six water oxygens. Thus Li+ cation and Cl- anion are 3:3 coordinated to water molecules. In Na(H2O)6 and K(H2O)6, the cations prefer being coordinated by four water oxygens because the larger size of these cations opens up the geometrical possibility of a higher coordination number than Li.



Figure 8 (a) The hydration energies and structures of LiCl(H2O)n (n=2-6) clusters in ReaxFF and DFT. (b) The hydration energies and structures of LiCl(H2O)6 clusters at different water configurations in ReaxFF and DFT. The colors of the atoms in the molecules represent element types – green (Cl), yellow (Li), red (oxygen) and white (hydrogen).


Page 16 of 49

Page 17 of 49 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


Figure 9 (a) The hydration energies and structures of NaCl(H2O)n (n=2-6) clusters in ReaxFF and DFT. (b) The hydration energies and structures of NaCl(H2O)6 clusters at different water configurations in ReaxFF and DFT. The colors of the atoms in the molecules represent element types – green (Cl), yellow (Na), red (oxygen) and white (hydrogen).

The potential energy profile for the bond dissociation of M-Cl bond (M = Li, Na, K) was also investigated in MCl(H2O)n clusters. In Figure 11, the bond dissociation energy in LiCl decreases from 40.9 kcal/mol to 4.7 kcal/mol as the number of water molecule in the solvation shell increases from n=1 to n=6, which is in agreement with DFT values (48.5 kcal/mol at n=1 and 7.3 kcal/mol at n=6). Also, the bond distance of Li-Cl at the energy minimum shifts from 2.1 Å (n = 1) to 3.9 Å (n = 6). Similarly, NaCl and KCl bond dissociations are easier when these species have a larger solvation shell – the bond dissociation energy at n=6 is 7.6 kcal/mol for NaCl and 2.6 kcal/mol for KCl in ReaxFF. According to the CsCl bond dissociation curve in Figure 6, the CsCl ionic bond is weaker than KCl at n=6 and can be dissociated into Cs and Cl ions with a lower bond dissociation energy about 2 kcal/mol.



Figure 10 (a) The hydration energies and structures of KCl(H2O)n (n=2-6) clusters in ReaxFF and DFT. (b) The hydration energies and structures of KCl(H2O)6 clusters at different water configurations in ReaxFF and DFT. The colors of the atoms in the molecules represent element types – green (Cl), yellow (K), red (oxygen) and white (hydrogen).

3.1.3. X/O/H force field (X = F, Cl, I). Figures 12-14 show the ReaxFF results for the binding energy of water to various X-[H2O]6 clusters – where X- = F-, Cl- or I- anion. We use the structure definitions according to Lee and co-authors.55 The binding energy was determined by subtracting the energy for water molecule and for X-[H2O]5 cluster from the energy for X-[H2O]6 cluster. The same X-[H2O]5 cluster configuration was used in all cases – as such, the data in Figures 12-14 depict the relative stability of the anion/water cluster as well as a water binding energy. We observe that ReaxFF gives a good reproduction of the relative cluster stability as well as the overall water binding energy.


Page 18 of 49

Page 19 of 49 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


Figure 11 Potential energy curve for the M-Cl (M = Li, Na, K) bond dissociation in MCl(H2O)n (n=1-6) clusters in ReaxFF and DFT. As the number of water molecules in the solvation shell increases, the M-Cl bond distances become longer and the strength of the bonds becomes weaker. The snapshot inserted in the graph represents the structure of the LiCl(H2O)6, NaCl(H2O)6 and KCl(H2O)6 clusters. The colors of the atoms in the molecules represent element types – green (Cl), yellow (alkali metal), red (oxygen) and white (hydrogen).

Figure 15 shows the ReaxFF and DFT energies for a proton transfer reaction between HX (X=F, Cl, I) and a 5-water cluster. We find that HI fully deprotonates – the equilibrium H-I distance is close to 2.3 Å, which is far longer than the gas phase H-I bond distance (1.61 Å), indicating that the H-I bond is fully dissociated. On the other extreme, the H-F bond equilibrium distance is close to 1.15 Å, indicating that the presence of the water solvation only marginally extends the H-F bond from its gas-phase value (1.05 Å). We observe good overall agreement between ReaxFF and DFT, indicating that both HI and HCl will rapidly undergo an ionic dissociation even in incomplete water dissolution, while HF will require a full solvation to fully 19 ACS Paragon Plus Environment


dissociate. This is strongly connected to the HF/HCl/HI gas phase dissociation energies, for which DFT predicts values of 135/104.2/71.2 kcal/mol, respectively. ReaxFF slightly underpredicts these gas-phase binding energies (132.8/96.0/63.6 kcal/mol, respectively), but reproduces the weakening of the H-X bond going form F → Cl → I, which is related to its dissociation behavior upon water solvation.

Figure 12 DFT and ReaxFF F-[H2O]6 cluster stabilities. Cluster definitions were taken from Lee et al.55

Figure 13 DFT and ReaxFF Cl-[H2O]6 cluster stabilities. Cluster definitions were taken from Lee et al.55 20 ACS Paragon Plus Environment

Page 20 of 49

Page 21 of 49 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


Figure 14 DFT and ReaxFF I-[H2O]6 cluster stabilities. Cluster definitions were taken from Lee et al.55

Figure 15 DFT (dashed line with open marker) and ReaxFF (solid line with filled marker) description of the ionic dissociation path for HX (X=F/Cl/I) in the presence of five explicit water molecules.

3.2. Ionization of aqueous electrolytes from ReaxFF MD simulations. The developed reactive force field was used to simulate the ionization of several common electrolyte molecules in water at the ambient temperatures of 300 K. Hydroxides and halide salts of alkali metals (M = Li, Na, K and Cs) and halide acids are known to be strong electrolytes, which are predominantly dissociated in water at standard temperature and pressure. Those systems were extensively 21 ACS Paragon Plus Environment


studied by experimental and theoretical methods56-59 and present convenient basis for force field validation. For ReaxFF simulations, a periodic box of 28.9  28.9  28.9 Å with 792 total initial electrolyte molecules (water or electrolyte) was used. The total density of the aqueous solutions at standard conditions was close to 1.00 kg/m3 in all simulations. To model an electrolyte solution, a number of water molecules were replaced with acid, base or salt molecules. For the acid and base studies, 10 molecules of respective electrolyte were introduced, which corresponded to the molal concentration of 0.71 mol/kg (Figure 16). For metal chloride studies, the electrolyte content varied from 10 to 72 molecules of salt, which corresponds to the molal concentration range from 0.71 to 5.55 mol/kg. The electrolyte systems used for MD simulations are listed in Table 2. The ReaxFF MD simulations were run for 100 - 500 ps at 300 K, which was sufficient to reach steady state conditions. The MD output provided speciation of the electrolyte species - including cationic and anionic complexes, molecules or clusters of molecules - and these results were compared (where possible) to the experimental observations on corresponding association and dissociation reactions.

Figure 16 Water-electrolyte system simulation box with 782 water molecules and 10 MOH molecules used for ReaxFF MD simulations.


Page 22 of 49

Page 23 of 49 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 2 Electrolyte systems used for the reactive force field validation. Base systems

Acid systems

Salt systems (M = Li, Na, K)

10 LiOH + 782 H2O

10 HCl + 782 H2O

10 MCl + 782 H2O (0.71 mol/kg)

10 NaOH + 782 H2O

10 HF + 782 H2O

21 MCl + 771 H2O (1.51 mol/kg)

10 KOH + 782 H2O

10 HI + 782 H2O

28 MCl + 764 H2O (2.03 mol/kg)

10 CsOH + 782 H2O

42 MCl + 750 H2O (3.11 mol/kg) 56 MCl + 736 H2O (4.22 mol/kg) 72 MCl + 720 H2O (5.55 mol/kg)

3.2.1. Ionization of alkali metal hydroxides in water. Ionization of metal hydroxide solution was studied by using both dissociation and association reactions - direct and reverse of the general process: MOH ⟺ M + + OH ―

(3)

In most cases good correspondence was achieved between dissociation and association runs over the simulation time. The electrolyte strength is expected to increase in the series LiOH - NaOH - KOH - CsOH based on thermodynamic basicity constants.56,

59

Hence, progressively more dissociation is

observed when moving down the alkali metal group from Li to Cs. The reactive force field developed in this work allowed this trend to be reproduced in ReaxFF MD simulations as shown in Figure 17. LiOH electrolyte was predominantly present in molecular (non-dissociated) form, existing as a mixture of 30% Li+ and 70% Li(OH)x in water. On the other hand, both NaOH and KOH exhibited partial ionization − predominantly dissociated with occasional association of ions up to 20–30%, and CsOH was completely dissociated with all Cs present in the ionic form. Dissociation simulation in the CsOH-H2O system also showed the fastest aqueous dissociation in water at 300 K.



Figure 17 Presence of M+ species in aqueous MOH solutions (System: 10 MOH + 782 H2O) at 300 K based on ReaxFF MD simulations. CsOH - green, KOH - black, NaOH – red and LiOH blue.

3.2.2 Ionization of halogenide acids in water. The ReaxFF MD simulation results for halogenide acid dissociation are shown in Figure 18. Out of the monoprotic halogenide acids, only HF exhibits partial dissociation behavior, with its dissociation constants values ranging between 10-3 and 10-4, according to the literature.58 That suggests that the fraction of F- ions in the HF solutions is expected to be less than 1% of total fluorine. The majority of HF remains in the molecular form at standard conditions and that fact justifies certain properties of this compound in practical chemistry. The ReaxFF MD simulations were performed in a small system with 10 HF molecules, and over the course of MD simulation for over 100 ps, on the average 9 out 10 molecules remained in the form of HF or H2F, and 1 out of 10 molecules exhibited dissociation into F- ion and proton. The distribution between HF and H2F+ species, which also bound a number of water molecules into the complex, dynamically varied over time, but the total ratio F/F(total) remained quite constant. Based on the average concentrations of species in the system throughout the simulation, the Ka(HF) value was determined at 7.7  10-3 which is slightly higher but close to the literature value of 6.6 10-4.58 Due to relatively small system size used in this validation, more accurate estimation of the Ka value could not be pursued, but the obtained result is consistent with expected behavior. Both HCl and HI acids are known as strong acids, which dissociate completely into proton and halide ion in water at standard conditions. That renders the practical attempt at estimating the 24 ACS Paragon Plus Environment

Page 24 of 49

Page 25 of 49 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


dissociation constant of those compounds unnecessary, since statistically negligible amounts of associated HClo or HIo species would produce very high and inaccurate values for the equilibrium constants. ReaxFF molecular dynamics simulations, however, allow approximation of the aqueous dissociation constant by taking into account random occurrence of HCl and HI species at certain moments during the simulation. Typically those occurrences are short-lived and the associated molecules readily break apart to form ions. Thus on the average, 93.5% of total HCl is present in Cl- form, which yields the Ka value at around 9.7. For HI, complete dissociation (100% I-) is observed most of the time, with Ka value undefined. This is completely consistent with the known chemistry of HI, which is the strongest acid of all hydrogen halides.60

Figure 18 Dissociation of halogenide acids in HX-H2O (X = F, Cl, I) system at 300 K based on ReaxFF MD simulations.

3.2.3 Ionization of alkali metal chlorides in water. Ionization in alkali metal chloride solutions was studied as a function of electrolyte concentration at 300 K. The systems studied are listed in Table 2. While the dilute solutions of NaCl and KCl are fully dissociated, more association starts to occur in more concentrated solutions, with eventual clustering of salt molecules to form particles or crystallites. The average concentrations of free metal ions in chloride solutions of various concentrations are shown in Figure 19. The fraction of the associated molecules rapidly increases from 1% to ~80% in NaCl and from 1% to ~50% as the concentration increases from 0.71 to 5.55 mol/kg. At the same time, the increased aggregation of MCl molecules occurs as shown by the growth of the average molecule size on the right-hand plot in Figure 19. Examples of such aggregates or 25 ACS Paragon Plus Environment


clusters formed in concentrated NaCl and KCl solutions are shown in Figure 20. More pronounced associative behavior is observed in NaCl compared to KCl, and correspondingly NaCl produces crystallites of larger maximum size compared to KCl. The clusters appear in the form of MClMCl chains and may be considered precursors of growing MCl crystals in the long run.

Figure 19 Comparison of NaCl and KCl association behavior at different concentrations based on ReaxFF MD simulations at 300 K: (left) percentages of free M+ ions (dissociated form) and (right) the average MCl cluster size (associated form).

Figure 20 MxCly clusters formed in 5.55 mol/kg NaCl (left) and KCl (right) solutions - snapshots taken from ReaxFF MD simulations.


Page 26 of 49

Page 27 of 49 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


3.3. Radial distribution functions. The radial distribution functions (RDF), g(r), defined as the probability of locating an atom at a distance 𝑟 from the origin, were calculated from the ReaxFF MD simulations in various electrolyte solutions. The RDF values were obtained for the ion– oxygen (water) and ion–hydrogen (water) pairs in aqueous solutions of metal hydroxides and hydrogen halogenide acids. The other set of RDF data was obtained for the aqueous solutions of alkali metal chlorides, with consideration of the concentration dependence. RDF is one of the characteristics that allowed validating the force field against DFT data or experiment. All data in this section have been obtained at 300 K. The same simulation box, as described above for electrolyte ionization studies, was used for RDF calculation. For the metal hydroxide (MOH) and hydrogen halogenide (HX) solutions, the system contained 10 molecules of electrolyte per 782 molecules of water, which results in the molal concentration of 0.71 mol/kg. For the metal chloride (MCl) solutions, two electrolyte concentrations - 0.71 mol/kg (10 salt molecules) and 4.22 mol/kg (56 salt molecules) were used to observe the concentration effect on RDF. The total density of the aqueous solutions at standard conditions was close to 1.00 kg/m3 in all simulations. The NVT MD simulations at 300 K were preceded by the conjugate gradient style energy minimization. The system was allowed to equilibrate for 50 ps in MOH and HX systems and 300-500 ps in MCl systems prior to performing RDF measurements. The MD simulations were performed with a 0.25 fs time step and 100 fs temperature damping constant. The Berendsen thermostat was used for temperature control. The RDF patterns obtained by ReaxFF MD simulations were compared for different MO(water) and MH(water) pairs in metal hydroxide solutions, which are shown in Figures 21a and 21b. The RDF patterns for anionO(water) and anionH(water) pairs are shown in Figures 22a and 22b. The tested cations exhibit the RDF distribution that is consistent with their ionization behavior. Li hydroxide is more prone to being associated with hydroxyl at this concentration as observed in Figure 17, and the narrow Li-O(water) RDF peak centered at 2.40 Å indicates a shift toward longer distances due to strong Li-OH association but still relatively strong interaction with surrounding water. NaOH is partially dissociated at this concentration and the first peak is positioned at 2.50 Å. At the same time KOH which is also partially dissociated in water has a broader peak at 2.75 Å, which indicates less strong interaction with the solvation shell and relatively frequent exchange between the solvation shells. CsOH has been shown to be 27 ACS Paragon Plus Environment


fully dissociated at this concentration and temperature condition, which is also reflected by much more distant positions of its broad first RDF peak at 3.35 Å (Figure 21a). The RDF patterns for the M-H pairs as well as more pronounced RDF peaks for Li-H and Na-H compared to K-H and Cs-H, show the progression of H mobility in the sequence Li-Na-K-Cs, which is a sign of stronger binding of water or hydroxyls to Li and Na compared to K and Cs (Figure 21b).

Figure 21 Radial distribution functions for (a) cation–oxygen (water) and (b) cation–hydrogen (water) pairs obtained from the ReaxFF simulations for MOH (M = Li, Na, K, Cs) solutions.


Page 28 of 49

Page 29 of 49 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


Figure 22 Radial distribution functions for (a) anion-oxygen (water) and (b) anion-hydrogen (water) pairs obtained from the ReaxFF simulations for HX (X = F, Cl, I) solutions. Based on the RDF peaks for anion-oxygen, the shortest radius was exhibited by F-, 2.45 Å, while Cl- showed a peak much farther from water dipoles, and I- peak is not very well expressed (Figure 22a). The RDF peaks of anion-H (water) pairs logically show much shorter distances between the anions and hydrogen compared to much longer M-H distances in hydroxide solutions due to repulsive forces between the positive end of a water dipole and a positively charged metal ion. The position of the RDF peak for F- at 1.45 Å possibly indicates a strong affinity of this anion to hydrogen and high probability of existing bond. The high peak of Cl- at 2.10 Å and less pronounced peak for I- at 2.35 Å show the likeliness of bond with H, but it is not as pronounced as for F- (Figure 22b). The RDF data determined from the ReaxFF simulations were compared to the values obtained by Car-Parrinello Molecular Dynamics (CPMD) and experimental data available in the literature for similar systems. The numerical data on RDF main peak positions for the studied systems were collected in Tables 3 and 4. Thus the position of the first peak for ion-water interactions matched well with previously published results.1, 17, 19 For the comparison of different studies, the magnitude of the g(r) peak is not the best matching criterion since it depends on system density, 29 ACS Paragon Plus Environment


which may be different in different studies and sometimes is not reported. For proper comparison, integrated RDF were plotted for some of the ions along with the ranges of values obtained by Bankura and co-authors1 via a number of different computational methods (BLYP, BLYP-D, HCTH, HCTH-D, PBE, PBE-D, TIP4P, SWM4-NDP) (Figure 23). The proximity of the integrated RDF range indicates that developed reactive force field parameters quite accurately describe ion - water interactions. The best agreement is observed for the K-H and K-O pairs. Table 3 The RDF first peak positions obtained by ReaxFF MD simulations in metal hydroxide (MOH, M = Li, K, Na, Cs) aqueous solutions in comparison with data from other sources. Ion - (H/O)water

r (Å)

Method

Source

Li – H

2.65

ReaxFF

this study

2.60

Hartree-Fock

Clementi et al.61

2.40

ReaxFF

this study

1.97

ab initio

Feller et al.62

2.15

ab initio

Glendening and Feller50

2,26

RDVH FF

Pethes 63

3.20

ReaxFF

this study

2.98

CPMD

Ikeda and Boero64

2.98

CPMD

Bankura et al.1

2.50

ReaxFF

this study

2.40

CPMD

Ikeda and Boero64

2.40

CPMD

Bankura et al.1

2.40

Polarizable FF

Tazi et al.17

2.45

MD-SPC/E

Lee and Rasaiah65

3.45

ReaxFF

this study

3.40

CPMD

Ikeda and Boero64

3.25

CPMD

Bankura et al.1

2.75

ReaxFF

this study

2.85

CPMD

Ikeda and Boero64

2.75

Polarizable FF

Tazi et al.17

2.85

CPMD

Bankura et al.1

2.80

MD-SPC/E

Lee and Rasaiah65

Cs – H

3.95

ReaxFF

this study

Cs - O

3.35

ReaxFF

this study

3.20

CPMD

Ikeda and Boero64

Li – O

Na – H

Na - O

K–H

K–O


Page 30 of 49

Page 31 of 49 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


3.20

Polarizable FF

Tazi et al.17

3.05

MD-SPC/E

Lee and Rasaiah65

Table 4 The RDF first peak positions obtained by ReaxFF MD simulations in hydrogen halogenide (HX, X = F, Cl, I) aqueous solutions in comparison with data from other sources. Ion - (H/O)water

r (Å)

Method

Source

Cl - H

2.05

ReaxFF

this study

2.10

CPMD

Bankura et al.1

3.05

ReaxFF

this study

3.10

CPMD

Bankura et al.1

3.20

MD-SPC/E

Lee and Rasaiah65

1.45

ReaxFF

this study

1.50

CPMD

Izvekov and Voth19

1.55

MD

Chaumont and Wipff66

1.95

QMCF MD

Kriptayakornupong et al.67

2.45

ReaxFF

this study

2.60

MD-SPC/E

Lee and Rasaiah65

2.35

ReaxFF

this study

2.60

CPMD

Heuft and Meijer68

3.60

ReaxFF

this study

3.60

MD-SPC/E

Lee and Rasaiah65

3.50

CPMD

Heuft and Meijer68

Cl - O

F-H

F-O I-H I-O

The RDF peaks obtained from ReaxFF simulations for metal chlorides - LiCl, NaCl, and KCl - were used to demonstrate the concentration dependence of electrolyte properties. As shown above, concentrated electrolyte solutions undergo less dissociation, and in this case cations more readily recombine with Cl-, and formation of clusters or nano-crystallites can occur. This tendency is reflected in RDF patterns plotted for two concentrations  0.71 and 4.22 mol/kg (Figures 24-26). For both KCl and NaCl solutions, it is seen that the first RDF peak for M-Cl pair increases with increasing concentration. The first peak of K-Cl pair which is almost negligible at 0.71 mol/kg increases and forms a dual peak at 3.65 Å and 5.05 Å, indicating the formation of KCl clusters. In NaCl, the RDF of Na-Cl pair exhibits several distinct peaks even at low concentration (0.71 mol/kg). From the integration of the RDF, each Na+ is surrounded by 0.090 Cl- ions and five water molecules. In addition, it was observed that only one Na-Cl ion pair 31 ACS Paragon Plus Environment


was occasionally associated during the MD simulation. These indicate that the structure of NaCl in solution is ordered, yet no NaCl clusters are indeed formed at 0.71 mol/kg. The magnitude of the RDF peaks is higher for NaCl than for KCl, which also agrees with more dissociative behavior of KCl observed in MD simulations. For the RDF of M-O pair, the first peak decreases slightly with the increase of the salt concentration, indicating a lower degree of solvation of the ion at higher concentration. The first solvation shell of Na+ has an average number of 4.7 water molecules at 0.71 mol/kg and 3.4 water molecules at 4.22 mol/kg, indicating that Na+ is partially solvated at higher concentration. Similarly, K+ is solvated by 6.0 water molecules at 0.71 mol/kg and 5.2 water molecules at 4.22 mol/kg. The RDF and ion coordination numbers at low concentration are well comparable to the DFT and experiments.6, 69, 70 The M-H RDF in NaCl shows a dual peak, which evidently corresponds to the two hydrogens of water in the ion hydration shell. In case of KCl, the M-H peak is much broader; suggesting that hydration shell of K+ is much less rigid.


Page 32 of 49

Page 33 of 49 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


Figure 23 Comparison of the ReaxFF integrated RDF with CPMD and classical MD simulation data1 for the same atom pairs. The black bars indicate the range of published values obtained by different computational methods (BLYP, BLYP-D, HCTH, HCTH-D, PBE, PBE-D, TIP4P, SWM4-NDP).



Figure 24 Radial distribution functions (black line) and coordination numbers (red line) obtained by ReaxFF in the 0.71 mol/kg and 4.22 mol/kg KCl(aq) solutions.


Page 34 of 49

Page 35 of 49 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


Figure 25 Radial distribution functions (black line) and coordination numbers (red line) obtained by ReaxFF in the 0.71 mol/kg and 4.22 mol/kg NaCl(aq) solutions.



Figure 26 Radial distribution functions (black line) and coordination numbers (red line) obtained by ReaxFF in the 0.71 mol/kg and 4.22 mol/kg LiCl(aq) solutions.

Behavior of LiCl appears to be different. First, the Li-O RDF shows the narrower first peak at a shorter distance compared to those of NaCl and KCl, confirming the strong water binding by Li+. Each Li+ is solvated by an average number of 5.5 water molecules at 0.71 mol/kg and 5.2 water molecules at 4.22 mol/kg, which is in the range of 4-6 water molecules estimated from experimental studies.71-73 Second, the Li-Cl peak is located much further away from the central metal ion, showing that the water shell essentially isolates Li from Cl under any conditions studied. This is not observed in NaCl or KCl solutions, where M-Cl peak appearing in concentrated solutions is close in position to the M-O peak, meaning both water and Cl- share the space around the central cation. The RDF at LiCl concentrations above 10.1 mol/kg in Figure 27 has a sharp peak at 2.25~2.35 Å, which is positioned at a slightly shorter radius than the peak of rock-salt lattice structure (2.49 Å), and is very close to the peak in sphalerite and wurtzite crystal 36 ACS Paragon Plus Environment

Page 36 of 49

Page 37 of 49 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


structures (2.37, 2.38 Å). Like NaCl and KCl, there is no pronounced concentration dependence of the RDF peak for Li-H pair between the two marginal concentrations tested. Behavior of Cl- with respect to water species is very similar in all three MCl solutions, as indicated by almost identical locations of the primary Cl-O and Cl-H RDF peaks. The RDF peaks for the M-M pairs help understand the structure of salt clusters that start to form in relatively concentrated solutions. It is best seen in the case of NaCl, where a distinct double peak indicates quasi-crystalline aggregation of NaCl molecules, with the first and second neighbor atoms at 4.75 Å and 7.15 Å, respectively. Interestingly, the average distance between the Na atoms in a NaCl crystal is around 4.78 Å, taking into account the diagonal Na-Na and NaCl-Na distances in the face-centered lattice. Such a close correspondence of the Na-Na distances infers crystalline nature of the NaCl clusters even at very early stages of electrolyte association. This is not clearly observed in the KCl case, where the K-K RDF peak is very broad and cannot be easily correlated to specific lattice parameters.

Figure 27 Radial distribution functions of Li-Cl pair obtained by ReaxFF at LiCl concentrations of 0.72 mol/kg, 4.28 mol/kg, 10.0 mol/kg and 15.0 mol/kg.

3.4. Diffusion constants. ReaxFF MD atom trajectories were used to calculate diffusion constants using mean square displacements (MSD) of oxygen or hydrogen atoms in water molecules as well as metal ions in electrolyte solutions. This information may be useful in understanding how the mobility of electrolyte species varies in different electrolytes and at different concentrations. The MSD method of diffusion constant calculation has been described previously.74 The mean square displacement of an atom is defined as follows: 37 ACS Paragon Plus Environment


MSD(m)  r (t )  r

2

1 n   r (m  i )  r (i ) n i 1

Page 38 of 49

2

(4)

where r is the position of an atom, t is the time, m + i = k is the total number of frames read in the simulation, m is the maximum number of points allowed for the MSD calculation (m was taken as k/2 in our calculations), n is the number of data points used for averaging, and i is the step counter. The diffusion constant is further obtained by Einstein formula: 1

𝐷𝑀𝐷 = 6𝑡〈|𝑟(𝑡) ― 𝑟|2〉

(5)

It has been reported that in MD simulations of condensed phase under periodic boundary conditions, long-range interactions, for instance, hydrodynamic interactions can lead to significant finite system-size effects and cause deviation of the calculated self-diffusion constants from the infinite-system limit.75-77 To take into account these systematic errors, we applied the system-size corrections to the calculation of diffusion constants, which is proposed by Yeh and Hummer.78 ξ𝑘𝐵T

(6)

𝐷 = 𝐷𝑀𝐷 + 6𝜋ηL

where D is the diffusion constant in the thermodynamic limit by adding the corrections, DMD is the diffusion constant computed from the MD simulations by Eq. (5),   2.837297 is a dimensionless constant for a periodic cubic simulation box, kB is the Boltzmann constant, T is the temperature in K,  is the shear viscosity of the system at temperature T and L is the side length of the simulation box. The shear viscosity is calculated using the Green-Kubo approach, the autocorrelation of the off-diagonal components of the stress tensor as follows: 1 V

〈(

𝑡

η = lim 2𝑡𝑘𝐵T ∫0𝑃𝛼𝛽(𝑡′)𝑑𝑡′ 𝑡→∞

)〉 2

(7)

V is the volume of the system and P is three off-diagonal components of the stress tensor. The simulation details and the computed shear viscosities of MOH (0.71 mol/kg) and NaCl and KCl solutions (0.71, 1.51, 2.03, 3.11 and 4.22 mol/kg) are summarized in the Supporting Information


Page 39 of 49 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


(Tables S1-S2). Diffusion constants (D) for alkali metal hydroxide solutions are shown in Table 5. Self-diffusion constants of water are provided as reference. The self-diffusion constant of pure water obtained from MSD analysis of oxygen atoms in ReaxFF simulations is slightly higher but reasonably compared to the experimental value of 0.246 Å2/ps.79 It is generally expected that the diffusion constant of metal ions will depend on the degree of dissociation of the electrolyte. Thus D(M) is found to increase from LiOH to CsOH, the trend that correlates with the partial dissociation of LiOH(aq) (30%) to the complete dissociation of CsOH(aq) (100%). The self-diffusion constants of water are shown to be affected to a various degree by present electrolytes. Table 5 Diffusion constants using Eq. (6) for water and cations in different MOH aqueous solutions (0.71 mol/kg) at 300 K determined from ReaxFF simulations. The D(O), D(H), and D(M) values were measured by the mean square displacement of oxygen, hydrogen, and metal atoms, respectively. Electrolyte

D(O), Å2/ps

D(H), Å2/ps

D(M), Å2/ps

H2O

0.2739

0.3009

-

LiOH(aq)

0.2865

0.2895

0.1395

NaOH(aq)

0.2635

0.2885

0.2045

KOH(aq)

0.2398

0.2618

0.2468

CsOH(aq)

0.2652

0.2902

0.2602

The diffusion constants as a function of electrolyte concentration in NaCl(aq) and KCl(aq) solutions are shown in Figure 28. The water self-diffusion constants gradually decrease with increasing concentration of electrolyte. Similar trend has been observed in the previous experiments and simulations.80,

81

There is an excellent agreement of the obtained water self-

diffusion constants with the experimental values in 2.0 mol/kg NaCl (D = 0.195 Å2/ps) reported by Easteal and Woolf82 and in 0.7 mol/kg NaCl (D = 0.213 Å2/ps) and 0.7 mol/kg KCl (D = 0.216 Å2/ps) based on the work by Tanaka.79 Similarly, the diffusion constant measured for alkali metal ions decreases with increasing electrolyte concentration (Figure 28b), which agrees with the experimentally observed trend.83 The diffusion constants for Cl- ions are underestimated by ReaxFF, being on the order of 0.05 Å2/ps compared to average of 0.15 Å2/ps experimentally 39 ACS Paragon Plus Environment


measured in the work of Vitagliano and Lyons.80 Close agreement of the ReaxFF data and the experiment for Na+ (Figure 29) provides a good validation for the used force field in this study.

Figure 28 Diffusion constants in NaCl(aq) and KCl(aq) solutions as function of electrolyte concentration: (a) self-diffusion of water, (b) diffusion of metal ions Na+ or K+.


Page 40 of 49

Page 41 of 49 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


Figure 29 The diffusion constant of Na+ in NaCl(aq) solution as a function of concentration at 300 K. Experimental data from Vitagliano and Lyons.80 Note that the unit of concentration is converted to molarity (mol/L) to compare with the experimental data.

4. CONCLUSIONS As a result of this study a new reactive force field was developed for the water-electrolyte system including alkali metal cations (Li+, Na+, K+, Cs+) and halogen anions (F-, Cl-, and I-), which would represent a variety of salt, acid, and base solutions present in both techno-industrial systems and natural environments. The applicability of this force field for molecular dynamics simulations is therefore very wide. ReaxFF has been shown to accurately reproduce the trend of increasing ionization (aqueous dissociation) in the series LiOH, NaOH, KOH, and CsOH, which is known from the existing experimental data and thermodynamic considerations. Furthermore, the performed MD simulations create an insight in the concentration dependence of alkali metal chloride solutions, up to the region of supersaturation, where MCl molecule dissociation becomes limited, and clustering of single molecules into clusters and crystallites is observed. Those formations in highly concentrated salt solutions clearly represent the initial stage of salt crystallization and can provide mechanistic information for crystal growth phenomena and colloidal stability in brines. The developed force field has been validated through calculation of radial distribution functions for the main ion pairs, and close agreement with the published DFT data has been achieved. Significant growth of the first RDF peak for the M-Cl pairs (M = Li, Na, K) upon 41 ACS Paragon Plus Environment


increasing the electrolyte concentration provides additional evidence for ion pair association in concentrated solutions, and that result is in agreement with the electrolyte ionization behavior described above. The MSD analysis of the ReaxFF simulation data provides the values of the diffusion constants for all elements involved in the electrolyte system. It is clearly seen that water selfdiffusion is significantly affected by the presence of electrolyte species, while the diffusion constants of alkali metal ions regularly drop with increasing concentration - another fact that is observed experimentally and readily confirmed in this study. At this stage, we did not find sufficient agreement between the simulated and experimental values of halogen diffusion constants, hence further confirmation will be needed. The ability to describe computationally the electrolyte ion diffusion in aqueous media using reactive MD simulation tool would be a significant step in understanding the mechanisms of ion transport in confined media - interfaces, clays, and two-dimensional structures, such as MXenes, which enclose electrolyte species in interlayer spaces, and are extensively investigated as new promising materials for energy storage. With the achieved validation, the electrolyte modeling capabilities provided by the ReaxFF open new opportunities to explore new systems and materials with lower experimental and analytical costs and to make effective predictions for properties under extreme environments, such as high concentrations, elevated temperature and pressure, and aggressive chemistry. This work has contributed to the advancement of the ReaxFF method to wider range of applications. ASSOCIATED CONTENT Supporting Information The Supporting Information is available free of charge on the ACS Publication website. Figures S1-S3 comparing the ReaxFF and QM equations of state for the condensed alkali metal hydroxide and alkali metal oxide phases. Shear viscosities computed from all MD simulations for MOH and MCl solutions. ReaxFF training set for electrolyte/water system. ReaxFF reactive force field parameters for electrolyte/water system. AUTHOR INFORMATION Corresponding Author 42 ACS Paragon Plus Environment

Page 42 of 49

Page 43 of 49 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


*E-mail: [email protected] Notes The authors declare no competing financial interest. ACKNOWLEDGEMENTS M.M., M.O., and H.N. (JAEA) were partially supported by JAEA-NIMS cooperative project organized by T. Yaita and F-Trace Project. The DFT computational contribution by ADF and BAND was done on supercomputers in CCSE, JAEA. We thank the Director of CCSE, H. Takemiya and all other CCSE staff members for their supports. M.M. (JAEA) is partially supported by a Grant-in-Aid for Scientific Research (B) (Grant No.16H04624) from JSPS. M.F., Y.S., W.Z., N.D. and ACTvD (Penn State) acknowledge support by the Laboratory Directed Research and Development (LDRD) program of Sandia National Laboratories. Sandia National Laboratories is a multi-mission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525. REFERENCES 1. Bankura, A.; Carnevale, V.; Klein, M. L. Hydration structure of salt solutions from ab initio molecular dynamics J. Chem. Phys. 2013, 138, 014501. 2. Rempe, S. B.; Pratt, L. R. The hydration number of Na(+) in liquid water Fluid Phase Equilib. 2001, 183, 121-132. 3. Cavallari, M.; Cavazzoni, C.; Ferrario, M. Structure of NaCl and KCl concentrated aqueous solutions by ab initio molecular dynamics Mol. Phys. 2004, 102, 959-966. 4. White, J. A.; Schwegler, E.; Galli, G.; Gygi, F. The solvation of Na+ in water: First-principles simulations J. Chem. Phys. 2000, 113, 4668-4673. 5. Ikeda, T.; Boero, M.; Terakura, K. Hydration of alkali ions from first principles molecular dynamics revisited J. Chem. Phys. 2007, 126, 034501. 6. Bucher, D.; Guidoni, L.; Carloni, P.; Rothlisberger, U. Coordination Numbers of K+ and Na+ Ions Inside the Selectivity Filter of the KcsA Potassium Channel: Insights from First Principles Molecular Dynamics Biophys. J. 2010, 98, L47-L49. 7. Rowley, C. N.; Roux, B. The Solvation Structure of Na+ and K+ in Liquid Water Determined from High Level ab Initio Molecular Dynamics Simulations J. Chem. Theory Comput. 2012, 8, 3526-3535. 8. Ramaniah, L. M.; Bernasconi, M.; Parrinello, M. Ab initio molecular-dynamics simulation of K+ solvation in water J. Chem. Phys. 1999, 111, 1587-1591.



9. van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A. ReaxFF: A reactive force field for hydrocarbons J. Phys. Chem. A 2001, 105, 9396-9409. 10. van Duin, A. C. T.; Strachan, A.; Stewman, S.; Zhang, Q. S.; Xu, X.; Goddard, W. A. ReaxFF(SiO) reactive force field for silicon and silicon oxide systems J. Phys. Chem. A 2003, 107, 3803-3811. 11. Sedlmeier, F.; Horinek, D.; Netz, R. R. Spatial Correlations of Density and Structural Fluctuations in Liquid Water: A Comparative Simulation Study J. Am. Chem. Soc. 2011, 133, 1391-1398. 12. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water J. Chem. Phys. 1983, 79, 926-935. 13. Vega, C.; Abascal, J. L. F. Simulating water with rigid non-polarizable models: a general perspective Phys. Chem. Chem. Phys. 2011, 13, 19663-19688. 14. Abascal, J. L. F.; Vega, C. Dipole-quadrupole force ratios determine the ability of potential models to describe the phase diagram of water Phys. Rev. Lett. 2007, 98, 237801. 15. Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials J. Phys. Chem. 1987, 91, 6269-6271. 16. Wang, J. M.; Cieplak, P.; Kollman, P. A. How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules? J. Comput. Chem. 2000, 21, 1049-1074. 17. Tazi, S.; Molina, J. J.; Rotenberg, B.; Turq, P.; Vuilleumier, R.; Salanne, M. A transferable ab initio based force field for aqueous ions J. Chem. Phys. 2012, 136, 114507. 18. Fyta, M.; Netz, R. R. Ionic force field optimization based on single-ion and ion-pair solvation properties: Going beyond standard mixing rules J. Chem. Phys. 2012, 136, 124103. 19. Izvekov, S.; Voth, G. A. Effective force field for liquid hydrogen fluoride from ab initio molecular dynamics simulation using the force-matching method J. Phys. Chem. B 2005, 109, 6573-6586. 20. Manzano, H.; Pellenq, R. J. M.; Ulm, F. J.; Buehler, M. J.; van Duin, A. C. T. Hydration of Calcium Oxide Surface Predicted by Reactive Force Field Molecular Dynamics Langmuir 2012, 28, 4187-4197. 21. Jeon, B.; Sankaranarayanan, S. K. R. S.; van Duin, A. C. T.; Ramanathan, S. Reactive Molecular Dynamics Study of Chloride Ion Interaction with Copper Oxide Surfaces in Aqueous Media ACS Appl. Mater. Inter. 2012, 4, 1225-1232. 22. Gale, J. D.; Raiteri, P.; van Duin, A. C. T. A reactive force field for aqueous-calcium carbonate systems Phys. Chem. Chem. Phys. 2011, 13, 16666-16679. 23. van Duin, A. C. T.; Bryantsev, V. S.; Diallo, M. S.; Goddard, W. A.; Rahaman, O.; Doren, D. J.; Raymand, D.; Hermansson, K. Development and Validation of a ReaxFF Reactive Force Field for Cu Cation/Water Interactions and Copper Metal/Metal Oxide/Metal Hydroxide Condensed Phases J. Phys. Chem. A 2010, 114, 9507-9514. 24. Rahaman, O.; van Duin, A. C. T.; Bryantsev, V. S.; Mueller, J. E.; Solares, S. D.; Goddard, W. A.; Doren, D. J. Development of a ReaxFF Reactive Force Field for Aqueous Chloride and Copper Chloride J. Phys. Chem. A 2010, 114, 3556-3568. 25. Rahaman, O.; van Duin, A. C. T.; Goddard, W. A.; Doren, D. J. Development of a ReaxFF Reactive Force Field for Glycine and Application to Solvent Effect and Tautomerization J. Phys. Chem. B 2011, 115, 249-261. 26. Monti, S.; Corozzi, A.; Fristrup, P.; Joshi, K. L.; Shin, Y. K.; Oelschlaeger, P.; van Duin, A. C. T.; Barone, V. Exploring the conformational and reactive dynamics of biomolecules in 44 ACS Paragon Plus Environment

Page 44 of 49

Page 45 of 49 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


solution using an extended version of the glycine reactive force field Phys. Chem. Chem. Phys. 2013, 15, 15062-15077. 27. Shin, Y. K.; Kwak, H.; Vasenkov, A. V.; Sengupta, D.; van Duin, A. C. T. Development of a ReaxFF Reactive Force Field for Fe/Cr/O/S and Application to Oxidation of Butane over a Pyrite-Covered Cr2O3 Catalyst ACS Catal. 2015, 5, 7226-7236. 28. Raymand, D.; van Duin, A. C. T.; Spangberg, D.; Goddard, W. A.; Hermansson, K. Water adsorption on stepped ZnO surfaces from MD simulation Surf. Sci. 2010, 604, 741-752. 29. Abolfath, R. M.; van Duin, A. C. T.; Brabec, T. Reactive Molecular Dynamics Study on the First Steps of DNA Damage by Free Hydroxyl Radicals J. Phys. Chem. A 2011, 115, 1104511049. 30. Bochevarov, A. D.; Harder, E.; Hughes, T. F.; Greenwood, J. R.; Braden, D. A.; Philipp, D. M.; Rinaldo, D.; Halls, M. D.; Zhang, J.; Friesner, R. A. Jaguar: A high-performance quantum chemistry software program with strengths in life and materials sciences Int. J. Quantum Chem. 2013, 113, 2110-2142. 31. BAND2012. SCM, Theoretical Chemistry, Vrije Universiteit: Amsterdam, The Netherlands, http://www.scm.com. 32. Guerra, C. F.; Snijders, J. G.; te Velde, G.; Baerends, E. J. Towards an order-N DFT method Theor. Chem. Acc. 1998, 99, 391-403. 33. te Velde, G.; Bickelhaupt, F. M.; Baerends, E. J.; Guerra, C. F.; Van Gisbergen, S. J. A.; Snijders, J. G.; Ziegler, T. Chemistry with ADF J. Comput. Chem. 2001, 22, 931-967. 34. Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple Phys. Rev. Lett. 1996, 77, 3865-3868. 35. 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. 36. Van Lenthe, E.; Baerends, E. J. Optimized slater-type basis sets for the elements 1-118 J. Comput. Chem. 2003, 24, 1142-1156. 37. Chong, D. P. Augmenting basis set for time-dependent density functional theory calculation of excitation energies: Slater-type orbitals for hydrogen to krypton Mol. Phys. 2005, 103, 749761. 38. Vanlenthe, E.; Baerends, E. J.; Snijders, J. G. Relativistic Regular 2-Component Hamiltonians J. Chem. Phys. 1993, 99, 4597-4610. 39. Vanlenthe, E.; Baerends, E. J.; Snijders, J. G. Relativistic Total-Energy Using Regular Approximations J. Chem. Phys. 1994, 101, 9783-9792. 40. Boerrigter, P. M.; Velde, G. T.; Baerends, E. J. 3-Dimensional Numerical-Integration for Electronic-Structure Calculations Int. J. Quantum Chem. 1988, 33, 87-113. 41. Velde, G. T.; Baerends, E. J. Numerical-Integration for Polyatomic Systems J. Comput. Phys. 1992, 99, 84-98. 42. Swart, M.; Bickelhaupt, F. M. QUILD: QUantum-regions interconnected by local descriptions J. Comput. Chem. 2008, 29, 724-734. 43. Ghosh, M. K.; Re, S.; Feig, M.; Sugita, Y.; Choi, C. H. Interionic Hydration Structures of NaCl in Aqueous Solution: A Combined Study of Quantum Mechanical Cluster Calculations and QM/EFP-MD Simulations J. Phys. Chem. B 2013, 117, 289-295. 44. Velde, G. T.; Baerends, E. J. Precise Density-Functional Method for Periodic Structures Phys. Rev. B 1991, 44, 7888-7903.



45. Wiesenekker, G.; Baerends, E. J. Quadratic Integration over the 3-Dimensional BrillouinZone J. Phys.: Condens. Matter 1991, 3, 6721-6742. 46. Kresse, G.; Furthmuller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set Phys. Rev. B 1996, 54, 11169-11186. 47. Blochl, P. E. Projector Augmented-Wave Method Phys. Rev. B 1994, 50, 17953-17979. 48. Gonzalez, B. S.; Hernandez-Rojas, J.; Wales, D. J. Global minima and energetics of Li+(H2O), and Ca2+(H2O)(n) clusters for n

Development of the ReaxFF Methodology for ElectrolyteâWater

Recommend Documents

Development of the ReaxFF Methodology for ElectrolyteâWater