Solvation Structure and Dynamics of Ammonium (NH4+) in Liquid

Among the hydrogen-bonded liquids, ammonia (NH3) is special because it is ... combined quantum mechanics/molecular mechanics (QM/MM) method. ... to el...
3 downloads 0 Views 298KB Size
Download PDF
J. Phys. Chem. B 2008, 112, 885-891

885

Solvation Structure and Dynamics of Ammonium (NH4+) in Liquid Ammonia Studied by HF/MM and B3LYP/MM Molecular Dynamics Simulations Anan Tongraar*,† and Supot Hannongbua‡ School of Chemistry, Institute of Science, Suranaree UniVersity of Technology, Nakhon Ratchasima 30000, Thailand, Department of Chemistry, Faculty of Science, Chulalongkorn UniVersity, Bangkok 10330, Thailand ReceiVed: August 2, 2007; In Final Form: October 4, 2007

The characteristics of NH4+ solvated in liquid ammonia have been investigated by means of combined HF/ MM and B3LYP/MM molecular dynamics simulations, in which the ion and its surrounding ammonia molecules were treated by HF and B3LYP methods, respectively, using the D95* basis set. For both HF/MM and B3LYP/ MM simulations, it is observed that four nearest-neighbor ammonia molecules directly hydrogen-bonded to each of the ammonium hydrogen atoms, forming a well-defined tetrahedral cage structure of the NH4+ solvate. Nevertheless, the solvation shell of NH4+ is rather flexible, in which several possible species of solvated NH4+ exist, ranging from 4- to 7-fold and from 3- to 6-fold coordinated complexes for the HF/MM and B3LYP/MM simulations, respectively. In terms of the dynamical details, i.e., the self-diffusion coefficients and the mean residence times of ammonia molecules surrounding the ion, the B3LYP/MM simulation shows slower dynamics of the solvated NH4+ when compared with the HF/MM results. With regard to the reported tendency of density functional theory (DFT) methods to predict overly rigid ion solvations as well as hydrogen bonds that are too short, the ab initio HF formalism has been demonstrated to be more reliable for providing a detailed description of this solvated ion.

1. Introduction Among the hydrogen-bonded liquids, ammonia (NH3) is special because it is known as one of the weakest hydrogenbonded liquids found in nature and is often used as a solvent in a variety of organic reactions. Analogously, its associated ions, in particular ammonium (NH4+), also play an important role in chemical and biological processes.1-4 According to existing detailed analysis of the structure, thermodynamics, and spectral data for NH4+ in liquid ammonia,5,6 it has been demonstrated that this ion forms a stable solvation shell with a coordination number of 4. In a subsequent study by Liu and Tuckerman,7 the solvation structure and dynamics of NH4+ in liquid ammonia were evaluated on the basis of the Car-Parrinello molecular dynamics (CP-MD) approach.8 One clear finding is that the NH4+ possessed a tight 4-fold solvation shell, comprised almost entirely of hydrogen bonds donated by the ion. However, this approach was applied to a relatively small system, i.e., a replicated cubic box made up of only 31 ammonia molecules plus one ion, and all interactions were evaluated by a simple BLYP functional with limited accuracy. In particular, for the description of hydrogen bonds between bulk ammonia molecules, the CP-MD simulation has reported somewhat different results for the characteristics of liquid ammonia when compared with the corresponding data obtained from an earlier CP-MD study.9 By means of quantum-mechanics-based simulations, an alternative approach to obtain detailed knowledge of condensedphase systems is to apply a so-called combined quantum mechanics/molecular mechanics (QM/MM) method. During the * Corresponding author. E-mail: [email protected]; fax: 006644-224185. † Suranaree University of Technology. ‡ Chulalongkorn University.

