Reverse Monte Carlo Simulation of a Heteronuclear Molecular Liquid

Lívia B. Pártay , Pál Jedlovszky and George Horvai. The Journal of Physical Chemistry C 2009 ... John S. Loring and W. Ronald Fawcett. The Journal ...
0 downloads 0 Views 1MB Size
Download PDF
5994

J . Phys. Chem. 1994, 98, 5994-6002

Reverse Monte Carlo Simulation of a Heteronuclear Molecular Liquid: Structural Study of Acetonitrile Tamas Radnai’ and PA1 Jedlovszky Central Research Institute for Chemistry of the Hungarian Academy of Sciences, Budapest, P. 0. Box 17, H-1525 Hungary Received: November 9, 1993; In Final Form: March 2, 1994’

The liquid structure of acetonitrile (AN) has been studied by the recently developed reverse Monte Carlo ( R M C ) simulation technique. This is the first attempt to apply the method to a heteroatomic molecular liquid. The basic simulation box contained 512 AN molecules. Each molecule was represented by three atomic sites, which were held together by “coordination constraints” allowing a realistic flexibility of the molecule. The generated configurations have been refined against experimental X-ray diffraction data. A Metropolis Monte Carlo simulation was also carried out and proved the reliability of the results of R M C . Partial pair correlation functions and angular and spatial distributions of the neighbors around a central molecule were calculated and analyzed in terms of predominant orientational configurations. A strong preference for antiparallel and a slight preference for parallel, head-to-tail, T- or L-shaped configurations were found, in qualitative accordance with previous results. These preferences however, decay rapidly after the first two to three neighbors. The average coordination number of the first shell resulted in five to six neighbors; moreover, different “orientational states” could be distinguished within the first coordination shell of molecules.

Introduction Acetonitrile (AN) has been subjected to several experimental studies aiming at the derivation of its intra- and intermolecular structures in the liquid state. The intramolecular structure of the A N molecules can be considered as well-known: and neutron diffraction studies3-’ yielded a set of intramolecular structural parameters which are in good agreement with each other within the limits of the experimental errors; also the interatomic distances and vibrational amplitudes agree well with those given in the literature for the free molecules in the gas state.2-4 However, as far as the intermolecular structure is concerned, the structural models derived from diffraction studies are much less conclusive and contradictions also arose. A general agreement seems to exist among solution chemists that AN forms a rather structureless liquid in the sense that the correlation between the neighboring molecules is rather loose and thus the local order decreases rapidly with increasing intermolecular distances. However, the crystallographic model for the liquid structure, applied by Kratochwill et al.,I reflects a longer-range ordering in the liquid. Here a correlation cluster with a diameter of about 13 %, is postulated, comprising eight neighbors aligned in parallel or antiparallel directions related to a central molecule, an ordering which seems to be exaggerated. The nearest-neighbor model of Radnai et a1.2 resulted in a good fit between the experimental X-ray structure function and the theoretical one with the simplest assumption of two nearest-neighbor molecules being aligned in antiparallel directions and positioned about 3.3 %, from thecentral molecule. Beyond the first neighbors, a random distribution of the molecules was assumed. These model constructions were based on the concept that the preferred orientation of the neighboring molecules is governed by the dipoledipole interactions, while the distance of the closest approach is determined by the repulsive part of van der Waals interactions. The neutron diffraction experiments did not help to refine these structural models but supplied a very important contribution to the understanding of the character of the orientational correlations between the molecules. The so-called molecular pair correlation function was determined from a combination of an X-ray and two neutron diffraction experiments,3-’ and the three leading terms of an invariant series expansion were then computed. They Abstract published in Aduance ACS Abstracts. May 1, 1994.

provided information about the distance dependence of the orientational correlation between the molecules in the liquid and suggested that the dipole4ipole and also the quadrupolequadrupole terms play significant roles in forming the liquid structure. Theoretical calculations were also devoted to the determination of the liquid structure of AN, and while many details have been explored, new contradictions have also arisen. Thus, Steinhauser et a1.6 used a perturbation theory with a computer-generated reference function of liquid CS2 in order to determine the same three expansion coefficients of orientational correlation as above, but the agreement can be qualified as only partial. Several attempts were made to apply integral equation theory to the determination of the predominant factors which form the liquid structure of AN. Thus, the site-site Ornstein-Zernike (SSOZ) equation has been solved for liquid acetonitrile with various approximations, e.g. assuming six hard-core atoms without explicitly taking into account their electric charges,’ or a polar hard dumbbell fluid with site-site mean spherical closure including electrostatic charges in the modelEand employing hypernettedchain (HNC) closure and without explicitly incorporated charges.9-IO All of these models described the A N molecule with six atomic sites and provided structure functions and partial paircorrelation functions in acceptable agreement with the experimental results and other theories. Nevertheless, the interpretation of the liquid structure again evoked contradictory opinions. Namely, the models in which no explicit charges were involved for the atomic sites could still reproduce the main features of the measured pair-correlation functions, and thus a conclusion was drawn that the dipole-dipole interactions do not have a major effect on the structure of liquid AN, and therefore the good overall agreement with the diffraction experiments is not relevant in deciding the role of the dipole-dipole interactions. The integral equations, however, could not give a clear explanation of this fact, neither by providing an unambiguous set of partial paircorrelation functions, nor by describing the orientational correlation. Molecular dynamics”-12 (MD) and Monte Carlo13 (MC) simulations also targeted the liquid structure of AN. The MD simulations applied a six-site model for the pair-potential functions, including both Coulombic charges and Lennard-Jones

0022-3654/94/2098-5994$04.50~0 0 1994 American Chemical Society

A Reverse Monte Carlo Simulation of Acetonitrile

