Parameterization of Reactive Force Field - ACS Publications

Feb 8, 2013 - Cherry L. Emerson Center for Scientific Computation and. ‡ ... recently, Hill and coworkers19−21 and other groups22−24 have report...
0 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/JPCA

Parameterization of Reactive Force Field: Dynamics of the [Nb6O19Hx](8−x)− Lindqvist Polyoxoanion in Bulk Water

Alexey L. Kaledin,† Adri C. T. van Duin,§ Craig L. Hill,‡ and Djamaladdin G. Musaev*,† †

Cherry L. Emerson Center for Scientific Computation and ‡Department of Chemistry, Emory University, Atlanta, Georgia 30322, United States § Department of Mechanical and Nuclear Engineering, Penn State University, University Park, Pennsylvania 16802, United States S Supporting Information *

ABSTRACT: We present results on parameterization of reactive force field [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] for investigating the properties of the [Nb6O19Hx](8−x)− Lindqvist polyoxoanion, x = 0−8, in water. Force-field parameters were fitted to an extensive data set consisting of structures and energetics obtained at the Perdew−Burke−Ernzerhof density functional level of theory. These parameters can reasonably describe pure water structure as well as water with an excess of H+ and OH− ions. Molecular dynamics simulations were performed on [Nb6O19Hx](8−x)−, x = 0−8, submerged in bulk water at 298 K. Analysis of the MD trajectories showed facile H atom transfer between the protonated polyoxoanion core and bulk water. The number of oxygen sites labeled with an H atom was found to vary depending on the pH of the solution. Detailed analysis shows that the total number of protons at bridging (terminal), η-O (μ2-O), sites ranges from 3(1) at pH 7, to 2(0) at pH 11, to 1(0) at pH 15. These findings closely reflect available experimental measurements.

I. INTRODUCTION Polyoxometalates (POMs) are a large and rapidly growing class of inorganic cluster compounds that are made up of d0 metal ions (such as W(VI), Mo(VI), V(V), Nb(V), and Ta(V)) bridged by oxide anions.1−4 They exist in several structural forms, among which the most common ones are the Lindqvist, Keggin, and Wells−Dawson structures. Lindqvist polyoxoanions are iso-polyoxometalates and contain a single type of metal cation bridged by oxide anions with a generic formula of MmOvp−. Keggin and Wells−Dawson structures contain one or more “heteroatoms”, X, in addition to the metal cations and oxide anions and are referred to as “heteropolyanions” (with generic formula of XxMmOvp−). POMs have a number of attributes that render them attractive for applications in catalysis, materials science, and other technologies.5−18 For example, most POMs are readily prepared in quantity by selfassembly in water from nontoxic and readily available precursors. In addition, almost limitless combinations of elements can be synthetically incorporated into POMs, making it possible to control (fine-tune) virtually all of the molecular properties that impact their potential utility as catalysts, medicines and other functional materials. These synthetically alterable properties of POMs include the redox potentials, acidities of the acid forms (H+ salts), polarities and solubilities (by varying the counter cations) as well as their size, shape, and charge. At the appropriate pH values, POMs are thermodynamically stable in water, and thus these species are uniquely compatible with the environmentally optimal oxidant (O2 or © XXXX American Chemical Society

air) and solvent (H2O) for any oxidation process. For example, recently, Hill and coworkers19−21 and other groups22−24 have reported several water oxidation catalysts (WOCs) based on POMs. The collective POM properties listed above and the potential importance of POMs in various areas of science and technology have stimulated a great deal of recent research directed at understanding the interplay of factors that dictate the structural, electronic, and dynamic properties of these compounds. For example, extensive experimental studies show that an important factor controlling the Lindqvist ion stability is the pH of the aqueous solution. As pointed out by Hegetschweiler and coworkers,25 the [HM6O19]7− complexes are dominant species at pH ∼13, whereas at the lower pH values further protonation takes place, and in the pH range ∼10−12, [H2M6O19]6− is expected to be dominant. At 8 ≤ pH ≤ 19, dilute solutions contain mainly the triply protonated [H3M6O19]5− species. Below pH 8 these complexes rearrange to different species. Recently, Casey and coworkers26 demonstrated that oxygenexchange between POMs and bulk water could be qualitatively different for the Lindqvist ions [Nb6O19Hx](8−x)− and [Ta6O19Hx](8−x)−. Furthermore, rates of this process at both Special Issue: Joel M. Bowman Festschrift Received: December 6, 2012 Revised: February 4, 2013

A

dx.doi.org/10.1021/jp312033p | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

[Nb6O19Hx](8−x)− for x = 8 both in gas phase and in bulk water. These calculations were done with VASP4.6.32,33 Our choice of the PBE functional was guided by proper description of metal-oxide systems and their interaction with water, as reported in the literature.34−36 The reasons for using the neutral system [Nb6O19H8] for assembling training sets for parameter fitting are dictated by (i) the much more efficient convergence of the electronic relaxation part of VASP for neutral systems than for negatively charged ions and (ii) the presence of counter-cations in real solutions (modeled by protons). In the VASP calculations, we used periodic cubic cells and chose the box size for POM in water so that at least two solvation shells of water would fit between the POM (∼8.5 Å across) and each of the six cell walls, which results in a 15 Å box. Making more room would improve the quality of the solvent bath but would also make the calculations prohibitively time-consuming. A much smaller box for pure water was needed because the smallest water cluster, (H2O)2, is about 3.5 Å across, and thus an 8 Å box was chosen. To test the actual effect of image-cells, we did a single-point dimer calculation (nonperiodic) (in these calculations we used a 6-31+G(d,p) basis for O and H and LanL2DZ ECP and basis for Nb) of two charge-neutral POMs separated to 15 Å, that is, two neighboring cells filled with dielectric medium modeled by PCM,37 as implemented in Gaussian 09.38 The resulting interaction energy between the monomers was found to be ∼0.6 kcal/mol (or equivalently a temperature perturbation of ∼0.6 K for the given cell size), a value small enough to be negligible. The energy cutoff of 300 eV (∼22 Ry) and the (1 × 1 × 1) Γ-point sampling scheme were applied in all VASP calculations. These settings resulted in ∼39 600 plane waves for the 15 Å cell. II.2. Reactive Force Field (ReaxFF). The potential energy function of ReaxFF consists of bonding (EB) and nonbonding (ENB) terms written as