past decades, the QM/MM MD technique has been proven to be a very reliable simulation scheme to elucidate microscopic details of various ionic solutions.10-21 For the study of hydrogenbonded liquids, the QM/MM simulations of HF/MM, B3LYP/ MM, and MP2/MM types have been performed for liquid water.22 Interestingly, it has been demonstrated that the HFand MP2-based simulations provided detailed descriptions of liquid water that agree well with the most recent experimental data, while the density functional theory (DFT) method, even with the B3LYP functional, has predicted improper dynamics due to overly rigid hydrogen bonds. On the basis of the HF/ MM scheme, however, the correct description of liquid water could only be achieved when the second-shell waters were included in the QM region.22 The failure of the HF method to reproduce the structural and dynamical properties of liquid water has been demonstrated. For example, a recent MD study by Todorova and co-workers23 showed that the use of HF approximations has led to strongly under-structured liquid water. In light of the previous simulation data for water, the reliability of the HF method for describing other liquids, such as liquid ammonia, might be doubtful. Recently, the HF/MM and B3LYP/ MM MD simulations have been carried out for liquid ammonia,24 supplying information that the hydrogen bonds in liquid ammonia are rather weak. The structural features and related dynamical properties of this liquid are suggested to be determined by the steric packing effects, rather than by the hydrogenbond interactions. In particular, the HF/MM and B3LYP/MM simulations have predicted properties of liquid ammonia in better agreement with experimental data, especially when compared with previous studies based on either an empirical or CP-MD approach.9,25 Since the results yielded from the HF/MM and B3LYP/MM simulations24 (cf. Figure 2) are almost comparable, the HF method, which consumes less time, was chosen for this study. On the other hand, it could be expected that the HF

10.1021/jp076173t CCC: $40.75 © 2008 American Chemical Society Published on Web 01/01/2008

886 J. Phys. Chem. B, Vol. 112, No. 3, 2008

Tongraar and Hannongbua anions,17,21 while the B3LYP method often gives poor results for several cases.17,18,21,22 Nevertheless, the B3LYP method was employed in the present work in order to test its adequacy for the description of this particular system. In fact, it is known that the HF scheme could produce an error due to the neglect of electron correlation effects, while the DFT methods, although including such effects to a certain (uncontrollable) extent, are often found to overestimate the correlation energy.17,18,21,22 On the other hand, a comparison of the HF calculations with the DFT results could be helpful to give a qualitative estimate of the possible influence of correlation effects. 2. Methods On the basis of the QM/MM MD technique,10-22 the system is partitioned into two parts, namely, QM and MM regions. The total interaction energy of the system is defined as

Etotal ) 〈ΨQM|H ˆ |ΨQM〉 + EMM + EQM-MM Figure 1. (a) Ni-N and (b) Ni-H RDFs and their corresponding integration numbers, where Ni denotes the N atom of NH4+ and other atoms belong to ammonia molecules.

where 〈ΨQM|H ˆ |ΨQM〉 refers to the interactions within the QM region, while EMM and EQM-MM represent the interactions within the MM and between the QM and MM regions, respectively. The QM region, a sphere which includes the NH4+ and its surrounding ammonia molecules, is treated quantum mechanically using HF and B3LYP methods. Since the computational expense for QM force calculations using large basis sets is significant, all quantum mechanical calculations were carried out using the double-ζ plus polarization (D95*) basis set,26 which is considered a suitable compromise between the quality of the simulation results and limits on CPU time. For the size of the QM region, a QM diameter of 8.8 Å was chosen, which includes the whole first solvation shell of NH4+ plus some next nearest-neighbor ammonia molecules. For the interactions within the MM and between the QM and MM regions, a flexible model,27,28 which describes inter- and intramolecular interactions, was employed for ammonia. The use of a flexible model is favorable, ensuring compatibility and a smooth transition when ammonia molecules move from the QM region with their full flexibility to the MM region, and vice versa. The pair potential functions for describing NH4+NH3 interactions were newly constructed using the aug-cc-pvtz basis set.29-31 With respect to the symmetric tetrahedral geometry of NH4+, 3556 HF and 3568 B3LYP interaction energy points for various NH4+-NH3 configurations, obtained from Gaussian 0332 calculations, were fitted to an analytical form of 5

∆ENH4+-NH3 ) Figure 2. (a) N-N, (b) N-H, and (c) H-H RDFs for atoms in liquid ammonia and their corresponding integration numbers, comparing the results obtained from the previous QM/MM simulations24 and NDIS experiment.39

