Investigation of the Local Structure in Sub and Supercritical Ammonia

Oct 28, 2010 - MCCCS Towhee: a tool for Monte Carlo molecular simulation. Marcus G. Martin. Molecular Simulation 2013 , 1-11 ...
1 downloads 0 Views 2MB Size
Download PDF
J. Phys. Chem. B 2010, 114, 15003–15010

15003

Investigation of the Local Structure in Sub and Supercritical Ammonia Using the Nearest Neighbor Approach: A Molecular Dynamics Analysis I. Vyalov,† M. Kiselev,‡ T. Tassaing,§ J. C. Soetens,§ and A. Idrissi*,† Laboratoire de Spectrochimie Infrarouge et Raman, Centre d’Etudes et de Recherches Lasers et Applications, UniVersite´ des Sciences et Technologies de Lille, 59655 VilleneuVe d’Ascq Cedex, France, Institute of Solution Chemistry of the RAS, Akademicheskaya st.1, 153045 IVanoVo, Russia, and Institut des Sciences Mole´culaires, UMR UniVersite´ Bordeaux I - CNRS n 5255, 351, Cours de la Libe´ration, 33405 Talence Cedex, France ReceiVed: June 20, 2010; ReVised Manuscript ReceiVed: October 1, 2010

Molecular dynamics simulations of ammonia were performed in the (N,P,T) ensemble along the isobar 135 bar and in the temperature range between 250 and 500 K that encompasses the sub and supercritical states of ammonia. Six simple interaction potential models (Lennard-Jones pair potential between the atomic sites, plus a Coulomb interaction between atomic partial charges) of ammonia reported in the literature were analyzed. Liquid-gas coexistence curve, critical temperature, and structural data (radial distribution functions) have been calculated for all models and compared with the corresponding experimental data. After choosing the appropriate potential model, we have investigated the local structure by analyzing the nearest neighbor radial, mutual orientation, and interaction energy distributions. The change in the local structure was traced back to the change of the nonlinear behavior (which is more pronounced at low temperatures) of the average distance between a reference ammonia molecule and its subsequent nearest neighbor. Our results suggest to use the position of the maximum in the fluctuation of the average distance to define the border of the first solvation shell (particularly at high temperature when the minimum of the radial distribution is not well-defined). Indeed, the effect of the temperature on the position of this maximum shows clearly that the spatial extent of the solvation shell increases with a concomitant decrease of the involved number of ammonia molecules. Furthermore, our results show that the signature of the hydrogen bonding is mainly observed for temperature below 300 K. This signature is quantified by a short distance contribution to the closest radial nearest neighbor distribution, by a strong mutual orientation (defined by the angles between the axis joining the nitrogen atoms and the molecular axes) and by a strong attractive character of the total interaction energy. I. Introduction Ammonia molecule (NH3) has widespread use in many industrial sectors and is one of the most highly produced chemicals in the world.1 Main uses are in industrial refrigeration systems, as well as in fertilizer, explosives, chemicals production,2,3 and more recently in the CO2 capture technology.4,5 By the same token, supercritical ammonia has been shown to provide advanced technological processes (such as separation, chemical reaction, material processing, waste treatment, plastic recycling, and others).6 Although liquid ammonia has been the subject of numerous investigations both using experimental X-rays,7-9 neutron scattering,10-12 infrared,13 Raman scattering,14-16 nonlinear transient spectroscopy17 and theoretical simulations,18-25 a detailed understanding of the molecular structure of supercritical ammonia is still lacking. We then believe that some special features of supercritical ammonia are relevant to the understanding of supercritical phenomena. Indeed, NH3 is a simple molecule where a subtle balance between hydrogen bonding and van der Waals interactions is expected to provide some peculiarities of its structure and dynamics in the supercritical domain. Understanding the local structure is one of the fundamental and essential problems for the successful use of * To whom correspondence should be addressed. E-mail: nacer.idrissi@ univ-lille1.fr. † Universite´ des Sciences et Technologies de Lille. ‡ Institute of Solution Chemistry of the RAS. § UMR Universite´ Bordeaux I - CNRS.

supercritical technology using ammonia.17,26 These considerations prompt us to characterize the change in the local structure of ammonia from the liquid state up to the supercritical fluid by means of molecular dynamics simulations. Indeed, computer simulations provide deep insight on the local structure and a successful use of this technique is based on two steps: the first one consists in the development of accurate potential model where the parameters are adjusted to available experimental data. The second step is to choose the adequate molecular simulation technique to estimate physical properties as efficiently as possible. Using computer simulations, the radial distribution function is a major means by which structure is usually described in fluids systems. However, this statistical property does not provide sufficient details on the nature of the local environment in terms of first-, second- or higher-order neighbors distributed around a central particle in fluids. In fact, even in cases where the dynamics differs by many orders of magnitude, the radial distribution function may reveal only minor differences. Indeed, the changes in the structure may occur at a level of nearest neighbors and averaging through the whole space volume reduces the effect of the change of thermodynamic conditions (pressure, temperature, etc.). Therefore to describe the change in the local structure, it is necessary to introduce statistical functions that are evaluated over small space volume. These statistical functions include the distribution of Voronoi polyhedra,27 the mutual arrangement of Delaunay simplices,28 and the neighborship approach.29-32