the terminal (η=O) and bridging [μ2-O] sites for both the M = Nb6 and Ta6 complexes increase by increasing the number of protons in the system (i.e., by reducing the pH of the solution). (See Scheme 1.) Scheme 1. Schematic Presentation of the [Nb6O19]8− Polyoxoaniona

a

There are 6 terminal (η-O), 12 bridging (μ2-O), and 1 central (μ6-O) oxygen.

Despite the active investigations of the Lindqvist polyoxoanions, including [Nb6O19Hx](8−x)− and [Ta6O19Hx](8−x)−, as well as Keggin and Wells−Dawson POMs, the factors controlling the stability and rearrangement pathways of these molecules in aqueous solutions still require a comprehensive investigation. The use of computation to solve these problems looks very promising. Unfortunately, application of the timeconsuming ab initio or density functional theory (DFT) approaches to understand the dynamical properties of such large systems in water is not practical. Therefore, development and validation of optimal (reasonably accurate and less timeconsuming) computational/theoretical approaches is highly desirable. The subject of the current paper is two-fold and includes: (1) development and validation of a reactive force field (ReaxFF) 27,28 for the Lindqvist polyoxoanion [Nb6O19Hx](8−x)− and (2) investigation of stability of [Nb6O19Hx](8−x)− in an explicit water solution using the density functional theory and newly developed parameters for ReaxFF. This article is a first step in the literature toward developing a comprehensive reactive force field capable of describing equilibrium solvation dynamics of the Lindqvist polyoxoanions.

E B = E bond + Eangle + Etorsion + Eover + Eunder + Econj + Epen + E H ‐ bond

E NB = Evdw + Ecoul

(1a) (1b)

The first three terms in eq 1a are energies of the bond between two atoms, valence angle between three atoms, and torsion angle between four atoms, respectively. The following four terms are over/under coordination energies, conjugation energy, and a penalty term for three-atom angle energy in cases of two double bonds sharing an atom (similar to the middle C in H2CCCH2). The last term is the O···H···O hydrogen bond interaction energy. The two terms in eq 1b are energies associated with van der Waals, and Coulomb interactions. To account for electron cloud repulsion at short interatomic distances, a shielded Coulomb potential is used qiqj Ecoul,ij = 3 (rij + 1/γij3)1/3 (2)

II. COMPUTATIONAL PROCEDURES We adopted a three-fold approach. At first, we use DFT to explore the potential energy surface of [Nb6O19Hx](8−x)− in the gas phase for x = 8. In the second step, we place the neutral molecules in bulk water and equilibrate at 298 K. In the final step, we use the equilibrated DFT data to perform parameter optimization for ReaxFF in the gas phase and in water. Following the parameter optimization, we perform equilibrium molecular dynamics simulations and analyze the results. Below, we describe the methods/approaches used. II.1. Details of DFT Calculations. The generalized gradient approximation Perdew−Burke−Ernzerhof (GGA-PBE) density functional29 with plane-wave basis and PAW pseudopotentials 30,31 was used for DFT dynamics simulations of

where γij is a parameter. The atomic charges qi are calculated on-the-fly using the electronegativity equalization method (EEM).39 The detailed description of all of the terms in eqs 1a and 1b has been given before.28 Here we briefly mention that ReaxFF does not use atom types and associated B

dx.doi.org/10.1021/jp312033p | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

Scheme 2. Presentation of the Three-Tier Parameter Optimization Process Useda

a

First, a single set of parameters describing pure water, protonated water, and deprotonated bulk water is obtained, ReaxFFW−H−OH. Second, the parameters describing a single molecule of [Nb6O19H8] (ReaxFFPOMGAS) in the gas phase are optimized. Then, the parameters ReaxFFPOMGAS describing a single molecule of [Nb6O19H8] are re-optimized, incorporating the ReaxFFW−H−OH parameter sets to include solvent effects yielding final parameters (ReaxFFPOMWATER).

