J. Phys. Chem. B 2009, 113, 4111–4118
4111
A Modified MSEVB Force Field for Protonated Water Clusters† R. Kumar, R. A. Christie, and K. D. Jordan* Department of Chemistry and Center for Molecular and Materials Simulations, UniVersity of Pittsburgh, Pittsburgh, PennsylVania 15260 ReceiVed: July 26, 2008; ReVised Manuscript ReceiVed: September 16, 2008
A multistate empirical valence bond (MSEVB) model for protonated water clusters, which incorporates the TIP4P water model, is presented. This model which is designated MSEVB4P represents a significant improvement over the original model of Voth et al. (J. Phys. Chem. B 1998, 102, 5547) which was based on the TIP3P water model and a smaller improvement over the recently introduced MSEVB3 model (J. Phys. Chem. B 2008, 112, 467) which is based on the SPC/Fw (J. Chem. Phys. 2006, 124, 024503) water model. 1. Introduction Protonated water clusters have attracted considerable interest in recent years.1-20 This interest has been motivated by fundamental questions concerning the relative importance of Eigen21 and Zundel22,23 structures of excess protons in water, by new vibrational spectroscopic measurements, and by growing recognition that the protonated water clusters are important in atmospheric processes, proton exchange membrane fuel cells, as well as in biological proton transfer.24-30 One of the most intriguing experimental observations is the occurrence of the n ) 21 cluster as a magic number in the mass spectrum.17-19,31,32 Vibrational predissociation spectroscopy of protonated water clusters formed in supersonic jet conditions has shed light on the structure of the magic number cluster and has provided evidence for Eigen, Zundel, and intermediate species for the excess proton.17,18,33 Ab initio electronic structure calculations on small protonated water clusters have been carried out by a number of groups,2,4,7,8,15,34,35 and DFT-based molecular dynamics simulations have been carried out for clusters as large as H+(H2O)22.19,20,36 However, such calculations are very demanding for large clusters or condensed phase systems which has led, naturally, to efforts to develop model potentials for describing excess protons in water.8,9,37-39 Because of the high mobility of excess protons in water, realistic models must allow for the proton to move between water monomers, which makes the development of force fields for excess protons in water more challenging than for neutral unionized water itself. One of the most promising approaches for describing proton delocalization is the empirical valence bond (EVB) method. As introduced by Warshel,40 the EVB method allowed for the delocalization of the proton over two states, namely, the “reactant” and the “product”. The two-state EVB concept was utilized by Voth et al. to model the excess proton in water.39 Subsequently,VuilleumierandBorgis,38,41 Vothandco-workers,42-44 and Brancato and Tuckerman45 extended the EVB model to allow for the delocalization of the proton over multiple water molecules. Voth and co-workers have used their multistate empirical valence bond (MSEVB) method to study proton transport in bulk water and in biological systems as well as to characterize excess protons at water/vapor interfaces.29,46,47 † Part of the special section “Aqueous Solutions and Their Interfaces”. * To whom correspondence should be addressed. E-mail:
[email protected].
Figure 1. Structures of the H+(H2O)n, n ) 2-4, clusters and definitions of the bond lengths used in the tables.
In spite of the many successes of the MSEVB approach, the relative energies for small H+(H2O)n clusters obtained using this approach can differ appreciably from those of MP2 calculations.15 The original MSEVB model used the TIP3P48 water potential to describe the water-water interactions. However, it has been found that the TIP3P model does not do a good job at describing the relative energies of different isomers of neutral water clusters.49 Ideally, one would like to couple the MSEVB method with a polarizable water model such as TTM2,50 AMOEBA,51 DPP,52 or CC-Pol.53 Here, we explore an “intermediate” strategy of coupling the MSEVB model with the TIP4P48 two-body water force field, which displaces the minus point charge from the O atom of a monomer to a socalled M site located on the rotational axis 0.15 Å from the O atom, in contrast to the TIP3P, SPC,54 or SPC/E55 models which locate the negative charge at the O atom. The displacement of the minus charge to the M site leads to improved electrostatics, and as a result, the TIP4P model has been found to be more successful at predicting the relative energies of different isomers of (H2O)n clusters than models employing three atom-centered charges.49,56,57 The paper is organized as follows. Section 2 gives a brief description of the MSEVB method and the parametrization of the new model (MSEVB4P). Section 3 compares for several protonated water clusters results obtained using the MSEVB4P
10.1021/jp8066475 CCC: $40.75 2009 American Chemical Society Published on Web 11/12/2008
4112 J. Phys. Chem. B, Vol. 113, No. 13, 2009
Kumar et al.
TABLE 1: Parameters for the MSEVB4P Modela function damped
parameter
value
Vconst ij Vex ij Vex ij Vex ij Vex
DOO P k DOO a
2.88 Å 1.05 30.0 Å-2 2.36 Å 0.0
a Parameters not specified in the table are the same as those in the MSEVB1 model.
TABLE 2: Bond Lengths (Å) for the Optimized Structures of H+(H2O)n, n ) 2-4 method geometrical MP2a MSEVB1 MSEVB3 MSEVB4P parameter
cluster H+(H2O)2
r1 r2 r3 r4 r5 r6 r7 r8 r9
H+(H2O)3 H+(H2O)4 (4a) H+ (H2O)4 (4b)
a
2.39 1.20 2.50 1.03 2.56 1.01 2.38 2.61 1.19
2.40 1.20 2.52 1.06 2.57 1.04 2.38 2.62 1.19
2.39 1.20 2.50 1.07 2.55 1.05 2.37 2.59 1.19
2.40 1.20 2.50 1.05 2.54 1.04 2.38 2.61 1.19
MP2 results from ref 8.
TABLE 3: Formation Energies (kcal/mol) of H+(H2O)n, n ) 2-4, Clusters method a
cluster
MP2
MSEVB1
MSEVB3
MSEVB4P
dimer trimer tetramer (4a) tetramer (4b) rmsdb
-34.40 -57.60 -77.40 -73.90 0.0
-33.83 -57.12 -79.57 -72.33 1.4
-33.05 -57.18 -79.7 -73.90 0.9
-34.79 -57.26 -78.97 -73.08 0.7
a MP2 results from ref 8. energies.
b
rmsd given relative to the MP2
TABLE 4: Formation Energies (kcal/mol) of Low-Energy Isomers of H+(H2O)5 method isomer
MP2a
MSEVB1
MSEVB3
MSEVB4P
I II III IV V rmsdb
-92.60 -91.72 -91.74 -88.47 -88.60 0.0
-90.96 -92.20 -92.16 -84.41 -85.56 2.4
-93.64 -93.96 -93.85 -87.67 -89.41 1.6
-92.47 -91.82 -91.82 -84.73 -87.47 1.7
a
From ref 15. b rmsd given relative to the MP2 energies.
model, the original MSEVB model (hereafter referred to as MSEVB1), and the recently introduced MSEVB3 model of Voth et al.58 which incorporates the SPC/Fw59 flexible three-site water model that has been found to be superior to the TIP3P model for describing bulk water properties. The performance of the various MSEVB procedures is assessed by comparing the relative energies of selected isomers of H+(H2O)n, n ) 5, 6, and 22, clusters obtained with these methods with those from MP2 or RIMP260 calculations. Our conclusions are presented in section 4. 2. Method In the multiconfigurational EVB models, valence bond states are constructed by associating the excess proton with various
Figure 2. Potential energy curves for proton tranfer in H5O2+ for ROO distances of (a) 2.4, (b) 2.6, and (c) 2.8 Å.
water molecules, giving hydronium species (H3O+), and the electronic wave function is represented as a linear combination of the valence bond states. The valence bond configurations are generated by identifying the “pivot hydronium” and allowing the O atoms associated with water molecules within the first three solvation shells (two solvation shells in the case of MSEVB1) of the “pivot” hydronium to form hydronium species.43,44,58 The solvation shells are determined using hydrogenbonding definitions involving the O · · · H distance (with a cutoff of 2.5 Å) and the O-H · · · O angle (with a cutoff of 130°). In addition, in the MSEVB4P and MSEVB3 procedures, individual O atoms are allowed to participate in more than one hydronium species, which is especially important when the O atoms of the water molecules in the cluster are involved in two acceptor hydrogen bonds. For example, for isomer I of H+(H2O)5 (Figure 3), six valence bond states are constructed using the new selection algorithm compared to five using the original one. This
MSEVB Force Field for Protonated Water Clusters
J. Phys. Chem. B, Vol. 113, No. 13, 2009 4113
Figure 3. Structures of five low-energy isomers of H+(H2O)5.
TABLE 5: Formation Energies (kcal/mol) of Low-Energy Isomers of H+(H2O)6 method isomer
MP2
MSEVB1
MSEVB3
MSEVB4P
I II III IV V VI VII VIII IX X XI rmsda
-105.33 -105.22 -105.71 -105.02 -105.80 -105.80 -105.49 -105.24 -104.05 -104.41 -101.79 0.0
-101.51 -101.51 -101.51 -101.34 -103.29 -104.83 -103.12 -104.12 -105.46 -104.27 -97.41 2.7
-106.56 -106.56 -106.51 -105.45 -106.89 -108.83 -107.32 -107.7 -109.18 -107.48 -102.55 2.2
-104.41 -104.59 -104.59 -105.55 -104.48 -105.77 -104.41 -106.26 -106.15 -103.48 -101.86 1.0
a
allows for the interaction between water monomers, and Vi,H3O+ accounts for the interaction of the H3O+ ion with the ith water monomer. In the MSEVB4P model, we retain the functional intra employed in the MSEVB1 model but forms of V Hintra + and V j 3O modify Vij and Vi,H3O+ to employ the charges and Lennard-Jones terms appropriate for the TIP4P water model. The Coulombic interactions of the hydronium ion with the water molecules are scaled by 0.85 (as opposed to 0.76 in the MSEVB1 model) to allow for the change in the dipole moment in the water monomer in going from the TIP3P to the TIP4P model. The motivation for the charge scaling is given in ref 42. The off-diagonal terms in the MSEVB Hamiltonian42,43 are represented as
Hij ) V const · f(ROO) · g(q) + V ijex · f(ROO) · g(q)
rmsd given relative to the MP2 energies.
extension appears to have been originally introduced by Wales et al. in their implementation of the MSEVB method into the OPTIM program.61 The Hamiltonian in the MSEVB method contains diagonal terms (Hii) corresponding to the intramolecular and intermolecular interactions in each valence bond state and off-diagonal terms (Hij) that couple the states. Clearly, valence bond states involving water molecules closer to the excess proton are more important than those involving water molecules far from the excess proton which, in turn, can be used to restrict the number of states used to describe the system. The diagonal terms are of the form42,43 Nwater
Hii ) V Hintra + + 3O
∑ j)1
(2)
where f(ROO) ) (a + P exp(-k(ROO - DOO)2)) × 1 0 0 (1 - tanh[β(ROO - ROO )]) + b exp[-a(ROO - rOO )] (3) 2
(
)
and
ROO - ROH 2
(4)
g(q) ) exp(-γq2)
(5)
q)
Nwater
V jintra +
∑ V ij + ∑ i