10.1021/jp108701t  2010 American Chemical Society Published on Web 10/28/2010

15004

J. Phys. Chem. B, Vol. 114, No. 46, 2010

The purpose of this paper is then to report molecular dynamics simulation investigations on the local structure in sub- and supercritical ammonia using the nearest neighbor approach. We started this work by the inspection of six simple interaction potential models of ammonia reported in the literature. The essential criterion that was adopted in the selection of these models is the simplicity of the expression of the intermolecular interaction potential. They are written in terms of a LennardJones pair potential between atomic sites, plus a Coulomb interaction between distributed partial charges. Liquid gascoexistence curve, critical temperature, densities, energies, and structural data (radial distribution functions) have been calculated for all models, and compared with the corresponding experimental data. These comparisons allow a classification of the different models in terms of their relative agreement with experimental data, particularly those properties of fluid ammonia near the critical point and under the supercritical conditions. After choosing the adequate potential model, we have investigated the local structure by analyzing the nearest neighbor radial distributions. The signature of the hydrogen bonding along the thermodynamic point considered here is discussed as well as the mutual orientation between ammonia molecules. Finally a conclusion is given. II. Simulation Methodology II.1. Gibbs Ensemble Monte Carlo. To evaluate the liquid-gas coexistence curve and critical parameters for the different potential models, we performed Gibbs Ensemble Monte Carlo simulations (GEMC). This technique allows one to obtain the liquid-gas coexistence curve directly from the Monte Carlo run. Practically, two simulation boxes corresponding to the liquid and to the gas phase are considered. Standard Monte Carlo moves that include the center of mass translations, rotation around randomly chosen axes, are supplemented by trial moves that swap molecules between boxes and change volume of boxes. They are performed in such a way that overall volume of the boxes and the total number of molecules are kept constant. After several Monte Carlo cycles, one obtains the density corresponding to the liquid state in one box and the density of the gas in the other one. One can find a comprehensive description of this method in the original paper by Panagiotopoulos33 or in general courses on molecular simulations.34 In this work, Monte Carlo simulations have been performed using the program Towhee.35 In each simulation step, a randomly chosen molecule has been randomly translated by no more than 0.5 Å and has been randomly rotated around a randomly chosenfixed axis by no more than 0.05 radians, the probability of performing of the box interchanging move was chosen equal to 0.01. The edge length of the basic simulation boxes L1 (liquid density) and L2 (gas density) have been set in accordance with the experimental densities of liquid and gaseous ammonia respectively.36 All interactions have been truncated to zero beyond the center-center cutoff distance of 12.5 Å. The longrange part of the electrostatic interaction has been accounted for using Ewald summation technique.37 The system has been equilibrated by performing 50 000 Monte Carlo cycles. The production runs last 50 000 Monte Carlo cycles each. II.2. Molecular Dynamics Simulations. After choosing the adequate potential model, molecular dynamics simulations of the NH3 system has been performed in the isothermal-isobaric (N,P,T) ensemble at 135 bar and for the following temperatures 250, 300, 350, 380, 390, 394, 398, 400, 402, 406, 410, 420, 440, 460, and 500 K. The cubic simulation boxes contain 864 molecules. Standard periodic boundary conditions have been applied. All the models used here are the combination of the

Vyalov et al. simple Lennard-Jones potentials with electrostatic potentials. The total potential thus writes

U)

∑ ∑ [ULJ(riRjβ) + UEL(riRjβ)] ) iR

jβ

∑∑ iR

jβ

[ (( ) ( ) ) 4εiRjβ

σiRjβ riRjβ

12

-

σiRjβ riRjβ

6

+

qiRqjβ riRjβ

]

(1)