method might provide a reasonable description of liquid ammonia, despite its failure for liquid water. In the present study, it is of particular interest to apply the high-level QM/MM simulation technique for elucidating the solvation structure and dynamics of NH4+ in liquid ammonia. Since the high correlated methods, even at the MP2 level of accuracy, are still far too time-consuming, the QM/MM simulations of HF/MM and B3LYP/MM types are considered as the applicable choices. The HF method has been well demonstrated in previous QM/MM studies,10-22 even for the treatment of

(1)

4

∑ ∑ i)1 j)1

[

Aij

rij5

+

Bij rij6

+ Cij exp(-Dijrij) +

]

qiqj rij

(2)

where A, B, C, and D are the fitting parameters (see Table 1), rij denotes the distances between the ith atom of NH4+ and the jth atom of the ammonia molecule, and q values are atomic net charges. In the present study, the charges on N and H of NH4+ were obtained from natural bond orbital (NBO) analysis33-35 of the corresponding HF and B3LYP calculations. They were set to -0.8048 and 0.4512 (HF), and -0.8176 and 0.4544 (B3LYP), respectively. The charges on N and H of the ammonia molecule were adopted from the flexible ammonia model27,28 as -0.8022 and 0.2674, respectively. According to the fact that an exchange of ammonia molecules between the QM and MM regions can take place frequently

Structure and Dynamics of NH4+ in Liquid Ammonia

J. Phys. Chem. B, Vol. 112, No. 3, 2008 887

TABLE 1: Optimized Parameters of the Analytical Pair Potential for the Interaction of NH4+ with NH3 (Interaction Energies in kcal‚mol-1 and Distances in Å) C (kcal mol-1)

D (Å-1)

-7216.4877 278.4976 503.0343 53.5965

HF-Based Method 11694.3117 6333.4144 268.7770 -1212.1552 -402.0212 -416.0600 -36.2258 11.8534

1.9796 2.2727 1.5931 1.4536

-7280.9939 304.3508 503.9488 53.3950

B3LYP-Based Method 11763.8137 6392.2237 242.8072 -1247.2617 -402.7317 -415.91934 -36.0786 11.5735

1.9827 2.2785 1.5963 1.4716

pair

A (kcal mol-1 Å5)

Ni-N Ni-H Hi-N Hi-H Ni-N Ni-H Hi-N Hi-H

B (kcal mol-1 Å6)

during the QM/MM simulations, the forces acting on each particle in the system were switched according to which region the ammonia molecule was entering or leaving the QM region and is defined as

Fi ) Sm(r)FQM + (1 - Sm(r))FMM

(3)

where FQM and FMM are the quantum mechanical and molecular mechanical forces, respectively, and Sm(r) is a smoothing function36 described by

Sm(r) ) 1, for r e r1 Sm(r) )

(r02 - r2)2(r02 + 2r2 - 3r12) , for r1 < r e r0 (r02 - r12)3 Sm(r) ) 0, for r > r0

(4)

where r1 and r0 are the distances characterizing the start and the end of the QM region, applied within an interval of 0.2 Å (i.e., for Ni-N distances of 4.4-4.6 Å, where Ni and N denote nitrogen atoms of NH4+ and of ammonia molecules, respectively) to ensure a continuous change of forces at the transition between QM and MM regions. All QM/MM MD simulations were performed in a canonical ensemble at 273 K with periodic boundary conditions. The system’s temperature was kept constant using the Berendsen algorithm.37 A periodic box, with a box length of 21.80 Å, contains one NH4+ and 255 ammonia molecules, corresponding to the experimental density of liquid ammonia of 0.6990 g cm-3. The reaction-field method38 was employed for the treatment of long-range interactions. The Newtonian equations of motions were treated by a general predictor-corrector algorithm. The time step size was set to 0.2 fs, which allows for the explicit movement of the hydrogen atoms of ammonia molecules. The system was initially equilibrated by performing preliminary HF/ MM and B3LYP/MM MD simulations (i.e., QM/MM simulations in which only the NH4+ was treated quantum mechanically using either the HF or B3LYP method, while the rest of the system was described by means of classical pair potentials) for 200 000 time steps. Then, the actual HF/MM and B3LYP/MM simulations were subsequently started with the system’s reequilibration for 100 000 time steps, followed by 120 000 and 110 000 time steps, respectively, to collect the configurations at every 10th step. 3. Results and Discussion 3.1. Structural Data. Structural details of the solvated NH4+ can be obtained from Ni-N and Ni-H radial distribution functions (RDFs) and their corresponding integration numbers,