water was defined by a three-center term, referred to as O···H···O. The remaining interactions, that is, van der Waals and Coulombic, were defined by atomic parameters of Nb, O, and H. All of the developed and above-defined parameters are provided in the SI. II.3. Training Set Generation and Parameter Fitting. DFT-optimized structures provide starting points for all equilibrium molecular dynamics simulations. A Nosé−Hoover thermostat42,43 with the mass parameter of 0.8 to 1.0 and 1 fs time step was used for generating canonical ensembles. The systems were typically pre-equilibrated for ∼2 ps at 298 K before dynamics were recorded. We propagated trajectories up to 10 ps to generate training sets for force-field parameter fitting. A three-tier approach for parameter optimization was used. First, the parameters responsible for describing all O−H interactions (bonded and nonbonded) were optimized. For this purpose, we generated a training set consisting of three parts: (a) 1 g/cm3 bulk water, modeled by 17 H2O molecules in an 8 Å periodic cube; (b) proton solvated in bulk water, modeled by an H+(H2O)16 system in an 8 Å periodic cube, and (c) hydroxide solvated in bulk water, that is, OH−(H2O)16 in an 8 Å periodic cube. Each of these systems was initially equilibrated at 298 K before a 10 ps molecular dynamics simulation was performed. This resulted in the force field for water/H+/OH−, called ReaxFFW−H−OH. Second, the parameters describing all Nb−O/H interactions were preoptimized based on a 10 ps 298K simulation of [Nb6O19H8] in the gas phase. The Nb6O19H8 structure we used for these simulations was chosen from a preliminary search for the most energetically stable distribution of eight protons at the 18 available oxygen sites. (See the SI for details.) Finally, these parameters were corrected for solvent effects by performing a fit of ReaxFF parameters for [Nb6O19H8] in bulk water, that is, a 15 Å cubic cell containing [Nb6O19H8] and 97 H2O molecules roughly corresponding to a 1 g/cm3 bulk water density. Scheme 2 demonstrates the overall process. In the last fitting procedure, all relevant Nb−O/H parameters were varied, whereas the O, H, and O−H bond and angle parameters were frozen. Although the freezing of the O−H parameters removes variational degrees of freedom and thus lessens the flexibility of

connectivities. These characteristics are determined by the nature of the elements and the geometry of the system. As a result, the bond energy, Ebond, becomes a function of the distance between any two bonded atoms through the introduction of bond order (cf. eqs 2a and 3). Because bond order can change in the course of a reaction or a molecular dynamics simulation, existing bonds can break, and new ones can form. This is based on the concepts put forward by Tersoff40 and later Brenner.41 For example, when a typical covalent bond between atoms “i” and “j” changes from triple to single, the bond order is calculated directly from interatomic distance rij, as follows BOij = exp[pbo,1 (rij/rσ ) pbo,2 ] + exp[pbo,3 (rij/rπ ) pbo,4 ] + exp[pbo,5 (rij/rππ ) pbo,6 ]

(2a)

where pbo,n and rσ, rπ, and rππ are parameters describing σ- and two π-bonds. The energy of the bond is expressed using a Morse-like function p

E bond, ij = −f (BOij )De exp[pbe,1 (1 − BOij be,2 )]

(3)

with two extra parameters pbe,n. The factor f(BOij) appearing at the front is used to modulate the bond order in cases of overcoordination. Its exact functional form can be found in previous publications.27,28 Similarly, Eangle and Etorsion are functions of three and four atom angles, respectively, with added dependence on bond order. For a more detailed description of the ReaxFF method, see Chenoweth et al.27 and van Duin et al.28 In the systems considered in the present work, we defined the following classes of interactions. The two-atom bonded interactions are Nb−O, Nb−H, and O−H. We did not explicitly define Nb−Nb, O−O, and H−H bonds because Nb2, O2, and H2 are extremely unlikely to be present as molecules in bulk water under normal conditions. The three-atom angular interactions considered in this work are H−O−H, H−O−Nb, O−Nb−O, and Nb−O−Nb. Four-atom torsional terms were not considered in the present model because they are not expected to be essential for the description of either the Nb6O19 frame or the bulk water. The energy of hydrogen bonding in C

dx.doi.org/10.1021/jp312033p | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

the ReaxFFPOMWATER fit, it retains a correct description of bulk water, that is, water molecules far from the solute, which makes the force field ReaxFFPOMWATER easily transferable to water solvation of other systems. The parameter fitting procedure is carried out by minimizing the energy function defined as F({C}n) = (1/M )

∑

b. Parameterization of the Nb−O/H Interactions for Nb6O19H8 in Bulk Water. As a preliminary step, at first, we optimize the reactive force-field parameters of the studied Lindqvist ion in the gas phase. Exploratory calculations on a charge-neutral Nb6O19H8 system showed that the lowest energy configuration contains protons at three η-O and five μ2-O sites. Therefore, this structure was used in all subsequent DFT molecular dynamics simulations. Following a short equilibration run, we recorded a 10 ps trajectory in gas phase at 298 K and used the 10000 geometries and corresponding energies to obtain Nb−O/H parameters in the gas phase. In Table 1, we compare key geometrical values obtained at the DFT and ReaxFF levels for [Nb6O19]8− in the gas phase. As

[V ReaxFF(R i) − V DFT(R i)]2 *Wi

i = 1, M

(4) ReaxFF

DFT

where V and V are the force-field and DFT potential energy values at geometry Ri. The summation runs over a set of geometries determined by the equilibrium molecular dynamics. The weights Wi can be used to place greater importance on the lower (or higher) energy configurations, although in practice we have used Wi = 1 throughout. The set of force-field parameters {C}n is varied until an acceptable threshold is reached. These sets are provided in the SI.

Table 1. Comparison of the Optimized Bond Lengths (in angstroms) and Bond Angles (in degrees) of the [Nb6O19]8− Lindqvist Polyoxoanion Computed at the DFT (PBE) and ReaxFF Levels of Theory