where riRjβ is the distance between a site R of molecule i and a site β of molecule j, qiR is the charge on site R of molecule i, qjβ is the charge on site β of molecule j. σ and ε are the LJ parameters. The simulations have been performed using the DL_POLY program.38 The temperature and pressure of the systems have been kept constant by means of the weak coupling algorithms of Berendsen,39 using the coupling parameter values of 0.1 ps (temperature) and 0.5 ps (pressure). All site-site interactions have been truncated to zero beyond the center-center cutoff distance of 12.5 Å. The long-range part of the electrostatic interactions has been accounted for using the Ewald summation method. The equations of motion have been integrated using the leapfrog algorithm, employing an integration time step of 2 fs. The systems have been equilibrated by generating 5 ns long trajectories for each. Then, the production stage lasts 800 ps (400 000 configurations). The coordinates of molecules were stored each 8 fs and were used to calculate the relevant statistical properties. II.3. Models of Ammonia. We have chosen to study six rigid potential models of ammonia proposed in the literature. The essential criterion we have adopted in the selection of the models is the simplicity of the expression for the intermolecular potential. The models we have selected were those written entirely in terms of a Lennard-Jones pair potential and Coulomb interaction between partial charges. Among the models that satisfy this criterion, we have chosen those proposed by Impey et al.40 and Ferrario et al.41 They were developed to predict structural properties of liquid ammonia and are then good for describing of the experimental radial distribution functions and structure factors. At the same time, they fail to predict the thermodynamic parameters of the ammonia closely to the critical point and along liquid-gas coexistence curve (CC) in general. The model proposed by Gao et al.42 was developed from pure quantum chemical considerations to take into account the bimolecular interactions and was reported to well predict the structure of the liquid ammonia and particularly the hydrogen bonds distribution. We have also used other recently developed potential models by Kristof et al.43 which parameters were adjusted to reproduce the liquid-gas coexistence curve by means of the Grand Canonical Ensemble Monte Carlo simulations. However no molecular dynamics simulations using these models were reported. In the remainder of the paper, we use the following acronyms for the various ammonia models: IK for the model proposed by Impey et al.40 FHMK for the model proposed by Ferrario et al,41 GXG for the model proposed by Gao et al42 and KRI1, KRI2, and KRI3 for the potential models proposed by Kristof et al.43 It is worthy to mention that in the five-sites models (IK, FHMK, and KRI2), one additional site is introduced along the molecular axis between the nitrogen atom and the hydrogen atoms plane. The parameters of the six potentials models of ammonia are given in Table 1. II.4. Methods for Studying the Local Structure. The nearest neighbor approach has been already shown to give interesting results in analyzing the local structure of hard sphere

Investigation of the Local Structure in Ammonia

J. Phys. Chem. B, Vol. 114, No. 46, 2010 15005

TABLE 1: Geometrical and Energetical Parameters of the Six Potential Models of Ammoniaa site

ε (kcal)

σ (Å)

q (e)

rNH (Å)

∠HNH

rNM (Å)

∠HNM

µ (D)

acronyms

N M H N M H N H N H N M H N H

0.278 0 0 0.278 0 0 0.210 0 0.348 0.35 0.368 0 0 0.338 0

3.4 0 0 3.4 0 0 3.36 0 3.4 0 3.37 0 0 3.39 0

0 -1.386 0.462 0 -1.455 0.485 -1.026 0.342 -1.050 0.350 -1.275 0.420 0.285 -1.035 0.345

1.0124

106.7

0.156

67.8

1.50

IK40

1.0127

110.96

0.156

67.8

1.50

FHMK41

1.0124

106.7

1.88

GXG42

1.0124

106.7

1.92

KRI162

1.0124

106.7

1.95

KRI262

1.0124

106.7

1.90

KRI362

0.191

67.8

a ∠HNH is the angle between the NH axes of ammonia while ∠HNM is the angle between the NH axis and that joining the N atom to the virtual site M; µ is the dipole moment.

systems,27,29,31 supercooled liquids,30,44 supercritical CO2,45,46 and binary mixtures.47,48 In this approach, the neighbors of a central atom are sorted by distance into the first neighbors, second neighbors, and so forth. Separate radial distribution functions, pRβ(n, r, T) may be defined for each set of nearest neighbors atoms β (indicated by n), and at distance r from the central atom R. T is the temperature. The parameters chosen for the analysis of these functions are the average distance 〈rR-β(n,T)〉 between a reference atom R and an atom β belonging to the nth neighbor distribution pRβ(n, r, T); the corresponding root-meansquare average distance ∆rR-β(n,T). These parameters are defined by the following equations

〈rR-β(n, T)〉 ) 2 〈rR-β (n, T)〉

)