Figure 3. Probability distributions of the coordination numbers, calculated up to the first minimum of the Ni-N RDFs.

as depicted in Figure 1. To reliably describe the ionic nature of NH4+, the corresponding atom-atom RDFs for liquid ammonia obtained at a similar QM/MM level of accuracy24 and from NDIS experiments39 were utilized for comparison (see Figure 2). As can be seen in Figure 1a, the Ni-N RDFs obtained from the HF/MM and B3LYP/MM simulations are significantly different, especially for the feature of the Ni-N peaks at the distances between 3.5 and 5.5 Å. The HF/MM simulation exhibits a sharp first Ni-N peak with a maximum at 2.95 Å, whereas the corresponding B3LYP/MM peak starts to be detected earlier, with a maximum at a shorter distance of 2.83 Å. A discrepancy is noticeable even between the B3LYP/MM results and the recent CP-MD study,8 which reported a more pronounced first Ni-N peak with a maximum at a longer distance of 3.0 Å. Thus, the observed difference between the HF/MM and B3LYP/MM simulations could be attributed to a consequence of the approximations and the parametrizations of the DFT methods, rather than to the neglect of the electron correlation effects at the HF level of theory. In both HF/MM and B3LYP/MM simulations, although a well-defined first solvation shell of NH4+ can be obtained (i.e., compared to the structural feature of the liquid ammonia, cf. Figure 2a), the first minimum of Ni-N RDFs is not distinctly separated from the bulk, suggesting that ammonia molecules in the solvation shell of NH4+ can easily exchange with other molecules in the bulk. In the present study, the minimum positions of the Ni-N RDFs obtained by the HF/MM and B3LYP/MM simulations may be roughly estimated as 3.7 and 3.3 Å, where integrations up to these minima yield average coordination numbers of 5.1 and 4.3, respectively. Figure 3 shows the probability distributions of the number of ammonia molecules surrounding the NH4+ ion, calculated up to the first minimum of the Ni-N RDFs. In the HF/MM simulation, a preferred coordination number of 5 (followed by 6, 4, and 7 in descending order) is observed for NH4+, whereas the B3LYP/ MM simulation favors a coordination number of 4 (followed by 5, 3, and 6 in descending order). The observed multiple coordination numbers, varying from 4 to 7 and 3 to 6 for the HF/MM and B3LYP/MM simulations, respectively, clearly indicate a high flexibility of the solvated NH4+. This supplies information that numerous ammonia molecules can mutually play a role in the hydrogen bond formation of NH4+ in liquid ammonia. In particular, for the B3LYP/MM simulation, it is observed that a small shift of the Ni-N RDF minimum will have significant impact on the number of surrounding ammonia molecules. For example, at the slightly larger Ni-N distance

888 J. Phys. Chem. B, Vol. 112, No. 3, 2008

Tongraar and Hannongbua

Figure 5. Distributions of the N-Ni-N angle, calculated within the first minimum of the Ni-N RDFs.

Figure 4. (a) Hi-N and (b) Hi-H RDFs and their corresponding integration numbers, where Hi denotes the H atoms of NH4+.