III. RESULTS AND DISCUSSION a. Parameterization of the O−H Interactions Using Bulk Water and H+/OH− in Water. The parameters obtained for pure water and water with an H+ and OH− ions were first validated against the available DFT data. We compared (see Figure 1) radial distribution functions for O−H, O−O, and H−

bonds and angles

DFT

ReaxFF

Nb-(μ2-O) Nb-(η-O) Nb-(μ6-O) Nb−Nb (η-O)-Nb-(μ2-O) Nb-(μ2-O)-Nb

2.01−2.05 1.83−1.86 2.43−2.53 3.47−3.54 104 120

2.00−2.07 1.83 2.06−2.88 3.41−3.57 103 113−121

in the case of bulk water, the overall structure appears to be well-reproduced by the force field. Most importantly, the distinction between the double Nb=(η-O) and single Nb-(μ2O) bonds is captured remarkably well. Accurate description of these two types of Nb−O bonds is important for future studies of oxygen exchange in bulk water. A point of concern arises for the rather asymmetric position of the central μ6-O produced by the fit. This artificial Jahn−Teller-like distortion in ReaxFF can be corrected by further refinement of the parameters using a data set containing symmetric geometries, for example, with H atoms at the six terminal (η-O) sites and two bridging (μ2-O) sites. However, the MD simulations presented below show that this asymmetry does not affect the stability of the Lindqvist ion in water. The parameters describing gas phase [Nb6O19]8− were subsequently refined using 8000 geometries obtained from a corresponding MD simulation of [Nb6O19H8] in bulk water at 298 K. In this final step of the fitting procedure, as was mentioned above, the O/H parameters were not varied. Figure 2 shows the error of the fit as compared with the reference DFT data. One can see that despite a significant deviation of ReaxFF from the DFT data there is a close correlation between the two. In other words, the force field is able to uniformly and effectively describe the low-energy as well as the high-energy structures, making it applicable to MD simulations up to ∼320 K. Higher temperature simulations could be problematic because the fit tends to underestimate the high-energy structures but does much better near the average. c. Stability of the Lindqvist Polyoxoanion in Bulk Water. Analysis of the DFT molecular dynamics data shows facile H transfer between the [Nb6O19] core and bulk water. Both the (η-O) and (μ2-O) sites are involved in H exchange with water. No significant structural changes to the Lindqvist core occur during the 10 ps simulations; that is, all Nb-[μ2O(H)] and Nb-[η-O(H)] bonds remain intact, whereas the

Figure 1. Radial distribution functions of pure water, singly protonated water, and singly deprotonated water obtained with newly developed ReaxFFW−H−OH parameters.

H pairs using a 10 ps MD simulation at 298 K. As seen in Figure 1, the main features presented in the DFT spectra are captured in the corresponding ReaxFF studies. Nevertheless, there are notable differences in O−H and H−H pair densities. The O−H pair density, appearing at 1 Å, is slightly broader at the ReaxFF level. However, the second O−H peak (spread over 1.5 to 2 Å in the DFT spectra) is sharper for ReaxFF. This can be explained by the fit’s apparent overestimation of the hydrogen bond and underestimation of the O−H bond in H2O. Similarly, the H−H pair density shows a less sharply defined intramolecular H−H distance and its extended overlap with the secondary, broad, H−H peak, which arises from intermolecular H2O···HOH dimer structures. This points to a weaker HOH bending force constant in the fit. D

dx.doi.org/10.1021/jp312033p | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

Figure 4. Radial distribution functions of [Nb6O19]8− in bulk water at 298 K. The periodic cell used in the simulation was charge neutral.

Figure 2. Statistical accuracy of the ReaxFF in comparison with the DFT data for the Nb6O19H8 system in a 15 Å water box. Deviation of the black dots from the diagonal red line constitutes the error. The potential energy values on the vertical and horizontal axes are relative to the thermal average of the DFT trajectory at 298 K. The rms of the fit is 9.4 kcal/mol. The largest error is 31.7 kcal/mol.

corresponding to the 24 Nb-(μ2-O) bridges. In the ReaxFF spectrum, the terminal/bridging Nb−O bond distinction appears blurred, reflecting the force field’s difficulty in describing the single Nb-[μ2-O(H)] bonds. This is obviously due to the effect of the solvent on the gas phase Nb−O parameters: in water, the distinction between the terminal and bridging oxygens tends to wash out. The long-range tail of the density, however, is reproduced very well. One can draw similar conclusions for the Nb−H pair density, which for DFT exhibits two broad but distinct peaks on 2.5 to 3 Å for Nb−H@(μ2-O) and 3 to 3.5 Å for Nb−H@(η-O) (where H atom is bound to (μ2-O) and (η-O)). ReaxFF, as in the case of Nb−O pair density, tends to merge the two peaks by partially washing out the distinction between the bridging and terminal oxygen centers. The main problem resulting from this is overestimation of the (η-O)-H bond and underestimation of the (μ2-O)-H bond, as shown in further analysis below. To test whether the ReaxFF correctly predicts the equilibrium of H transfer between bulk water and the Lindqvist ion, we integrated the normalized (ηO)-H and (μ2-O)-H bond pair densities, weighted by the spherical volume element r2dr. We further scaled each curve by the corresponding site multiplicity factor, that is, 6 for (η-O)