∫0∞ dr'r'pRβ(n,r'T) ∫0

∞

2

dr'r' pRβ(n, r', T)

(2)

(3)

and 2 ∆rR-β(n, T) ) (〈rR-β (n, T)〉 - 〈rR-β(n, T)〉2)1/2

Figure 1. Pictorial representation of the orientation between a reference ammonia molecule and its first nearest neighbor. The angles θΑ is that between the molecular axis and that joining the nitrogen atoms of the two molecules. The angles θΒ is that between the neighbor molecular axis and that joining the nitrogen atoms of the two molecules.

(4)

The root-mean-square average distance, ∆rR-β(n,T), characterizes the width of the distribution; a large value of ∆rR-β(n,T) reflects in an average sense, a fluctuating character of the corresponding nearest neighbor. The behavior of these parameters as a function of the temperature allows one to quantify its effect on the radial distribution of each neighboring distribution. In the same manner, we can also have access to the total interaction energy distribution PUTot(n,T) between a reference molecule and its nearest neighbors. The electrostatic (PUEl(n,T)) and Lennard-Jones (PULJ(n,T)) interactions energies were also calculated to get more insight into the contribution of the different component of the intermolecular interaction to the resultant local structure. Information on the orientation between a probe molecule and another belonging to a class of neighboring molecules, have been obtained by the calculation of the orientation distribution, P(θA,θB,n,T) as follows: for a distance r between the central atom R and an atom β belonging to a class n of neighboring atoms, we calculated the distribution of the angles θA and θB. These angles are defined as follows: for a two nearest neighbor ammonia molecules A and B, θA and θB are the angles between the axis joining the two nitrogen

atoms and the molecular axis of molecules A and B, respectively (see Figure 1). III. Results and Discussion III.1. Structure. There are several experimental data available from X-ray diffraction7 and neutron scattering49 on ammonia and it is possible to establish the validity of a chosen potential model by comparing the calculated radial distribution function (rdf) and the experimental one. We compare in Figure 2a-c the calculated partial rdfs gNN(r), gNH(r), and gHH(r) with experimental ones measured at 273 K. It can be seen that all models have reasonable agreement with experiment. The shoulder in the gNN(r) at 2.7 Å is not reproduced by any model considered in the present study. But to our knowledge, no one could reproduce this feature even using such sophisticated methods as DFT, QM/MM, and path-integral molecular dynamics.19,50 One should also mention that there is no contribution at this distance in the obtained gNN(r) from the X-ray experiment. The first solvation shell defined by the first minimum of gNN(r) contains 12-13 ammonia molecules. The peaks at distance less than 2 Å in the experimental NH and HH rdfs, correspond to the intramolecular contribution and of course cannot be seen from the MD simulations with rigid potential model. gNH(r) rdf has a noticeable shoulder at 2.25 Å that was considered to be a fingerprint of the hydrogen bonding in liquid ammonia.49 This peak is more resolved for IK and FHMK potential models. Integration of this shoulder gives less than one hydrogen bond per one ammonia molecule.49 III.2. Liquid-Gas Coexistence Curve. To qualify one of the potential models for further analysis, we have also compared

15006

J. Phys. Chem. B, Vol. 114, No. 46, 2010

Vyalov et al.

Fl - Fg ) B(T - Tc)β

(6)

In these equations, A and B are constants, β is the critical exponent, and is equal to 0.32, Fl and Fg are the liquid and gas densities, respectively. The comparison between the experimental and calculated liquid-gas CC is drawn in Figure 3. It is obvious to notice that the KRI1, KRI2, KRI3, and IK models reproduce quite well the highest density side of the liquid state as well as the lowest density of the gas state. The FHMK and GXG potential models fail to reproduce the density of the liquid state of ammonia. As the temperature increases along the liquid-gas CC, the KRI1, KRI2, and KRI3 models give density values which are close to the experimental ones. Furthermore, the calculated critical temperature and the critical density are compared with the corresponding experimental values in Table 2. As one can see the IK, GXG, and FHMK models predict critical temperatures that are far from the experimental one compared with those predicted by KRI1, KRI2, and KRI3 potential models. A good agreement was found between the experimental and calculated values using the KRI3 for quantities such as density,36 thermal expansion36 and diffusion coefficient.51 We then chose to carry the analysis of the local structure of ammonia along the isobar 135 bar using the KRI3 potential model. III.3. Nearest Neighbor Radial Distributions. First we start by analyzing the connection between the gNN(r) and gNH(r) rdfs and the corresponding nearest neighbor radial distributions. The nearest neighbor radial distributions pNN(n, r, T ) 250 K) and pNH(n, r, T ) 250 K) are shown in Figure 4, panels a and b, respectively. One can see here that up to the first minimum in the NN and NH rdfs, these functions are well reproduced by the sum of the nearest neighbor radial distributions. The advantage of this approach can already be shown here. Indeed, as shown in Figure 5, the shoulder in the gNH(r) rdf at 2.25 Å is associated with the radial distribution contribution of the first nearest hydrogen atom to the nitrogen atom. Furthermore, the pNH(n ) 1, r, T ) 250 K) distribution may be resolved in two contributions at distances d1 ) 2.40 Å and d2 ) 2.91 Å. Using DFT calculations implemented in Spartan package,52 we opti-

Figure 2. Comparison between the calculated radial distribution functions (top) gNN(r), (middle) gNH(r), and (bottom) gHH(r) with experimental ones measured at 273 K.49

the calculated and experimental liquid-gas coexistence curves (CC). We have then performed Monte Carlo simulations in the Grand Canonical Ensemble using the Towhee program package.35 To calculate the critical density, Fc, and the critical temperature, Tc, the obtained T-F coexistence curves were fitted to the law of rectilinear diameters (see eq 5) and to the scaling law (see eq 6)

Fl + Fg ) Fc + A(T - Tc) 2

Figure 3. Comparison between the experimental36 and the calculated liquid-gas coexistence curves for the six potential models used in our simulations.

TABLE 2: Comparison between the Calculated Critical Temperature, Tc and Critical Density, Gc and the Corresponding Experimental Values IK

(5)

FHMK GXG KRI1 KRI2 KRI3 experiment36