of 4.0 Å, the most frequently observed number of ammonia molecules increases from 4 to 7. In Figure 1b, the Ni-H RDFs obtained from both HF/MM and B3LYP/MM simulations provide detailed information in good accord with the recognizable “structure-making” ability of the NH4+ ion. The characteristics of hydrogen bonds between NH4+ and ammonia can be interpreted through the Hi-N and Hi-H RDFs, as shown in Figure 4. The HF/MM and B3LYP/MM simulations reveal first Hi-N peaks with maxima at 1.95 and 1.82 Å, respectively. These peaks are identified as being indicative of the hydrogen bonds between the hydrogen atoms of NH4+ and their nearest-neighbor ammonia molecules. In comparison with the N-H RDF of liquid ammonia (Figure 2b), e.g., in terms of shape and peak height, it is apparent that NH4+ can certainly form hydrogen bonds with its surrounding ammonia molecules. With respect to the Hi-N RDFs of both HF/MM and B3LYP/ MM simulations, integrations up to the corresponding first Hi-N minimum yield about 1.1 and 1.0 nitrogen atoms being bonded to each hydrogen atom of NH4+. These data suggest that the first solvation shell of NH4+ is rather dominated by hydrogenbond donation to its neighboring ammonia molecules. In this work, the first Hi-N peak of the B3LYP/MM simulation starts to be detected at a distance about 0.25 Å shorter than that observed in the HF/MM simulation, implying that the inclusion of electron correlation reflects in stronger ion-ammonia hydrogen-bond interactions. Nevertheless, according to the recent CP-MD study,8 which reported a tight solvation shell structure of the NH4+ solvate together with a distinct first Hi-N peak, the overestimation of electron correlations by the DFT methods is apparent.17,22,40 With regard to this point, geometry optimizations of the NH4+(NH3)4 complex were carried out at the HF, BLYP, B3LYP, and CCSD levels of accuracy using the D95* basis set. According to the optimized NH4+(NH3)4 complex, the observed Hi-N distances of 2.01 (HF), 1.89 (BLYP), 1.89 (B3LYP), and 1.95 Å (CCSD) clearly demonstrate the overestimation of hydrogen-bond strength by the DFT methods. Figure 5 shows the probability distributions of the N-Ni-N angle, calculated from the subset of configurations within the first minimum of the Ni-N RDFs. The B3LYP/MM simulation clearly predicts a tetrahedral cage structure of the solvated NH4+, as indicated by the pronounced large peaks between 80 and

Figure 6. Distributions of the angle R, calculated within the first minimum of the Hi-N RDFs.