(LJ) type potentials, and it was concluded that even if the effect of the dipole-dipole interactions on the structure functions is minor, their role is predominant for the extension of the shortrange order in the structure, for the orientational correlations and for the dynamics (e.g. reorientational motions). The same qualitative conclusion was also emphasized by Jorgensen et al.13 from a M C simulation when a three-site model was applied with and without electric charges incorporated in the potential parameters, and the removal of the charges from the atomic sites led to drastic changes in the local structure. In order to make an attempt to clarify more details about the liquid structure of AN and to determine the extent in which diffraction experiments can provide detailed information on the local structure, the recently developed reverse Monte Carlo method (RMC) has been applied in the present work. Since, according to the concept of the RMC method, a very good agreement between the experimental structure function and the simulated one is a priori ensured, the obtainable sets of A N molecules can be considered as an extremely good model for the evaluation of the structure. Indeed, as it could be shown earlier by other examples1&15 of atomic systems, the agreement between the experimental and RMC structure functions is usually much better than in the cases where theoretical results are derived from integral equation methods or computer simulations. Moreover, the RMC method does not use potentials, and therefore there is no artificial separation between the Coulombic and non-Coulombic parts of interactions hidden in the selection of the potential, and as a consequence, the correlation functions can be studied on the configuration sets without pre-assumptions. We also note that, according to our knowledge, there are only two examples of the application of the RMC method to the structure of molecular (homonuclear, diatomic) liquids.16J7Our attempt is the first one where more than two different atoms are involved in the simulated molecules, and therefore several methodical details have also to be discussed here. Experimental Methods and Details of Simulations X-ray Diffraction Experiment. Technical details of the X-ray diffraction experiment on liquid acetonitrile and the method of the evaluation of the measured intensity have been published previously.2 However, the reported experimental structure function and radial distribution corresponded to a six-site model of AN; Le. the hydrogen atoms were explicitly taken into account as independent scattering units. In addition, all functions were normalized to the coherent scattered intensity of one nitrogen atom. In order to apply a three-site model, all structure functions were recalculated with the atomic scattering factors for N , C, and Me units, where Me represents a methyl group. All these scattering factors were calculated from analytical formulas, and the parameters necessary for the analytical functions were taken from the International Tables of X-ray Crystallography for N and C atoms1*and from Narten for the methyl group.19 Therefore, the experimental structure function was defined as

where k is the scattering variable, Iab(k)is the observed intensity corrected for background, absorption, polarization, and Compton scattering, and converted to atomic units, xi is the atomic fraction of the atom or atomic group i, andfi(k) is the related coherent scattering factor. The k scattering variable is defined as

where A is the wavelength of the incident radiation and B is the scattering angle. The M ( k ) modification function had the form of

The Journal of Physical Chemistry, Vol. 98, No. 23, 1994 5995

where /3 is a damping factor selected as @ = 0.01 A2. The corresponding radial distributions were calculated by Fourier transform of the structure function in the usual way.2 TheReverseMonte Carlo (RMC) Method. The RMC method for atomic liquids has been previously described in detail. The procedure is similar to the ordinary (Metropolis) Monte Carlo method (MC),20with the significant difference that instead of the minimization of the potential energy of the system the difference between the measured and the calculated structure functions or radial distribution functions of the system is minimized here. A peculiarity of the method is that several experimental functions, determined by different diffraction experiments, can be simultaneously employed within one simulation process and thus a more complete structural information can be extracted. In addition, any kind of previously fixed constraints imposed on the structure can also be incorporated into the system; this offers a facility to test a great variety of structural models. For a monoatomic liquid a t the beginning of the procedure an initial configuration of N points, each of them representing an atom, is placed in a cubic box of edge length L,which is selected in correspondence with the experimental number density. The usual periodic boundary condition, in which the simulation box is surrounded by images of itself, is applied. During the RMC procedure the partial pair-correlation functions gij(r)’s are calculated for each pair of i andj type atoms. Fourier transforms of gij(r)’s provide the partial structure functions Aij(k)’s: Ai,(k) = 1

+ *J”r(gij(r) k o

- 1) sin(kr) dr

(4)

wherepis theatomicnumber density. Thetotalstructurefunction is a linear combination of the partials, multiplied by the same M(k) modification function as defined in eq 3:

Finally the quantity x2is calculated from the difference between the simulated total structure functions kH(k)and the experimental one, kHE(k), in the following way: ‘k