Tc (K) 347 375 346 390 399 400 Fc (g cm-3) 0.225 0.240 0.235 0.225 0.236 0.218

405.4 0.225

Investigation of the Local Structure in Ammonia

J. Phys. Chem. B, Vol. 114, No. 46, 2010 15007

Figure 6. Average values of the interaction energy between a reference ammonia molecule and its first nearest neighbor, as a function of temperature.

Figure 4. Nearest neighbor radial distributions for n ) 1-54 (top) gNN(r), and pNN(n, r, T ) 250 K) and (bottom) gNH(r) and pNH(n, r, T 54 54 ) 250 K). The sum ∑i)1 pNN(n,r,T) and ∑i)1 pNH(n,r,T) are also presented.

Figure 5. First nearest neighbor radial distributions pNH(n ) 1, r, T) calculated along the isobar 50 bar and at various temperatures in the range between 259 and 500 K. The arrows indicate the position of the main contributions at distance d1 and d2 (see text).

mized the dimer geometry of two ammonia molecules under the constraints that the intermolecular distance between the nitrogen atoms of the two molecules is equal to 〈rN-N(n ) 1, T ) 250 K)〉 ) 3.28 Å and that the intermolecular distance between nitrogen and one of the hydrogen atoms is equal to d1. The optimized geometry clearly shows that the intermolecular

distance between the nitrogen and another hydrogen is equal to 2.92 Å, which may be associated with the contribution at d2 in pNH(n ) 1, r, T ) 250 K). Figure 5 shows that as the temperature increases, the contribution at short distance, d1 is drastically reduced indicating a weakening of hydrogen bond and the pNH(n ) 1, r, T) distributions are then characterized by the occurrence of a large peak centered at larger distance (higher than 3 Å). If we rely on the definition that the hydrogen bonding between an hydrogen atom and X atom is considered as weak when the H · · · X distance is typically between 2 and 3 Å,53 our results indicate that at higher temperature than 300 K the hydrogen bonding in ammonia may be considered as weak. To reinforce this conclusion, we have calculated the total interaction energy distributions, pUTot(n ) 1, T) as well as the electrostatic pUEl(n ) 1, T) and Lennard-Jones pULJ(n ) 1, T) contributions between a reference ammonia molecule and its first nearest neighbor. The corresponding average values are given in Figure 6. These results show that the weakening of the hydrogen bonding may be correlated with the balance between the rate of decrease of both the repulsive character of the LJ interactions and the attractive character of the electrostatic interactions. These results also suggest that to analyze the signature of hydrogen bonding in ammonia using techniques (such as vibration spectroscopy, NMR, etc.) that probe the local structure one should carry the experiments in the range of temperature below 300 K. At low temperature, the well-defined first minimum of the gNN(r) rdf defines the spatial extension of the first solvation shell and the integration of gNN(r) at this position gives the number of molecules (coordination number) in the first solvation shell. As the temperature increases, the position of the first minimum is not clearly defined (see Figure 7). The nearest neighbor approach allows to overcome this difficulty. Indeed the superimposition of the gNN(r) and ∆rN-N(n,T ) 250 K) (as a function of 〈rN-N(n, T ) 250 K)〉) shows that the maximum of ∆rN-N(n, T ) 250 K) occurs in the region of the first minimum of gNN(r) (see Figure 8). It is very interesting to note that such behavior was observed for other liquid systems such as water, carbon dioxide, acetone.45-47 Therefore, the position of the maximum of ∆rN-N(n, T ) 250 K) provides a useful technique to find the border of the first solvation shell in fluids. We have plotted in Figure 9 the behavior of the average distance to each neighbor 〈rN-N(n,T)〉 and in Figure 10 the fluctuations of distance that are available for this neighbor

15008

J. Phys. Chem. B, Vol. 114, No. 46, 2010

Vyalov et al.

Figure 10. Fluctuation of the average distance ∆rR-β(n,T) along the isobar 50 bar and in the range of temperature between 250 and 500 K.

Figure 7. Radial distribution functions gNN(r) and coordination number (right axis) at various temperatures along the isobar 135 bar.

Figure 8. Radial distribution function gNN(r) and the fluctuations of distance of the nearest neighbors ∆rN-N(n, T ) 250 K).

Figure 9. Average distance, 〈rN-N(n,T)〉, between a reference ammonia molecule and its nearest neighbors, n, along the isobar 50 bar and in the range of temperature between 250 and 500 K. The derivative of 〈rN-N(n,T)〉 with respect to n is given in the inset.

∆rN-N(n,T). It can be seen that the increase of the average distance has a nonlinear character which is more pronounced at low temperatures. This feature is clearly seen in the behavior