position of the (μ6-O) center fluctuates moderately inside the core. Repeating the dynamics using the ReaxFF POM WATER optimized parameters over the same time window as the DFT simulation, that is, 10 ps, yields very similar behavior. Figure 3 shows three snapshots along the trajectory: at the beginning, at 5 ps, and at 10 ps. The main feature is the continuous exchange of H+ particles between the oxygen sites of the POM and the surrounding water. Figure 4 shows the radial distribution functions of Nb−Nb, Nb−O, and Nb−H pairs. The agreement between ReaxFF and DFT is reasonable for all three pairs. The agreement is closest for the Nb−Nb and Nb−O pairs, understandably, because ReaxFFPOMWATER accurately reproduces all of the relevant bond distances (cf. Table 1) and the heavier (Nb, O) atoms fluctuate little during the dynamics. For the Nb−O pair, the DFT density shows a small, sharp peak at 1.8 Å, which corresponds to the six Nb-(η-O) bonds, and a tall and broad peak at ∼2 Å,

Figure 3. Snapshots along the ReaxFF trajectory of [Nb6O19]8− in bulk water at 298 K. The dashed lines are used to mark hydrogen bonds and shared protons. The 15 Å cubic periodic cell used in the simulation was charge-neutral. E

dx.doi.org/10.1021/jp312033p | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

and 12 for (μ2-O). The resulting curves are plotted in Figure 5. The window of 0.95 < r < 1.02 Å should provide a clear

Table 2. Average Number of H Atoms at the Oxygen Sites of the Lindqvist Ion in the [Nb6O19H8](H+)97‑n(OH−)97 {i.e. [Nb6O19H8](H2O)97} System at 298 K, Obtained with ReaxFF Simulationsa n

cell charge

pH

0 2 4 6 8

0 −2 −4 −6 −8

7 9 11 13 15

@ (η-O) 1.85 0.13 0.01 0.77 0.25

(0.92) (0.06) (0.01) (0.38) (0.13)

@ (μ2-O) 0.86 2.47 0.13 1.03 0.35

(1.72) (4.94) (0.26) (2.06) (0.70)

a

Sphere of radius 1.025 Å was used to define an O−H bond. The values in parentheses are scaled by 0.5 for (η-O) and by 2 for (μ2-O) (see the text for more details).

Figure 5. Average number of H atoms found within a given distance of the (μ2-O) and (η-O) centers.

indication of the formation of a typical O−H bond. One can see from this Figure that ReaxFF overestimates the proton binding at terminal (η-O) sites by a factor of 2 while underestimating the binding at the bridging (μ2-O) sites, also by a factor of about two. In the former case, (η-O) is undercoordinated and prefers to bind a proton or to closely share it with a water molecule, as in the classical HO−···H+···OH− structural motif. The discrepancy in the latter case is attributable to the fact that (μ2-O) is similar to the oxygen atom in H2O and prefers to share a proton with another water molecule rather than fully bind it. Because all O···H interactions, regardless of whether the O atom is bound to Nb or is in bulk water, are treated with the parameters optimized for water (ReaxFFW−H−OH), this behavior is not unexpected, and further adjustments in the Nb−O parameters should be able to improve the description of the (μ2-O) and (ηO) sites. In the applications below, we will use a correction factor of two for the average number of hydrogens at the bridging oxygen sites and 0.5 for the terminal oxygen sites. d. Effect of the Environment on Proton Affinity. Following the calibrations reported in the previous sections, we investigated the effect of pH on the hydrogen exchange dynamics. The pH of the cell is related to the number of “unpaired” OH− or H+ ions, which is essentially the net charge of the cell. Thus, pH 7 is defined as the charge neutral cell, such as [Nb6O19H8](H2O)97. Removing nH atoms and setting the charge of the cell to −n leads to pH 7 + n. Short-time dynamics (10 ps simulation time) produced closely converged averages, and these results are summarized in Table 2 and Figure 6. Analysis of the pH effect on protonation was found to be greatly facilitated by additional fitting of the scattered MD data to a function that decays to zero at high pH values, that is, lack of free protons, and that has a ceiling at low pH values (large excess of protons) with the value of the corresponding number of maximum available oxygen sites. An obvious choice is a twoparameter hyperbolic tan, for example, G(1-tanh[A(pH-B)])/2, where G is the site multiplicity. The graph shows (see Figure 6) that the total number of protons varies from 4 at pH 7 to 1 at

Figure 6. The average number of H atoms on [Nb6O19] in liquid water obtained with ReaxFF MD simulation at 298 K, taken from the scaled data in Table 2. The function used to fit the data is G(1tanh[A(pH-B)])/2, where G is a site multiplicity factor (6 for η-O, 12 for μ2-O). A is 0.7/0.067, and B is 5.78/−1.52 for η-O/μ2-O, respectively.

pH 15. We find negligible proton density at the terminal sites for pH >8, whereas the density at the bridging sites is measurable up to at least pH 15. These findings closely reflect the descriptions of the [Nb6O19Hx](8−x)− Lindqvist ions made by Hegetschweiler and coworkers25 and also agree qualitatively with experiments of Balogh et al., which indicate that at 12 < pH < 14.5 the oxygen exchange with bulk water is significantly faster at the bridging sites.26 This suggests that at the high values of pH protons label predominantly bridging oxygens, thereby weakening the Nb-(μ2-O) bond and lowering the activation barrier for exchange reaction.

■

CONCLUSIONS The first steps toward developing a full, reactive force field for describing solvation of the [Nb6O19]8− Lindqvist polyoxoanion in water have been completed. An optimized set of ReaxFF parameters is built on the main principle that pure liquid water as well as protonated/deprotonated liquid water should be described with reasonable accuracy, as compared with the available density functional data. The parameters describing the polyoxoanion as a single molecule are then systematically readjusted in the presence of liquid water to include the solvent effects. After validation against reference density functional data, the final set of parameters is used to perform a set of exploratory MD simulations on the effect of water pH on protonation equilibrium of the Lindqvist ion at 298 K. These preliminary findings reveal that in basic solutions (pH >7) protons prefer to F