(HE(kJ - H(kx))’

x2=x A= 1

(6) U2

where kx means the value of k for the Xth discrete point of the structure function, nk is the number of points in the k space, and u is the estimated error of the experiment. Here we supposed that only one experimental structure function is employed. Generalization of eq 6 for more experimental functions issimple.15 From a given configuration of the atoms a new configuration is generated by a random displacement of a randomly selected particle, and xncw2 is calculated for this new configuration in the same way as x2 was. If xncw2 I x2, the new configuration is accepted, otherwise it is accepted only with a probability of p = exp(-(Xnew2- x2)/2). The role of this probability in the MC method is to provide sampling according to a given distribution. However, in R M C it is applied in order to simply prevent the configuration from being trapped in a local minimum.15 The process is repeated until x2starts to oscillate, indicating that structural equilibrium is reached. If the procedure is further continued, a set of configurations consistent with the experimental data, within the given error u, can be generated. Although there is no theoretical reason why exactly the same method cannot be used for modeling molecular liquids, to our

5996

The Journal of Physical Chemistry, Vol. 98, No. 23, 1994

knowledge, studies on homonuclear systems only as liquid halogens, nitrogen, sulphur, and phosporus were reported so far.16J7 Thereason for this neglect of morecomplicated molecular systems might be that when rigid molecular models are used, the molecular part of the structure function can hardly be reproduced, and thus an overall fit of only rather poor quality is obtained. In order to overcome this problem, the concept of flexible molecules was introduced in the following way. Instead of whole molecules weconsidered atoms as the smallest units in the system, and in each simulation step we moved single atoms rather than whole molecules. (The methyl group is called here an atom for simplicity). During the simulation we set up a number of coordination constraints in order to confine the intramolecular distances into definite ranges and, in such a way, prevented the atoms from forming molecules of unrealistic geometry. In the present work the coordination of atoms of type A by atoms of type B was determined by two constraints. The first one madeeach atom A zero-coordinated by atoms B within the sphere of radius R 1 ,while thesecond made it one-coordinated within Rz. The two constraints together caused each A atom to have exactly one B neighbor between distances R1 and Rz. For the sake of secure separation of intra- and intermolecular distances, one might use a third constrant on coordination of atoms of type A by B as well, which makes atom A one-coordinated within the sphere of radius R3 > Rz. Because this coordinating B atom is also closer than R2, the introduction of this third constraint means that between RZand R3 there is not any B atom, so this constraint can act as an intermolecular cutoff distance. However, due to the special distance distribution of the system of liquid acetonitrile molecules, we did not need to use this third constraint for any pairs of atoms. The constraints between neighboring atoms allow the bond lengths to vary within a limited range, and therefore the bond angles can also vary within a limited range only. A drawback ofthe method is that theseconstraintsaredefined by coordination numbers exclusively, and no guarantee is given that the coordinating atom is the one which belongs to the same molecule during the whole simulation run. Consequently, it becomes very important to start from a configuration in which the atoms are grouped corresponding to realistic molecular geometries and to fully satisfy all of the constraints. Since these constraints cannot be violated during the simulation, there is no possibility for coordinating atoms to be exchanged and then realistic molecular geometries to always be provided, and those atoms which form a molecule a t the beginning of the simulation remain unchanged during the whole process. Simulation Details for Three-Site Liquid Acetonitrile. The RMC simulation has been run with 512 acetonitrile molecules placed in the basic box. During the RMC procedure the experimental X-ray structure function was fitted. Since each acetonitrile molecule was considered as consisting of three ’atomic” sites of N , C, and Me, the total number of atomic sites in the box was 1536. Representation of the sites and a local coordinate system for the A N molecule are shown in Figure 1. The number density of the system was adjusted to the experimental density and was fixed at p = 0.034 16 sites/A3, corresponding to a half-box length of 17.78 A. The distances of allowed closest intermolecular approach for N-N, C-C, and MeMe pair dmin(‘cutoff values”) were fixed together with the constraints imposed on the intramolecular distances (the closest allowed approach R1 and the maximum separation distances R2 for each pair within a molecule). The cutoff and constraint distance values are given in Table 1. During the fitting procedure the maximum distance for an instant move of an atomic site was fixed a t 0.1 A. With these parameters a relatively fast convergence could be achieved, with an acceptance ratio of successful moves to the total number of movesof about 1:3. Thestartingconfigurationfor thesimulation

Radnai and Jedlovszky 0.08

I

2 1

..

0.06

I

3.04

*

1

I

I I *

I

. ’

I

***.‘

0.00

1200

136.0

1460

156.0

166.0

170.0

186.0

a (deg)

F w e 1. DistributionP(a)of theintramolecularangle aof theacetonitrile molecules in the RMC simulation run. The insert shows a schematic view of a three-siteacetonitrilemolecule, the molecular coordinateframe, and the definition of the angle a.

TABLE 1: Constraints Imposed on the Three Sites C, N, and Me in Terms of RI, Rs,and d h in Order To Define the Intramolecular Structure of the Acetonitrile Molecule and the Closest Intermolecular Approach. For Definitions See Text site A

site B

N N N

N

C C Me

RI(&

&(A)

1.12 2.57

1.18 2.67

1.44

1S O

(A)

&in

1.50

C Me C Me Me

1.90 3.00

TABLE 2 Potential Parameters Applied for a Three-Site MC Simulation of Liquid Acetonitrile. All Listed Parameters Have Their Conventional Meanings and Are Explained in the Literature” t

site q (e)

u (A)

N

3.2 3.65 3.775

C Me

-0.43 0.28 0.15

A

B

(kcal mol-]) (kcal mol-’ Al2) (kcal mol-’ A b ) 0.1698 0.15 0.2068

783 064.3 3 354 792.2 6 928 097.6

729.29 1418.76 2393.93

was generated with a short M C simulation with suitably adjusted potential parameters. When a starting configuration in equilibrium was obtained, the RMC simulation started. After about half a million moves, when a structural equilibrium was reached again, 200 sample configurations were collected and evaluated from the next 300 000 steps; thus the statistical error on the resulted functions was held within acceptable limits. Monte Carlo (MC) Simulation. For comparison, a conventional Monte Carlo simulation on the same system was run. The conditions and potential parameters were essentially the same as those published by Jorgensen et al.IO An NVT ensemble was used. The number of molecules, the box length, and the number density agreed with that of the RMC simulation. The (rigid) acetonitrile molecule was again represented by three atomic sites and kept fully linear with interatomic distances rN4 = 1.157 A and rC-Me = 1.458 A. The pair potentials included Coulombic terms as well as Lennard-Jones terms. For the detailed mathematical formulas of potentials we refer to the literature.10 The potential parameters are listed in Table 2. The interaction energies were calculated up to a limiting distance of 10.5 A, affecting about 55 molecules in the neighborhood of a central one. The correction energy for the cutoff was about 3.5% of the total energy. The temperature of the system was 25 O C . From a random configuration, after a preliminary equilibration of the system during the first 1 million attempted moves, 100 independent configurations were collected for further evaluation, selected at each 10 000 moves, thus spanning an additional 1 million total

A Reverse Monte Carlo Simulation of Acetonitrile

The Journal of Physical Chemistry, Vol. 98, No. 23, 1994 5997

OI

100

0.00 2.00

-lOi, , , '

-1 5

00

20

40

1 .oo

,

,

,

60

80

100

, 120

~

140

k (A.')

Figure 2. Experimentalstructure function kHE(k)obtained by the X-ray diffraction method (full line) and the best fitting structure function kH( k ) (stars) determined by the R M C procedure for liquid acetonitrile.

0.00 2.00 1.00

0.00

-

2.r4.60

moves. Theacceptanceratioofmoves wasabout 2:3. Theallowed upper limits were 0.1 A for a move and 10' for a rotation.

Results and Discussion Intramolecular Structure. An examination of the results of previous diffraction experiments2 shows that the average structure of an acetonitrile molecule is linear and that the intramolecular distances fall within the ranges of rN4 = 1.14-1.15 A, rN-Me = 1.46-1.48 A, and R c - M= ~ 2.61-2.62 A. The peak positions of the total X-ray experimental g(r) function correspond to these distances, and results of all reported geometrical models are in good agreement with these parameters. Due to the RMC procedure described previously, the total g(r)function is consistent with the experimental data and therefore the average structure has to coincide with the experimental findings. However, due to the present definition of the acetonitrile molecules uia independent sites held together by constraints, the individual acetonitrile molecules can deviate from their average linear shape in any instantaneous structure. This fact is characterized by the distribution of the N-C-Me angles denoted by a. Since the MC and thus the RMC method cannot account for the dynamical properties, no results for frequencies of intramolecular vibrations can be directly obtained, and thus the effect of the introduced elasticity in the molecules cannot be directly evaluated. Nevertheless, a rough estimation of the reliability of the deviation from the linearity has been carried out in the following way. Figure 1 shows the distribution P(a)of angle a. The maximum is located at about 168', which means that the most frequent deviation from linearity is about 12'. The figure also shows that the probability of extremely large deviations from linearity is very small, practically zero for values larger than 35' (a1145'). If, at a rough approximation, weconsider the potential parameters used in the M C simulation to be valid also for the intramolecular interactions, a relative energy change of about 1 kcal/mol can be calculated, associated with the bending of a molecule from a linear shape to a = 168'. Parallel ab initio calculations for a free acetonitrile molecule,21 with MP2/6-31G** basis set and full geometry optimization, resulted in a rather flexible molecule, with a bending potential value of about 2.5 kcal/mol for a = 168O, which corresponds to a bending frequency at about 200300 cm-1. In fact, from force field calculations based on spectroscopic calculations, the bending mode frequencies resulted in 364 cm-1 for the gas phase, with a -17-cm-I shift in the liquid phase,22 which proves the reliability of the ab initio result, and which, in turn, produces a little less flexible molecule than that which results from the coordination constraints in our RMC calculation. This showsqualitatively that (i) the molecular model in our RMC calculations is realistic in the sense that the flexibility is not exaggerated at all and (ii) the assumption for fully linear molecules used for our M C calculation and for the previous molecular models in the literature1-l3needs appropriate correction. Comparison of RMC and X-ray Results. Figure 2 shows the agreement between the experimental structure function deter-

6.00 8.00 l0.bO 2.00 4.00

6.dO

8.00 lO.bO

r(A)

Figure 3. Partial pair-correlation functionsg&) as computed from the R M C procedure for all six ij pairs of atomic sites in the three-site