of the derivative of 〈rN-N(n,T)〉 with respect to n (see the inset of Figure 9). For instance, the derivative of 〈rN-N(n,T)〉 has a maximum for the 11th nearest neighbor and the minimum for the third nearest neighbor. This behavior reflects the fact that the rate of increase of the average nearest neighbor distance with respect to the number of the neighbor depends on the local number density of the molecules in this region. For small n values (close nearest neighbors), the decrease of this rate (particularly the position of its minimum) is associated with high density regions and then with the dominance of the attractive component of the total interaction energy and the increase of this rate indicates the opposite situation. As the temperature increases along the isobar 135 bar, the minimum in the rate of increase of 〈rN-N(n,T)〉 disappears completely near the critical temperature (see the inset of Figure 9) indicating that the total interaction energy becomes less attractive. Furthermore, for higher temperatures, the rate of increase has a plateau shape indicating that these nearest neighbors are experiencing a homogeneous density distribution. At large n values, the rate of increase is almost the same for all the temperatures suggesting that a reference ammonia molecule sees these neighboring molecules randomly distributed. An important role in the description of the nearest neighborhood of the reference molecule is given by the fluctuations of the average distances that describe the deviations of the distance between the reference molecule and its nearest neighbor from its average value. The effect of temperature on this parameter is given in Figure 10. The position of the maximum which defines the spatial extent of the first solvation shell shifts to occur for close nearest neighbors. The number of molecules in the first solvation shell then decreases from 11 to 7 near the critical temperature and down to only 5 for the highest temperatures considered here. The spatial extent of the first solvation increases from 4.7 to 6.5 Å. The large gap in the behavior of ∆rR-β(n,T) near the critical temperature is associated with the maximum of the thermal expansion. III.4. Nearest Neighbor Orientational Distributions. Previous findings42 emphasize the importance of a detailed description of mutual orientation of the dimer of ammonia molecules to check the validity of the intermolecular potential model. We have used the nearest neighbor approach to analyze the mutual orientation distribution P(θA,θB,n,T) between a reference ammonia molecule and its nearest neighbor. The angles θA and θB are defined in Figure 1. These distributions are depicted in Figure 11 for T ) 250, 380, 410 and 500 K. The peaks positions of P(θA,θB, n ) 1, T ) 250 K) are at the values of θA ) 40, 115° and θB ) 115, 40° (the average distance 〈rN-N(n ) 1, T ) 250 K)〉 ) 3.26 Å). The equilibrium geometry of the ammonia dimer has been the subject of a considerable body of experimental and theoretical studies and is then a subject of controversy. Indeed, using

Investigation of the Local Structure in Ammonia

J. Phys. Chem. B, Vol. 114, No. 46, 2010 15009

Figure 11. Nearest neighbor relative orientational distribution. The angles θA and θB are defined in Figure 1. For comparison purposes the color scaling has been preserved.

electric deflection method, Odutola et al.54 found the values θA ) 50.6° and θB ) 66.5°; using high resolution microwave spectroscopy, Nelson et al55 found the average values of θA ) 48.6° and θB ) 115.5° for a distance between nitrogen atoms of two ammonia molecules rNN equal to 3.34 Å. Using quantum calculations, Olthof et al.56 found the values θA ) 45° and θB ) 95°; Liu et al.57 found the values θA ) 13° and θB ) 102° and rNN ) 3.34 Å; Lee et al.58 found the values θA ) 36° and θB ) 95° and rNN ) 3.24 Å. Beu et al.59 report the values θA ) 27° and θB ) 95° for a distance rNN ) 3.26 Å. Other values are compiled in the paper by Nelson et al.55,60 Although the comparison between our results and those reported in the literature cited before is not easy since in our case the mutual orientation is environment dependent (many body interactions), our results seem to correlate with the experimental data reported by Nelson et al.55,60 The effect of temperature on the mutual orientation between a reference ammonia molecule and the first nearest neighbor is illustrated in Figure 11. The behavior of P(θA,θB,n ) 1,T) clearly indicates that the mutual orientation defined by the angles θA ) 40, 115° and θB ) 115, 40° remains the most probable although the associated probability decreases with increasing the temperature. IV. Conclusion In this paper, we report the results of a number of molecular dynamics and Monte Carlo simulations of sub and supercritical ammonia system along the isobar 135 bar. In the first step, we promote one of six potential models proposed in the literature on the basis of a favorable comparison of the calculated and experimental data such as the liquid-gas coexistence curve, the partial structure factors and the density. The analysis of the local structure along the isobar 135 bar was carried out using the nearest neighbor approach. The change of the local structure was traced back to the change in the behavior of both the