dx.doi.org/10.1021/jp312033p | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

(9) Weinstock, I. A.; Barbuzzi, E. M. G.; Wemple, M. W.; Cowan, J. J.; Reiner, R. S.; Sonnen, D. M.; Heintz, R. A.; Bond, J. S.; Hill, C. L. Equilibrating Metal-Oxide Cluster Ensembles for Oxidation Reactions Using Oxygen in Water. Nature 2001, 414, 191−195. (10) Neumann, R.; Dahan, M. A Ruthenium-Substituted Polyoxometalate as an Inorganic Dioxygenase for Activation of Molecular Oxygen. Nature 1997, 388, 353−355. (11) Nishiyama, Y.; Nakagawa, Y.; Mizuno, N. High Turnover Numbers for the Catalytic Selective Epoxidation of Alkenes With 1 Atm of Molecular Oxygen. Angew. Chem., Int. Ed. 2001, 40, 3639− 3641. (12) Kozhevnikov, I. V. Catalysis by Polyoxometalates; Wiley: Chichester, England, 2002. (13) Neumann, R. Polyoxometalate Complexes in Organic Oxidation Chemistry. Prog. Inorg. Chem. 1998, 47, 317−370. (14) Santoni, M.-P.; Pal, A. K.; Hanan, G. S.; Proust, A.; Hasenknopf, B. Discrete Covalent Organic-Inorganic Hybrids: Terpyridine Functionalized Polyoxometalates Obtained by a Modular Strategy and Their Metal Complexation. Inorg. Chem. 2011, 50, 6737−6745. (15) Yamase, T.; Prokop, P. V. Photochemical Studies of Alkylammonium Molybdates. Part 11. Photochemical Formation of Tire-Shaped Molybdenum Blues: Topology of a Defect Anion, [Mo142O432H28(H2O)58]12−. Angew. Chem., Int. Ed. 2002, 41, 466−469. (16) (a) Long, D.-L.; Tsunashima, R.; Cronin, L. Polyoxometalates: Building Blocks for Functional Nanoscale Systems. Angew. Chem. Int. Ed 2010, 49, 1736−1758. (b) Long, D.-L.; Burkholder, E.; Cronin, L. Polyoxometalate Clusters, Nanostructures and Materials: From Self Assembly to Designer Materials and Devices. Chem. Soc. Rev. 2007, 36, 105−121. (17) (a) Ibrahim, M.; Lan, Y. H.; Bassil, B. S.; Xiang, Y. X.; Suchopar, A.; Powell, A. K.; Kortz, U. Hexadecacobalt(II)-Containing Polyoxometalate-Based Single-Molecule Magnet. Angew. Chem., Int. Ed. 2011, 50, 4708−4711. (b) Kortz, U.; Müller, A. Introduction: a Special Issue Dedicated to Michael T. Pope. J. Cluster Sci. 2006, 17, 139−141. (18) Umena, Y.; Kawakami, K.; Shen, J.-R.; Kamiya, N. Crystal Structure of Oxygen-Evolving Photosystem II at a Resolution of 1.9 Å. Nature 2011, 473, 55−60. (19) (a) Geletii, Y. V.; Botar, B.; Kögerler, P.; Hillesheim, D. A.; Musaev, D. G.; Hill, C. L. An All-inorganic, Stable, and Highly Active Tetraruthenium Homogeneous Catalyst for Water Oxidation. Angew. Chem., Int. Ed. 2008, 47, 3896−3899. (b) Geletii, Y. V.; Huang, Z.; Hou, Y.; Musaev, D. G.; Lian, T.; Hill, C. L. Homogeneous LightDriven Water Oxidation Catalyzed by a Tetraruthenium Complex with All Inorganic Ligands. J. Am. Chem. Soc. 2009, 131, 7522−7523. (c) Geletii, Y. V.; Besson, C.; Hou, Y.; Yin, Q.; Musaev, D. G.; Quinonero, D.; Cao, R.; Hardcastle, K. I.; Proust, A.; Kögerler, P.; Hill, C. L. Structural, Physicochemical, and Reactivity Properties of an AllInorganic, Highly Active Tetraruthenium Homogeneous Catalyst for Water Oxidation. J. Am. Chem. Soc. 2009, 131, 17360−17370. (d) Kuznetsov, A. E.; Geletii, Y. V.; Hill, C. L.; Morokuma, K.; Musaev, D. G. Dioxygen and Water Activation Processes on Multi-RuSubstituted Polyoxometalates: Comparison with the “Blue-Dimer” Water Oxidation Catalyst. J. Am. Chem. Soc. 2009, 131, 6844−6854. (e) Quiñonero, D.; Kaledin, A. L.; Kuznetsov, A. E.; Geletii, Y. V.; Besson, C.; Hill, C. L.; Musaev, D. G. Computational Studies of the Geometry and Electronic Structure of an All-Inorganic and Homogeneous Tetra-Ru-polyoxotungstate Catalyst for Water Oxidation and Its Four Subsequent One-Electron Oxidized Forms. J. Phys. Chem. A 2010, 114, 535−542. (f) Besson, C.; Huang, Z.; Geletii, Y. V.; Lense, S.; Hardcastle, K. I.; Musaev, D. G.; Lian, T.; Proust, A.; Hill, C. L. Cs9[(γ-PW10O36)2Ru4O5(OH)(H2O)4], a New All-Inorganic, Soluble Catalyst for the Efficient Visible-Light-Driven Oxidation of Water. Chem. Commun. 2010, 46, 2784−2786. (20) Yin, Q.; Tan, J. M.; Besson, C.; Geletii, Y. V.; Musaev, D. G.; Kuznetsov, A. E.; Luo, Z.; Hardcastle, K. I.; Hill, C. L. A Fast Soluble Carbon-Free Molecular Water Oxidation Catalyst Based on Abundant Metals. Science 2010, 328, 342−345. (21) (a) Huang, Z.; Luo, Z.; Geletii, Y. V.; Vickers, J. M.; Yin, Q.; Wu, D.; Hou, Y.; Ding, Y.; Song, J.; Musaev, D. G.; Hill, C. L.; Lian, T.