acetonitrile molecule (stars). The lower and upper full lines show the contribution of the nearest-neighbor molecules and the nearest 10 neighboring molecules, respectively, to the partial pcfs. mined by the X-ray diffraction method and that obtained by the R M C fitting procedure. In order to demonstrate the excellent agreement of the two curves, we selected the structure function forms as defined by eqs 1 and 5 rather than the usual H(k) form which damps rapidly and does not show the smaller details of the functions. The quality of the agreement reaches the level of the earlier R M C simulations performed for atomic liquids; it is comparable with the best fits ever reached by fitting simple geometrical models to the structure functions by the least-squares method; and it is better than in most cases in which results from traditional MD or M C simulations are compared with experimental data. All characteristic features of the experimentalcurve (positions, periods, and amplitudes of the maxima) are fully reproduced, and the RMC structure function oscillates around the experimental one, producing a difference curve which cannot be related to the liquid structure but probably to the differences in the maximum distances up to which the radial distribution functions have been calculated (10 A for the experimental and 17 A for the RMC result). Partial Pair Correlation Functions. All partial (site-site) pair correlation functions gij(r)are shown in Figure 3, as computed on the basis of the set of configurations collected from the RMC run. All pcfs exhibit clear features in forms of peaks, shoulders, and minima, which are slightly disturbed by minor ripples caused by statistical uncertainties. The features of the curves can be compared to previous results obtained from integral equation calculations10 and from M D and MC~imulations.~~J3 As far as thelatter twostudiesareconcerned, a detailed comparison of the relevant pcfs could not reveal any significant differences in the overall shapes of the curves obtained from either a six-site molecular dynamicsll or a three-site Monte Carlo study." On the other hand, visible discrepancies can be observed when the pcfs from the simulations are compared to thoseobtained by solving the SSOZ integral equations with HNC closure.1° The situation is further complicated when attempts are made to distinguish between curves when site charges are explicitly included or omitted in the integral equation study, e.g. for g " ( r ) and gMeMc(r), both being remarkably different from each other and from the MD result as well. At this level it is difficult to decide whether the forms of the potentials, the assumptions of the theory, or structural consequences of the explicitly included charges are responsible for these discrepancies. Previous attempts to investigate the effect of the charges proved their dramatic influence on the pair correlation functions, but only when the fully charged potentials were compared to those in which the charges have been completely omitted.1°J2J3 The non-Coulombic part of the pair potentials, however, has been kept intact, which makes the comparison physically unrealistic.

5998

The Journal of Physical Chemistry, Vol. 98, No. 23, 1994