130°. In the HF/MM simulation, a slight distortion from the ideal tetrahedral arrangement can be understood due to the interference between the first four ammonia molecules that are directly hydrogen bonded to NH4+ and the other nearestneighbor (i.e., the fifth) ammonia molecules positioned close to the ion. As a consequence, a discrepancy between the structural data obtained from the HF/MM and B3LYP/MM methods could be expected to reflect in different dynamics of the solvated NH4+. Further information on the formation of hydrogen bonds between NH4+ and ammonia can be obtained from the distributions of the angle R, calculated within the first minimum of the Hi-N RDFs, as shown in Figure 6. In both HF/MM and B3LYP/MM simulations, it is observed that the Ni-Hi-N hydrogen bonds are nearly linear, with peak maxima around 160-165°. Figure 7 shows the dipole-orientated arrangement of nearest-neighbor ammonia molecules toward NH4+. The orientations of ammonia molecules are described in terms of the distributions of angle θ, as defined by the N-Hi axis and the dipole vector of the ammonia molecule. On the basis of both HF/MM and B3LYP/MM simulations, a clear dipole-oriented arrangement of the first-shell ammonia molecules toward the ion is observed, as indicated by the pronounced peaks around 150-170°. By the QM/MM MD scheme, since the charges of all particles located within the QM region are varied dynamically during the simulation process, the polarization effects, at least within the solvation shell of NH4+, are automatically included. In both HF/MM and B3LYP/MM simulations, it is observed that firstshell ammonia molecules are slightly polarized by NH4+, leading to an increase of NBO charges of the atoms of ammonia molecules by about 6-8% (i.e., N atoms become more negative

Structure and Dynamics of NH4+ in Liquid Ammonia

J. Phys. Chem. B, Vol. 112, No. 3, 2008 889

Figure 7. Distributions of the angle θ between the N-Hi vector and the vector of ammonia’s dipole moment, calculated within the first minimum of the Hi-N RDFs.

and H atoms more positive). In addition to the correct description on the polarization effects, QM calculations can be used to obtain an estimate of charge transfer between solute and solvent. This makes the QM treatment even more significant, since the transfer of charge between the solvent molecules and the ion is reliably corrected according to the structural changes of the solvated ion. In the HF/MM and B3LYP/MM simulations, the average NBO charges on NH4+ are +0.871 and +0.853, respectively. 3.2. Dynamical Data. 3.2.1. Intramolecular Geometry. The geometrical arrangements of NH4+ and of ammonia molecules in the vicinity of the ion are described in terms of distributions of Ni-Hi and N-H bond lengths, together with the Hi-Ni-Hi and H-N-H angles, as shown in Figures 8 and 9, respectively. In Figure 9, the details with respect to the intramolecular geometry of liquid ammonia obtained from a compatible QM/ MM simulation24 were also plotted for comparison. According to Figure 8, both HF/MM and B3LYP/MM simulations reveal a substantial change in the local structure of NH4+, with halfheight widths of about 1.02 ( 0.04 and 1.05 ( 0.05 Å for the distributions of Ni-Hi bonds and of about 108 ( 4 and 108 ( 7° for the distributions of Hi-Ni-Hi angles, respectively. In the B3LYP/MM simulation, the longer Ni-Hi bond length and broader Hi-Ni-Hi angle, when compared with the HF/MM results, can be regarded as a consequence of the observed stronger ion-ammonia hydrogen-bonding interactions (cf. Figures 1 and 4). For the geometrical arrangements of firstshell ammonia molecules (Figure 9), both HF/MM and B3LYP/ MM simulations show a similar tendency as in the case of NH4+, with half-height widths of about 1.00 ( 0.05 and 1.02 ( 0.06 Å for the distributions of N-H bonds and about 111.5 ( 6 and 111 ( 6° for the distributions of H-N-H angles, respectively. In comparison with the data for liquid ammonia, the HF/MM and B3LYP/MM results clearly indicate a significant influence of NH4+ on the local structure of its surrounding ammonia molecules. 3.2.2. Translational Motions and Exchange Processes of Ammonia Molecules in the SolVation Shell of NH4+. The selfdiffusion coefficients (D) for NH4+ and for ammonia molecules in the bulk and in the solvation shell of the ion are calculated from the center-of-mass velocity autocorrelation function (VACF), Cv(t), using the Green-Kubo relation:41

1 D ) lim 3 tf∞

∫0t Cv(t)dt

(5)

Figure 8. Distributions of (a) Ni-Hi bond length and (b) Hi-Ni-Hi angle of NH4+.

All of the calculated D values are summarized in Table 2. In both HF/MM and B3LYP/MM simulations, the D values for ammonia molecules in the vicinity of NH4+ are significantly less than that of liquid ammonia.24 This result clearly demonstrates the “structure-making” ability of the NH4+ ion. Compared with the HF/MM results, the smaller D value in the B3LYP/ MM simulation corresponds to the stronger ion-ammonia hydrogen-bonded interactions demonstrated earlier. According to the data shown in Table 2, the D value for liquid ammonia derived from the HF method is in poorer agreement with experiment than those of the B3LYP/MM and CP-MD results.7 In this context, the observed larger D value in the HF/MM simulation can be attributed to the effects of electron correlation, which are neglected in the HF treatment. With respect to the flexibility of the NH4+ solvate, several ligand exchange processes are observed when the Ni-N distances are plotted during the HF/MM and B3LYP/MM simulations, as shown in Figures 10 and 11, respectively. For better understanding of the ligand exchange processes, the snapshots showing the 4-, 5-, and 6-fold coordinated complexes are also given as an inset. In order to estimate the rate of ligand exchange processes at NH4+, the mean residence times (MRT) of the ligands were evaluated using the direct method43 with t* values of 0.0 and 0.5 ps. The time parameter t* has been defined as the minimum duration of the ligand’s displacement from its original coordination shell. In general, the MRT data obtained with t* ) 0.0 ps are often used as an estimation of hydrogen-

890 J. Phys. Chem. B, Vol. 112, No. 3, 2008

Tongraar and Hannongbua

Figure 10. Time dependence of Ni-N distances, selecting only the first 15 ps of the HF/MM simulation.

Figure 11. Time dependence of Ni-N distances, selecting only the first 15 ps of the B3LYP/MM simulation. Figure 9. Distributions of (a) N-H bond length and (b) H-N-H angle of ammonia molecules in the bulk and in the first solvation shell of NH4+.

TABLE 2: Self-Diffusion Coefficients (D) of Ammonia Molecules in the Bulk and in the Solvation Shell of NH4+ self-diffusion coefficient (×

10-5

2 -1

cm s )

method

NH4+

first-shell ammonia

liquid ammonia

HF/MM MD B3LYP/MM MD CP-MD experiment

9.67 7.54 9.9b

9.70 7.97

11.34a 9.96a 8.8-11.0b 8.75c

a

Reference 24. b Reference 7. c Reference 42.

bond lifetimes, whereas the data obtained with t* ) 0.5 ps can be considered as a good measure for exchange processes.43 All the calculated MRT values are summarized in Table 3. In comparison with the MRT data for liquid ammonia, both HF/ MM and B3LYP/MM simulations clearly demonstrate the “structure-making” ability of the NH4+ in liquid ammonia. In comparison with the HF/MM results, the B3LYP/MM simulation predicts a significantly slower exchange rate of ligands at the ion. It is known that the possible weaknesses of the B3LYP scheme could include the incompleteness of the kinetic energy term, the self-interaction error, and the parametrizations of the B3LYP method, which did not contain any hydrogen-bond systems. In this context, it is worth noting that the DFT method for describing such a particular system must be employed with caution and that the HF/MM scheme could be regarded as a

TABLE 3: MRTs of Ammonia Molecules in the Bulk and in the Solvation Shell of NH4+ t* ) 0.0 ps ion/solute

CN

tsim

N0.0 ex

+

5.1 4.3 13.6 13.3

24.0 22.0 18.0 14.0

340 168 744 483

NH4

liq. NH3 a

t* ) 0.5 ps

0.0 τNH 3

N0.5 ex

0.5 τNH 3

method

0.38 0.56 0.33 0.38

67 33 151 112

1.83 2.87 1.62 1.66

HF/MM B3LYP/MM HF/MMa B3LYP/MMa

Reference 24.

more reliable simulation technique to provide a detailed description of this solvated ion. 4. Conclusion In this study, HF/MM and B3LYP/MM MD simulations have been performed to investigate the solvation structure and dynamics of NH4+ in liquid ammonia. On the basis of both HF/MM and B3LYP/MM simulations, a well-defined tetrahedral cage structure of the NH4+ solvate is observed, but the shell is quite flexible. Consequently, several species of the solvated NH4+ appear, including 4- to 7-fold and 3- to 6-fold coordinated complexes for the HF/MM and B3LYP/MM simulations, respectively. With respect to the detailed analysis of the ligand exchange process and the MRTs of ammonia molecules surrounding the ion, the behavior related to the “structuremaking” ability of this ion in liquid ammonia is clearly confirmed. However, the overestimation of hydrogen-bond strength by the B3LYP method is recognizable. Following the

Structure and Dynamics of NH4+ in Liquid Ammonia reported tendency that the DFT methods are often found to overestimate the hydrogen bonding, the present HF/MM data are considered to be more reliable. When computational facilities become more feasible, a better description of the structure and dynamics of the solvated NH4+ can be achieved, i.e., by a substantial increase of the QM size in conjunction with the use of more accurate ab initio correlated methods. Acknowledgment. Financial support for this work by the Thailand Research Fund (TRF) is gratefully acknowledged (Project RTA4980006 and Project BRG4880010). The authors thank the National Electronics and Computer Technology Center (NECTEC) for providing computational facilities. Additional support (for computer equipment) by the Suranaree University of Technology (SUT) is also acknowledged. References and Notes (1) Westaway, K. C.; Bourns, A. N. Can. J. Chem. 1972, 50, 2332. (2) Bunnett, J. F.; Connor, D. S.; O’Reilley, K. J. J. Org. Chem. 1979, 44, 4197. (3) Rykowski, A. Pol. J. Chem. 1983, 57, 467. (4) Zoltewicz, J. A.; Oestrich, T. M. J. Org. Chem. 1991, 56, 2805. (5) Tang, I. N.; Castleman, A. W. J. Chem. Phys. 1975, 62, 4576. (6) Meot-Ner, M.; Speller, C. V. J. Phys. Chem. 1986, 90, 6616. (7) Liu, Y.; Tuckerman, M. E. J. Phys. Chem. B 2001, 105, 6598. (8) Car, R.; Parrinello, M. Phys. ReV. Lett. 1985, 55, 2471. (9) Diraison, M.; Martyna, G. J.; Tuckerman, M. E. J. Chem. Phys. 1999, 111, 1096. (10) Kerdcharoen, T.; Liedl, K. R.; Rode, B. M. Chem. Phys. 1996, 211, 313. (11) Tongraar, A.; Liedl, K. R.; Rode, B. M. J. Phys. Chem. A 1997, 101, 6299. (12) Tongraar, A.; Liedl, K. R.; Rode, B. M. J. Phys. Chem. A 1998, 102, 10340. (13) Schwenk, C. F.; Loeffler, H. H.; Rode, B. M. J. Chem. Phys. 2001, 115, 10808. (14) Tongraar, A.; Sagarik, K.; Rode, B. M. J. Phys. Chem. B 2001, 105, 10559. (15) Tongraar, A.; Sagarik, K.; Rode, B. M. Phys. Chem. Chem. Phys. 2002, 4, 628. (16) Tongraar, A.; Rode, B. M. Phys. Chem. Chem. Phys. 2003, 5, 357. (17) Rode, B. M.; Schwenk. C. F.; Tongraar, A. J. Mol. Liq. 2004, 110, 105. (18) Armunanto, R.; Schwenk, C. F.; Tran, H. T.; Rode, B. M. J. Am. Chem. Soc. 2004, 126, 2582. (19) Intharathep, P.; Tongraar, A.; Sagarik, K. J. Comput. Chem. 2005, 26, 1329. (20) Intharathep, P.; Tongraar, A.; Sagarik, K. J. Comput. Chem. 2006, 27, 1723. (21) Tongraar, A.; Tangkawanwanit, P.; Rode, B. M. J. Phys. Chem. A 2006, 110, 12918.

J. Phys. Chem. B, Vol. 112, No. 3, 2008 891 (22) Xenides, D.; Randolf, B. R.; Rode, B. M. J. Chem. Phys. 2005, 122, 174506. (23) Todorova, T.; Seitsonen, A. P.; Hutter, J.; Kuo, I. W.; Mundy, C. J. J. Phys. Chem. B 2006, 110, 3685. (24) Tongraar, A.; Kerdcharoen, T.; Hannongbua, S. J. Phys. Chem. A 2006, 110, 4924. (25) Boese, A. D.; Chandra, A.; Martin, J. M. L.; Marx, D. J. Chem. Phys. 2003, 119, 5965. (26) Dunning, T. H., Jr.; Hay, P. J. Modern Theoretical Chemistry; H. F. Schaefer, Ed.; Plenum III: New York, 1976. (27) Hannongbua, S. J. Chem. Phys. 2000, 113, 4707. (28) Hannongbua, S. Aust. J. Chem. 1991, 44, 447. (29) Dunning, T. H. J. Chem. Phys. 1989, 90, 1007. (30) Kendall, R. A.; Dunning, T. H.; Harrison, R. J. J. Chem. Phys. 1992, 96, 6769. (31) Woon, D. E.; Dunning, T. H. J. Chem. Phys. 1993, 98, 1358. (32) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Peterson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision D.1; Gaussian, Inc.: Wallingford, CT, 2005. (33) Carpenter, J. E.; Weinhold, F. J. Mol. Struct. (THEOCHEM) 1988, 169, 41. (34) Reed, A. E.; Weinstock, R. B.; Weinhold, F. J. Chem. Phys. 1985, 83, 735. (35) Reed, A. E.; Curtiss, L. A.; Weinhold, F. Chem. ReV. 1988, 88, 899. (36) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187. (37) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Phys. Chem. 1984, 81, 3684. (38) Adams, D. J.; Adams, E. H.; Hills, G. J. Mol. Phys. 1979, 38, 387. (39) Ricci, M. A.; Nardone, M.; Ricci, F. P.; Andreani, C.; Soper, A. K. J. Chem. Phys. 1995, 102, 7650. (40) Rode, B. M.; Schwenk, C. F.; Hofer, T. S.; Randolf, B. R. Coord. Chem. ReV. 2005, 249, 2993. (41) McQuarrie, D. A. Statistical Mechanics; Harper & Row: New York, 1976. (42) O’Reilly, D. E.; Peterson, E. M.; Scheie, C. E. J. Chem. Phys. 1973, 58, 4072. (43) Hofer, T. S.; Tran, H. T.; Schwenk, C. F.; Rode, B. M. J. Comput. Chem. 2004, 25, 211.