average distance between a reference ammonia molecule and its subsequent nearest neighbors, and the corresponding fluctuation. Our results show that the average distance has a nonlinear behavior with an occurrence of a minimum and maximum in the rate of increase of the average distance. The minimum coincides with the position of the third nearest neighbor and was interpreted as indication of high density region and a dominance of attractive character of the interaction potential energy. The maximum coincides with the position of the eleventh nearest neighbor that was associated with a low density region and a reduced contribution of the attractive contribution to the total interaction energy. Furthermore, we suggest to use the position of the maximum in the fluctuation of the average distance to define the spatial extent of the first solvation shell, showing then that this spatial extent increases and that the number of the involved ammonia molecules is reduced. At low temperature (far from the critical one), hydrogen bonding was inferred from the occurrence of short distance peak in the first nearest neighbor radial distribution between nitrogen and hydrogen atoms. The position of this peak coincides with that observed in the gNH(r) radial distribution function. As this peak disappears with increasing the temperature (higher than 300 K) in favor of high distance peak, this suggests that the hydrogen bonding in ammonia becomes weak. This result is in accordance with Buback conjecture that hydrogen bond in ammonia is weak at and above room temperature and, in particular, under supercritical conditions.14,61 Furthermore our results provide unambiguous information on the change of the mutual orientation between a reference ammonia molecule and its first nearest neighbor, upon the change in the thermodynamic points considered in this work. Indeed, our results clearly show that the angles between each molecular axis and the axis joining the two nitrogen atoms are equal respectively to 40 and 115° for a distance of 3.26 Å between the two nitrogen atoms. When increasing the temper-

15010

J. Phys. Chem. B, Vol. 114, No. 46, 2010