We note here that the influence of different point charges incorporated into the potential models on the structure functions has also been analyzed for a M D simulation of liquid methanol,23 with similar conclusions. A further remark has to be made here. All pair-correlation functions obtained in the present study from the conventional MC simulation are in very good agreement, within the limits of statistical noise, with those published by Jorgensen et al.,I3 and therefore we omitted them from the discussion of the structure at this level. When the partial pcfs of Figure 3, determined by the RMC method, are compared to the literature results, a good overall reproduction can be observed. Indeed, the maxima and minima, andmost ofthe peak positionsand heights,arethesameasreported by previous authors. This demonstrates that the RMC method is able to produce structural information of reliability at least comparable with those obtained by other theoretical methods. Some minor details, however, need further investigation. As a starting point, the gN-Mc(r)function exhibits the sharpest peak of all curves at about 3.8 A, and these features were also found in the simulation studies. The position, width, and height of this peak of the RMC curve is in good agreement with those from simulation studies, and this is also valid for the minimum and second maximum, except that the first peak is more asymmetric and the second peak is slightly lower and broader. On the contrary, the SSOZ study with q = 0 charge did not reproduce the height of the first peak. This feature of the gpJMc(r) function has led the previous authors to the conclusion of a preferred head-to-tail configuration of the neighboring molecules or the antiparallel alignment, and the orientations were further explained by the dipole-dipole interactions or, more generally, by the effect of the charges. The RMC study reproduced the result, but without a direct involvement of the electric charges in the structural picture. The C-C, C-Me, and N-C functions are reported to be rather similar to each other, computed either from simulation or from integral equation studies. In fact, the present result does not represent an exception. The first peak of gcc(r) is somewhat broader and lower than those reported earlier; the shoulder on the right hand side of the first peak of gCMe(r)is more pronounced, and the first peak of gNc(r) is shifted down by about 0.2 A, and it is narrower than those in the case of M C simulation, but it is hard to draw any significant conclusion on the impact of these differences on the liquid structure. Major discrepancies have been reported among the shapes of curves determined by different methods for N-N and Me-Me pcfs. The mutual orientations of the molecules are strongly reflected in these distributions of like atoms a t end positions in the molecules. The behavior of the pcfs in Figure 3 is interesting: the position and width of the first peak of g”(r) reproduces those determined from integral equations with q = 0 (but its peak height is increased by about 0.5 units), which contradicts any other model where the charges are explicitly included. On the contrary, the gMeMc(r) function from R M C is very similar in its shape to what was obtained from the MD (or MC) studies. The first peak is shifted down by about 0.5 A, which is a significant change in the peak position, but it is still about 0.3 higher than that obtained from the integral equation study with q = 0. These features show again that the role of the charges cannot be simplified to a yes-or-no question, and that the orientation of neighboring molecules can play a more complicated role than analyzed previously. Coordination Numbers. A literature survey reveals that the reported estimates for the average number of coordinating neighbors around a central AN molecule varies between 2 and 13. There are obvious definition problems in the background of such a discrepancy. The coordination number CAB of atom A coordinated by B atoms is most frequently defined by an integral

Radnai and Jedlovszky

TABLE 3: Coordination Numbers for Each ij Pairs of Atoms Derived by (I) Integration up to the First Minimum r h of the gdr) Function and (11) Integration up to 5.1 A, the Minimum Position of the m c ( r )Function definition definition site i sitej r d n (A) I I1 N N N

N C

C

C

C

Me Me

Me

Me

4.4 4.2 5.1 6.1 6.2 6.2

3.2 2.1 5.9 10.6 11.2 11.5

5.2 5.6 5.9 5.8 6.2 5.9

over the pair-correlation function gAB(r)

(7) hence referred to as definition I, where X B is the atomic fraction of atom B,rl might be the first non-zero point of the function, and the upper limit r2 is usually the first minimum position. Another definition for thecoordination number can be introduced when one fixes a certain given distance for r2 and each pcf is integrated up to this common upper limit (definition 11). For illustration see Table 3, where a set of coordination numbers is shown by using the above two definitions. A surprisingly large variation in the CAB values can be observed from 2.7 for N-N to 11.5 for Me-Me pairs. These variations are certainly connected with the different sizes of the atomic sites and also point to the role of different possible orientations in the formation of the first coordination shell. The lower limit suggests that a maximum number of 2.7 neighbors can approach the N atom of the central molecule with the N atoms. On the contrary, the relatively bulky Me groups are scattered around a central Me group, producing a broad first peak in the gMeMe(r) curve and, consequently, a much higher coordination number, too. If the coordination shell occupies a broader space, the probability of finding differently oriented molecules increases. In other words, the different coordination numbers suggest that the neighboring molecules can exist in different coordination states: nearest neighbors in close approach with preferred orientations and other, more or less randomly oriented, neighbors at larger distances. However, when an upper integration limit of 5.1 A, corresponding to the minimum position of the gNMe(r),is chosen, the obtained CAB values fall in the close range of 5.2-6.2 for each pcf. This value is in good agreement with Jorgensen’s result.13 The latter fact shows that, within a sphere of approximate radius of 5 A, the orientational correlations damp rapidly and a similar average coordination number results for each atom pair. In order to further clarify the interpretation of the coordination numbers, another approach was also introduced. Instead of integrating up to a certain distance and obtaining the average number of neighboring molecules, the distribution of the first n neighbors was determined (n = 1,2, ..., 10)by direct evaluation of thepositiondistributionsoftheatoms. Theresulteddistribution curves are shown for each atom pair in Figure 3 for n = 1 and 10, and in an enlarged scale for n = 1-3 and 10 for gNMe(r) in Figure 4. In order to derive these distributions, a local coordinate system was defined in which the neighbors are located. In our approach, the local coordinate system was defined by thecentral acetonitrile molecule as it can be seen in the insert in Figure 1. The origin of the coordinate system coincides with the N atom of the central molecule, and the positive half of the Z axis was directed toward the Me group. Since we had planar molecules rather than linear ones, the molecular plane was used for the definition of the X Z plane in such a way that X coordinate of the C atom was always positive. In the case of the M C simulation with fully linear molecules, the X Z plane was defined by a random selection of the X axis, which is allowed by the cylindrical symmetry of the

A Reverse Monte Carlo Simulation of Acetonitrile

The Journal of Physical Chemistry, Vol. 98, No. 23, 1994 5999

a

b

3.00 h v i

c, 2.50

2.00

2.00

ii

1.50

1.00

-

0.50

-

0.00 2.00

4.60

6.60 8.60 lO.bO 12.00

r(A)

r(A)

Figure 4. Partial pair-correlation functions g,j(r) as computed from the RMC procedure for N-Me pairs (stars) and contributions computed for the nearest n neighbors of the central molecule (full lines). The neighbors are selected according to the location of the centers of the molecules (a) or