label bridging rather than terminal oxygen sites, in agreement with experimental inferences of kinetic data as well as density functional values for proton affinity in gas phase. Detailed calculations show that the total number of bridging + terminal protons found on the ion ranges from 3 + 1 at pH 7, to 2 + 0 at pH 11, to 1 + 0 at pH 15. The results of these simulations are quite encouraging for follow-up investigations, such as, (a) refinement of the current set of Nb−O/H parameters to yield a more accurate description of protonation equilibrium, (b) extension of the reference density functional data set and reoptimization of the parameters for simulations of oxygen exchange reactions, and (c) optimization of parameters for the [Ta6O19]8− Lindqvist ion and comparison with the experimental properties for [Nb6O19]8−, as reported.

■

ASSOCIATED CONTENT

* Supporting Information S

Cartesian coordinates of the lowest energy configuration of [Nb6O19H8] optimized in a 15 Å cubic cell at the DFT (PBE) level of theory; Cartesian coordinates of the initial geometry for DFT (PBE) MD simulations of [Nb6O19H8] in bulk water using a 15 Å cubic cell; and parameter set ReaxFFPOMWATER. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The present research is supported by grant DE-FG0203ER15461 from the U.S. Department of Energy. The authors gratefully acknowledge NSF instrumentation MRI-R2 grant (CHE-0958205) and the use of the resources of the Cherry Emerson Center for Scientific Computation.

■

REFERENCES

(1) Pope, M. T. Heteropoly and Isopoly Oxometalates; Springer-Verlag, Berlin, 1983. (2) Day, V. W.; Klemperer, W. G. Metal Oxide Chemistry in Solution: The Early Transition Metal Polyoxoanions. Science 1985, 228, 533−541. (3) Pope, M. T. In Comprehensive Coordination Chemistry; Wilkinson, G., Gillard, R. D., McCleverty, J. A., Eds.; Pergamon Press: New York, 1987; Vol. 3, Chapter 38. (4) Pope, M. T.; Müller, A. Polyoxometalate Chemistry: An Old Field with New Dimensions in Several Disciplines. Angew. Chem., Int. Ed. Engl. 1991, 30, 34−48. (5) Klemperer, W. G.; Wall, C. G. Polyoxoanion Chemistry Moves toward the Future: From Solids and Solutions to Surfaces. Chem. Rev. 1998, 98, 297−306. (6) Kozhevnikov, I. V. Catalysis by Heteropoly Acids and Multicomponent Polyoxometalates in Liquid-Phase Reactions. Chem. Rev. 1998, 98, 171−198. (7) Lv, H. J.; Geletii, Y. V.; Zhao, C. C.; Vickers, J. W.; Zhu, G.; Luo, Z.; Song, J.; Lian, T.; Musaev, D. G.; Hill, C. L. Polyoxometalate Water Oxidation Catalysts and the Production of Green Fuel. Chem. Soc. Rev. 2012, 41, 7572−7589. (8) Hill, C. L. Introduction: Polyoxometalates Multicomponent Molecular Vehicles to probe Fundamental Issues and Practical Problems. Chem. Rev. 1998, 1, 1−2 and references therein. G

dx.doi.org/10.1021/jp312033p | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

(36) Alimohammadi, M.; Fichthorn, K. A. A Force Field for the Interaction of Water with TiO2 Surfaces. J. Phys. Chem. C 2011, 115, 24206−24214. (37) Cammi, R.; Tomasi, J. Remarks on the Use of the Apparent Surface Charges (ASC) Methods in Solvation Problems: Iterative Versus Matrix-Inversion Procedures and the Renormalization of the Apparent Charges. J. Comput. Chem. 1995, 16, 1449−1458. (38) Frisch, M. J., et al. Gaussian 09, revision A.02; Gaussian, Inc.: Wallingford, CT, 2009. (39) Janssens, G. O. A.; Baekelandt, B. G.; Toufar, H.; Mortier, W. J.; Schoonheydt, R. A. Comparison of Cluster and Infinite Crystal Calculations on Zeolites with the Electronegativity Equalization Method (EEM). J. Phys. Chem. 1995, 99, 3251−3258. (40) Tersoff, J. Empirical Interatomic Potential for Silicon with Improved Elastic Properties. Phys. Rev. B 1988, 38, 9902−9905. (41) Brenner, D. Empirical Potential for Hydrocarbons for Use in Simulating the Chemical Vapor Deposition of Diamond Films. Phys. Rev. B 1990, 42, 9458−9471. (42) Nosé, S. A Unified Formulation of the Constant-Temperature Molecular-Dynamics Methods. J. Chem. Phys. 1984, 81, 511−519. (43) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695−1697.