ature, this mutual orientation remains the most probable, however, the associated probability decreases drastically. Acknowledgment. The authors would like to thank C. Gusinde for reading the manuscript. The Institut du De´veloppement et des Ressources en Informatique Scientifique (IDRIS), the Centre de Ressources informatiques (CRI) de l’universite´ de Lille, and the centre de ressource Informatique de HauteNormandie (CRIHAN) are thankfully acknowledged for the CPU time allocation. The Centre d’Etudes et de Recherches Lasers et Applications is supported by the Ministe`re Charge´ de La Recherche, the Re´gion Nord/Pas de Calais and Les Fonds Europe´en de De´veloppement Economique des Re´gions. References and Notes (1) Appl, M. Ammonia. In Ullmann’s Encyclopedia of Industrial Chemistry; Wiley-VCH: Weinheim, 2006. (2) Wada, M.; Heux, L.; Isogai, A.; Nishiyama, Y.; Chanzy, H.; Sugiyama, J. Macromolecules 2001, 34, 1237. (3) Desmoulins-Krawiec, S.; Aymonier, C.; Loppinet-Serani, A.; Weill, F.; Gorsse, S.; Etourneau, J.; Cansell, F. J. Mater. Chem. 2004, 14, 228. (4) Darde, V.; Thomsen, K.; van Well, W. J. M.; Stenby, E. H. Energy Procedia 2009, 1, 1035. (5) Martin, A.; Kant, M.; Jackstell, R.; Klein, H.; Beller, M. Chem. Ing. Tech. 2007, 79, 891. (6) Darr, J. A.; Poliakoff, M. Chem. ReV. 1999, 99, 495. (7) Narten, A. H. J. Chem. Phys. 1977, 66, 3117. (8) Sette, F.; Ruocco, G.; Cunsolo, A.; Masciovecchio, C.; Monaco, G.; Verbeni, R. Phys. ReV. Lett. 2000, 84, 4136. (9) Giura, P.; Angelini, R.; Datchi, F.; Ruocco, G.; Sette, F. J. Chem. Phys. 2007, 127, 084508. (10) Bausenwein, T.; Bertagnolli, H.; David, A.; Goller, K.; Zweier, H.; Todheide, K.; Chieux, P. J. Chem. Phys. 1994, 101, 672. (11) Chieux, P.; Bertagnolli, H. J. Phys. Chem. 1984, 88, 3726. (12) Thompson, H.; Wasse, J. C.; Skipper, N. T.; Hayama, S.; Bowron, D. T.; Soper, A. K. J. Am. Chem. Soc. 2003, 125, 2572. (13) Buback, M.; Franck, E. U. J. Chim. Phys. 1975, 72, 601. (14) Buback, M.; Schulz, K. R. J. Phys. Chem. 1976, 80, 2478. (15) Schwartz, M.; Wang, C. H. J. Chem. Phys. 1973, 59, 5258. (16) Bradley, M.; Zerda, T. W.; Jonas, J. J. Chem. Phys. 1985, 82, 4007. (17) Schafer, T.; Schwarzer, D.; Lindner, J.; Vohringer, P. J. Chem. Phys. 2008, 128, 064502. (18) Nelson, J. D. D.; Fraser, G. T.; Klemperer, W. J. Chem. Phys. 1985, 83, 6201. (19) Boese, A. D.; Chandra, A.; Martin, J. M. L.; Marx, D. J. Chem. Phys. 2003, 119, 5965. (20) Tongraar, A.; Kerdcharoen, T.; Hannongbua, S. J. Phys. Chem. A 2006, 110, 4924. (21) Kiselev, M.; Kerdcharoen, T.; Hannongbua, S.; Heinzinger, K. Chem. Phys. Lett. 2000, 327, 425. (22) Janeiro-Barral, P. E.; Mella, M. J. Phys. Chem. A 2006, 110, 11244. (23) Spirko, V. J. Mol. Spectrosc. 1983, 101, 30. (24) Spirko, V.; Kraemer, W. P. J. Mol. Spectrosc. 1989, 133, 331. (25) Rajamaki, T.; Kallay, M. l.; Noga, J.; Valiron, P.; Halonen, L. Mol. Phys. 2004, 102, 2297. (26) Kiselev, M.; Kerdcharoen, T.; Hannongbua, S.; Heinzinger, K. Chem. Phys. Lett. 2000, 327.

Vyalov et al. (27) Lavrik, N. L.; Voloshin, V. P. J. Chem. Phys. 2001, 114, 9489. (28) Medvedev, N. N.; Geiger, A.; Brostow, W. J. Chem. Phys. 1990, 93, 8337. (29) Mazur, S. J. Chem. Phys. 1992, 97, 9276. (30) Keyes, T. J. Chem. Phys. 1999, 110, 1097. (31) Bhattacharjee, B. Phys. ReV. E 2003, 67, 041208. (32) Simon, S. H.; Dobrosavljevic, V.; Stratt, R. M. J. Chem. Phys. 1991, 94, 7360. (33) Panagiotopoulos, A. Z. Mol. Phys. 1987, 61, 813. (34) Frenkel, D.; Smith, B. Understanding Molecular simulation; Academic Press, Inc.: Orlando, 2001. (35) MCCCS Towhee; http://towhee.sourceforge.net (accessed 2010). (36) Lemmon, W.; Mclinden; M. O. Friend, D. G. Thermophysical Properties of Fluid Systems; National Institute of Standards and Technology: Gaithersburg, MD, 2005; p 20899 (http://webbook.nist.gov). (37) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, 1987. (38) The DL_POLY Molecular Simulation Package; www.ccp5.ac.uk/ DL_POLY/ for serial and parallel molecular dynamics simulations (accessed 2010). (39) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (40) Impey, R. W.; Klein, L. M. Chem. Phys. Lett. 1984, 104, 579. (41) Ferrario, M.; Haughney, M.; McDonald, I. R.; Klein, M. L. J. Chem. Phys. 1990, 93, 5156. (42) Gao, J.; Xia, X.; George, T. F. J. Phys. Chem. 1993, 97, 9241. (43) Kristof, T.; Vorholz, J.; Liszi, J.; Rumpf, B.; Maurer, G. Mol. Phys. 1999, 97, 1129. (44) Saitta, A. M.; Strassle, T.; Rousse, G.; Hamel, G.; Klotz, S.; Nelmes, R. J.; Loveday, J. S. J. Chem. Phys. 2004, 121, 8430. (45) Idrissi, A.; Damay, P.; Kiselev, M. Chem. Phys. 2007, 332, 139. (46) Idrissi, A.; Vyalov, I.; Damay, P.; Frolov, A.; Oparin, R.; Kiselev, M. J. Phys. Chem. B 2009, 113, 15820. (47) Idrissi, A.; Gerard, M.; Damay, P.; Kiselev, M.; Puhovsky, Y.; Cinar, E.; Lagant, P.; Vergoten, G. J. Phys. Chem. B 2010, 114, 4731. (48) Sur, P.; Bhattacharjee, B. J. Chem. Sci. 2009, 121, 929. (49) Ricci, M. A.; Nardone, M.; Ricci, F. P.; Andreani, C.; Soper, A. K. J. Chem. Phys. 1995, 102, 7650. (50) Diraison, M.; Martyna, G. J.; Tuckerman, M. E. J. Chem. Phys. 1999, 111, 1096. (51) O’Reilly, D. E.; Peterson, E. M.; Scheie, C. E. J. Chem. Phys. 1973, 58, 4072. (52) Wavefunction, Inc.; http://www.wavefun.com (accessed 2010). (53) Jeffrey, G. A. An Introduction to Hydrogen Bonding; Oxford University: New York, 1997. (54) Odutola, J. A.; Dyke, T. R.; Howard, B. J.; Muenter, J. S. J. Chem. Phys. 1979, 70, 4884. (55) Nelson, J. D. D.; Fraser, G. T.; Klemperer, W. J. Chem. Phys. 1985, 83, 6201. (56) Olthof, E. H. T.; van der Avoird, A.; Wormer, P. E. S. J. Chem. Phys. 1994, 101, 8430. (57) Liu, S.-Y.; Dykstra, C. E.; Kolenbrander, K.; Lisy, J. M. J. Chem. Phys. 1986, 85, 2077. (58) Lee, J. S.; Park, S. Y. J. Chem. Phys. 2000, 112, 230. (59) Beu, T. A.; Buck, U. J. Chem. Phys. 2001, 114, 7848. (60) Nelson, D. D.; Fraser, G. T.; Klemperer, W. Science 1987, 238, 1670. (61) Buback, M. Ber. Bunsen-Ges. Phys. Chem. 1974, 78, 1230. (62) Kristof, T.; Vorholz, J.; Liszi, J.; Rumpf, B.; Maurer, G. Mol. Phys. 1999, 97, 1129.

JP108701T