according tothe location of the nearest atoms (b). molecule. After the local coordinate system was defined, all molecular positions and orientations were determined relating to it. For the collection of nearest-neighbors' statistics, two alternate definitions were used: In definition 1 the molecules were distinguished and sorted according to their center+enter distances (for practical purposes we considered the positions of C atoms as molecular centers). In definition 2 the distances of the closest atom pairs between the central and each neighboring molecule were searched and the molecules were sorted according to these distances. It is obvious that both definitions give their preferences to certain extreme orientations of molecules: e.g. if parallel or antiparallel alignments of the neighbors are preferred, the closest approach is achieved by their centers, while in T- or L-shaped configurations the end-atoms of the molecules can better approach each other. Thus, the definition used for the selection of the neighbors can sort them, to a certain extent, also according to their preferred orientations. This is demonstrated in Figure 4 by the different shapes of the neighbor distributions from n = 1 ton = 3. In the case of definition 1 (Figure 4a), the nearest neighbors form an increasingly higher and broader distribution. The shapes of the gi;(r) curves are similar up to n = 10, with a shift in the peak position from about 3.8 to 4.0 A, due to the fact that with increasing n the number of molecules with more distant centers increases in the statistics. On the contrary, definition 2 produces a split peak for neighbors n = 1-3 (Figure 4b). For n = 1 and 2, the left side of the peak is higher, while for n = 3 the tendency changes and by n = 10, the right side of the peak grows up to the main peak of the gNMe(r) function. This fact illustrates that, on average, the first two neighbors are often oriented differently from the further ones. In this way, two different orientational states in the coordinating molecules can roughly be distinguished. The interesting feature is the gradually emerging right shoulder in the neighbors' distribution. This represents, on one hand, a fraction of molecules which are located at definitely larger distances than the first one to two neighbors, and also beyond the minimum of the gNMc(r) curve. On the other hand, the shoulder a t the right side is partly due to the consequences of definitions used for the selection of the neighbors. Since each atom pair between the central and the neighboring molecules contributes to the statistics, in the case of certain orientations the short r,j distance may be coupled with a much longer rji distance. The most extreme example is a headto-tail configuration where, if the nearest neighbor's Me group is located 3.8 8, from the N atom of the central molecule, corresponding to the maximum position of the gNMe(r) function, then a much longer rMeN = 6.4-6.5 A distance can also be expected to appear. The shape of the shoulder suggests that the probability of having such head-to-tail configurations must be small but supports the idea that a largevariety of orientations with different orientational angles appear with increasing n.

-

0.06

.. . . . ZMC MC

Q

0.02

I1 \ 1

'

0 00 -1.0

-0 5

0.0

.. .. ,.............* ...* . *..

1 .o

0.5

cos$

Figure 5. DistributionsP(cos 19) of the cosine of the angle 19 as derived from RMC (full lines) and MC (stars) simulationsfor liquid acetonitrile, when the 1 or 10nearest neighborsare sampled around a central molecule. The angle 19 represents the direction of the neighbor's Z axis related to the central one's. The neighbors were selected by using definition 1 (see the text). Angular Distributions and Orientations of Molecules. The mutual orientation of neighboring molecules can be characterized, in a first approximation, by the distribution of an angle between the characteristic vectors of the two molecules. In terms of our definition of the local coordinate system we have selected the angle 6 between the 2 axes of the molecules, which should be good approximation for the dipole moment's direction since the shapes of the molecules are not far from linear. Figure 5 compares the distributions P(cos 6 ) derived from both the traditional M C and the RMC simulations. Two distribution curves are presented for both cases; for n = 1 and n = 10, definition 1 (i.e. based on the center-center distances of the molecules) has been used for the selection of the neighbors. In the case of randomly oriented molecules the cosine statistics should result in a constant value, independent of angle. On the contrary, the parallel, antiparallel, and T- or L-shaped orientations should result in peaks a t 6 = 0, 180, and 90°, corresponding to cos t9 = 1, -1, and 0, respectively. The first overall observation is that the main features of the cosine distributions for the two methods are very similar to each other. For n = 1 the distributions show a strong preference for cos 9 = -1, which clearly demonstrates that the nearest neighbors tend to be oriented in an antiparallel form. Other bumps for possible spots of slight preferences can also be observed, e.g. at about cos 19 = -0.7, -0.4, and 0.52, corresponding to 135, 115, and 60°,and especially at about cos 6 = 1, which corresponds to the parallel alignment of the neighbors. When n = 10, the distribution is almost constant (except a minor peak at cos 6 = -1, the consequence of the first neighbors' predominating peak), which proves that the preferential orientation is lost very quickly when the averaging is extended to several neighboring molecules.

Radnai and Jedlovszky

6000 The Journal of Physical Chemistry, Vol. 98, No. 23, 1994 b

a

0 01

m1

0 00 0 00

30

I

I

50 00

100 00

150 00

0 0'

0 00 0