Efficient Light-Driven Carbon-Free Cobalt-Based Molecular Catalyst For Water Oxidation. J. Am. Chem. Soc. 2011, 133, 2068−2071. (b) Zhu, G.; Geletii, Y. V.; Kögerler, P.; Schilder, H.; Song, J.; Lense, S.; Zhao, C.; Hardcastle, K. I.; Musaev, D. G.; Hill, C. L. Water Oxidation Catalyzed by a New Tetracobalt-Substituted Polyoxometalate Complex: [{Co4(μ-OH)(H2O)3}(Si2W19O70)]11−. Dalton Trans. 2012, 41, 2084−2090. (22) (a) Sartorel, A.; Carraro, M.; Scorrano, G.; Zorzi, R. D.; Geremia, S.; McDaniel, N. D.; Bernhard, S.; Bonchio, M. Polyoxometalate Embedding of a Tetraruthenium(IV)-oxo-core by Template-Directed Metalation of [γ-SiW10O36]8−: A Totally Inorganic Oxygen-Evolving Catalyst. J. Am. Chem. Soc. 2008, 130, 5006−5007. (b) Sartorel, A.; Miro, P.; Salvadori, E.; Romain, S.; Carraro, M.; Scorrano, G.; Valentin, M. D.; Llobet, A.; Bo, C.; Bonchio, M. Water Oxidation at a Tetraruthenate Core Stabilized by Polyoxometalate Ligands: Experimental and Computational Evidence To Trace the Competent Intermediates. J. Am. Chem. Soc. 2009, 131, 16051−16053. (c) Natali, M.; Orlandi, M.; Berardi, S.; Campagna, S.; Bonchio, M.; Sartorel, A.; Scandola, F. Photoinduced Water Oxidation by a Tetraruthenium Polyoxometalate Catalyst: Ion-pairing and Primary Processes with Ru(bpy)32+ Photosensitizer. Inorg. Chem. 2012, 51, 7324−7331. (d) Sartorel, A.; Carraro, M.; Toma, F. M.; Prato, M.; Bonchio, M. Shaping the Beating Heart of Artificial Photosynthesis: Oxygenic Metal Oxide Nano-Clusters. Energy & Env. Sci. 2012, 5, 5592−5603. (23) Car, P.-E.; Guttentag, M.; Baldridge, K. K.; Albertoa, R.; Patzke, G. R. Synthesis and Characterization of Open and Sandwich-Type Polyoxometalates Reveals Visible-Light-Driven Water Oxidation via POM-Photosensitizer Complexes. Green Chem. 2012, 14, 1680−1688. (24) Stracke, J. J.; Finke, R. G. Electrocatalytic Water Oxidation Beginning with the Cobalt Polyoxometalate [Co4(H2O)2(PW9O34)2]10−: Identification of Heterogeneous CoOx as the Dominant Catalyst. J. Am. Chem. Soc. 2011, 133, 14872−14875. (25) Hegetschweiler, K.; Finn, R. C.; Rarig, R. S., Jr.; Sander, J.; Steinhauser, S.; Worle, M.; Zubieta, J. 1,3,5-Triamino-1,3,5-TrideoxyCis-Inositol, a Ligand with a Remarkable Versatility for Metal Ions. 10. Surface Complexation of [Nb6O19]8− with NiII: Solvothermal Synthesis and X-ray Structural Characterization of two Novel Heterometallic Ni-Nb-Polyoxometalates. Inorg. Chim. Acta 2002, 337, 39−47. (26) Balogh, E.; Anderson, T. M.; Rustad, J. R.; Nyman, M.; Casey, W. H. Rates of Oxygen-Isotope Exchange between Sites in the [HxTa6O19](8−x)‑(aq) Lindqvist Ion and Aqueous Solutions: Comparisons to [HxNb6O19](8−x)‑(aq). Inorg. Chem. 2007, 46, 7032−7039. (27) Chenoweth, K.; van Duin, A. C. T.; Goddard, W. A. ReaxFF Reactive Force Field for Molecular Dynamics Simulations of Hydrocarbon Oxidation. J. Phys. Chem. A 2008, 112, 1040−1053. (28) 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. (29) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (30) Blöchl, P. E. Projector Augmented-Wave Method. Phys. Rev. B 1994, 50, 17953−17979. (31) Kresse, G.; Joubert, D. From Ultrasoft Pseudopotentials to the Projector Augmented-Wave Method. Phys. Rev. B 1999, 59, 1758− 1775. (32) The Vienna Ab initio Simulation Package (VASP) is a computer program for atomic scale materials modeling. http://www.vasp.at. (33) Kresse, G.; Furthmüller, J. Efficient Iterative Schemes for ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B 1996, 54, 11169−11179. (34) Liu, H.; Tan, K. Adsorption of Water on Single-Walled TiO2 Nanotube: a DFT Investigation. Comput. Theor. Chem. 2012, 991, 98− 101. (35) Zhao, Z.; Li, Z.; Zou, Z. A Theoretical Study of Water Absorption and Decomposition on the Low-Index Stoichiometric Anatase TiO2 Surfaces. J. Phys. Chem. C 2012, 116, 7430−7441. H

dx.doi.org/10.1021/jp312033p | J. Phys. Chem. A XXXX, XXX, XXX−XXX