11726
J. Phys. Chem. B 2008, 112, 11726–11733
Al(III) Hydration Revisited. An ab Initio Quantum Mechanical Charge Field Molecular Dynamics Study Thomas S. Hofer, Bernhard R. Randolf, and Bernd M. Rode* Theoretical Chemistry DiVision, Institute of General, Inorganic and Theoretical Chemistry, UniVersity of Innsbruck, Innrain 52a, A-6020 Innsbruck, Austria ReceiVed: March 27, 2008; ReVised Manuscript ReceiVed: May 17, 2008
To assess the novel quantum mechanical charge field (QMCF) molecular dynamics (MD) approach, two simulations of hydrated Al(III) have been carried out, as this system proved to be a well-suited test case for hybrid ab initio/molecular mechanics simulations. Two different population analysis schemes according to Mulliken and Lo¨wdin have been applied to evaluate the atomic charges in the QM region. It is shown that the QMCF MD approach yields a substantially improved description of the system and that, due to the fact that solute-solvent potentials can be renounced, the QMCF MD framework is a more convenient approach to investigate solvated systems compared to conventional ab initio QM/MM MD approaches. Introduction In recent publications the novel quantum mechanical charge field molecular dynamics (QMCF MD) approach has been introduced,1,2 a QM/MM MD scheme utilizing electrostatic embedding and a dynamical fluctuating charge field to describe the interactions between the QM and MM regions. One of the main advantages of this framework is the renounceability of solute-solvent potentials due to a large radius of the QM region. The construction of these potential functions is a timeconsuming and tedious task, and the associated accuracy is limited. Investigations utilizing the QMCF methodology have already shown its superiority over conventional QM/MM schemes. The only required potential functions in this approach are those for solvent species located in the MM region, while full migration to and from the QM region is possible. To further assess the capabilities of this method compared to the conventional QM/MM MD scheme, QMCF MD simulations of hydrated Al(III) have been carried out. This particular ion induces a considerable polarization in its surroundings resulting from its high charge and small ionic radius (0.54 Å,3 q/r ≈ 5.5 e/Å). Due to the low number of electrons, the computational effort is not as demanding as for other strong polarizing agents such as transition-metal ions. Furthermore, the influence of special relativity can be considered negligible in the case of this light atom. Several properties of this system can be reliably deduced after a few picoseconds of simulation time due to the rigidity of the aluminum(III) hydrate, and structural and dynamical data have been well investigated by experiment. Hence, it can be concluded that this system is a well-suited candidate for a comparative methodical study. Aluminum is employed in industry as a basic material for numerous tasks, and hence, exposure to Al(III) is quite common. It is believed to exhibit no toxicity except in very high doses.4 On the other hand, it is assumed that Al(III) ions are associated with biological defects related to Alzheimer’s disease.4 Besides the high polarizing effects, its amphoteric character makes Al(III) a particular case for solution chemistry. Its high Lewis * To whom correspondence should be addressed. E-mail: Bernd.M.Rode@ uibk.ac.at. Phone: +43-512-507-5160. Fax: +43-512-507-2714.
acidity is utilized in various chemical processes such as Friedel-Crafts-type reactions. The Al(III) ion has been recently studied employing a conventional QM/MM MD approach.5,6 It was concluded that the inclusion of the first and second hydration shells into the QM region was mandatory to obtain a reliable description of this hydrate. The simulation data were in good agreement with experimental results as far as available. As a two-shell treatment is inherent in QMCF MD simulations, its computational effort is similar to that of the aforementioned QM/MM MD study. It was, therefore, of special interest whether the QMCF MD framework with its advanced consideration of charges is able to yield improved results compared to the conventional QM/ MM MD scheme and, at the same time, offers a more convenient approach as solute-solvent interaction potentials are not required. Investigations of similar systems have shown that the QMCF scheme offers a straightforward route for the investigation of even more complex systems.1,7 Thus, the treatment of composite ions such as TiO+ 8 and Hg22+ 9 and oxo ions10–12 such as ClO4-, SO42-, and PO43- have greatly benefited from the advantages of this novel approach. Methods Similar to all QM/MM methodologies,13–16 the chemical system is divided into two parts: The first contains the chemically most relevant region and is treated by a suitable quantum mechanical technique. The remaining part of the system is treated by means of classical mechanics. This procedure combines the accuracy of quantum theory with the affordability of classical mechanics. While the theoretical treatment of interactions within the regions is straightforward, the coupling of the subsystems is by no means a trivial task. In conventional QM/MM schemes molecular mechanics potentials are required to account for the interactions between all QM and MM particles, which is often referred to as potential embedding. The construction of these potential functions is a delicate, time-consuming, and sometimes difficult task. However, it was observed in QM/MM MD simulations of hydrated ions that the non-Coulombic interactions of these potentials become negligibly small after a distance of a few angstroms and that the only remaining contribution results from the
10.1021/jp802663h CCC: $40.75 2008 American Chemical Society Published on Web 08/26/2008
Al(III) Hydration Revisited
J. Phys. Chem. B, Vol. 112, No. 37, 2008 11727
Figure 1. Partitioning of the simulation box in the QMCF approach.
Coulombic interactions. As long as the solute (in this case an ion) resides at the center of the QM region, the distance remains large and only Coulombic interactions need to be evaluated. Solvent particles close to the center, e.g., first-shell ligands of an ion, can in many cases be treated in the same way, whereas atoms located close to the QM/MM border do no fulfill this condition; non-Coulombics need to be applied as the interatomic distances are low. As interaction potentials are available due to the MM description of the solvent, this does not pose further difficulties. Thus, solvent particles have to be distinguished depending on their location in the QM region. This leads to a partition of the QM region as displayed in Figure 1. Due to the short range of the non-Coulombic parts of the water-water potential the six first-shell ligands of the Al(III) ion were included in the core zone as well. The interaction between particles in the core region with atoms in the MM zone is accounted for only via Coulombics. The layer region on the other hand contains the entire second hydration shell. Besides Coulombic forces, also non-Coulombic forces are utilized to describe the interaction between particles in this layer zone and the MM region. As the exchange rate of water molecules in the first shell of Al(III) is low, no exchange of first-shell water molecules was expected; thus, the partitioning of this system will be maintained. According to this definition, the forces in the respective subregions are defined as
Fcore ) FQM J J + M
Flayer ) FQM J J +
∑
M
N1+N2
∑ FBJH IJ + ∑ I)1
MM qQM J qI
I)1
rIJ2
∑
MM qQM J qI
I)1
) FMM J
M
I)1
rIJ2
∑ FMM-nC IJ
(2)
I)1
N2
+
HCF ) HHF + V′ M
V′ )
q
∑ riJJ
(4)
J)1
M
+
MM qQM I qJ
rIJ2
(1)
discussed earlier the non-Coulombic forces FMM-nC between QM IJ particles located in the layer region (N2) and MM particles (M) have to be included via potential embedding. The force acting on an MM particle consists of three contributions, namely, the interactions with all other MM atoms (M) calculated on the basis of classical potential functions FMM IJ , the non-Coulombic interactions FIJMM-nc with all layer atoms (N2), and the charge field contributions with all QM atoms (N1 + N2). According to these expressions, the particles obtain proper forces and the conservation of the total linear momentum is ensured, as the sum of forces equals zero. All Coulombic interactions between QM and MM atoms are obtained utilizing partial charges assigned to the QM atoms via population analysis. This assignment is performed in every step of the simulation. It has been noticed from a variety of QM/ MM simulations that the assignment of fixed partial charges could result in inconsistencies near the border region between the QM and MM zones. As the quantum mechanical treatment accounts for polarization and charge transfer effects, in contrast to the majority of MM schemes, discrepancies occur whenever a particle moves from the QM to the MM region and vice versa. For example, the formal charge of +3.0 assigned to Al(III) in connection with fixed, unpolarized charges assigned to the atoms of the first-shell molecules appears to be not realistic, as the formal charge of the ion is depolarized by its chemical surroundings and the charge of the ligand atoms is altered as well. This can be accounted for by performing a population analysis incorporated in every step of the simulation to obtain a more reliable description of the Coulombic interactions. The re-evaluation of the partial charges in every subsequent step accounts for all structural changes of the system and delivers a more flexible representation of the electron distribution than fixed charges even if they are derived by means of a more complex framework17 to include polarization effects in an averaged way. To obtain a more reliable description of the quantum mechanical region, the point charges of all MM particles (according to the respective MM model) are embedded in the Hamiltonian as a perturbation potential. This technique is routinely applied in various methods to polarize the QM region.18–20 Otherwise the treatment of the QM region would occur in an artificial vacuum environment which does not account for the surrounding bulk, and hence, artificial surface effects may result. As the positions of the MM point charges vary along the simulation, a fluctuating field of charges is obtained:
∑ FMM-nC IJ
Furthermore, a special treatment is required whenever solvent molecules leave and enter the QM region. To ensure a continuous change of forces, a smoothing function S(r) is applied between the radii r0 (5.0 Å) and r1 (4.8 Å):
(3)
I)1
I*J
Fcore corresponds to the quantum mechanical force acting on a J particle J in the core zone, Flayer to the forces acting on a particle J J located in the solvation layer, and FMM to the forces acting J on a particle J in the MM region. The forces acting on all atoms in the core region (N1) consist of the QM forces FQM plus the J charge field interaction according to Coulomb’s law. As
MM Fsmooth ) S(r)(Flayer - FMM J J J ) + FJ
(5)
Fsmooth is the force of a particle J in the smoothing region, Flayer J corresponds to the force of the particle as a layer atom, and FMM corresponds to the force resulting if the particle were J already included in the MM region. The smoothing factor S(r) is derived on the basis of the distance between the QM center and the center of mass of the respective molecules:
11728 J. Phys. Chem. B, Vol. 112, No. 37, 2008
S(r) ) 1 S(r) )
for
r e r1
(r02 - r2)2(r02 + 2r2 - 3r12)
S(r) ) 0
(r02 - r12)3 for
Hofer et al.
for
r1 < r e r0
TABLE 1: Range of Partial Charges Obtained in Single-Point Calculations of Various Al(III)-Water Clusters Obtained from Different Population Analysis Schemes
(6)
r > r0
r is the distance of a given solvent molecule’s center of mass from the center of the QM region, r0 the radius of the layer region, and r1 the inner border of the smoothing region. A thickness of 0.2 Å is usually employed as this value was found to be the minimum distance to ensure smooth transitions of solvent molecules. An important consideration for QM/MM studies is the selection of suitable basis sets, the respective quantum mechanical level, and the associated computational effort. To be consistent with the previous QM/MM MD study, the same basis sets were employed.5,6 Efficiency considerations are of particular importance related to the affordable level of accuracy. The expected number of particles included in the QM region was ∼20. The associated average computing time per MD step should be lower than 2.5 min (executing the quantum mechanical computations in parallel), thereby leading to a net computation time of 3.5 months for a 10 ps trajectory. As the parallelization is in general limited to about four to eight processors, the only affordable quantum chemical methods are those of the single-determinant level, i.e., density functional theory (DFT) or ab initio Hartree-Fock. Test calculations on various levels of theory (HF, CC-D, QCISD, MP4-SDQ, MP2, B3PW91, B3P86, and B3LYP) have shown5 that the ab initio HF-SCF method is preferred over all mentioned density functional methods. This finding is not surprising as the Al(III)-water interaction is strongly dominated by electrostatic interactions. The known shortcomings of the present implementations of density functional theory such as self-interactions and the erroneous treatment of kinetic energies could be responsible for the observed discrepancies. Detailed re-evaluation and refitting of the parameters of the respective functionals could result in an improved description of this particular system. Besides the applied level of theory, the choice of the population analysis scheme is most crucial in the case of a QMCF MD study. Different approaches are available in the Turbomole 5.9 package,21,22 namely, population analysis according to Mulliken23,24 and Lo¨wdin25 and the natural bond orbital (NBO) scheme26 derived by Reed, Weinstock, and Weinhold. The fourth scheme derives the partial charges according to the respective occupation numbers.27 Test calculations of random configurations including the Al(III) ion and first plus second hydration shell embedded in a set of point charges taken from the trajectory of a previous two-shell QM/ MM MD study5 were carried out (cf. Table 1). The partial charges assigned to the MM particles employed for embedding were taken in accordance with the BJH-CF2 water model,28,29 which was utilized to describe the MM water molecules throughout the simulation. The charges are -0.65966 and 0.32983 unit charges for oxygen and hydrogen, respectively. The different population analysis schemes show some significant differences (cf. Table 1). NBO and Mulliken show similar partial charges for oxygens and hydrogens, but the charge of the Al(III) ion is depolarized in the Mulliken case. The atomic partial charges of the first-shell ligands are amplified, indicating the influence of polarizaton. This effect is not as pronounced in the case of Lo¨wdin population analysis, and an actual distinction between first- and second-shell ligands is hardly possible. Furthermore, the partial charge of Al(III) of 2.10
atom Al
NBO
Mulliken
Lo¨wdin 2.1
occupation
2.95
2.55
O H
-1.15 to -1.2 0.5 to 0.6
First Shell -1.1 to -1.2 -0.65 to -0.7 0.55 to 0.65 0.35 to 0.40
2.95 -0.75 to -0.8 0.35 to 0.4
O H
-1.0 to -1.05 0.45 to 0.55
Second Shell -0.95 to -1.05 -0.6 to -0.7 0.45 to 0.55 0.30 to 0.35
-0.5 to -0.55 0.25 to 0.30
appears to be rather low. In the case of the charge assignment according to the occupation numbers, a considerable difference between first- and second-shell ligands can be noticed. However, the charge of Al(III) seems to be not affected, and the partial charges of the second-shell ligands appear to be too low compared to those of the BJH-CF2 model. On the basis of these preliminary test calculations, the Mulliken and Lo¨wdin population analyses were chosen to derive partial charges in the QM region. As the Mulliken and NBO charges are similar, an explicit simulation with the latter method appeared renounceable. The population analysis based on occupation number yielded too low charges compared to the MM model and was not considered, therefore, for a separate simulation as well. Furthermore, it is important to notice that this method requires a significant computing time (up to 15 s), whereas the computation is rapidly executed with the other population analysis schemes. All parameters of the simulation protocol were chosen identical to those of the previous QM/MM MD study5 of hydrated Al(III). The density of the system was fixed at 0.997 g/cm3, which is the density of pure water at the simulation temperature. A cubic box with a 24.6 Å side length was employed containing 1 Al(III) ion and 499 water molecules. The simulation was carried out in an NVT ensemble employing the Berendsen thermostat30 (bath relaxation time 0.1 ps) to maintain a temperature of 298.15 K. The reaction field method31 was employed to correct the cutoff of long-range electrostatic interactions. The Newtonian equations of motion were treated by a predictor-corrector algorithm of the Adams-Bashforth family. Since the QM treatment as well as the BJH-CF2 water model28,29 allows explicit hydrogen movements, a time step of 0.2 fs was chosen. Dunning DZP basis sets were applied for aluminum,32 hydrogen, and oxygen.33 The QMCF MD simulations were started from the final configuration of the previous QM/MM MD simulation,5 including the Al(III) ion and its first and second hydration shells into the QM region (QM radius of 5.0 Å). Although it would be desirable to include even the third shell in the quantum mechanical treatment, the associated increase of the computational effort (the third shell contains more molecules than the first plus second shell) leads to a considerable prolongation of the computing time of up to a few years. As the QM/MM methodology has been explicitly designed for such scenarios, the discussed partitioning of the system is justified. Important charge transfer, polarization, and many-body effects close to the ion are treated on the basis of quantum mechanics, and the improved coupling realized in the QMCF framework also accounts for these effects. The development of advanced computational facilities as well as the ongoing improvement of quantum mechanical computations will enable three-shell treatments of hydrated ionic systems in the near future. These data will help to interpret and assess the results presented in this
Al(III) Hydration Revisited
J. Phys. Chem. B, Vol. 112, No. 37, 2008 11729
Figure 2. (a) Al-O and (b) Al-H radial distribution functions and their running integration numbers for QMCF MD simulations with Mulliken (solid line) and Lo¨wdin (dashed line) charges and a conventional QM/MM MD simulation (dotted line) of Al(III) in aqueous solution.
TABLE 2: Maxima (rM, Å) and Minima (rm, Å) of the Al(III)-O Radial Distribution Function (gAl(III)-O) and Average Coordination Numbers (CNs) of the Respective Shells rM1 QMCF (Mulliken) QMCF (Lo¨wdin) two-shell QM/MM5 exptl value45–47 Al(III) MM MD38 AlCl3 MM MD36 Al(III) CP MD39 AlCl3 CP MD40 Al(III) CP MD41
1.88 1.88 1.86 1.87-1.90 1.9 2.0 1.92 1.92 1.93
rm1
2.6 2.5
rM2 4.15 4.15 4.10 3.98-4.15 4.0 4.3/4.5 4.09 4.09 4.08
study, but it is expected that only minor deviations occur, leading to a fine-tuning of the results. After 3 ps of equilibration 100 000 steps were sampled, resulting in a simulation time of 20 ps. The total computing time of one simulation amounted to six months on two AMD Opteron 280 processors (4 CPU cores). Besides radial, angular, and coordination number distribution functions, vibrational frequencies and the evaluation of ligand exchange dynamics were employed to obtain information about structural and dynamical properties of the aluminum(III) hydrate and to compare the QMCF and the conventional QM/MM MD simulations. Results and Discussion The Al-O radial distribution functions (RDF; Figure 2a) obtained from the QMCF MD simulations agree well with the result of the previous investigations including two hydration shells into the QM region. Three distinct hydration layers can be deduced; their characteristic data are listed in Table 2. The most pronounced difference is the significant reduction of the first-shell peak’s intensity when the QMCF approach is applied. While the peak is located at 1.9 Å in the case of all three simulations, the intensity decreases from 32.0 in the case of the conventional QM/MM MD simulation to 26.9 and 25.8 in the case of the QMCF MD simulation employing Mulliken and Lo¨wdin population analysis, respectively. As no exchange of first-shell ligands took place in any of the simulations as expected, the first-shell coordination number remains fixed at 6. The half-widths of the first-shell peaks determined as 0.116, 0.142, and 0.146 Å increase, indicating that the QMCF scheme predicts more flexibility for the hydration layer, and as the width of the first-shell peak is an indicator of vibrational motion, differences in the frequencies of the vibrational modes are expected. The second-shell structure is significantly altered on changing from the QM/MM to the QMCF methodology. In addition to a
rm2
rM3
rm3
4.85 5.05 4.75
6.25 6.6 6.25
7.20 7.25 7.15
4.6 ∼5.3
∼6.2 ∼6.8
∼7.2 ∼7.5
CNav, 1 6.0 6.0 6.0 6.0 6.0 6.0 6 6 6
CN
av, 2
12.8 13.0 12.2 12-14 ∼14 ∼19 11.2 12 12
CN
av, 3
31.6 28.8 37.3
decrease of the peak’s intensity (from 3.5 to 3.1 and 2.6), an increase of the second-shell distance is observed. While the peak maximum is found at 4.1 Å in the case of the QM/MM simulation, application of the QMCF scheme yields secondshell distances of 4.15 Å. The differences of the positions of the minimum between the second and third shells are even more pronounced. While the QM/MM and QMCF (Mulliken) simulations show a rather similar location of the minimum at 4.7 and 4.8 Å, a value 5.1 Å is obtained from the simulation employing Lo¨wdin population analysis. The average second-shell coordination number and the peak’s half-width have increased slightly. The CNs are 12.2, 12.8, and 13.0; the half-widths were determined as 0.46, 0.50, and 0.56 Å in the case of the QM/ MM and the QMCF MD simulations, respectively. One should also take note of the deformation of the second-shell peak observed in the conventional QM/MM treatment. The steepness of the slope near the border region is considerably high, whereas the peaks resulting from the QMCF treatment are symmetric. A close investigation of the third shell (see the inset in Figure 2a) reveals some major differences as well. In the case of the conventional QM/MM framework the third-shell peak rises rapidly after the minimum. The peak is very broad, and the associated coordination number of 37.4 is higher compared to the results of the QMCF framework, namely, 31.6 and 28.8 for Mulliken and Lo¨wdin charges, respectively. On the other hand, the third-shell peak observed from the QMCF MD simulation utilizing Lo¨wdin population analysis shows besides its rather low population a significantly larger distance of 6.6 Å compared to the peak maximum of 6.25 Å resulting from the other two simulations. The minima of the third-shell peak are found at similar distances, however, namely, at 7.15 (QM/MM), 7.20 (QMCF Mulliken), and 7.25 Å (QMCF Lo¨wdin). The comparison of the coordination number distributions (CNDs; cf. Figure 3) further demonstrates the differences observed in the simulations. While the distribution of the secondshell coordination number obtained from the conventional QM/
11730 J. Phys. Chem. B, Vol. 112, No. 37, 2008
Figure 3. Coordination number distributions obtained from QMCF MD simulations with Mulliken and Lo¨wdin charges and a conventional QM/MM MD simulation of Al(III) in aqueous solution.
MM simulations is similar to the QMCF (Mulliken) scheme, the result for the third shell is clearly overestimated. The CND of the QMCF (Lo¨wdin) simulation on the other hand indicates a more labile structure, resulting in lowered intensities. An increased exchange rate of the solvent molecules between the second and third shells compared to that of the other simulations is expected. Comparison with available experimental and theoretical data leads to the conclusion that the QMCF MD simulation utilizing Mulliken population analysis yields the most reliable structural data for the first and second hydration shells. No experimental data are available for the third hydration shell. In fact, the existence of a third shell has been controversial in the literature; however, recent experimental investigations utilizing infrared spectroscopy in HDO34 and dielectric spectroscopy35 confirmed the existence of preferentially oriented water molecules beyond the second hydration shell. Most classical simulations tend to overestimate the ion-solvent interactions and hence lead in many cases to an inaccurate description of the system.5,36,37 The hydrated ion model utilized by Martinez and co-workers38 yields data in agreement with experimental estimations. Nevertheless, clear indications of third-shell structures are obtained by most classical simulations. The Al(III) ion in aqueous solution has been investigated by Car-Parrinello (CP)-type simulations as well. Amira and coworkers39 carried out two CP MD simulations utilizing different exchange-correlation functionals in a box containing 32 D2O molecules. Although structural parameters such as bond lengths and coordination numbers were predicted in agreement with experimental data, sensitive parameters such as the OD stretching frequency showed deviations due to “shortcomings of the BLYP and PBE functionals”.39 It was concluded that a number of empirical correction factors are required to compensate for shortcomings related to the exchange-correlation functionals, the fictious electron masses, and anharmonicity effects. Although it was postulated that box size effects should have no influence on the properties of the hydrate, the low number of particles does not allow the formation of a complete hydrate, i.e., a system containing three hydration shells and a substantial amount of bulk molecules. The low second-shell coordination number reported as 11.239 underlines this conclusion. Similar structural properties have been obtained from a CP MD study of AlCl3 in 64 D2O molecules performed by Ikeda and co-workers.40 A recent CP study41 indicated that an increase of the number of particles (64 and 128 D2O molecules, respectively) is
Hofer et al. mandatory to yield a reasonable description of the transition region between the second shell and the bulk. Consequently, the increased number of particles has an impact on the results for the hydrate structure as well as on rate constants for migration to and from the second hydration shell, and it was concluded “to reliably predict the behavior in these regions, sizes larger than 64 [solvent molecules] may be required”. Although this study41 indicated the presence of dipole-oriented molecules beyond the second shell, no data could be extracted on the third hydration shell. Some differences are also visible in the Al-H radial distribution functions (cf. Figure 2b). Whereas the first-shell peaks located at 2.55 Å differ only in the intensity, significant changes for the second and third shells can be deduced. Similar to the case of the Al-O RDF, the second-shell peak is shifted to larger distances (from 4.7 to 4.8 Å) upon application of the QMCF framework. The peak heights resulting from the latter simulations were found at 2.15 (Mulliken charges) and 1.57 (Lo¨wdin charges), resulting in coordination numbers of 22.6 and 18.8, respectively, compared to a value of 19.8 deduced from the conventional QM/MM MD treatment. The most significant differences, however, can be identified in the transition region between the second and third shells, coinciding with the QM/MM transition region. Similar to the case of the Al-O RDF, the conventional QM/MM MD framework yields a deformed second-shell peak and a weakly defined minimum between the second and third shells, and the third-shell peak forms a broad plateau rather than a distinct peak. The description of this region is significantly improved when the QMCF ansatz is utilized. These data demonstrate the sensibility of QM/MM simulations on the coupling of the subregions. In the presented studies the partial charges of the QM atoms defined by population analysis or by assignment of fixed charges result in substantially different Coulombic forces between QM and MM atoms, which apparently have a high impact on particles close to the QM/ MM interface. The clearly ambiguous behavior of the conventional QM/ MM framework can be explained by a “pressure artifact” resulting from the assignment of fixed group charges according to the MM model on the QM atoms. Along with a high charge of +3.0 assigned to the Al(III) ion, the charges of the water molecules do not account for polarization and many-body effects; thus, the central charge of the ion is not shielded properly. As a consequence water molecules (in particular the oxygen atoms) located in the MM region are attracted too strongly, leading to a significant “pressure” of the MM molecules on the QM region, which leads to the observed high population of the third shell, the lowered second-shell distance and deformations of the second and third shell peaks, which are especially pronounced in the Al-H RDF. Although the QM region is capable of accounting for polarization, charge-transfer, and many-body effects, the associated relaxation of the hydration structure is not possible due to the ambiguous behavior of the MM molecules. In the case of less rigid systems these inconsistencies lead to artifacts near the QM/MM transition region. Particles leaving the QM region due to an adequate treatment (i.e., the inclusion of polarization, charge-transfer, and many-body effects) are forced to migrate back into the QM zone by the inaccurate, i.e., too strong Coulombic forces acting on MM atoms. As a consequence an accumulation of particles in the transition region occurs, which results in small peaks in the radial distribution functions near the QM/MM border.
Al(III) Hydration Revisited The QMCF scheme on the other hand assigns varying partial charges to the QM atoms on the basis of population analysis which reflect the respective polarization as well as many-body interactions and take into account all geometrical changes. As a consequence the charge on the ion is smaller than the formal oxidation number, while the charges on oxygen and hydrogen atoms are modulated according to the polarization and charge transfer induced by the ion. Compatibility of the charges obtained from the quantum mechanical treatment with those of the respective force field model utilized to represent the solvent in the MM region has to be ensured when the population analysis scheme is selected. Further considerations involve the selected basis sets, the respective quantum mechanical level, and the number of QM particles, which is associated with the size of the quantum mechanical region. Sample calculations have to be performed in advance to select the best method for a particular study as has been discussed in the Methods. The structural data obtained from the RDF have indicated that the description of the QM/MM transition region and the third hydration shell is less satisfactory in the case of the QMCF simulation utilizing Lo¨wdin population analysis. The results obtained from its radial distribution functions appear less accurate compared to those from the simulation with Mulliken charges, which, despite some general criticism found in the literature, proved to be the most adequate choice in this particular study. The intensities of the minima deduced from the RDF plots as well as the distribution of the coordination numbers indicate differences in the solvent exchange rate between the second and third hydration layers. Therefore, the mean ligand residence times of second-shell ligands were determined according to direct measurements of all occurring exchange events, accepting only ligand displacements lasting longer than the hydrogen bond lifetime of 0.5 ps as a valid exchange event.50 Migrations below this threshold are considered as short time fluctuations and discarded. The MRT values obtained from the different simulations are 26.4 (QM/MM), 17.7 (QMCF Mulliken), and 13.3 ps (QMCF Lo¨wdin). The observed decrease is proportional to the decrease of the effective charge of the Al(III) ion according to the different charge models. The QM/MM framework has a fixed charge of 3.0; the average charges from the simulations with Mulliken and Lo¨wdin population analysis are 2.5 and 2.05 unit charges, respectively. Although the difference of the MRT appears significant, the logarithmic scale of the water exchange rates has to be considered. The range of first-shell ligand exchange rates of cationic hydrates covers more than 20 orders of magnitude ranging from of Ir(III) to Cs(I).51 Thus, the values are considered in good agreement with each other, and the influence of the different Coulombic forces has to be considered minor. Consideration of the MRTs and the structural data determined from the RDFs leads to the conclusion that an upper and lower bound for the properties of this systems can be given on the basis of the conventional QM/MM and the QMCF (Lo¨wdin) framework resulting from either too strong or too weak charge field interactions. The MRT value deduced from the investigation with Mulliken charges of 18 ps is halfway between the data obtained from the other two simulations and therefore appears to be the best estimate. Accordingly, the structural description appears most reliable in the case of the QMCF (Mulliken) study as well. Despite the requirement of a suitable QM/MM coupling framework, the main target of this type of study is the solute located in the center of the QM region. In this particular case
J. Phys. Chem. B, Vol. 112, No. 37, 2008 11731
Figure 4. Angular distribution functions of the first hydration shell of Al(III) in aqueous solution obtained from QMCF simulations with Mulliken (solid line) and Lo¨wdin (dashed line) charges and a conventional QM/MM MD simulation (dotted line).
the interaction of the strongly polarizing Al(III) ion with the solvent is the target of the study. Besides the structural description via radial distributions, the characterization of the first hydration shell via angular distribution functions (ADFs) enables a detailed insight into hydration phenomena. Figure 4 depicts the normalized O-Al-O, Al-O-H, θ and tilt angle distributions as well as the tilt-angle distributions within the first shell (for definitions see the sketches in Figure 4). The O-Al-O angles (Figure 4a) reflect the octahedral arrangement of the six first-shell ligands. The peaks at 90° show only small deviations, the QM/MM scheme displays the lowest (8.4) and the QMCF (Lo¨wdin) the highest intensity (10.2), and the peak of the QMCF (Mulliken) study has a height of 9.2. The intensity at 180° displays stronger differences. Although both QMCF simulations yield similar peak intensities of ∼75 and ∼80, the corresponding value deduced from the QM/MM scheme is approximately 105. As the broadness of all peaks is inversely proportional to the peak height, the latter is a direct measure of the flexibility of the first-shell ligands. Accordingly, it can be concluded that the QMCF simulations predict a higher degree of flexibility than the conventional QM/MM schemes. The distribution of the Al-O-H angle (Figure 4b) peaks at 125° for all simulations; however, the respective intensities are different. The order and the differences are similar to those in the case of the 90° peak of the O-Al-O angles. Apparently, angular flexibility of the first-shell ligands is reduced when the QMCF scheme is applied. A similar conclusion can be drawn on the basis of the θ angle distribution (cf. Figure 4d). This effect appears to be more pronounced in the simulation utilizing the Lo¨wdin charge scheme. The tilt-angle distributions (cf. Figure 4c) are almost identical, the QMCF (Mulliken) simula-
11732 J. Phys. Chem. B, Vol. 112, No. 37, 2008 tions yielding a slightly lower intensity, whereas the peak is slightly displaced in the QMCF (Lo¨wdin) framework. Despite these minor changes visible in the various ADFs, all structural features are retained when the QMCF methodology is applied, thus leading to the conclusion that this methodology leads to results similar to those of the conventional QM/MM framework. Still the QMCF framework has the benefit that no potential functions between Al(III) and H2O need to be applied. In view of the tedious construction process of potentials required for the QM/MM study and the associated lower accuracy resulting in artifacts visible in the radial distribution functions, it can be concluded that the QMCF framework offers a substantial improvement for simulation studies of solutes. This advantage has already been utilized in simulation studies of the hydration properties of more complex systems with nonspherical potential distribution such as Pd(II) and Pt(II)1,7 or composite species such as Hg22+,9 SO42-,11,12 and PO43-.10 As the values are only slightly modified, comparison to experimental data does not allow an assessment of whether the accuracy of the results is indeed increased by the QMCF framework. Comparative simulations including all three hydration shells could provide valuable information to determine which of the simulations yields the best results. As the third shell contains more solvent molecules than the first plus second shell, the computational effort associated with the QM treatment is not yet manageable, however. The vibrational motions are sensitive indicators for the accuracy of a simulation study as they are proportional to the second derivative of the energy, and therefore, the octahedral modes of the vibrations were analyzed. Velocity autocorrelation functions (VACFs) are utilized to derive spectra from velocity data obtained in the simulations. A correlation length of 2 ps (2000 configurations, time step between configurations 1 ps) was utilized to analyze the VACFs over the whole trajectory after an approximative normal coordinate analysis42 had been carried out. The corresponding spectra were obtained by Fourier transformation of the time-averaged VACFs and are compared with the data of the conventional QM/MM methodology in Figure 5; the corresponding data are given in Table 3. It is important to note that no scaling has been applied to the spectra. The vertical lines in Figure 5 indicate the experimental values. The Eg, T2g, and the higher T1u frequencies are in excellent agreement with the experimental value if the QMCF MD approach is applied. A1g differs by about 6.7%, the lower T1u frequency by 14.5%. In the case of all six frequencies, the results obtained with the QMCF scheme are at least equal and in many cases improved compared to the conventional QM/MM framework. The shapes of the peaks are improved, and side peaks visible in the spectra resulting from the QM/MM framework have vanished or at least significantly decreased. Another interesting fact is that the influence of the population analysis scheme on the spectra is almost negligible. This indicates that the Coulombic forces between QM and MM atoms are too small to influence the vibrational motions significantly. All important contributions are, therefore, treated at the quantum mechanical level. The differences of some vibrational motions from the experimental results can be explained by the missing treatment of electron correlation in the Hartree-Fock method. Furthermore, it should be considered that the system is not a true octahedron, which will also lead to deviations in the frequencies. On the other hand, the presence of counterions in the experimental investigations influences the dielectric properties, and hence, the experimental environment differs from that of the system utilized in this simulation study.
Hofer et al.
Figure 5. Unscaled power spectra of the Al(III)-O octahedral modes obtained from QMCF simulations with Mulliken (solid line) and Lo¨wdin (dashed line) charges and a conventional QM/MM MD simulation (dotted line). The experimental values are indicated by vertical lines.
TABLE 3: Unscaled Peak Maxima (cm Octahedral Modesa QMCF (Mulliken) QMCF (Lo¨wdin) 2 shell QM/MM6 experimental data48 experimental data49 a
A1g
Eg
560 561 583 525 525
428 416 440 438 438
-1)
606 595 640 598 598
of the Al(III)-O T1u
T2g
(357) (357) (360)
325 325 336 332 332
(309, 268)
Side peaks/shoulders are given in parentheses.
As mentioned earlier the scaling factor of 0.89 for HartreeFock43,44 has not been applied. This empirical correction factor corrects not only for the missing treatment of electron correlation, but also for the typical vacuum environment and the influence of anharmonicity. The latter contributions are explicitly included in QMCF MD studies as well as the influence of bulk on the hydrate. Cluster computations at different levels of theory have indicated that the influence of electron correlation is minor,5 which explains the accuracy of the spectra obtained via application of the QMCF approach. In particular, the improved QM/MM coupling and the polarization of the QM region via inclusion of MM point charges into the QM regions greatly enhance the quality of the results compared to those of the conventional QM/MM framework. Conclusion On the basis of the structural as well as dynamical data obtained from the MD simulations of Al(III) in aqueous solution presented here, it can be concluded that the QMCF framework is in several aspects a significant improvement of the conventional QM/MM methodology. One important aspect in the execution of QMCF MD simulations is the choice of a suitable population analysis scheme. In this study Mulliken partial charges led to an excellent description of the system, whereas the application of Lo¨wdin charges yielded a too weakly
Al(III) Hydration Revisited structured hydrate, in particular concerning the third hydration shell. Test calculations at the respective computational level are recommendable to find the best point charge scheme for any particular application. Dynamical data obtained from the simulation are in good agreement with the experimental data. Again the QMCF simulation proved superior to the conventional QM/MM scheme. It is important to note that the vibrational spectra do not require the application of empirical scaling factors to be comparable with the experimental data. The increased accuracy of the novel simulation method in combination with the nonrequirement of solvent-solute potential functions proves the QMCF MD methodology a most valuable tool for the treatment of solvated systems with a very wide spectrum of applications. Future computational equipment will enable the inclusion of even the third hydration shell, but it is expected that, besides a significant increase of the computational effort, only minor changes of the presented data will be observed. Acknowledgment. Financial support for this work from the Austrian Science Foundation (FWF) and the Hypo Tirol Bank is gratefully acknowledged. References and Notes (1) Rode, B. M.; Hofer, T. S.; Randolf, B. R.; Schwenk, C. F.; Xenides, D.; Vchirawongkwin, V. Theor. Chim. Acta 2006, 115, 77. (2) Rode, B. M.; Hofer, T. S.; Randolf, B. R.; Strangret, J.; Pribil, A. B.; Vchirawongkwin, V. In Trends and PerspectiVes in Modern Computational Science; Maroulis, T. S. G. Ed.; International Science Publishers (VSP): Leiden, The Netherlands, 2006; p 441. (3) Shannon, R. D. Acta Crystallogr. 1976, 751, A32. (4) da Silva, J. J. R. F.; Williams, R. J. P. The Biological Chemistry of Life; Clarendon Press: Oxford, U.K., 1991. (5) Hofer, T. S.; Randolf, B. R.; Rode, B. M. Phys. Chem. Chem. Phys. 2005, 7, 1382. (6) Hofer, T. S.; Randolf, B. R.; Rode, B. M. Chem. Phys. Lett. 2006, 422, 492. (7) Shah, S. A. A.; Hofer, T. S.; Fatmi, M. Q.; Randolf, B. R.; Rode, B. M. Chem. Phys. Lett. 2007, 445, 193. (8) Fatmi, M. Q.; Hofer, T. S.; Randolf, B. R.; Rode, B. M. J. Comput. Chem. 2007, 28, 1704. (9) Hofer, T. S.; Randolf, B. R.; Rode, B. M. Chem. Phys. 2008, 349, 210–218. (10) Pribil, A. B.; Hofer, T. S.; Vchirawongkwin, V.; Randolf, B. R.; Rode, B. M. Chem. Phys., in press. (11) Vchirawongkwin, V.; Persson, I.; Rode, B. M. J. Phys. Chem. B 2007, 111, 4150. (12) Vchirawongkwin, V.; Rode, B. M. Chem. Phys. Lett. 2007, 443, 152. (13) Warshel, A.; Levitt, M. J. Mol. Biol. 1976, 103, 227. (14) Field, M. J.; Bash, P. A.; Karplus, M. J. Comput. Chem. 1990, 11 (6), 700–733.
J. Phys. Chem. B, Vol. 112, No. 37, 2008 11733 (15) Gao, J. J. Am. Chem. Soc. 1993, 115, 2930–2935. (16) Bakowies, D.; Thiel, W. J. Phys. Chem. 1996, 100 (25), 10580– 10594. (17) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. J. Chem. Phys. 1993, 97, 10269. (18) Lin, H.; Truhlar, D. G. Theor. Chim. Acta 2007, 117, 185. (19) Laio, A.; VandeVondele, J.; Rothlisberger, U. J. Chem. Phys. 2002, 116, 6941. (20) Voloshina, E.; Gaston, N.; Paulus, B. J. Chem. Phys. 2007, 126, 134115. (21) Ahlrichs, R.; Ba¨r, M.; Ha¨ser, M.; Horn, H.; Ko¨lmel, C. Chem. Phys. Lett. 1989, 162 (3), 165. (22) Brode, S.; Horn, H.; Ehrig, M.; Moldrup, D.; Rice, J. E.; Ahlrichs, R. J. Comput. Chem. 1993, 14 (10), 1142–1148. (23) Mulliken, R. S. J. Chem. Phys. 1955, 23 (10), 1833–1840. (24) Mulliken, R. S. J. Chem. Phys. 1955, 23 (10), 1841–1846. (25) Lo¨wdin, O. P. J. Chem. Phys. 1950, 18, 365. (26) Reed, A. E.; Weinstock, R. B.; Weinhold, F. J. Chem. Phys. 1985, 83, 735. (27) Heinzmann, R.; Ahlrichs, R. Theor. Chim. Acta 1976, 42, 33. (28) Stillinger, F. H.; Rahman, A. J. Chem. Phys. 1978, 68 (2), 666. (29) Bopp, P.; Jansco´, G.; Heinzinger, K. Chem. Phys. Lett. 1983, 98 (2), 129. (30) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Phys. Chem. 1984, 81, 3684–3690. (31) Adams, D. J.; Adams, E. M.; Hills, G. J. Mol. Phys. 1979, 38 (2), 387–400. (32) Dunning, T. H., Jr.; Hay, P. J. In Modern Theoretical Chemistry: Methods of Electronic Structure Theory; Schaefer, H. F., III Ed.; Plenum Press: New York, 1977; Vol. 3. (33) Dunning, T. H., Jr J. Chem. Phys. 1970, 53, 2823–2833. (34) Strangret, J.; Gampe, T. J. Phys. Chem. A 2002, 106, 5393–5402. (35) Schro¨dle, S.; Rudolph, W. W.; Hefter, G.; Buchner, R. Geochim. Cosmochim. Acta 2007, 71, 5287. (36) Lauenstein, A.; Hermansson, K.; Lindgren, J.; Probst, M.; Bopp, P. A. Int. J. Quantum Chem. 2000, 80, 892. (37) Spångberg, D.; Hermansson, K. J. Chem. Phys. 2004, 120, 4829. (38) Martinez, J. M.; Pappalardo, R. R.; Marcos, E. S. J. Am. Chem. Soc. 1999, 121, 3175. (39) Amira, S.; Spångberg, D.; Hermansson, K. J. Chem. Phys. 2006, 124, 104501. (40) Ikeda, T.; Hirata, M.; Kimura, T. J. Chem. Phys. 2003, 124, 104501. (41) Bylaska, E. J.; Valiev, M.; Rustad, J. R.; Weare, J. H. J. Chem. Phys. 2007, 126, 104505. (42) Bopp, P. Chem. Phys. 1986, 106, 205. (43) Scott, A. P.; Radom, L. J. Phys. Chem. 1996, 100, 16502–16513. (44) DeFrees, D. J.; McLean, A. D. J. Chem. Phys. 1985, 82 (1), 333– 341. (45) Ohtaki, H.; Radnai, T. Chem. ReV. 1993, 93 (3), 1157–1204. (46) Bol, W.; Welzen, T. Chem. Phys. Lett. 1977, 46, 189. (47) Johansson, G. AdV. Inorg. Chem. 1992, 39, 159. (48) Rudolph, W. W.; Mason, R.; Pye, C. C. Phys. Chem. Chem. Phys. 2000, 2, 5030. (49) Mink, C. N. J.; Hajba, L.; Sandstro¨m, M.; Goggin, P. L. J. Mol. Struct. 2003, 141, 661–662. (50) Hofer, T. S.; Tran, H. T.; Schwenk, C. F.; Rode, B. M. J. Comput. Chem. 2004, 25, 211–214. (51) Helm, L.; Merbach, A. E. Coord. Chem. ReV. 1999, 187, 151.
JP802663H