fl( d e 9 1 fi( d e 9 1 Figure 6. Distributions P(9) of the angle I9 between the neighboring molecules' Z axes and the central one's, for the first n neighbors ( n = 1, 2, 3, ..., 10). The neighbors were sampled by using definition 1 (a) and definition 2 (b). It is important to emphasize that both the traditional MC and the RMC methods lead to very similar distribution curves and qualitative conclusions as above. However, the MC seems to slightly underestimate the preferred specific orientation when compared to the RMC result, and more preference is given to the random orientation. This difference, however, is not significant: if we integrate the distribution curves from cos 9 = -1 to cos 9 = -0.9 and -0.85 (approximately within the half-widths of the RMC and MC peaks, respectively), we find that 22.0 and 23.4% of the first neighbor molecules fall within this range, Le. when the maximum deviation from the antiparallel alignment is about 30'. If the integration limit corresponds to an angle range of 40°, the percentage values increase to 37.4 (RMC) and 32.4 (MC). Similar integration of the small peak located near cos 9 = 1 until the local minimum at cos 9 = 0.83-0.85 reveals that about 5.1% (RMC) and 4.7% (MC) of the molecules are aligned parallel to the central molecules with an uncertainty of about 30'. These data show that the antiparallel alignment is much more predominant for the first neighbors than the parallel, and this is confirmed both by the RMC and by the traditional MC simulations. Figure 6 shows the effect of the different definitions of the neighbors on the weighting of the preferential orientation. For thesakeof better visualization, the angular distributions areshown instead of the cosine distributions of 9. Here the random orientation should correspond to a sine curve. Parallel alignment should cause a local maximum near ,'O antiparallel alignment at 1 80°, while a T- or L-shaped configuration should be detected in the form of a local peak or bump superimposed on the broad maximum of the sine curve at 90 'C. When definition 1 is used (Figure 6a), the features discussed in relation to thecorresponding cosine distributions can be clearly detected: when n = 1, strong preference is given for the antiparallel alignment; when n ranges from 1 to 3, the predominant antiparallel orientation decreases gradually, and it is converted to an almost random orientation for n = 10. When definition 2 is used (the neighbors are sorted according to the distances of the nearest atoms, Figure 6b), smaller preference is given to the antiparallel orientation, and the fraction of intermediate configurations increases. Moreover, this feature does not changedrastically when going from n = 1 ton = 10. This means that there is no clear indication of the relative preference of T- or L-shaped molecular orientations, except a minor peak near 90'. The effect can easily be explained if one considers that, for the antiparallel alignment, in order to reach the closest approach, the centers of two neighboring molecules cannot move too far from theXYplane (because of the steric hindrance caused by the larger Me groups to the smaller N atoms), but T- or L-shaped configurations can be realized at various distances from theXYplane, and thus the chance of having more T- or L-shaped configurations at most distant locations from the central molecule increases. It is worth noting that for n = 10 the angular distributions in Figure 6a,b become very similar, which shows

that the selection criteria set in the two definitions play a significant role only for the nearest neighbors. Based on the angular distributions alone, it is very difficult to analyze the other possible orientations, as e.g. head-to-tail configurations (along the 2 axis), because they cannot be clearly distinguished from the parallel orientation. Moreover, in the liquid state, the deviation from the favorable alignment along the same axis can be large when the intermolecular distances are larger. Spatial Distributions of the Neighboring Molecules. Another approach to reveal the coordination structure of the nearestneighbor molecules is through the analysis of the spatial distribution of the molecules around a central one. The same coordinate system and the same definitions can be used for the selection of the neighbors as above. Three-dimensional pictures were constructed for the two-dimensional projections of the position distributions of selected atoms of n neighboring molecules, which are presented in Figures 7-9. Figure 7 shows projections onto the XY plane, Le. perpendicular to the main axis of the central molecule. Although some asymmetries can be observed when n = 1 and definition 1 is used for both the N atoms and Me groups projected, the cylindrical symmetry is obvious in each case. Drawings of greater variety and different shapes can be obtained when projections onto the XZ plane are generated for both definition 1 (Figure 8 ) and definition 2 (Figure 9). The most characteristic feature of the figures for n = 1 is the appearance of two main peaks below and above the origin (Le. the N atom in the center) along the X axis, irrespective of which atom has been projected and which definition has been used. This split, however, is a trivial consequence of the excluded volume of the central molecule. The asymmetry in the two peaks can also be explained as a consequence of the model introduced, because a larger volume is excluded in the positiveXaxis direction where the C atom is located when the molecule deviates from linearity. The third artifact is a shift of the distributions to the positive 2 axis direction. The structure of the coordination shell is reflected in some minor, but still visible, effects on the shapes of the projections. When definition I is applied and n = 1 (Figure 8a-c), the projection of each atom results in roughly similar distributions, which again can be explained by the predominance of antiparallel orientations. However, the non-zero value of each distribution at its origin witnesses that a portion of the first neighbors are either T-oriented or, more probably, shifted along the 2 axis. When definition 2 is applied and n = 1 (Figure 9a-c), significant contributions appear at positive and negative directions along the Zaxis even when theN atoms are projected. Since thedistribution is broad enough, it can prove the presenceof either shifted parallel or antiparallel neighbors, or head-to-tail configurations as well. The filling up of the gap above the origin, when the Me groups are projected, can be explained by the predominance of the antiparallel orientation with Me atoms of the neighboring

A Reverse Monte Carlo Simulation of Acetonitrile

The Journal of Physical Chemistry, Vol. 98, No. 23, 1994 6001 2

t

X

a

b

d

C

Figure 7. Three-dimensional drawings of the projections of the position distributions of atom i of the n nearest-neighboracetonitrile molecules around a central one onto the XYplane of the central molecule’s local coordinatesystem (shown in the insertion of Figure l), calculated from the RMC simulation of liquid acetonitrile: (a) i = N, n = 1, definition 1; (b) i = Me, n = 1, definition I ; (c) i = Me, n = 1, definition 2; (d) i = Me, n = 10, definition 1.

a

N

b

d

e

Figure 8. Three-dimensional drawings of the projections of position distributions of atoms N, C, and Me of the n = 1 (a+) or n = 10 (d-f) nearest-neighboracetonitrile molecules around a central one onto theXZ plane of the central molecule’s local coordinate system (shown in the insertion of Figure l), calculated from the RMC simulation by using definition 1.

molecules located below and above the N atom of the central molecule in the origin. When n = 10 neighbors are projected (Figures 8d-fand 9d-f), the distributions broaden and tend again to exhibit a cylindrical symmetry. Apart from the characteristics caused by thedefinition of the coordinate system, only smaller contributions can be assigned to certain favorable positions. The two definitions lead again to similar distributions. All of these facts show the loss of the predominant orientation and the gradual change to completely random orientations.

Summary and Conclusions The reverse Monte Carlo method has been applied to determine the microscopic structure of liquid acetonitrile. Our results show that the appropriate modification of the original form of the method can provide a reproduction of the experimental (X-ray)

a

b

N

d

e

Figure 9. Three-dimensional drawings of the projections of position distributions of atoms N, C, and Me of the n = 1 (a+) or n = 10 (d-f) nearest-neighboracetonitrilemolecules around a central one onto theXZ plane of the central molecule’s local coordinate system (shown in the insertion of Figure l), calculated from the RMC simulation by using definition 2.

structure function as good as that for atomicliquids. An additional gain of the coordination constraints method is that instantaneous deviations from linearity of the acetonitrile molecules can be simulated and thus the intramolecular geometry becomes more realistic. The microscopic ensembles generated by the R M C method provide a reliable set of molecular configurations, from which many properties of the coordination shells can be deduced. The extremely good agreement between the experimental and the computer-generated structure functions authorized us to analyze the partial pair-correlation functions, angular distributions, and spatial distribution of molecules in terms of a “best”, although evidently not unique, microscopic model of the structure. The R M S result helped to check the reliability of an auxiliary traditional MC simulation, and an overall good agreement was produced between the structure functions and partial pcf‘s.

6002

The Journal of Physical Chemistry, Vol, 98, No. 23, 1994

Detailed analysis of the liquid structure showed that a reasonable average coordination number of 5-6 can be derived from the first coordination shell if the integration is extended up to about 5 A, almost twice the length of an A N molecule. The coordination number is in good agreement with an earlier M C simulation result;'3 it is higher than that reported previously on the basis of a simplified geometrical model fitted to an X-ray diffraction experiment2 and lower than those reported from an X-ray diffraction analysis' and an MD simulation.'' Beyond the discussed differences in the applied methods and the definition for the coordination number itself, this discrepancy points to the fact that, in a relatively less ordered molecular liquid like AN, it is difficult to define unambiguously the boundary of the first coordination shell. The partial pair-correlation functions presented here suggest, and the angular and spatial distributions confirm, that the predominant orientation of the nearest-neighbor molecules is antiparallel, in accordance with previous studies. Indications for parallel aligned neighbors, which are presumably shifted along and/or aside the Z axis of the molecule, are found, but we agree with Bohm et al." that their amount must not be significant. Also, the perpendicular (T- or L-shaped) alignment is not significant as far as the first two to three neighbors are concerned, in agreement with the conclusions of Steinhauser et a1.6 We propose two essential conclusions from our R M C results for the orientational correlations between the molecules. The first is that the predominance of the antiparallel orientation of the neighbors decays rapidly. In this sense the two or three nearest neighbors are in a different coordination state than the further ones still in the same first coordination shell. The second conclusion is that while the predominance of the antiparallel orientation is lost, noother (e.g. perpendicular) orientation seems to become predominant. In terms of the expansion coefficients of the generalized pair correlation functions characteristic of the orientational correlations, this means that our result predicts that the relating coefficients should damp quickly, like those determined by the perturbation theory: while those derived from the diffraction experiments are probably much too structured.3 This statement needs verification, and therefore the analysis of the expansion coefficients derived from the RMC result is in progress at our laboratory. A question of major importance in the description of the liquid structure of diatomic and linear triatomic molecules was whether the electrostatic interactions play a decisive role in the formation of the (static) structure or not. According to the present results, our opinion is that the major factor in the formation of the structure is the molecules' size and shape. This is relevant both for the total structure function, the partial pcf s, and for the orientational correlations as far as they are accessible. However, two remarks have to be made: (i) the molecular shapes as used in our representation comprise the interatomic interactions implicitly,

Radnai and Jedlovszky of both the van der Waals type and the electrostatic ones, since all intramolecular dimensions were adjusted to fit the experiments directly. Artificial separation of the different kinds of interactions is not conclusive enough to decide which has the major effect on the structure. (ii) Even if the agreement between experimental structure functions and those obtained from theory is visibly good, the smaller deviations become important when the reliability of the applied potentials have to be checked; from this point of view the R M C method is able to produce superb results for a further check. Especially when integral equation theories are to be checked, the R M C method can play an important role in the research of molecular liquids.

Acknowledgment. The authors are indebted to Dr. L. Pusztai (ELTE University, Budapest) for his valuable initial discussions on the topic and to Drs. M. Howe and R. L. McGreevy (Oxford University, U.K.) for the R M C code provided. We are also grateful to Prof. Veszprbmi (Technical University, Budapest) for his ab initio calculations and remarks on the flexibility of the A N molecules. Discussions with Dr. G. Ptilinkls and I. Bak6 (C.R.I.C, Budapest) are also acknowledged. The research was partly supported by the National Scientific Research Fund (OTKA), Project No. 1087 (459/91). References and Notes (1) Kratochwill, A.; Weidner, J. U.; Zimmermann, H. Ber. Bunsen-Ges. Phys. Chem. 1973. 77.408. Radnai, T.; Itoh, S.; Ohtaki, H. Bull. Chem. Soc. Jpn. 1988, 61,

Bertagnolli, H.; Zeidler, M. D.Mol. Phys. 1978, 35, 177. Bertagnolli, H.; Chieux, P.; Zeidler, M. D. Mol. Phys. 1976,32,759. Bertagnolli, H.; Chieux, P.; Zeidler, M. D.Mol. Phys. 1976, 32, Steinhauser, 0.;Bertagnolli, H. Chem. Phys. Lett. 1981, 78, 5 5 5 . Hsu, C. S.;Chandler, D. Mol. Phys. 1978, 36, 215. Morriss, G. P.; Monson, P. A. Mol. Phys. 1983, 48, 181. Fraser, K. J.; Morriss, G. P.; Dum, L. A. Mol. Phys. 1986,57,1233. Fraser, K. J.; Dunn, L. A.; Morriss, G. P. Mol. Phys. 1987,61,775. BGhm, J. B.; McDonald, I. R.; Madden, P. A. Mol. Phys. 1983,49, Ohba, T.; Ikawa, S. Mol. Phys. 1991, 73, 999. Jorgensen, W. L.; Briggs, J. M. Mol. Phys. 1988, 63, 547. McGreevy, R. L.; Pusztai, L. Mol. Simul. 1988, 1, 359. Howe, M. A. Physica B. 1989, 160, 170. Howe, M. A. Mol. Phys. 1990, 69, 161. Howe, M. A.; McGreevy, R. L.; Pusztai, L.; Borzsik, I. Phys. Chem.

Liq. 1993, 25. 205. (1 8) International Tablesfor X-ray Crystallography; The Kynoch Press: Birmingham, England, 1974; Vol. 4, p 99. (19) Narten, A. H. J. Chem. Phys. 1979, 70,299. (20) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. M.; Teller, A. H.; Teller, E. J. Chem. Phys. 1953, 21, 1087. (21) VeszprCmi, T. Private communication. (22) Jdkli, Gy.; Koritsinzky, T.; Jancs6, G. Ber. Bunsen-Ges. Phys. Chem.

1981, 85. 773. (23) PBlinMs, G.;Tamura, Y.;Spohr, E.; Heinzinger, K. 2.Naturforsch. 1988,43a, 43.