摘要:

GENERAL DISCUSSIONS OF THE FARADAY SOCIETY 317 Date 1964 1964 1965 1965 1966 1966 1967 1967 1968 1968 1969 1969 1970 1970 1971 1971 1972 1972 1973 1973 1974 1974 1975 1975 1976 1977 1977 1978 1978 1979 Subject Volume Chemical Reactions in the Atmosphere Dislocations in Solids The Kinetics of Proton Transfer Processes Intermolecular Forces The Role of the Adsorbed State in Heterogeneous Catalysis Colloid Stability in Aqueous and Non-Aqueous Media The Structure and Properties of Liquids Molecular Dynamics of the Chemical Reactions of Gases Electrode Reactions of Organic Compounds Homogeneous Catalysis with Special Reference to Hydrogenation and Bonding in Metallo-Organic Compounds Motions in Molecular Crystals Polymer Solutions The Vitreous State Electrical Conduction in Organic Solids Surface Chemistry of Oxides Reactions of Small Molecules in Excited States The Photoelectron Spectroscopy of Molecules Molecular Beam Scattering Intermediates in Electrochemical Reactions Gels and Gelling Processes Photo-effects in Adsorbed Species Physical Adsorption in Condensed Phases Electron Spectroscopy of Solids and Surfaces Precipitation Potential Energy Surfaces Radiation Effects in Liquids and Solids Ion-Ion and Ion-Solvent Interactions Colloid Stability Structure and Motion in Molecular Liquids Oxidation For current availability of Discussion volumes, see back cover.37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66GENERAL DISCUSSIONS OF THE FARADAY SOCIETY 317 Date 1964 1964 1965 1965 1966 1966 1967 1967 1968 1968 1969 1969 1970 1970 1971 1971 1972 1972 1973 1973 1974 1974 1975 1975 1976 1977 1977 1978 1978 1979 Subject Volume Chemical Reactions in the Atmosphere Dislocations in Solids The Kinetics of Proton Transfer Processes Intermolecular Forces The Role of the Adsorbed State in Heterogeneous Catalysis Colloid Stability in Aqueous and Non-Aqueous Media The Structure and Properties of Liquids Molecular Dynamics of the Chemical Reactions of Gases Electrode Reactions of Organic Compounds Homogeneous Catalysis with Special Reference to Hydrogenation and Bonding in Metallo-Organic Compounds Motions in Molecular Crystals Polymer Solutions The Vitreous State Electrical Conduction in Organic Solids Surface Chemistry of Oxides Reactions of Small Molecules in Excited States The Photoelectron Spectroscopy of Molecules Molecular Beam Scattering Intermediates in Electrochemical Reactions Gels and Gelling Processes Photo-effects in Adsorbed Species Physical Adsorption in Condensed Phases Electron Spectroscopy of Solids and Surfaces Precipitation Potential Energy Surfaces Radiation Effects in Liquids and Solids Ion-Ion and Ion-Solvent Interactions Colloid Stability Structure and Motion in Molecular Liquids Oxidation For current availability of Discussion volumes, see back cover.37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

ISSN:0301-7249    

DOI:10.1039/DC97866FX001

出版商:RSC    

年代:1978

数据来源: RSC

摘要:

GENERAL DISCUSSIONS OF THE FARADAY SOCIETY 317 Date 1964 1964 1965 1965 1966 1966 1967 1967 1968 1968 1969 1969 1970 1970 1971 1971 1972 1972 1973 1973 1974 1974 1975 1975 1976 1977 1977 1978 1978 1979 Subject Volume Chemical Reactions in the Atmosphere Dislocations in Solids The Kinetics of Proton Transfer Processes Intermolecular Forces The Role of the Adsorbed State in Heterogeneous Catalysis Colloid Stability in Aqueous and Non-Aqueous Media The Structure and Properties of Liquids Molecular Dynamics of the Chemical Reactions of Gases Electrode Reactions of Organic Compounds Homogeneous Catalysis with Special Reference to Hydrogenation and Bonding in Metallo-Organic Compounds Motions in Molecular Crystals Polymer Solutions The Vitreous State Electrical Conduction in Organic Solids Surface Chemistry of Oxides Reactions of Small Molecules in Excited States The Photoelectron Spectroscopy of Molecules Molecular Beam Scattering Intermediates in Electrochemical Reactions Gels and Gelling Processes Photo-effects in Adsorbed Species Physical Adsorption in Condensed Phases Electron Spectroscopy of Solids and Surfaces Precipitation Potential Energy Surfaces Radiation Effects in Liquids and Solids Ion-Ion and Ion-Solvent Interactions Colloid Stability Structure and Motion in Molecular Liquids Oxidation For current availability of Discussion volumes, see back cover.37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66GENERAL DISCUSSIONS OF THE FARADAY SOCIETY 317 Date 1964 1964 1965 1965 1966 1966 1967 1967 1968 1968 1969 1969 1970 1970 1971 1971 1972 1972 1973 1973 1974 1974 1975 1975 1976 1977 1977 1978 1978 1979 Subject Volume Chemical Reactions in the Atmosphere Dislocations in Solids The Kinetics of Proton Transfer Processes Intermolecular Forces The Role of the Adsorbed State in Heterogeneous Catalysis Colloid Stability in Aqueous and Non-Aqueous Media The Structure and Properties of Liquids Molecular Dynamics of the Chemical Reactions of Gases Electrode Reactions of Organic Compounds Homogeneous Catalysis with Special Reference to Hydrogenation and Bonding in Metallo-Organic Compounds Motions in Molecular Crystals Polymer Solutions The Vitreous State Electrical Conduction in Organic Solids Surface Chemistry of Oxides Reactions of Small Molecules in Excited States The Photoelectron Spectroscopy of Molecules Molecular Beam Scattering Intermediates in Electrochemical Reactions Gels and Gelling Processes Photo-effects in Adsorbed Species Physical Adsorption in Condensed Phases Electron Spectroscopy of Solids and Surfaces Precipitation Potential Energy Surfaces Radiation Effects in Liquids and Solids Ion-Ion and Ion-Solvent Interactions Colloid Stability Structure and Motion in Molecular Liquids Oxidation For current availability of Discussion volumes, see back cover.37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

ISSN:0301-7249    

DOI:10.1039/DC97866BX003

出版商:RSC    

年代:1978

数据来源: RSC

摘要:

Structure and Dynamics of Diatomic Molecular Fluids BY PETER A. EGELSTAFF Physics Department, University of Guelph, Guelph, Ontario N l G 2W I , Canada Our basic knowledge about molecular fluids is in a much less satisfactory state than that for atomic fluids, and the historical development of these two topics is reviewed. Reasons for the slow develop- ment of our understanding of molecular fluids are discussed. The significance of several kinds of ex- perimental data is considered, and the importance of improving the quality of measurable correlation functions is emphasised. It is concluded that a significant underdeveloped experimental field is that of high precision gas diffraction at low to moderate densities and various temperatures, and some preliminary results are given to illustrate this point.The current status of work on fluid nitrogen and hydrogen chloride is reviewed, and it is concluded that the properties of nitrogen are better understood on the basis of recent models than those of HCl. In spite of the insoluble features of the molecular problem, a great improvement should be possible by the patient application of modern techniques. It is an honour and privilege to be asked by the Faraday Division of the Chemical Society to present the Spiers Memorial Lecture at this General Discussion on “ Struc- ture and Motion in Molecular Liquids ”. But at this time an extra responsibility falls on an introductory speaker because the subject has reached a turning point. Conse- quently the initial part of my lecture will be a philosophical overview from an experi- mentalist’s standpoint, in which the fundamental problems facing him will be discussed.Then I shall try to see how the present base might be extended over the next few years and will conclude with a brief review of two types of diatomic molecular fluids. HISTORICAL INTRODUCTION The study of fluids is divided historically into two parts: the static and the dynamic properties. The modern treatment of the static properties followed the measurement of the equation of state of carbon dioxide by Andrews in 1869 and the general explana- tion of the three states of matter provided by the van der Waals equation of state in 1873. The dynamics of fluids really began with the work of Clausius, Maxwell and Boltzmann on the kinetic theory of gases, and their explanation of the experimental behaviour of the transport coefficients of dilute gases.At the same time experimental studies of Brownian motion were carried out, and a theoretical treatment was given by Einstein in 1905. The underlying behaviour of the atoms and molecules in the fluid was being described by interaction potentials. van der Waals’ potential consisted of a hard core and a long range attractive tail of unknown origin. Later potentials due to Sutherland (1 893) or Mie (1 903), for example, were more realistic but progress in using them was difficult due to mathematical difficulties. In the 1920s the work of Enskog and collaborators extended the kinetic theory of gases to higher densities so that transport coefficients at moderate fluid densities could be discussed.But a proper understanding of the interactions between the atoms and molecules was delayed until the development of quantum mechanics. By the thirties the properties of gases were8 SPIERS MEMORIAL LECTURE understood in both a general way and semi-quantitative way, but liquids were still regarded as being a problem area. However, experimental techniques had improved and equilibrium and transport properties of many fluids over a large part of the equa- tion of state were becoming available: Bridgman, in particular, extended the range studied. In addition Zernike and Prins used X-ray diffraction to examine the microscopic structure and determine g ( r ) experimentally. The work of Kirkwood and Yvon in developing the pair theory of fluids, and the density expansions of Mayer and Mayer, together established the framework for many developments in the theory of liquids.Static properties could then be discussed within the framework of formal statistical mechanical results and through a growing number of approximate integral equations for the pair Correlation function. In parallel the statistical mechanics of non-equili brium phenomena was developed, and general application of time correlation functions was accepted. Finally the van der Waals picture was restated as a modern perturbation theory using the hard sphere system as a reference. The combination of computer simulation techniques, modern versions of old theories and advanced experimental techniques have established a realistic foundation for the physics of fluids. From the overall point of view we have seen an initial attempt starting about 100 years ago to understand the properties of dense gases in some detail, and after its success an extension to liquids.In the past twenty years this has led to a fairly satisfactory general understanding of fluids at all densities, especially for atomic fluids. To obtain a fully quantitative treatment may require a renewed attempt to understand gases before the study of liquids can be taken further. We examine this question in the case of atomic fluids. BRIEF REVIEW OF ATOMIC FLUIDS The theoretical results and the results of the computer simulation work have been used together to develop a theory of fluids, and for its proper application to real atomic fluids we need to know the pair potential with high precision and to have some information about the short ranged three-body potential and the three-body correla- tion function.To obtain data of this kind, one probably has to go back to the study of gases. Data on the scattering of atoms by atoms, on the energy levels of the noble gas dimers, on second virial coefficients and gas viscosities have provided enough infor- mation to obtain good pair There are practical difficulties connected with the fact that this method is indirect, since the parameters in a model are optimised. This limits our firm knowledge about the pair potential, but well-depths and size parameters are believed to be known to z 1 % due to errors in the primary data and there is probably an equal uncertainty due to the indirect method. This method uses diffraction data [i.e., the structure factor S(q)], and although experiment- ally difficult it is slowly being made practicable.The expansion of the pair correlation function, g(r), as a power series in the density may be written : Fortunately there is in addition a direct method of obtaining the potential. S-'(q) = (1 + [g(r) - 11 e'er dr1-l = 1 - p / f ( r ) ei+r dr 1 - p2 11 eiqWr dr ds { f ( r ) f ( s ) f ( r - s) where f ( r ) = e-Bu2(') - 1, and ut(r), u3(r, s) are the pair and three body potentials, respectively. This equation is analogous to the virial series through which P,V,TP. A . EGELSTAFF 9 FIG. ].-Structure factor for krypton gas at 297 K and a density of 0.316 x lozz atoms ~ m - ~ .A smooth line has been drawn through the experimental points which were measured at steps of 0.025 A-' from q = 0.15 to 4.0 A-', and were determined absolutely with an uncertainty of 0.005. The dashed line is a Monte Carlo simulation for pair interactions only, and the dotted curves are inter- polations between q = 0 and the data. data are frequently interpreted. The coefficient of p is the fourier transform of exp - pu(r), and so u(r) can, in principle, be measured directly. For such an interpreta- tion it is necessary to determine the complete fluid structure factor for a large number of state points along a single isotherm, each structure factor being determined to an accuracy of z+%. Preliminary results of this kind have been obtained for krypton gas at room temperature by Teistma and Egelstaff2 in the density range (0.02-0.6) x atoms ~ m - ~ , and an example is given in fig.1 compared with a Monte Carlo calculation using the pair potential only. The difference between the two curves may -30 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 density I lo2* atoms FIG. 2.-Comparison between the direct correlation function at q = 1.33 k' as a function of p and data predicted from eqn ( 1 ) using the previously known pair and three body potentia1s.j 0 , experimental data; - - -, pair potential only; -, pair and Axilrod-Teller potential.10 SPIERS MEMORIAL LECTURE be due to errors in the pair potential or to the three body terms. These effects can be p for each value of q, since only two terms from the series (1) are found to be needed for q > 0.4.For q = 1.33 A-', where the coefficient of p in (1) is nearly zero; such data are given in fig. 2. It is evident that both the slope and the intercept at p = 0 are not predicted correctly, so that both effects are important. The slope is related to the three-body potential and hence models of three-body forces may be tested. In particular, it may be shown that the Axilrod-Teller triple dipole term does not account for the slope in fig. 2. In addition it may become tech- nically possible to scatter atoms from dimers, so providing a further check on models of three-body forces. The evidence in either case is " circumstantial " since insufficient variables are measured to give the complete function. But for a higher order term this may be acceptable.With accurate knowledge of potentials, the density dependence of the pair correlation function can be predicted quantitatively and compared to ex- periments, and this checks our knowledge of higher order function^.^ In the case of dynamic phenomena the situation is not so satisfactory but progress has been made experimentally through more extensive measurements of transport coefficients and through inelastic neutron scattering experiments which measure the Van Hove space-time correlation function. At this stage of development, computer simulation of dynamic behaviour has yielded many properties which cannot be meas- ured directly. There have been a number of papers on the kinetic theory of dense fluids which attempt to derive space-time correlation functions from first principles, and the results compared with experimental data and computer simulations with moderate success.This favourable situation might be further improved by new types of scattering experiments which will become possible, when new sources of radiation (e.g., syn- chrotron radiation) are generally available. In the case of mixtures of different kinds of atoms the situation is also reasonably satisfactory. From atom-atom scattering experiments one can obtain the pair potential for unlike atoms. In principle this same information can be obtained from the partial diffraction data. The technique of isotopic substitution with neutron diffraction experiments on favourable combinations of atoms, can yield all the partial structure factors [sab(q)], i.e., separated by plotting the direct correlation function, c(q) = p-"l - S-' (d1, against where a, b denote two species of atoms.is exp - puab(r). The cross section measured in a diffraction experiment is : At low density the leading term in &(r) where c,, b, are the concentration and scattering amplitude of species, a, respectively. Several measurements are required in which the values of b,, bb are varied by isotopic substitution, so that &b(q) can be obtained algebraically. However if p - 0, the concentration can be used as a variable. In practice such experiments would be very time-consuming and so have not yet been carried out. Thus we can in general expect both for pure and mixed atomic fluids to be able to (essentially) completely solve our problem: that is to describe all the properties of the system in a fully quantitative manner, using the correct underlying information and the theoretical results.P.A . EGELSTAFF 11 It is probably fair to say that a successful study of fluids requires two distinct phases. There is an initial phase involving crude potentials and physical models where one starts with low density gases and moves in understanding towards high density liquids. In the case of atomic liquids this phase is complete. The second phase, beginning again with low density gases, is in progress but this time the inter- atomic potentials are measured with great care, and via more sophisticated theories and precise measurements of correlation functions the properties of the dense fluid will eventually be fully understood. THE MOLECULAR PROBLEM “ Even the clearest and most perfect circumstantial evidence is likely to be at fault, after all, and therefore ought to be received with great caution.” - Mark Twain, Pudd’nhead Wilson’s Calendar.We contrast a fluid composed of atoms with one of rigid molecules, and show The reason is that we have moved from a soluble problem to an insoluble one. FIG. 3.-Ideal molecule-molecule scattering experiment in which the orientation (C) of the molecules in the beams and target are assumed to be known throughout. This idealistic experiment method is required if u(R1; 52,; 52,) is to be measured in full. The following information is collected: incident beam alignment, incident beam direction, incident beam velocity, target alignment, target molecule movement (if any), scattered beam alignment, scattered beam direction, scattered beam velocity.simply that the number of variables has escalated beyond our capacity to measure them. In the case of a diatomic molecule there is a natural axis along the line joining the centres of the two atoms and in addition to measuring the position of the molecule we must measure the direction of this axis. To measure the angular depend- ent pair potential would require an ideal molecule-molecule scattering experiment of the type shown in fig. 3. If we refer to a particular direction for the axis as a polarisa- tion (c) of the molecule, then we require a polarized beam of mono-energetic molecules to fall on a polarized target, and we have to study the number of molecules scattered at a given angle with a given polarization.Obviously this experiment is beyond the realm of possibility at the present time. However the quantities that can be measured, such as cross-sections for molecule-molecule scattering, properties of dimers and12 SPIERS MEMORIAL LECTURE second virial coefficients, can be used in an indirect way to test models and optimize their parameters. The diffraction experiments do not contain enough information because the results are averages over the orientations of the molecules. An ideal experiment would involve a rigid molecule containing a nearly infinite number of different atoms, so that in principle a nearly infinite number of partial structure factors could be measured. The results of this hypothetical experiment could be used to deduce some spherical harmonic components of the general pair correlation function but not all of them.These data are related to the correlation function g(R12; 0,; SZ,), or the normalised probability of finding molecule 1 at R, in orientation SZ, and molecule 2 at R2 in orientation a,. It is often expressed as a sum of spherical harmonics : several differ- ent systems have been used, namely the polar axis taken along Rlzs or fixed relative to laboratory space' or taken along the scattering vector 4,9910 and the latter choice advantageous in the present application. State selected scattering data will be useful in the future. Thus we write: 1 dR eiqmR [g(R; 0,; Q,) - 11 = h(g; 0,; 0 , ) = 2 hl1lzrnk) YlI,ffll(Ql) YlZ.ffll(Q2)~ (4) Il,[a,m Since the partial structure factor h,b(q) is: it can be shown that; where the direction of the molecular axis is from 6 to a, r,, r b are the distances from the (arbitrary) molecular centre to atoms a and 6, respectively, i, j label two molecules, j , is a spherical Bessel function and S~~""(q) is the partial structure factor for the molecular structure.In principle if we have a nearly infinite number of h,b(q) we might solve to find the harmonics hl,,,o(q), but this set does not include the harmonics for m # 0. Harye and Stell" say that this corresponds to leaving out the effect of " twisting forces ". (The same problem would apply to the atom-atom model dis- cussed in the next section.) In practice for homonuclear diatomic molecules (e.g., N,) there is only one hub, for heteronuclear molecules (e.g., HC1) there are three and for a trinuclear asymmetric molecule (e.g., HCN) there are six.Thus eqn (6) is not usually soluble, and the in- formation content is least for the homonuclear diatomic case. Once again the data can be used to test models and optimise parameters only. Thus our knowledge about the pair interaction and the three-body effects in molecu- lar fluids is much less satisfactory than in the case of atomic fluids. This field is more dependent on models and less dependent on underlying ideas, so that it is in a less de- veloped state. The basic information about molecular interactions is obtained from " circumstantial evidence " and there is a parallel between the ways of studying the three-body potential for atomic fluids and those for the pair potential of molecular fluids.AVAILABLE DATA A N D MODEL POTENTIALS Experimental data may be divided into four classes: the macroscopic and the Macro- microscopic data and both cases include static and dynamic properties.P . A . EGELSTAFF 13 scopic data such as the thermodynamic properties of the system, the dielectric constant for a polar fluid and similar properties may be expressed as volume integrals over certain microscopic correlation functions. For example, the isothermal compressibility ( x T ) (which would include the P, V,T data) is: p k T X T = + phab(q)] = ti%O L1 -k phaa(4)1 (7) and the dielectric constant spacing (d) is: for a non-polarizable rigid diatomic of internuclear where p is the dipole moment and ha&) is defined in eqn (5).In piinciple all the macroscopic data can be represented by single points on the appropriate microscopic functions, and (apart from the practical value of P,V,T data for example in allowing a state to be identified) the virtues of measuring macroscopic quantities lie in the pre- cision with which they can be determined and the large number of states which can be covered. However, in some cases, data of this type have been available for some years and have been extensively analysed (see for example Hirschfelder, Curtis and Bird).I3 It seems likely that the weaknesses in the study of molecular fluids will not be resolved by better macroscopic measurements, because these data cannot be used to define the potential unless its analytic form is predetermined.But if a microscopic correlation function could be measured to the same high precision and at as many states as its macroscopic counterpart, these data would contain new information of a potentially useful kind. In eqn (7) and (8) the relevant quantity is the partial structure factor ha&) and if measured many new fourier components, in addition to q = 0, would become available. This represents new and potentially more powerful (ak- though incomplete) information. These functions may be measured in suitable cases by the isotopic substitution technique of neutron diffraction or in others by combin- ing X-ray, neutron and electron diffraction data, but technical difficulties have ham- pered work on ha&) for many years. We shall describe later a new attempt to obtain extensive high quality data on these functions.At this point it should be noted that, by analogy with the atomic case discussed at eqn (l), data on gases as a function of density would be more informative than liquid data. The macroscopic dynamic quantities, such as the coefficients of diffusion and viscosity, may be expressed as time integrals over microscopic correlation functions : these functions include velocity correlation functions, stress correlation functions, orientational correlation functions, the van Hove space-time correlation functions, and so on. However in most cases microscopic dynamic functions cannot be meas- ured with high precision, and at present the macroscopic quantities are usually avail- able with a precision which may be ten times superior to that of the microscopic ones. In addition many more state points are covered by the macroscopic experi- ments. Thus (historically) information about intermolecular potentials has been deduced by fitting models to transport coefficients.For example, Monchick and assumed that HCl molecules interacted according to a Stockmeyer potential and obtained its parameters by fitting to the gas viscosity. However in the writer’s opinion this method can give only a rough idea of the real potential unless its ana- lytical form has been predetermined. Unfortunately this cannot be done at present in a convincing way. The van Hove and the orientational correlation functions are not related in a simple way to the potential.Usually a model is used which contains some parameters or14 SPIERS MEMORIAL LECTURE functions which may be adjusted to fit the experimental spectra. Frequently these quantities can be written in terms of spectral frequency moments, that is static quantities which may be expressed in terms of the potentials and the static correlation functions. As a simple example we consider the dynamic (self) structure factor observed in incoherent inelastic neutron scattering and the depolarized rotational Raman spectrum of a homonuclear diatomic fluid. The ratios of the fourth to the second frequency moments for these spectra are (in the classical limit) : 32 incoherent neutron depolarized Raman - - 48 kT ( r 2 ) 3 I + 2 1 k T where M , I are the molecular mass and moment of inertia, respectively. By combining these results the mean square torque (T~), and the mean square force (F2) may be deduced in principle.However, the spectra must be determined with high precision well into the wings for this purpose. The values of ( T ~ ) and ( F 2 ) may be predicted from pair potential models using standard formulae. The fourth moment of the total dynamic correlation function, S(q, co), is a com- plicated function of q, and can be written in terms of a pair potential model. It is a good quantity to use for comparison of theory with experiment, but is hard to measure accurately. In both the latter and former cases, it is better to calculate the full fre- quency spectra by computer simulation which are then compared directly with the observations.Recent pqpers by Hills" connect the i.r. spectrum with the van Hove correlation functions and the harmonics of the potential. Again these results would be hard to use except as part of a consistent test procedure. For these reasons it seems desirable to employ the static structure factor data as the most convenient func- tions from which to work back towards intermolecular forces. We may write a general model for the pair potential in this form: where a, b label the atoms on different molecules. The first term in this equation is an atom-atom interaction term and the second term represents a multipolar expansion about the centre of the molecule. The models used, for the first term, range from centre to centre hard spheres or Lennard-Jones potentials through hard shapes and atom-atom Lennard-Jones potentials.For the multipolar terms for linear molecules some of the dipole-dipole, quadrupole-quadrupole, quadrupole-dipole terms are usually employed with measured dipole and quadrupole moments. In addition for some models charges are added near the atomic sites and these terms are then in- cluded in the atom-atom part of the potential with the charge distribution arranged to simulate the multipolar terms. Sometimes these potentials are simplified to represent just the shape of the molecule, giving the so called " hard object " model. A particu- lar example is provided by the simulation of the molecular shape by a set of fused hard spheres, since the nearest neighbour correlations are often dominated by the molecular shape and so can be represented by a model of this kind.Frequently more than one model has been used in the literature to fit the same data. It is not clear whether the L-J form is accurate, how the multipolar terms should behave at small values of RI2 and so on. In addition to these difficulties the polarizability of the molecule can lead to collec- tive effects. Wertheim16 has discussed the effect of polarizability on the free energy of a fluid of polar polarizable hard spheres. For an example in which the strongly A proper analytic form for u(RI2, Ql, Q2) is not really known.P . A . EGELSTAFF 15 polar molecules are polarizable in a direction parallel to the dipole moment he finds a six-fold increase in the free energy of the dense liquid. Such effects will vary con- siderably from case to case and much more work is required.In the case of an atom-atom potential, the set of U,b(r) can be determined uniquely from the measured set of partial structure factors. This method is analogous to the method suggested for mixtures, but requires a more complicated relationship l7 between &(r) and uab(r). Of course if the real potential is not an atom-atom poten- tial then the results would depend upon the state of the fluid. It may become possible to determine effective partial pair potentials in this way and to analyse their variation with density and temperature to deduce more information about the real potential, Unfortunately precise measured correlation functions are not available in most cases. although the structures of liquid and gaseous nitrogen have been compared.'* In evaluating properties from a potential two approaches have been used. First there is the method of computer simulation where a number of molecules, usually between 100 and 500, are placed in a box and the positions and motions of these molecules followed on a computer.Within the limitations of the budget and the size of the computer, almost any property may be evaluated in this way. Secondly various theoretical methods have been used, for example the RISM equation for pair correlations, and the thermodynamic perturbation theory. These simulation and theoretical methods are fully covered in other reviews and will not be discussed further here. However in many cases the comparison of accurate measurements with accurate predictions will require quantum corrections to the latter.The first order (h) quantum correction of the van Hove correlation function is the leading term of its imaginary part and it is readily included by ensuring that the frequency spectra satisfy detailed ba1an~e.I~ The h2 correction is the leading correction to the real part of G(r, t ) , and therefore is the important correction to g ( r ) . A number of papers at this conference are concerned with improvements to theoretical predictions, and the experimentalist can assume that, in principle, the con- sequences of a theoretical potential model can be calculated reliably in the long run. We conclude therefore that an important aim for the experimentalist should be to improve both the precision of microscopic observations and the coverage of different states to the same level as now available for macroscopic ones.This will be a diffi- cult and time-consuming task, although the information content of accessible data can be improved substantially by this approach. The experimentalist may well ask " what is the best starting point? " In the writer's opinion the answer is, the partial (static) correlation functions of low density diatomic molecular gases.* For cases such as N2 or O2 there are some data available, but for AB molecules there is none, apart from a preliminary result described later. A new diffractometer for gas samples will be described in the next section. At liquid densities several recent publications,20-22 as well as papers at this conference, report results on partial structure factors which illustrate significant improvements over earlier work.In the balance of this review we will discuss the properties of two diatomic molecular fluids, namely nitrogen and hydrogen chloride. NEUTRON DIFFRACTOMETER FOR GAS SAMPLES Neutron diffraction is a useful technique for gas diffraction work, because neutrons readily penetrate the walls of pressdre vessels and the isotopic substitution method allows a number of different diffraction patterns to be measured. For each isotopic * Pioneering work in this field was done 30 years ago by N. Z. Alcock and D. G. Hurst, Phys. Reu., 1949, 75, 1609; 1951 , 83, 1100.16 SPIERS MEMORIAL LECTURE sample the diffraction pattern at a very low density is measured and subtracted from that at the density of interest, so that partial structure factors, h&(q), may be extracted algebraically.This step removes not only the S:F""(q) term but various other unwanted effects, such as incoherent scattering and the inelasticity correction. It is necessary, however, to design and build a specialised instrument if many states are to be covered with high precision in a reasonable time. We describe here an instrument -recently completed at the N.R.U. reactor by the University of Guelph neutron scat- tering group. The most intense part of ha&) occurs at q < 4 A-' for many gases, and it was de- cided to operate at 2.4 A and scattering angles up to 120" to cover this range in detail. e \ Y Q FIG. 4.-Layout of the University of Guelph neutron diffractometer at the N.R.U.reactor, Chalk River. The items lettered are: ( a ) beam shield, (6) graphite monochromator, (c) graphite filter, ( d ) alu- minium monochromator (rotatable), ( e ) sample environment, (f) counter shields, (g) counters. Fig. 4 is a layout of the instrument: the reactor beam is reflected twice to separate the " wanted " from the " unwanted " radiation. The first monochromator (graphite) reflects both 2.4 and 1.2 A neutrons intensely, and the second order 1.2 A neutrons are removed with a graphite filter (using a different arrangement for the second monochromator it is possible to select the 1.2 A beam so doubling the q range, but this facility is not used normally). The beam falling on the sample has AA = 0.05 A, A0 = 0.6", an area 5 cm wide by 8 cm high and an intensity of -4 x lo5 neutrons cm-2, but the beam area and intensity used in any experiment are determined by a series of diaphragms.A total of 18 detectors are placed at 2.2 m from the sample in three banks: four small (1.3 cm diameter) detectors are used at angles from 2" up- wards, six larger ones over the central range and eight 45 cm high x 5 cm wide detec- tors for angles >60". All detectors are recorded automatically and the three rangesP . A . EGELSTAFF 17 are covered simultaneously so reducing the length of time to record one pattern. Usually the detectors are stepped along in 0.5" intervals so that each detector yields a complete record of its section of the pattern, and the results for several detectors of the same kind are averaged to improve counting statistics.Since the gas pressure vessels are cylindrical there is no error in simultaneously recording different sections of the pattern. The krypton data shown in fig. 1 were taken with this instrument, and after averaging the points in each 0.025 A-1 interval the statistical error was 0.1 5%. This is less than, or of the order of, the error arising from uncertainties in the experimental corrections, giving a combined error of -0.5%. This instrument may also be used for inelastic scattering, by the time-of-flight method, since the second crystal may be rotated at 120 Hz. Thus the advantages of an instrument of this kind lie in the optimisation of the design for gas diffraction which permits a good quality result to be obtained after 30 h of counting. Because the instrument is dedicated to this work long experiments can be undertaken and systematic errors controlled so that it is possible to work along an isotherm in density steps of ~ 1 0 % .Normally the recording of data on back- grounds, empty vessels, reference materials such as vanadium and so on will involve about another 10 scans. Since the density range below the critical density requires z 10 density steps, this range may be covered in "30 days of reactor time for an atomic gas or z 75 days for a molecular gas such as HC1 where four different isotopic samples are used. (The whole experiment may take twice as long due to setting up time and reactor shutdowns.) This procedure should then be repeated on other isotherms; however, a program of this kind is limited by specimen construction time and data analysis time.The University of Guelph group is presently studying krypton and xenon at room temperature, nitrogen at various temperatures and hydrogen chloride at 80 "C. Preliminary results for nitrogen and hydrogen chloride will be given in the next two sect ions . This experience suggests that it is now technically possible to measure h,b(q) at many points along a number of isotherms. FLUID NITROGEN The nitrogen molecule is often treated as a rigid dumb-bell shaped body. Since the internuclear distance d < the molecular size, it is not far from spherical. Moreover, its quadrupole moment is known, and the polarizability is not high so that many- body forces are assumed to be small.Thus fair guesses at models for the inter- molecular potential may be made and we can refine the model parameters by optimis- ing the fit of predicted properties to experimental data. If a wide range of data can be fitted, then at least we have an empirically successful model and possibly are close to the real potential. Significant work of this kind has been done for nitrogen and will be reviewed below. The thermodynamic data and liquid structure factor have been calculated by the molecular dynamics simulation method with an atom-atom pair potential of the Lennard-Jones (12-6) type, by Barojas, Levesque and Q ~ e n t r e c ~ ~ and by Cheung and Powles 24a using slightly different parameters. Also Powles and Gubbins 24b have calculated the second virial coefficients for both sets of parameters.In all cases a classical method was used and no higher order forces were included: a successful fit was obtained to the available experimental data. In a later paper Cheung and Powles 24c added a quadrupole-quadrupole term to their potential which gave slightly better predictions. Table 1 lists the parameters of the two potentials, and C. P.18 SPIERS MEMORIAL LECTURE TABLE 1 .-PARAMETERS USED FOR NITROGEN ATOM-ATOM POTENTIAL Barojas et al. 3.341 44.0 1.100 (B. L. Q .) Cheung and Powles 3.31 37.3 1.090 (C.P.) obtain a reasonable fit to a wider range of data than B.L.Q. and so claim that their potential is significantly better. The intersite distance, d, is not forced by C.P. to be the same as the internuclear distance, but since it is found to be nearly the same this is regarded as supporting evidence for the atom-atom potential.Additional work of this kind has been done by Quentrec and B r ~ t , ~ ~ Streett and Tildesley2'j and by Mac- Rury et ~ 1 . ~ ~ I q1A-l FIG. 5.--Structure factor for nitrogen gasz8 at 297 K and a density of 0.473 x molecules cm-3 compared with the prediction^^^ of the C.P. (-) and B.L.Q. (- - -) potentials. The structure factor for the low density gas at room temperature and at 200 K has been measured recently by Sullivan and Egelstaff,28 and predictions for the two potentials of table 1 have been made by Johnson29 on the basis of the RISM method. One of these results is shown in fig. 5 , where the predicted and measured quantities are compared.There is an important discrepancy in the region q < 1.6 A-' and a smaller one near q - 4 A-'. If these effects are interpreted according to eqn (6) it is necessary to compare the behaviour of the coefficientsj,,j,, as a function of q. In this way we see that for q < 1.6 A-' the zero-order term hooo(q) is probably in error, while for q - 4 A-' the terms involving 1J2 = 1 or 2 will be in error. Both C.P. and B.L.Q. claim that their models fit the liquid structure and in this case the data of fig. 5 would imply that the model parameters vary with state conditions. However, neither C.P. nor B.L.Q. made quantum corrections to their calculations and the comparison with experiment did not extend below 0.8 A-'. Thus this conclusion may be premature. In the gas calculations quantum corrections were unimportant, but this is not so for accurate work on liquid nitrogen.The ti2 correction to the density independent part of g(12) is well known13 namely: g(12,p + 0) = [exp - /3u(12)] [l - W(12)]P . A . EGELSTAFF 2.4 2.2 2.0 1.8 1.6 1.4 - L Y cn 1.2 1.0 0.8 0.6 0.4 0.2 19 - - - - - - - - - - - - where and 412) = u(R12, R,, 0,) is the pair potential. The appropriate expressions for a higher density are not known for molecular fluids, but are for atomic Numeri- cal results for the hard sphere fluid have been given by Gib~on,~' and if we treat N2 as20 SPIERS MEMORIAL LECTURE where the superscript " c " means classical part. In this case, if Sc(q -w 0) is small ( z 0.1 for liquid N,) the correction will be small. From the work of Thompson et al.32 we may evaluate the leading correction term, i.e., p,/ W(12)g(12) dR d o , dR2 (12) and for the C.P.potential with liquid nitrogen at 66 K the result is 0.34. By compari- son the same quantity for neon at 35 K is 0.40. In the latter case this term is increased by (MkT)-l and reduced by the loss of the rotational term: these two effects nearly cancel, making nitrogen and neon similar. The complete result could be calculated \ 1 I I I I 2 4 6 8 1 0 time/ 10-13s FIG. 7 . 4 = 2 reorientational correlation function for liquid nitrogen from ref. (24c). (---) experi- mental result from depolarized Raman spectrum, (- - -)calculated result for C.P. potential with added quadrupole term. if molecular dynamics or Monte Carlo data on the density derivatives of g ( r ) became available.However the size of the known term (0.34) suggests that the full correction will be significant at higher Q in agreement with the conclusions drawn from fig. 6. This may modify the optimum values of the potential parameters. A successful potential should predict dynamic properties also, and Cheung and Powles compare their predictions with reorientational correlation functions deduced from depolarized Raman data. Fig. 7 shows their comparisons for liquid nitrogen (the only data available at the time) both with and without the quadrupole term. Although the latter makes some improvement there is a discrepancy for times -1 ps. De Santis et al.33 measured the depolarized Raman spectrum of nitrogen gas at room temperature. A rotational line spectrum is seen at low pressures but the lines merge into an envelope at about 100 atmospheres.De Santis et al. use a moment analysisP. A . EGELSTAFF 21 to deduce the intermolecular mean square torque [( T ~ ) - eqn (9)] and show that there is an almost linear variation with pressure. It would be useful to compare these results with computer simulation data for the C.P. potential. Another quantity which it is useful to calculate (but requiring more computer time) is the generalised van Hove correlation function34 namely : where i, j run over the molecules in the fluid. angle average over this function or : The dynamic structure factor measured in neutron inelastic scattering is a weighted and F(q, ni) = (2 b, eiQ.rz)/2 b,, where a runs over the atoms in the ith molecule, b, is the scattering amplitude and r,“ is a vector joining the centre of the ith molecule and atom a.This function has been measured by Carneiro and M c T a g ~ e ~ ~ for liquid nitrogen at 66 K and by Hawkins and Egel~taff~~ for nitrogen gas at 297 K. Weiss and Levesque3’ have calculated S(q, w) for liquid nitrogen by molecular dynamics using the B.L.Q. (table 1) potential and a state near to, but not the same as, that of Carneiro and McTague’s experiment. These simulation results disagree in magnitude with the experimental data and this may be due to the difference in the two states. There are also some spectral shape differences at low q which may be due to the poor potential employed, but Weiss and Levesque attribute this to weaknesses in the ex- periment.a a More simulation work is required in this area. FLUID HYDROGEN CHLORIDE Hydrogen chloride differs from nitrogen in possessing a dipole moment, having a lower amount of inertia, being more polarizable and having a different shape. Be- cause of these differences, the anisotropy of the intermolecular potential is larger, the many-body forces are probably larger, the rotational quantum effects are larger, and there are three partial structure factors rather than one. Thus the properties of this fluid are in many ways more interesting than those of nitrogen, but are much more difficult to measure and analyse. Monchick and Mason14 have used a model potential of the Stockmeyer type [that is two terms in u(R, R,, a 2 ) ; a central Lennard-Jones (12-6) potential plus a di- pole-dipole term] in order to calculate the gas viscosity, and the parameters in their potential were then optimised to fit the experimental viscosities.At the time this work was done it was difficult to evaluate the collision integrals, so approximate methods were used. As commented above, this method can yield different results for different model potentials, and a re-evaluation and optimization of the potential parameters would be worthwhile using more terms in eqn (10) and better methods of calculation. Some thermodynamic properties of hydrogen chloride have been calculated by Calado et aZ.38 who used a central Lennard-Jones (12-6) potential plus dipole-dipole, quadrupole-quadrupole and the two dipole-quadrupole terms [that is five terms in eqn (lo)].They obtained a satisfactory fit to the saturated liquid density and the vapour pressure, but did not fit all the P, V, T data available.39 Their potential para-22 SPIERS MEMORIAL LECTURE TABLE 2.-POTENTIAL PARAMETERS FOR HYDROGEN CHLORIDE (&lk)lK p x 10l8 Q x 1026/e.s.u. 018, Monchick and Mason l4 3.36 328 1.08 0 Calado et aZ.38 3.641 191.4 1.07 3.80 meters are shown in table 2 in comparison with those of Monchick and Mason. There is a substantial difference between the two results which may lie partly in the different data employed and partly in the (simplified) analysis used by Monchick and Mason. However, more work will be needed before the potential for HC1 is known even as well as that for nitrogen. 0 1 2 3 4 w A-1 FIG.8.--Intermolecular structure factors4' for three different isotopic samples of HC1 gas. (a) mixture of HCI and DC1 such that c H ~ H + c D b D = 0, (6) HCI, (c) DCI. In the case of HCl there are three different partial structure factors and these have been measured in a preliminary experiment by Soper and Egelstaff40 for HCl gas at 80 "C and 100 atm near the principal peak at q - 1.8 A-' (the diffractometer used was a primitive version of that shown in fig. 4 but using a wavelength of 1.09 A). The results are shown in fig. 8 : the quantity D,(q) = S(q) - S,(q) has been plotted where S,(q) is the result of the low pressure experiment. The partial h,,(q) deduced from these data are given in fig. 9, and using CI as a centre an approximate prediction for the h,,,(q) term is given also.Similar results are obtained for either of the poten-P. A . EGELSTAFF 23 tials in table 2. It can be seen that there is a substantial discrepancy, particularly for the hydrogen-hydrogen correlation function which reflects the higher harmonics in eqn (6). A rough estimate suggests that this difference could not be given correctly by the above potentials. It is well known, for example Rank et aL4I and Lavor et aZ.,42 that dimers can be found in HCl gas under suitable conditions. In principle the correct potential should predict the proper formation of the dimer; predictions of +0.4 +0.3 +0.2 + 0.1 0 0- -0.1 -0.2 G & 2 -0.3 c4 Y -a c .c c $ 0 V s" -0.1 s -0.2 v1 - b -0.3 1 2 3 4 I I I I FIG. 9.-Data of fig. 8 separated algebraically into the three intermolecular partial structure factors.(a) H-H, (b) CI-CI, (c) CI-H; (-) experimental; (- - -) calculated hooo(q) term for potential of Calado et aL3* The experimental uncertainty is large, but the relative magnitudes and signs are significant. Curve ( a ) is the least accurate. these potentials could be further tested against the dimer abundance. If this predic- tion is not satisfactory the potential is unlikely to give a reasonable representation of the partial structure factors either. There is a large rotational quantum correction to classical calculations of correlation functions, but it seems likely that this would broaden the H-H correlation function, and so (in this case) is not a reasonable ex- planation of the difference in fig. 9. The rotational behaviour of free and hindered molecules was studied in 1939 by West43 using i.r.techniques. He showed that sharp rotational lines observed for low pressure states along the saturated vapour curve disappeared at a pressure of z 30 atm. A simplified model of pressure broadening due to Gray44 is able to pre-24 S P I E R S MEMORIAL LECTURE dict this effect at least semi-quantitatively. Such results are interesting since they are probably sensitive to the anisotropic part of the potential in the neighbourhood of the distance of closest approach of two molecules. The same effect may be seen in the inelastic scattering of neutrons by HCI, since the observed spectra are related mainly to the van Hove self correlation function for the hydrogen nucleus.Lurie4' has calculated the spectra expected for the free molecule and has shown that spectra for q - 0.1 to 0.5 A-' would be required. Since neutron scattering intensities are readily measured and interpreted it may be useful to study this phenomenon in this way too. meV 100 50 30 20 10 7 5 4 meV FIG. 10.-Inelastic neutron spectrum46 for scattering of 4 A neutrons at 90" from (a) liquid HF at 219 K (m.p. + 27 K) and (b) liquid HC1 at 183 K (m.p. + 24 K). At higher densities the rotational lines merge together to give a broad envelope, and a discussion through rotational correlation functions may be used. Boutin and S a f f ~ r d ~ ~ measured the spectrum of neutrons scattered inelastically by liquid HCl at 183 K. Agrawal, Yip and Gordon4' successfully predicted this spectrum using the I = 1 rotational correlation function which they deduced from the data of West.43 Both the experiment and the analysis were simplified because of the preliminary nature of these investigations.Nevertheless there is a sharp contrast to the spectra observed for HF and HCI as shown in fig. 10. This comparison suggests that the HCl molecules do not form strongly bound groups in the liquid even though many body forces may be quite important. Molecular dynamics simulations of HCl fluids to give these dynamical correlation functions would be useful, and would stimulate further experimentai work. McDonald and Klein 48 report work on liquid HF which is in qualitative agreement with the HF spectrum of fig. 10.P . A . EGELSTAFF 25 SUMMARY The study of real molecular fluids, in contrast to model fluids, is in an unsatisfac- tory state.We do not know how to measure the things we wish to know, and the approximate models are already quite complicated yet do not fit all the data. Investi- gations are generally quite time-consuming and sometimes end up with only a small advance. Even for simple diatomic molecules there is no acceptable way of measuring the full anisotropic pair potential, and the role of many body effects is not understood. It is difficult to make reliable predictions for a wide range of static and dynamic properties ; some thermodynamic properties can be predicted successfully but when dynamic properties or detailed studies of correlation functions a x considered the results are poor. The best that can be hoped for at the present time is to justify em- pirical potentials and to use them to try to discover the significance of many body forces.Two opposing examples, nitrogen and hydrogen chloride, were reviewed to illustrate the range of our problems and the general state of the field. It is suggested that if experimentalists became devoted to measuring partial structure factors and dynamic correlation functions with as high as possible a precision and for many states (especi- ally of the gas), the situation depicted in this review could be improved. Hopefully the precision might approach that with which the macroscopic data have been meas- ured leading to a significant improvement in model potentials. Computer simulation methods do not at present cover this field adequately, some properties are not calcu- lated properly and others take too long to compute. Another major effort would be required to obtain quantum corrected computer simulation results with pair and many body forces included for all the properties of interest.A comparison of these two sets of data would then provide the best forseeable test of our " circumstantial " models of these fluids. (a) J. A. Barker, R. A. Fischer and R. 0. Watts, Mol. Phys., 1971, 21, 657; (b) R. Aziz, J. Chem. Phys., 1976, 65, 490; (c) G. C . Maitland and W. A. Wakeham, Mol. Phys., 1978, 35, 1429. (a) J. A, Barker, R. 0. Watts, J. K. Lee, T. P. Shafer and Y. T. Lee, J. Chem. Phys., 1974,61, 3081 ; (6) R. A. Aziz and W. L. Kocay, Mol. Phys., 1975,30, 857. K. E. Gubbins, C.G. Gray and P. A. Egelstaff, Mol. Phys., 1978, 35, 315. S. Yip, Varenna Summer SchooI on New Directions in Physical Acoustics, ed. D. Sette (1976), p. 55. U. Buck, F. Huisker, H. Pauly and J. Schleusener, J. Chem. Phys., 1978, 68, 3334. This may be proved by an extension of the method used by R. L. Henderson, Phys. Letters, 1975, 49A, 197. R. L. Henderson and C. G. Gray, Canad. J. Phys., 1978,56, 571. C. G. Gray, 1978, unpublished; I am indebted to Dr. Gray for an extensive discussion of this topic. * A. Teitsma and P. A. Egelstaff, 1978, unpublished. a W. A. Steele, J. Chem. Phys., 1963, 39, 3197. lo J. S. H0ye and G. Stell, J. Chem. Phys., 1977, 66, 795. l2 D. Chandler, J. Chem. Phys., 1977, 66, 1 1 3. l3 J. 0. Hirschfelder, C. F. Curtis and R. B. Bird, Molecular Theory ofcases and Liquids (J.Wiley, l4 L. Monchick and E. A. Mason, J . Chem. Phys., 1961,35,1676. l5 B. P. Hills, Mol. Phys., 1978,35,793; B. P. Hills and P. A. Madden, Mol. Phys., 1978,35,807. l6 M. Wertheim, Theory of Polar FluEds-V, 1978, to be published in Mol. Phys. l7 B. M. Ladanyi and D. Chandler, J. Chem. Phys., 1975,62,4308. l8 P. A. Egelstaff, D. Litchinsky, R. McPherson and L. Hahn, Mol. Phys, 1978, 36,445. 'O G. W. Stanton, J. H. Clarke and J. C. Dore, Mol. Phys., 1977, 34, 823. 21 A, H. Narten, R. Agrawal and S. I. Sandler, Mol. Phys., 1978, 35, 1077. New York, 1954). P. Schofield, Phys. Rev. Letters, 1960, 4, 239.26 SPIERS MEMORIAL LECTURE 22 G. Palinkas, E. Kalman and P. Kovacs, Mol. Phys., 1977,34, 525. 23 J. Barojas, D. Levesque and B. Quentrec, Phys. Rev. A, 1973,7, 1092. 24 (a) P. Cheung and J. G. Powles, Mol. Phys., 1975,30,921; (b) J. G. Powles and K. E. Gubbins, Chem. Phys. Letters, 1976,38,405; (c) P. Cheung and J. G. Powles, Mol. Phys., 1976,32, 1383. 25 B. Quentrec and C. Brot, Phys. Rev. A, 1975, 12, 272. 26 W. B. Street and D. Tildesley, Proc. Roy. SOC. A , 1977, 355, 239. 27 T. B. MacRury, W. A. Steele and B. J. Berne, J. Chem. Phys., 1976,64, 1288. 28 J. Sullivan and P. A. Egelstaff, to be published. 29 E. Johnson, to be published. 30 W. G. Gibson, Mol. Phys., 1975, 30, 1 . 32 S. M. Thompson, D. J. Tildesley and W. B. Street, Mol. Phys., 1976, 32, 71 1 . 33 A. De Santis, M. Sampoli, P. Morales and G. Signorelli, MoI. Phys., 1978,35, 1125. 34 P. A. Egelstaff, C. G. Gray, K. E. Gubbins and K. C. Mo, J. Stat. Phys., 1975, 13, 315. 35 K. Carneiro and J. P. McTague, Phys. Rev. A, 1975,11, 1744. 36 R. K. Hawkins and P. A. Egelstaff, Mol. Phys., 1975, 29, 1639. 37 J. J. Weiss and D. Levesque, Phys. Rev. A, 1976, 13, 450. 38 J. G. C. Calado, C. G. Gray, K. E. Gubbins, A. M. F. Palavra, V. A. M. Soares, L. A. K. Staveley and C. H. Twu, J.C.S. Faraday I, 1978, 74, 893. 39 E. U. Franck, M. Brose and K. Mangold, in Progress in International Research on Thermo- dynamic and Transport Properties, ed. J. F. Masi and D. H. Tsai (Academic Press, N.Y., 1962), p. 159. W. G. Gibson, Mol. Phys., 1975, 30, 13. 40 A. K. Soper and P. A. Egelstaff, to be published. 41 D. H. Rank, B. S. Rao and T. A. Wiggins, J. Chem. Phys., 1962,37,2511. 42 M. Lavor, J. P. Houdeau and C. Haensler, Canad. J. Phys., 1978,56, 334. 43 W. West, J. Chem. Phys., 1939, 7, 795. 44 C. G. Gray, J. Chem. Phys., 1974,61,418. 45 N. A. Lurie, J. Chem. Phys., 1967, 46, 352. 46 H. Boutin and G. J. Safford, in Inelastic Scattering of Neutrons (I.A.E.A., Vienna, 1965), vol. 47 R. Agrawal, S. Yip and R. Gordon, Chem. Phys. Letters, 1968,2, 584. 48 I. R. McDonald and M. L. Klein, paper at this Discussion, p. 48. 11.

ISSN:0301-7249    

DOI:10.1039/DC9786600007

出版商:RSC    

年代:1978

数据来源: RSC

摘要:

Fluids of Hard-core Triatomic Molecules Computer Simulations of Polyatomic Molecules. Part 4.-Monte Carlo Simulations of Hard-core Triatomics BY W. B. STREETT School of Chemical Engineering, Cornell University, Ithaca, New York 14853, U.S.A. AND D. J. TILDESLEY Science Research Laboratory, United States Military Academy, West Point, New York 10996, U.S.A. Received 9th May, 1978 Monte Carlo simulations have been carried out for systems of 500 symmetric, linear triatomic molecules composed of fused hard spheres. Orientational ordering in these systems has been exam- ined by calculating two types of distribution functions : (1) angle-dependent pair distribution func- tions, g(rI2co1co2), expanded in spherical harmonics; and (2) site-site distribution functions, gaB(raB), which describe the distribution of atom centres.Pressures calculated from a virial-like equation have been compared with several theoretical estimates based on modifications of the Carnahan-Starling equa- tion. Structure factors, H&), have been calculated by Fourier transforming the gaB(raB), and the results compared with recent X-ray and neutron scattering data for CS2. The calculated gaB(raB) have been compared with values predicted by the reference interaction site model (RISM) equation. 1. INTRODUCTION Computer simulations of hard-core models provide a cornerstone for the develop- ment of theories for the equilibrium properties of molecular liquids. They play the same role in this field that studies of the hard-sphere system played in increasing our understanding of simple 1iquids.l The simulation results provide data on the thermo- dynamic properties and distribution functions of simple model systems which are useful reference systems for perturbation theories.2 They provide stringent tests of statistical mechanical theories of molecular fluids3 and help to elucidate the role played by the molecular shape in determining the structure of such liquid^.^ This is the fourth in a series of papers in which we have examined the properties of molecular liquids, (I,5 11,6 111,’).In this work we have constructed a linear model of three fused hard spheres. The bond length (I), the reduced density [p* = (6/n)pVm, where p = N/V and V, is the volume of the molecule] and the diameters of the two outside atoms (a,) are fixed and we have varied the diameter of the inner atom (a,).Three typical models are shown in fig. 1. The first of these three models (i) corresponds to a hard-core model for liquid carbon disulphide (CS,) at 298 K and 1 atmosphere.8 The other two models have not been chosen to model real liquids, but to examine the effect of changes in molecular geometry on the structure and thermodynamic properties of the fluid. Section 2 con-28 FLUIDS OF HARD-CORE TRIATOMIC MOLECULES tains a brief description of the simulation method. In section 3 we compare the calculated site-site distribution functions (s.d.f.s) of model (i) with those of the RISM integral equation for the same model, and we calculate the X-ray and neutron scatter- ing curves to compare with available experimental data for CS2. In section 4 we calculate the second virial coefficients of the models and use these to make estimates of the pressure of the fluids.lo These estimates of the pressure are compared with the results from the simulation.In section 5 we demonstrate the usefulness of the (ii) (iii) FIG. 1.-Hard-core models of three linear triatomic molecules. The parameters of the models are I* = 0.897, p* = 0.897 and (i) a A / a B = 0.857, (ii) aA/aB = 1.0, (iii) oA/oB = 1.2. In these hard- core models no overlap between molecules is possible. spherical harmonic coefficients in the expansion of the total distribution function, g(r,,co,co,), for calculating properties of molecular fluids." 2. MONTE CARL0 CALCULATION The Monte Carlo simulation has been carried out for systems of 500 molecules using the program described by Streett and Tilde~ley.~ The initial configuration in each run was a face-centred cubic lattice.The simulation proceeded until the root mean square displacement of particles from the original lattice sites exceeded 0.75 oB and continued for 5 x lo6 configurations during which contributions to the ensemble averages were obtained. The site-site distribution functions, ga&), have been calcu- lated from histograms based on 1.7 x lo6 trials; a trial consists of the selection of a particular %-site as the origin, followed by counting the number of sites lying in shells of thickness 6r = 0 . 0 2 5 ~ ~ centred at the origin. In addition we have calculated the following glllam(r12) coefficients in the spherical harmonic expansion of g(r12co,co2) : 1112m = 000, 200, 220, 221, 222, 400, 420, 421, 422, 440, 441, 442, 443, 444, 600, 620, 640, 660, 800, 820, 840, 860 and 880.The maximum error in s.d.f.s is ~ 0 . 5 % and the error in the spherical harmonic coefficients ranges between 1 and 5%, [for a dis- cussion of the error in the coefficients see section 4(b) of ref, (6)].W. B . STREETT AND D . J . TILDESLEY 2 . 5 - 2.0 29 r0 0 0 - 0 0 3 . SITE-SITE DISTRIBUTION FUNCTIONS 30 The site-site distribution functions in a molecular liquid are defined by gadr> = (1/p2)(NN - 1)4rXri - r)> (3.1) where rp is the position vector of site a on molecule i. gQB(r) is the normalized proba- bility that site /3 on molecule 2 is a distance r from site a on molecule 1 at the origin.a D+ 0 I I I I 1 2 3 4 1 /a, FIG. 2 (a)-The three site-site distribution functions for model (i) determined by computer simulation. C marks the position of the cusp in each distribution function; 0, g B B ; 0, g A B and A, g A A . (b) Difference between the simulation results and the RISM integral equation results for the site-site distribution functions of model (i). The simulation results are accurate to 3 % so that differences between the two results that lie outside k0.03 can be attributed to approximations in the RISM The three distinct s.d.f.s for model (i), g B B , g B A and g A A are shown in fig. 2(a). The cusps in these functions are marked C and the origin and location of these features have been discussed elsewhere.’ The comparisons of the computer-generated s.d.f.s equation and its solution. (--) g E B , (- - -1 g B A and (- - - * -) g A A .30 FLUIDS OF HARD-CORE TRIATOMIC MOLECULES with those calculated from the RISM integral equation are presented as a difference plot in fig.2 (b). Since the linear triatomic molecules discussed in this paper are the most asymmetric models for which this comparison has been made, these simulation results provide a most stringent test of the closure approximation in RISM. The g B B s.d.f. (corresponding to a sulphur-sulphur distribution in CS2), produced by the RISM theory has an error of z 20 % at contact, but is elsewhere in good agree- ment with the simulation results. If a fluid of polyatomic molecules is viewed as a mixture of its constituent atoms, then the central integral equation of the RISM approximation is equivalent to the Ornstein-Zernike equation of the mixture.12 The closure relationships used to solve the equation are similar to those used in the Percus-Yevick (PY) approximation for systems of hard spheres,13 and we can think of the RISM equation as a PY-like approxi- mation for molecular liquids.However for hard sphere fluids, the contact values of the distribution functions in the PY approximation were always too low,l4 whereas in this case the RISM contact value is too high. The sign and size of this difference at contact is a sensitive function of density and elongation. We discern no pattern of differences related to molecular geometry. In the case of the g B A and g A A s.d.f.s the simulation results contain secondary peaks on the leading side of the first peak which are not present in the RISM solution.These subdued peaks have been noted before (see section 2.3 of the paper by Hsu et aZ.”); they become more prominent the higher the asymmetry of the model and the higher the packing fraction of the fluid. They are indicative of the failure of the RISM equation to handle the fine detail of the short range angular ordering in the liquid. The s.d.f.s for these models can be Fourier transformed to give X-ray and neutron scattering factors. We follow the notation of Sandler and Nartens who show that the distinct part of the total structure factor for CS,, Hd(k), is given to a good approxima- tion by: X-ray: neutron : Hd(k) = 0.03acc(k) + 0.26acs(k) + 0.71ass(k), Hd(k) = 0.29&k) + 0.50acs(k) + 0.21U&).In our model the partial structure factors a,p(k) are given by (3.4) (3.5) 1 aJ a&) = (4rrp/k)[/ hap(r) sin (kr)rdr + d,qjl(ka,p) , hap(r) = gap(r) - 1 where andj,(x) is the ith order spherical Bessel function. The upper limits in the integra- tions are replaced by 40, and the transforms performed using Filon’s method.16 Assuming the diameter of the sulphur atom (aB) to be 3.5 A, the s.d.f.s for model (i) were transformed to calculate the Hd(k) functions for the X-ray and neutron scattering. These are compared with the experimental data for CS2 in fig. 3 and 4. The overall agreement between X-ray scattering for the model and experiment (fig. 3) is remarkable. In the region below 1 A-’ the results agree to within 10% at k > 0.5 A-’ and to with- in 207< at k = 0 A-’.The Hd(k) are subject to errors for k < 0.7 A-’ because the model s.d.f.s are truncated at r = 40,. There is a disagreement between model and experiment at k > 3 A-‘ where we expect the hard-core model to oscillate more strongly than the real &(k) for CS,; however, the disagreement is slight. The cal-W. B . STREETT AND D. J . TILDESLEY 31 culated Hd(k) have been found to be different from the values obtained from the RISM equation which can be found in fig. 5 of Sandler and Narten.* The differences between the &(k) obtained from the RISM s.d.f.s for model (i) and the Hd(k) from the simulation of the same model can only be explained by the differences in their r- space precursors. There are two independent sets of experimental data available.The results of Suzuki and Egel- staff” disagree with the model calculations in the region of the first peak and for k > 6 A-1. In a personal communication Egelstaff has pointed out that the rescaling of their data within this range results in a substantially improved agreement with the simulation results. The curvature in the experimental results at large k suggests an incomplete The comparison with neutron scattering data is shown in fig. 4. These authors state that their data contain a systematic error of 10%. 0 1 2 3 I 5 6 7 k / A - ’ FIG. 3.-Distinct part of the X-ray scattering from CS2, H,(k). The fourier transform of the simula- tion data from model (i), ____ . Experimental data of Sandler and Narten, 0.reduction of the raw scattering data. The more recent results of Gibson and Dore18 are only preliminary experimental results, but for k > 1 A-’ they are in excellent agreement with the model calculations. The comparisons shown in fig. 3 and 4 have been made from tables of experimental Hd(k) functions. The conclusion is that model (i) is an excellent hard core model for CS2. In the past, evidence of this kind has been interpreted to mean that the liquid structure is determined mainly by short range repulsive force^.^*^*'^ At first sight it is tempting to draw the same conclusion here; however, this may be an oversimplification. In the case of spherical molecules it is straightforward to separate the intermolecular poten- tial into attractive and repulsive parts ; l9 however, for polyatomic molecules the potential is a multidimensional function of distance and angle and the positions and motions of molecules are influenced by torques as well as forces.In our opinion ir has not been clearly demonst rated that these complicated intermolecular potentials can be divided into repulsive and attractive “ parts ” in any simple and consistent way. It becomes difficult to retain notions of shape and repulsive force when one realises that two molecules interact through a field that involves anisotropic overlap,32 FLUIDS OF HARD-CORE TRIATOMIC MOLECULES dispersion, and induced and permanent multipole interactions. We can only observe that this hard-core model effectively reproduces those features of the true intermolecu- lar potential that vary rapidly with separation and orientation, and that these features are primarily responsible for determining the structure of liquid CS2.FIG. 4.-Distinct part of the neutron scattering from CS2, If&). Fourier transform of the simulation data from model (i), ~- . Experimental data of Suzuki and Egelstaff," m. Experimental data of Gibson and Dore,'" 0. 4. VIRIAL PRESSURE The virial pressure of a fluid of hard-core triatomic molecules is given by20 pP/p = 1 + ( 2 ~ ~ / 3 ) { 4 a ~ g ~ B ( a B ) t ~ B ( a , > + 4aiAgAA(aBA)ziA(oBA) + + akiA(aA)TiA(aA)l ( 4 - 1 ) where P is the pressure, p = l/kT, oBA = (UB + 6 A ) / 2 , gai(aaB) is the contact value of the ap s.d.f. and t$(aaB) is the contact value of the function20 TaB(raB) = <r12( arafi/ 8r12)u~u~> shell rap to raB + drab' (4.2) The distance r12 is between the centres of mass of two triatomic molecules, (it is this distance which is scaled in the normal derivation of the virial equation).rap is the site-site distance and the partial derivative is at fixed relative orientation. In the simulation we have obtained values of gafi(rafi)Tas(rafi) in spherical shells of thickness dr = 0.025aB, centred on atom a. We have extrapolated to contact by fitting log [ g a b ( r a b ) ~ a s ( r a b ) ] to a cubic polynomial in rap. The results for the pressure are given in table 1 and they are estimated to be accurate to within 5%. B* is reduced by a factor of four times the molecular volume ( Vm) and is given by In the low density limit eqn (4.1) gives the second virial coefficient B.B* = (2n/3 Vm)[~iXeoo0(~d + [l*el~~(ad/2/3ll + dA{eOOO(aBA) + [z*e100(aBA>/22/3]> + d{eOOO(aA)/4>1, (4.3)W . B . STREETT AND D . J . TILDESLEY 33 where I* = Z/a is the reduced bond length and Y;m(@) is the complex conjugate of the normalized spherical harmonic. The integral, eqn (4.4), is performed by fixing two sites X and Y at contact, oxy, and integrating over all relative orientations. The orientations are defined with respect to the line XY = BB, BA or AA, as shown for a linear triatomic in fig. 5 . These integrals for FIG. 5.--Three coordinate axes systems for the integrations in eqn (4.3). Two atoms are fixed at contact and with this constraint the molecules take on all values of Ox, O2 and qXy.The integrand contributes 1 if the molecules do not overlap and zero if they overlap. the second viral coefficient have been performed numerically using a three dimensional Simpson’s rule quadrature, and the results are shown in table 1. Two methods of TABLE 1 .-THERMODYNAMIC PROPERTIES FOR THE THREE HARD MODELS. (SEE TEXT FOR A FULL DISCUSSION.) model B:,,,, PPIP AAINkT Pk TKT MC BN R BN R BN R (i) 1.2102 14.84 14.36 14.43 5.706 5.733 0.01641 0.01631 (ii) 1.1532 12.84 12.67 12.65 5.085 5.078 0.01897 0.01901 (iii) 1.1137 12.88 12.19 12.61 4.910 5.067 0.01984 0.0190634 FLUIDS OF HARD-CORE TRIATOMIC MOLECULES estimating the pressure of a fluid of hard-core molecules have recently been pro- posed.21v22 They both rely on a hard convex-body equation of state of the form,23 PPIP = [1 + - 2)q + (1 - + (u2/3)1272 - (r2/9)q31/(1 - q ) 3 , (4.5) where the packing fraction q = Vmp and u is the shape factor.The validity of this type of equation is discussed elsewhere23 and the problem remains to fit the shape fac- tor for the models discussed here. For a convex body, within the approximations of scaled particle theory24 v u = RS/Vm, (4.6) where V, is the molecular volume, S is the surface area and I? is (1/4n) x the mean curvature The hard interaction site models used here and in previous papers5*' are generally concave bodies. The volume and surface area of these bodies are well defined, and following Boublik and Nezbeda2l (BN) we chose a corresponding convex body to calculate k. The choice of that body is arbitrary and for the three models here we have used a spherocylinder of elongation I and diameter oB.We now approximate the shape factor using eqn (4.6). The pressure calculated from this shape factor using eqn (4.5) is listed under BN in table 1. Rigby22 (R) has suggested that the shape factor can be calculated from the exact second virial coefficient of the model. From eqn (4.5) we note B" = (u + 1)/4. (4.7) The exact second virial coefficient calculated from eqn (4.3) is used to calculate u which is then used in eqn (4.5) to calculate the pressure listed under R in table 1. Once we have obtained a suitable prescription for u the other thermodynamic properties of the liquid are readily calculated. The residual free energy, AA, is given by AAINkT = [{(u2/9) + U)V - ~q']/(l - q)' - (1 -- (u'/9)} In (1 - q) (4.8) and the isothermal compressibility I C ~ by pkTIcT = (1 - ~ ) ~ / ( l + (2u - 2 ) ~ + (u2 - 2u + 1)q2 - (4u2/9)q3 -+ (u2/9)q4).(4.9) These properties calculated using uBN and aR are listed in table 1 for the three models. The reduced second virial coefficients calculated from eqn (4.3) agree with the reduced second virial coefficients calculated from the geometrical shape factor and eqn (4.7) to within 2%. The pressures calculated from eqn (4.5) using aR and uBN agree with the Monte Carlo results to within the estimated error of the simulation. We can recommend both methods (BN, R) for calculating the pressure, but we cannot determine which is the more accurate because of the difficulty of extrapolating the simulation results to contact in eqn (4.1).5 . SPHERICAL HARMONIC COEFFICIENTS OF THE DISTRIBUTION FUNCTION The spherical harmonic coefficients, glllsm(r12), of the total molecular distribution function, g(r,,co,cu,), can be used to calculate a number of interesting properties of molecular liquids.W . B . STREETT AND D . J . TILDESLEY -0.1 35 - I I I 1 I The isothermal compressibility of the liquid K~ is given by W PkTKT = 1 -I- 4ZP 1 -kooo(312> - 1>r:2 d312. ( 5 4 0 The integrand at large r12 (the upper section of fig. 6) contains considerable scatter and the value of the isothermal compressibility obtained depends critically on where the integral is truncated. The values of p k T ~ , for various truncations of the integral are shown in the lower section of fig. 6 and compared with the value obtained from eqn t (4.9).calculation of from eqn (5.1).26 polarisability which is related to G2. It can be seen that this truncation problem is an important source of error in the A measure of the overall effect of orientation is the effective anisotropic optical P2 is the second-order Legendre polynomial and t,bij is the angle between the axes of molecule i and moleculej. It has been shown that2’ In table 2, values of G2 calculated from eqn (5.3) are compared with those calculated from the ensemble average according to eqn (5.2). Our conclusion is that accurate tabulated functions of the spherical harmonics can be used to calculate G2. The value of G2 calculated from eqn (5.2) is shown in fig. 7 as a function of the number of con- figurations generated for the three models.It can be seen that accurate values of G2 can only be obtained from Monte Carlo runs larger than those usually carried out for homogeneous liquids.36 FLUIDS OF HARD-CORE TRIATOMIC MOLECULES TABLE 2.-Gz VALUES FOR THE THREE MODELS FROM THE ENSEMBLE AVERAGE OF P z ( C 0 S $ij) AND FROM INTEGRALS OVER THE SPHERICAL HARMONIC COEFFICIENTS gZzm(Y12). model G2 integration over g 2 2 m (P2(c0s #u)> (9 0.34 (ii) 0.26 (iii) 0.42 0.30 0.23 0.43 It is possible to consider these models as reference systems in the calculation of the thermodynamic properties of fluids interacting with more realistic intermolecular potentials. For a symmetric linear triatomic molecule the largest contributor to the 0.751 0.5 - 0.25 - I / 1 I I I & 1 2 3 4 5 0 conf i gurat i ons x 1 o6 l + FIG.7.-G2 factor as a function of the number of Monte Carlo configurations. Runs of at least two million configurations with 500 particles are required to obtain reliable estimates of the property. -0- model (i), -El- model (ii) and -A- model (iii). perturbation potential (the difference between the real and reference system potentials) is likely to be the quadrupole-quadrupole interaction. The residual free energy (dA) is written for the non-spherical core/quadrupole system as LfA = LfAo + dA1 + . . . dAo is the residual free energy of the hard-core system and AAl is the residual free energy due to the quadrupole interaction. (5.4) Sandler2* has shown thatW . B . STREETT AND D . J . TLLDESLEY 37 where Q is the quadrupole moment. As an example of the relative importance of the zero and first order term in this perturbation we consider the case of liquid CS, at 298 K and one atmosphere. If we use model (i) as the reference system, AAo/NkT = 5.733.We take the value of the quadrupole moment of CS, from a quantum mechanical calcu- l a t i ~ n , ~ as 2.4515 x Using this value in eqn ( 5 . 9 , AAJNkT = 0.0152, which is two orders of magnitude smaller than the zero-order term. Since carbon disulphide has a small quadrupole moment we would expect the spherical harmonic coefficients of a hard-core plus point quadrupole model to be qualitatively similar to those of the simple hard-core model, and that eqn (5.4) truncated after two terms is a good estimate of the free energy of the model.C m2. I I u 2 3 4 521 4e O ‘ i FIG. 8.-Probability of the parallel orientation(& = O2 = x/2, q12 = 0) as a function of intermolecular separation. At the position of the first two maxima, molecule two has been rotated from 0 to n and the probability shown in the two insets in the figure. This calculation has been performed for model (i). The spherical harmonic coefficients themselves can be recombined to give the total distribution function g(r12cu1w2) according to the equation g ( r 1 2 ~ 1 ~ 2 ) = 4~ 2 g l l l z m (r12) Yl,rn(~l> Y~z-m(~2)- (5.6) lil2rn The convergence of this series and the usefulness of the harmonic coefficients as an efficient method of storing all the information concerning the orientational structure of the liquid have been n ~ t e d .~ . ~ ~ In this work we have fixed the relative orientations of two molecules and examined the probability of this particular orientation as a function of the separation r12. This conditional distribution gWIw,(rl2) has been calculated with 23 terms in the series. In fig. 8 we have plotted gu,u,(rll) for the parallel orientation (0, = 8, = 4 2 , q12 = 0). The insets show the dependence on the azimuthal angle q12, with O1 = 8, = n/2, at distances corresponding to the two peaks in this function. It is readily shown that if Nu,,, is the number of pairs of molecules in the liquid with a particular orientation in a spherical volume of radius rc, then (5.7)38 FLUIDS OF HARD-CORE TRIATOMIC MOLECULES We have calculated N,,,, for a large number of relative orientations out to rc = 4.0 oB where the orientational ordering is random.No overall orientation is preferred in the liquid; that is, NuIu, is practically independent of o1 and 02. Thus the second- ary peaks in the g B A and g B B s.d.f.s noted in section 3 cannot be explained by calculat- ing site-site distances in pairs with preferred orientations, since such orientations do not exist. These intriguing features in the s.d.f.s remain unexplained. We are indebted to I. Gibson and J. Dore of the University of Kent for making available their unpublished scattering data and to the Academic Computer Center of the United States Military Academy, West Point, New York for a grant of com- puting time. We would like to thank Dr I. Aviram, of the Nuclear Research Centre, Negev, for discussions on the virial pressure. This work was supported by a grant from the engineering division of the National Science Foundation U.S.A.For useful reviews see: (a) Physics of Simple liquids, ed. H. N. V. Temperley, J. S. Rowlinson and G. S. Rushbrooke (John Wiley and Sons, New York, 1968); (6) Statistical Mechanics, ed. K. Singer (Special Periodical Rep., The Chemical Society, London), 1973, vol. I. K. C. Mo and K. E. Gubbins, J. Chem. Phys., 1975, 63, 1490. D. Chandler, C. S. Hsu and W. B. Streett, J. Chem. Phys., 1977,66, 5231. L. J. Lowden and D. Chandler, J. Chem. Phys., 1973,59,6587. W. B. Streett and D. J. Tildesley, Proc. Roy. SOC. A., 1976, 348, 485. W. B. Streett and D. J. Tildesley, Proc. Roy. SOC. A, 1977, 355, 239. ' W. B. Streett and D. J. Tildesley, J. Chem. Phys., 1978, 68, 1275. S. I. Sandler and A. H. Narten, Mol. Phys., 1976, 32, 1543. B. M. Ladanyi and D. Chandler, J. Chem. Phys., 1975,62,4308. lo I. Nezbeda and T. Boublik, Czech J. Phys., 1977, B27, 1071. l1 W. A. Steele, J. Chem. Phys., 1963, 39, 3197. l2 J. S. Hoye and G. Stell, preprint of unpublished material. l3 M. S. Wertheim, Phys. Rev. Letters, 1963, 10, 321 ; and section 2 of ref. (4). l4 L. Verlet and J. J. Weiss, Phys. Rev. A, 1972, 5, 939. l6 L. N. G. Filon, Proc. Roy SOC. Edinburgh, 1928, A49, 38. l7 K. Suzuki and P. A. Egelstaff, Canad. J . Phys., 1974, 52, 241. la I. Gibson and J. Dore, preprint of unpublished data. l9 D. Chandler and J. D. Weeks, Phys. Rev. Letters, 1970, 25, 149. 2o I. Aviram, D. J. Tildesley and W. B. Streett, Mol. Phys., 1977, 34, 881. 21 T. Boublik and I. Nezeda, Chem. Phys. Letters, 1977, 46, 315. 22 M. Rigby, Mol. Phys., 1976, 32, 575. 23 T. Boublik, J. Chem. Phys., 1975, 63,4084. 24 R. M. Gibbons, Mol. Phys., 1969, 17, 81. 25 T. Kihara, Rev. Mod. Phys., 1953, 25, 831. 26 R. 0. Watts, J. Chem. Phys., 1969, 50, 4122. 27 W. B. Streett and K. E. Gubbins, Ann. Rev. Phys. Chem., 1977, 28, 373. 2a S. I. Sandler, Mol. Phys., 1974, 28, 1207. 29 C. R. Fisher and P. J. Kemmy, Mol. Phys., 1971,22, 1133. 30 Site-site distribution functions for fluids of hard-core linear triatomic molecules are given in Supplementary Publication No. SUP22598 (24 pp). See Notice to Authors in J.C.S. Faraday, 1977, Index Issue. C. S. Hsu, D. Chandler and L. J. Lowden, Chem. Phys., 1976, 14, 213.

ISSN:0301-7249    

DOI:10.1039/DC9786600027

出版商:RSC    

年代:1978

数据来源: RSC

摘要:

X-ray Diffraction Study of Liquid Neopentane and Tertiary Butyl Alcohol 7 BY A. H. NARTEN Chemistry Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37830, U.S.A. AND S. I. SANDLER AND T. A. RENSI Department of Chemical Engineering, University of Delaware, Newark, Delaware 1971 1, U.S.A. Received 5th May, 1978 The X-ray scattering functions for neopentane (CMe,) at - 17, 25 and 150 "C, and tertiary butyl alcohol (CMe30H) at 26 "C have been measured. These data have been analysed to obtain intra- molecular structures, and by Fourier inversion to obtain an average nearest neighbour hydrogen- bonded oxygen-oxygen separation distance of 2.75 A in CMe30H and methyl-methyl separation distances of 4.0 and 3.73 A in CMe4 and CMe,OH, respectively. The data are also studied using the reference interaction site model (RISM).By fitting a five-centre fused hard-sphere model to the neopentane data, we obtain estimates for the Me-Me, Me-C and C-C correlation functions and the molecular coordination number. Also, by comparing RISM predictions with the results of a Monte Carlo study, we obtain some insight into the accuracy of RISM theory for fused sphere molecules. Finally, we find that it is not possible to fit the tertiary butyl alcohol data to a five-centre reference interaction site model. X-ray diffraction is a direct probe of the short-range structure of liquids. In the past, have used data obtained from such experiments to broaden our under- standing of the liquid state and, together with Monte Carlo computer simulations, test molecular interaction potential models and the reference interaction site model (RISM) developed by Chandler and In this paper we report the results of our re- cent experimental studies of neopentane (CMe,) at - 17, 25 and 150 "C, and tertiary butyl alcohol (CMe,OH) at 26 "C.The analysis of these data using Fourier trans- formation and RISM theory is presentcd, as is a test of RISM using Monte Carlo computer simulation. While this communication is the first in our present programme to study alkanes and alcohols, it is a continuation of our previous work on molecular liquids. Con- sequently, the details of the experimental and theoretical methods will not be repeated. Instead, only a summary description of the experiments and their analysis will be given. EXPERIMENT A N D DATA REDUCTION The X-ray measurements were made at Oak Ridge National Laboratory using reflection geometry6,' and monochromatic MoK, radiation (A = 0.7107 A).$ The range of scattering angles considered correspond to the interval 0.6 A-' I k i 16 A-' for CMe, and 0.4 A-' I k S 16 A-' for CMe30H, where k == (47c/A) sin 0.The counting strategy and procedures for correcting the raw data have all been given elsewhere. In fig. 1, the distinct structure functions, H,(k), determined from experiment for CMe, at - 17,25 and 150 "C are displayed, t Research sponsored by the U.S. Department of Energy and the National Science Foundation. I A = 10-*0m.40 LIQUID NEOPENTANE AND BUTYL ALCOHOL and in fig. 2 &(k) for CMe30H at 26 "c is given. The distinct structure function &(k) is defined by Hd(k) M(k)[S(k)'- Ss(k) - S,(k)] (1) H(k) - Hs(k) - Hm(k) where S(k) is the static coherent scattering function derived from experiment and normalized to the scattering from the atoms in a single molecule, Ss(k) = C f i * ( k ) , at large k (where this and the following sums are over all scattering centres in a moiecule).Herefa(k) is the X-ray scattering factor for site CI; for the work reported here, the methyl and hydroxyl groups were a I I I I I I I 'I 0.4 0.2 - - 0 - -0.2 +0.2 - 5 0 - 2 -0.2 +0.2 0 -0.2 - -0.4 - - C) - L ) - FIG. 1 .-Distinct structure function H&) for neopentane at (a) - 17, (b) 25 and (c) 150 "C. The circles are derived from experiment and the solid lines from RISM theory as described in the text.treated as single scattering (and force) centres. Thus, for example, neopentane was considered to be composed of four methyl groups and a carbon atom. The X-ray scattering factors for the methyl and hydroxyl groups were computed as described in ref. (8). The other functions in eqn (1) are defined as M(k) = [Xfal-' a and Sm(k) = molecular = C C fa(k)fs(k) exp (- I&kZ/2)j0(kRap) scattering a @ f a function (3) where Zap is the root-mean-square deviation of the local and instantaneous site-site separation distances from the mean value Rap, and jo(x) = x-l sin x. Since, in general, the function Hd(k) decays to zero much faster than U,(k), H(k) z H,(k) at large values of k. By a least squares fitting of eqn (3) to the large k portion of the experi- mental data, we have found that in CMe, the average value of over the temperature range studied was 1.55 8, and Ic--M~ was 0.05 A.For CMe,OH, Rc-M~, = 1.54 A, IC-Me = 0.05 A,A . H . NARTEN, S . I . SANDLER A N D T . A . RENSI 41 RC-oH = 1.43 8, and Ic-oH = 0.05 A. in tetrahedral locations around the central carbon atom. In both cases, the groups (methyl and hydroxyl) were 0.6 0.4 0.2 - Y - t o 2 -0.2 -0.6 -0.41 I I I 1 I I I I I 1 k/A-' FIG. 2.-Distinct structure function &(k) for tertiary butyl alcohol at 26 "C. The circles are derived from experiment, the dashed line from RISM for the van der Waals model of CMe30H, and the solid line from RISM with oMe = 3.73 A, oOH = 2.75 A and oc = 3.2 A. RESULTS DATA ANALYSIS The main focus of our work is the determination of the intermolecular, rather than intramolecular, structures of molecular liquids ; thus, the remainder of this paper will deal with the analysis of the distinct structure function, Hd(k).In fig. 1 we see that as the temperature increases (and density decreases), the peaks in Hd(k) decrease in magnitude, as expected. By comparing fig. 1 and 2, we see that the replacement of one methyl group in CMe, by a hydroxyl group to produce CMe,OH results in several significant changes in H,(k). In particular, an inner peak in Hd(k) for CMe,OH occurs near k - 0.7 k', which is characteristic of hydrogen-bonded liquids and has been reported for other alcohol^.^*'^ Also, the dominant peak in Hd(k) is of lower magnitude than in CMe, at a similar temperature. Finally, the secondary peaks in Hd(k) for CMe, and CMe,OH are of similar magnitude, but out of phase.To proceed further, note that the intermolecular site-site correlation function gaa(r) is defined so that pgaa(r) dv is the number of p sites in a volume element dv a distance Y from an cc site on another molecule in a fluid of bulk number density p. Then, Hd(k) can be written as H d k ) = M(k) Z: Z: .fa(k)fp(k>aadk) (4) a a where the partial structure function aap(k) is defined by42 LIQUID NEOPENTANE AND BUTYL ALCOHOL 1.2 1.0 0.8 h L Y CS" 0.6 0.4 0.2 Fourier inversion of the function Hd(k) derived from experiment yields a correla- tion function - - - - - - r / A FIG. 3.-Correlation function C&) for CMe4 at (a) -17, (b) 25 and (c) 150 "C. The circles are ob- tained from Fourier inversion of the experimental data, and the solid line from RISM theory.9. oo i 0 'rr' 0 0 0 0 0 0 0 0 0 0° .r#' %. 0 . 0 4= =O: 0' I I 1 I I I I I I 2 3 4 5 6 7 0 9 10 I t 12 l / A FIG. 4.-Correlation function Gd(r) for CMe30H at 26 "C obtained from Fourier inversion of the experimental data.A . H . NARTEN, S . I . SANDLER AND T . A . RENSI 43 which, from eqn (4) and (9, is seen to be a weighted super-position of site-site corre- lation functions modified by products of site scattering factors. This function is mainly of use in identifying the peak locations of the dominant atom-atom correla- tions, and thus for obtaining estimates of the size parameters for the statistical mech- anical models which follow. Fig. 3 and 4 give the correlation functions Gd(r) for neopentane at the three temperatures studied, and tertiary butyl alcohol, respectively.The most striking feature of these plots is the presence of a strong peak at 2.75 A for CMe,OH which is absent in CMe,; this peak is due to hydrogen bonded hydroxyl groups. Also, there is a shift at the methyl-methyl peak clearly evident at 4.0 A in the - 17 "C CMe, study to 3.73 A in CMe,OH. Finally, it is also of interest to note that the small k methyl-methyl peak evident in the lowest temperature CMe, correla- tion function disappears in the higher temperature (lower density) experiments, though the remainder of the Gd(r) functions are quite similar. Further information on the intermolecular structure can be obtained using interaction potential models and statistical mechanics.The procedure used was to model the CMe, molecule by four hard spheres of equal diameter fused in tetrahedral locations to a central sphere (of smaller diameter) with = 1.55 A. The site- site correlation functions for this model were then computed using RISM with the methyl group and carbon atom sphere diameters (oMe and a,-) adjusted to give the best agreement between theory and experiment for Hd(k) at -17 "C in the vicinity of the first peak. The computed values of Hd(k) for the - 17 "C data appear in fig. 1; the values of oMe = 4.0 A [in agreement with our expectations based on Gd(r)], a, = 3.2 A and oC-Me = $(aMe + a,) = 3.6 A were used in the calculations shown. We have also included, in fig. 1, the RISM predictions for Hd(k) for CMe, at 25 and 150 "C using the same molecular diameters as above (that is, in these cases we did not search for a model which gave a best fit of RISM predictions and experiment).Fig. 5 gives the RISM predictions of the methyl-methyl, methyl-carbon and carbon- carbon correlation functions at each of the three temperature (densities) studied. Also, we have obtained estimates of the numbers of nearest neighbour molecules sur- rounding each CMe, molecule using the relation N,, = 4np g,,(r)r2dr 1:'" where rmin is the value of the first minimum of the central carbon-carbon (centre of mass-centre of mass) correlation function. We find that the number of CMe, nearest neighbours is 12.3, 12.2 and 10.3 in first coordination shells which extend out to 8.25, 8.40 and 9.45 A, respectively, at - 17, 25 and 150 "C.The fact that the RISM predictions for our CMe, model fit the first peak in Hd(k) is no assurance that we have obtained good estimates for all the site-site correlation functions. Indeed, since neither the amplitude nor phase of the predicted H&) agree with experiment for k > 2 A-', this is unlikely. This disagreement at large k suggests that the site-site correlation functions predicted from RISM may not be indicative of the correlations in the real fluid at small intermolecular separations. Further evidence for this comes from comparing, in fig. 3, the correlation function Gd(r) derived from experiment with that obtained using eqn (4)-(6) and the RISM site-site correlation functions. There the first (methyl-methyl) peak in the predicted Gd(r) is both too high and too steep, which suggests that methyl group interactions are softer than our hard sphere representation.Succeeding peaks are in better agreement with the experimentally derived Gd(r) suggesting that, at larger values of r, the RTSM site-site correlation functions are approximately correct. We have, in44 LIQUID NEOPENTANE AND BUTYL ALCOHOL fig. 6, displayed the predicted Gd(r), together with a representative collection of experimental points, and the terms and fic--.(r) = (2n2pr)-' kM(k) fc(k)acc(k) dk La so defined that Gd(r) = 1 t hM+Me(r) + /?Me-C(r) + hcJr). In this way, we see the dominant influence of the methyl-methyl correlation function in determining Gd(r). I / A 3 4 5 6 7 8 9 10 II 12 I ' 1 1 I I 1 I 1 0 I I I I 1 1 I I I I I I I I 2t 0 *I- 1 I I I I I I -3 4 5 6 7 8 9 1011 12 I / A FIG.5.-Methyl-methyl, methyl-carbon and carbon-carbon correlation functions for CMe4 at -17 "C predicted from RSIM, as described in text. -- 3 g M e - M e ; ---, gMe-C; -- -3 gc-c. (6) T = 25" C, p = 0.004 884 A-3; (c) T = 150 "C, p = 0.003 039 8, - '. ( a ) T = -17 "C, /I = 0.005 277 The RlSM theory for fused hard sphere molecules is only an approximate repre- sentation of real liquids; two questions arise concerning its use. First is the question of the accuracy of the theory for a complicated fused hard sphere interaction model.A . H. NARTEN, S. I . SANDLER AND T . A . RENSI 45 I .4 r I .2 I .o 0.8 Me -C Me -Me - -0.4 -0.G t 1 , -0.8 I I , I I I I I I I r / A FIG. 6.-Correlation functions Gd(r) -, for CMe4 at - 17 "C, obtained from experiment and RISM theory, together with the functions hag which demonstrate the contribution that each of the predicted site-site correlation functions makes to Gd(r).0, experimental values; - - -, h ~ ~ - ~ ~ ( r ) ; - - -, hMe-C(r); * * * 7 hC-C(r)* 4 5 6 7 8 9 10 II r / A FIG. 7.-Comparison of the site-site correlation functions in liquid carbon disulphide obtained from Monte Carlo computer simulation and RISM.46 LIQUID NEOPENTANE A N D BUTYL ALCOHOL The second is the validity of representing any real molecule by a fused hard sphere model. Partially to answer the first of these questions, we have developedZb*" a Monte Carlo computer simulation program and are currently testing RISM theory for the fused hard sphere neopentane model.We can, however, report the results of a study of a fused hard sphere model of carbon disulphide, a liquid which was the subject of an earlier diffraction study. With the exception of a change in the inter- action model, and increasing the number of molecules to 256, the computer program is essentially the one described in our earlier work on The site-site correlations obtained from the simulation for CS2, presented in fig. 7, are in excellent agreement with the RISM predictions with the exception of the region near the first peak in the carbon-carbon (or centre-of-mass) correlation function. This suggests that RISM theory, as has been pointed out by others, is reasonably accurate for hard molecules. In the application of RISM theory to CMe,, the hard sphere diameters for the methyl groups were close to their van der Waals diameters." Consequently, our first attempt at fitting the distinct structure function H,(k) for CMe,OH used RISM theory, the bond lengths determined from our experi- ments, and the van der Waals diameter for the hydroxyl group, uOH = 3.0 A, together with uMe = 4.0 A, oC == 3.2 A and additive diameters.The values of H,(k) obtained appear in fig. 2. There are two characteristics of the calculations based on this van der Waals model which do not agree with experiment (especially, the Gd(r) correla- tion function mentioned earlier.) First, the experimental hydrogen bonded oxygen- oxygen correlation occurs at 2.75, not 3.0 A as predicted from the model. Next, what appears to be the methyl-methyl peak in Gd(r) for CMe,OH occur at 3.7 P\, rather than at 4.0 A as for both the present model and neopentane.Therefore, another RISM calculation with additive diameters was made for a model in which a, = 3.73 A, aOH = 2.75 A and a, = 3.2 A (the calculations are only slightly sensitive to the value of ac). The results of these calculations also appear in fig. 2. While this last calculation is in better agreement with the experimental data than the van der Waals model for CMe,OH, the overall agreement is still very poor and, in particular, much worse than the agreement between theory and experiment found for carbon disulphide or neopentane. Furthermore, no realistic change in the parameters would substantially improve this situation.Next, consider the CMe,OH results. DISCUSSION The work reported here is similar to our previous studies of liquid carbon disul- phide and bromine, and is part of a continuing investigation into the liquid state behaviour of real molecules. It is useful, therefore, to compare the results here with those found previously. Neopentane, like carbon disulphide, has no strong electro- static multipole moments, so that its structure is dominated by repulsive forces and RISM theory should work well. Here, however, the RISM prediction for the experi- mental distinct structure function is not as accurate as was the case for CS,. In particular, beyond the first peak, the predictions for Hd(k) are of larger magnitude and out of phase with the experimental results. The generally good agreement of the computer simulation studies of the site-site correlation functions by us and others l 3 with RISM theory suggests that the discrepancy between theory and experiment, for this molecule, is not due to the approximate nature of the theory, but rather to the inaccuracy of the five-site fused hard sphere model.Consequently, we are presently exploring the use of models in which the hydrogens in CMe, are considered explicitly, rather than lumped as part of a methyl group.A . H . NARTEN, S . I . SANDLER A N D T . A . RENSI 47 Another failing of the present calculations is that the temperature (and density) dependence of the distinct structure function for CMe, is not predicted correctly. This, presumably, is because our model for CMe4 consists of fused spheres of fixed diameter.For a softer, more realistic potential (for example, Lennard-Jones force centres), the effective hard sphere diameter would decrease with increasing temperature, which is precisely the change that is needed to bring our high temperature predictions into better agreement with experiment. Studies are presently underway in using softer interaction models in both our Monte Carlo simulation work and RISM theory.14 The most serious failure of our calculations occurs in the study of tertiary butyl alcohol. In this case, the fused sphere model which results in the best fit between experiment and RISM calculations uses hard sphere diameters which are considerably smaller than the van der Waals diameter (similar observations have been made when attempting to use RISM for strongly quadrupolar liquid brominell and polar liquid antimony trichloride.) Furthermore, here, as with bromine and antimony trichloride, the agreement between calculation and experiment is quite poor.We presently plan a Monte Carlo study of tertiary butyl alcohol in order to obtain better information about intermolecular potentials and liquid state structure. S. I. Sandler and A. H. Narten, Mol. Phys., 1976, 32, 1543. S. I. Sandler and A. H. Narten, Mol. Phys., 1978, 35, 1087. L. J. Lowden and D. Chandler, J. Chem. Phys., 1974,61,5228. C . S. Hsu, D. Chandler and L. J. Lowden, Chem. Phys., 1976,14,213. D. Chandler, Mol. Phys., 1976, 31, 1213. A. H. Narten and H. A. Levy, Water: A Comprehensive Treatise, ed. F. Franks (Plenum, London, 1974), vol. I, p. 31 1. H. A. Levy, M. D. Danford and A. H. Narten, Data Collection and Evaluation with an X-ray Diflractometer Designed for the Study of Liquid Structure, 1966, ORNL-3960. Available upon request from Laboratory Records, Oak Ridge National Laboratory, Oak Ridge, Tenn. (a) W. C. Pierce and D. P. McMillan, J. Amer. Chem. SOC., 1938,60,779; (b) G . G . Harvey, J. Chem. Phys., 1939, 7, 878. * (a) A. H. Narten, R. Agrawal and S. I. Sandler, Mol. Phys., 1978, 35, 1077; (b) R. Agrawal, * A. H. Narten, J. Chem. Phys., 1977, 67,2102. lo D. L. Wertz and R. K. Kruh, J. Chem. Phys., 1967,47, 1967. l1 R. Agrawal, M.S. Thesis (Univ. of Delaware, 1977). l2 See, for example, L. Pauling, The Nature of the Chemical Bond (Cornell University Press, l3 D. Chandler, C. S. Hsu and W. B. Streett, J. Chem. Phys., 1977,66,5231. l4 E. Johnson, personal communication. l5 E. Johnson, A. H. Narten and W. Thiessen, this Discussion p. 287. Ithaca, 1948), p. 187ff.

ISSN:0301-7249    

DOI:10.1039/DC9786600039

出版商:RSC    

年代:1978

数据来源: RSC

摘要:

Molecular Dynamics Studies of Hydrogen-bonded Liquids BY IAN R. MCDONALD Department of Chemistry, Royal Holloway College, Egham, Surrey TW20 OEX AND MICHAEL L. KLEIN Chemistry Division, National Research Council of Canada, Ottawa, Canada KIA OR6 Received 26th April, 1978 A review is given of results obtained from computer simulation of the hydrogen bonded liquids HF, H20 and NH3. An assessment is made of the relative merits of potentials obtained by fitting to quantum mechanical calculations, gas phase data and solid state properties. The pioneering work of Rahman and Stillinger1-3 has shown how the method of molecular dynamics (MD) can be used in assessing the extent to which experimental data on real molecular liquids (in their case, liquid H20) can be interpreted on the basis of relatively simple intermolecular force models.The purpose of the prcsent paper is to bring together the results of similar calculations for other hydrogen bonded systems, namely HF and NH,. Like the models adopted for HzO, the poten- tials we have used are obvious oversimplifications. In particular, many body effects, including polarization, appear only in some average sense. The intermolecular potentials we shall discuss have been derived from various sources. Information is most straightforwardly obtained from gas phase data, especi- ally from structural studies of dimers and measurement of the temperature dependence of the second virial c o e f f i ~ i e n t . ~ ~ ~ Quantum mechanical calculations on dimers can be a useful g ~ i d e , ~ ' ~ if carried out with sufficient care, and measurement of solid state properties such as the sublimation energy, crystal structure and optically active (Raman and infrared) frequencies are also of considerable value, though for hydrogen bonded solids this possibility has been largely ignored9 (see, however, the review by Schuster).'O The models we use in the present work are all of the same basic type,6*1'*'2 consisting of short range atom-atom potentials plus point charges chosen to rcproduce the known lower order electric moments of the isolated molecules. Our MD calculations have all been carried out on rather small systems, 216 mole- cules in the case of HF and H 2 0 and 108 in the case of NH3.A typical computation lasts a few thousand time steps of ~ 0 . 2 x s, so that a given MD " experiment " covers a total of about 10 ps.However, when interest is focused solely on static properties, the proton can be given an artificially increased mass and a much larger time step used. Mostly we have worked with rigid molecules, though the central force models of Stillinger,3*1' which allow for molecular vibrations, possess a number of theoretical and computational advantages, since formally the liquid can be treated as a two-component mixture. The major computational problem lies in taking proper account of the long range of the coulomb potential. The optimum choice of bound- ary conditions is a matter of much controversy13-16 which we do not have the space to enter into here. It is sufficient to say that we have employed an Ewald transforma-I . R .MCDONALD AND M . L . K L E I N 49 tion,” which though not without its own disadvantages appears to avoid the grosser errors of other methods of summation. At tbe same time it is clear, given the size of system we have used, that orientational correlations are certain to be distorted. For the present, therefore, we omit any discussion of dielectric properties. RESULTS WATER There now exists a number of intermolecular potentials which give good results for the liquid state properties of HzO, notably the ST2 model of Stillinger and Rahman,2 the latest central force model of the same workers3 and the CI model of Lie and Cle- menti.” The success of the Lie-Clementi potential is puzzling, since although it represents a numerical fit to the quantum mechanical potential energy surface for (H20)2, it yields poor results for the only dimcr property against which it has been tested, namely the second virial coefficient of steani.18 Since the ST2 and central force models also fare badly in reproducing the second virial coefficient, the question raises itself of how profitable it is to proceed in the opposite sense, that is to say by simulating the liquid with the help of potentials obtained (as in the early work of Rowlinson)19 from gas phase data.Such a potential has recently been proposed by Watts6 Based in part on fitting to the second virial coefficient, it also accounts satis- factorily for data coming from molecular beam scattering experiments. However, as can be seen from fig. 1, the results obtained for the radial distribution functions of the *t I I P I I I I I I I I I r l i 1 2 3 .4 5 6 7 8 FIG. 1 .-Radial distribution functions for H20. The solid curves are our MD results using the gas phase potential of Watts6 at 0°C and at a density of 1 g cmP3. The circles are the Monte Carlo results of Lie et a/.” for their CI model at 25°C. The dashed curves are the distribution of molecular scattering centres derived by N a ~ t e n ~ ~ from his X-ray data at 4°C.50 MOLECULAR DYNAMICS OF H-BONDED LIQUIDS liquid are disappointingly poor.*O While the first peak in goo is moderately well reproduced, the second peak is much too weak and displaced to a separation which is much too large. The other feature of note is the fact that the first peak in g H H is rather broad, with a wing extending down to a proton separation comparable with that found in the H20 molecule.Similar results were obtained in the first Monte Carlo calculations of Clemepti and coworkers12 and, indeed, there seem21 to be two classes of H 2 0 potentials: those which do, and those which do not yield a satisfac- tory second peak in goo (and hence a well-defined hydrogen bonded structure). A review of the available calculations for H20 shows that those models which give poor secondary structure in go, all have relatively weak proton-proton repulsions at the relevant interatomic separations, and this may well be the source of the difficulty. No set of pair potentials is able to reconcile the experimental data on gas and liquid, and though it is by no means obvious that such a model necessarily exists, it is also clear that the detailed form of the intermolecular potential in H20 remains an open question.AMMONIA The success of the Rahman-Stillinger central force models in the case of H20, together with the attendant computational simplifications, encouraged us to adopt a similar approach to the problem of NH3. Our first calculations22 were based on a potential due to Stillinger," designed to yield the correct monomer and dimer geo- metry, in which the non-electrostatic interactions are represented (in units of kcal mol-I and A) by the expressions 2738 400 3087.2 R12 R6 (1) V,,(R) = ~ - - (2) 9.7 1 1 R8 V N H ( R ) = - - 9.14 exp [-25(R - 1.0124)2] - 1 + exp [lO(R - 2.4)] - 7.61 exp [-25(R - 1.6)2] (3) 10 1 + exp [20(R - 2)] VHH(R) = to which must be added the coulombic interaction between charges Q = +0.267e placed on the protons and -3Q placed on the nitrogen atom.This choice of Q I I I I I I I I I 0 1 2 3 4 5 6 7 8 r l A2 PI 5 h 1 0 I . R . MCDONALD AND M . L . KLEIN n 0 1 2 3 4 5 6 7 8 r t A :; - - - I I I I ;; - - - I J 51 0 1 2 3 4 5 6 7 8 r I; FIG. 2.-Radial distribution functions for NH3. The dashed curves are MD resultsZ2 for the Stil- linger central force model" at -78°C and a density of 0.731 g ~ r n - ~ . The solid lines are our MD results at -3 "C and a density of 0.64 g for the modified Stillinger model described in the text. The dots show the distribution of molecular scattering centres derived by N a ~ t e n ~ ~ from his X-ray data at 4°C. ensures that the known gas phase dipole moment (1.47 D) is correctly reproduced. Some of our results for conditions close to the triple point are shown in fig.3. The most interesting features of the calculated distribution functions are the asymmetry of52 MOLECULAR DYNAMICS OF H-BONDED LIQUIDS the main peak in g" (which is possibly linked to the pronounced shoulder seen in the experimental results) 23 and the well-defined hydrogen bond peak in gNH. Narten's X-ray results 23 for the distribution of molecular scattering centres at 4 "C are shown as the dotted curve in the graph of g". The discrepancy between computation and experiment is considerable and due only in part to the difference in state conditions: a later calculation under the conditions of Narten's experiment yielded a peak height in g" of ~ 3 .2 , which remains much too large. Moreover, a high proportion of the molecules were found to deform irreversibly in an unphysical fashion. Attempts have been made to improve the model by varying the parameters in eqn (1) to (3), and we thank Dr. F. H. Stillinger for helpful correspondence on this subject. After much trial and error we arrived at the set of central force potentials which give the results shown by solid lines in fig. 2. The modified Stillinger (MS) potential has the form 4000000 2700 R12 R6 V"(R) = - - (4) ( 5 ) 9.71 1 1 + exp [5(R - 2.6)] VNH(R) = -jp- - 9.14 exp [-25(R - 1.0124)2] - with other terms remaining unchanged. The main effect of these modifications is to improve the agreement with experiment in the positions of the peaks in g", but at the expense of greatly weakening the hydrogen bond peak in gNH.Clearly it is difficult to achieve a liquid structure in which there is a significant degree of hydrogen bonding but in which there is not also an artificially enhanced ordering in the distribution of nitrogen atoms. In subsequent work we have largely abandoned the use of central force models and worked instead with rigid molecules. Following the programme already pur- sued for water, we have explored the properties of a number of potential models based 3 2 U u_ Lh 1 0I . R . MCDONALD AND M. L . KLEIN I 53 I 3 2 I LL -8 1 0 3 2 I I h 1 0 0 -- 2 4 6 8 r / % 12 3 a n I 4 0 2 4 r / A (c 1 6 8 '2 a 3 T I 4 1 FIG. 3.-Radial distribution functions for HF. density of 1.002 g The solid curves are our MD results for 5 "C and a using the three-centre model.27 The dashed curves give our MD results at 0 "C and a density of I .002 g cm-3 for the modified Stillinger model described in the text.54 MOLECULAR DYNAMICS OF H-BONDED LIQUIDS in part6 on fitting to the second virial coefficient of gaseous NH3.The results parallel those already described for H 2 0 to the extent that although the first peak in g" has approximately the correct height and position, the long range order is very poorly reproduced, and g H H at small separations suffers from the same defect as in Watts' model for H20. This suggests once more that the H-H potentials derived from gas phase data are not sufficiently repulsive at the relevant distances in the hydrogen bonded liquid.Solid NH3 is the simplest hydrogen bonded solid and lattice dynamics calculations are, therefore, relatively ~traightforward.~~ This offers the possibility of deriving potentials by fitting to solid-state data, and the situation here is much more en- couraging. The major computational problem is associated with the fact that in the low-temperature solid it is possible to distinguish between those N-H interactions which correspond to hydrogen bonds and those which do not. When the molecules are free to rotate, as in the liquid, it is necessary to introduce a rather artificial switching function for the N-H potential. For the models we have studied this gives rise to a double minimum in VNH, that corresponding to the hydrogen bond distance being much deeper than the other.Very recently we have found that a model of this type yields a g" in good agreement with the results of N a ~ t e n , ~ ~ even to having a shoulder on the main peak, and also gives rise to a hydrogen bonded peak in g N H which is much more pronounced than in the case of the MS potential. The computed internal energy (-4.7 kcal mol-') is also in good agreement with the experimental value. This work is continuing. HYDROGEN FLUORIDE For our initial studies on liquid H F we again used a central force model due to Stillinger," in which charges of &0.41e are assigned to the atoms, yielding a gas phase dipole moment of 1.82 D. The original model was characterised by an unrealistically large dispersion term in the F-F interaction.Changes were therefore required and the modified (MS) potential finally used (in kcal mol-' and lacking the coulombic terms) was of the form VFF(R) = 7 - exp [-4(R - 2.9)2] (6) 30 000 3 87 R VFH(R) = -+- - 18.59 exp [-25(R - 0.917)2] - 0.7 exp [-lO(R - 3)2] 5 - 1 + exp [5(R - 2)1 (7) - 21.65 exp [-45(R - 1.412)2] 20 1 + exp [20(R - 1.95)] VHdR) = +0.7 exp [-12(R - 2.7)2]. (8) The main input information used in arriving at the potentials of eqn (6) to (8) is the gas phase geometry of the monomer and the known structure of the " bent " dimer. The calculated radial distribution functions for the model liquid at 0 "C are shown as dashed curves in fig. 3. The most notable feature is the fact that gFF has a split main peak, the first component of which is due to hydrogen bonded fluorine atoms and the second arises from geometrical packing requirements consistent with the chosen density.There are, unfortunately, no experimental data with which to com- pare this interesting result. The function g F H has a strong hydrogen bond peak and the overall structure ing,, and g H H and the running coordination numbers are consistentI . R . MCDONALD AND M . L . KLEIN 55 with dimer formation. Detailed investigation of configurations generated in the MD run yields no clear evidence for the existence of rings; the structure is best described as a tangled network, with some rather large voids. The mean of the two computed self-diffusion coefficients (for fluorine atoms and protons) is 7.5 x cm2 s-l, which is only slightly larger than the experimental value of O'Reilly.25 The velocity autocorrelation functions for the two types of atom both decay rather rapidly, being indistinguishable from the noise level by about 0.2 ps.The Fourier transforms of these functions are shown in fig. 4; the proton spectrum is dominated by a peak around w/cm-' 50 0 1000 i l l l ~ l i l l ~ l l l l l l 1 ~ ' l l l ~ l l l l ~ l l 50 100 150 \ w /meV 0 70 20 30 40 ( w /2 A ) /THz FIG. 4.-Fourier cosine transform of the atomic velocity autocorrelation functions for liquid HF at 0 "C and a density of 1.002 g cm-j for the modified Stillinger model discussed in the text. 70 meV, while the fluorine atom spectrum has a pronounced shoulder at 30 meV. The simulated self motion of the protons is therefore in remarkably good agreement with experiment,26 since the incoherent spectrum of inelastically scattered neutrons is also dominated by a peak around 70 meV, a feature associated with " hydrogen bond oscillation ".Overall, therefore, the MS potential gives a good account of the few experimental facts that are available. The success of Lie-Clementi calcu1ationsl2 for H20 make it natural to follow a similar procedure in the case of HF, since reliable quantum mechanical calculations are again available. The extensive Hartree-Fock calculations of Yarkony et aL7 on rigid (HF)2 dimers were therefore fitted27 to an empirical model of the following type. Each molecule consists of three centres of interaction, with charges +0.68e assigned to the fluorine atom and proton (separated by 0.912 A) and a charge of - 1.36e placed 0.165 In this way the correct monomer dipole from the F atom.56 MOLECULAR DYNAMICS OF H - B O W E D LIQUIDS and quadrupole moments are reproduced.Short range Born-Mayer repulsions were used for the H-H and F-F interactions and a Morse function for the hydrogen bonded F-H interaction. The parameters characterising these interactions were adjusted to give a reasonable fit to the theoretical calculations, whilst a1 the same time requiring that the model give the sequence “ linear ” -w “ bent ” - “ ring ” for the most stable configuration of the dimer as the F-F distance is decreased. A correction to the F-F interaction was also added to account for the dispersion forces, the contribution from which is not included in the Hartfree-Fock calculations.The details are given el~ewhere.~’ Since the modified Stillinger and three-centre models are constrained in roughly the same way (i.e., to produce the known dimer structure), it is not surprising that they yield broadly similar results for the liquid, which contains once again a high proportion of dimers. However, as is the case in H20, the model based on quantum calculations gives rise to a liquid which is rather less structured. Some results from an MD calculation for the three-centre model are shown as the solid lines in fig. 3. We see again a splitting in the main peak of g,, and the number of atoms contributing to each component (as measuring by the running coordination number) is much the same as in the modified Stillinger model. The diffusion coefficient is now higher (8.2 x lo-’ cm’ s-’), however, and the internal energy a little lower (-6.2 kcal mol-’, compared with -5.6 kcal mol-I).CONCLUSIONS Our calculations have been far from exhaustive and much remains to be done, but a number of relatively clearcut conclusions have emerged. For H20 therc scems to be no immediate hope of reconciling the properties of gas and liquid states on the basis of a single intermolecular force model. More generally, the work on both H,O and NH, suggests that potentials derived from gas phase data are not sufficiently specific to yield good liquid structures, though they do give reasonably good ncarest- neighbour peaks in g,,(X = N, 0) and satisfactory internal energies. Clearly the problem is linked to the fact that gas phase data are rather insensitive to details of the anisotropic parts of the intermolecular potential.Solid state studies are in principle a valuable source of valuable information on the anisotropy of the potential but have so far not been exploited in a systematic way. Our most successful work to date has been that on HF, where we have relied heavily on accurate quantum mechanical calculations. It seems that NH, does not represent a simpler problem, but the ex- perience with H20 and NH, suggests that the task of constructing a satisfactory poten- tial would be greatly helped if comparable quantum mechanical results were available. In the case of HF progress is hampered by the lack of experimental data, particularly on the structure of the liquid. In summary, the best available potentials account reasonably well for the growth of hydrogen bonding in the sequence NH, + H,O --+ HF, and on balance the outlook for further work along these lines is encouraging.This work has been supported in part by the S.C.R. and by NATO. We are glad to acknowledge the hospitality of the Laboratoire de Physique ThCorique des Liquides, Universitk Pierre et Marie Curie (Paris) and of CECAM (Orsay). A. Rahman and F. H. Stillinger, J. Chem. Phys., 1971, 55, 4213. F. H. Stillinger and A. Rahman, J. Chem. Phys., 1972, 57, 1281 ; 1974, 60, 1545. A. Rahman, F. H. Stillinger and H. L. Lemberg, J. Chem. Phys., 1975,63,5223; F. H. Stillinger and A. Rahman, J. Chem. Phys., 1978, 68, 666.I . R . M C D O N A L D A N D M . L . KLEIN 57 T.R. Dyke, B. J. Howard and W. Klemperer, J. Chem. Phys., 1972, 56, 2442; T. R. Dyke K. M. Mack and J. S. Muenter, J. Chem. Phys., 1977,66,498. D. J. Evans and R. 0. Watts, Mol. Phys., 1974,28, 1233. R. 0. Watts, Chem. Phys., 1977, 26, 367; G. Duquette, T. H. Ellis, G. Scoles, R. 0. Watts and M. L. Klein, J. Chem. Phys., l978, 68, 2544. ’ G. H. F. Diercksen and W. P. Kraemer, Chem. Phys. Letters, 1970, 6, 419; D. R. Yarkony, S, V. O’Neil, H. F. Schaeffer 111, C. P. Baskin and C. F. Bender, J. Chem. Phys., 1974,60, 855. 0. Matsuoka, E. Clementi and M. Yoshimine, J. Chem. Phys., 1976,64,1351; G.H. F. Diercks, W. P. Kraemer and B. 0. Roos, Theor. Chim. Acta., 1975,36,249. F. H. Stillinger, Adu. Chem. Phys., 1975, 31, 1 ; Phil. Trans. B, 1977, 278, 97; E. Clementi, Liquid Water Structure : Lecture Notes in Chemistry, (Springer-Verlag, Berlin, 1976), vol.2. lo P. Schuster, in The Hydrogen Bond, ed. P. Schuster, E. Zundel and C. Sandorfy (North Holland, Amsterdam, 1976); A. Karpfen and P. Schuster, Chem. Phys. Letters, 1976,44,459. F. H. Stillinger, Israel J. Chem., 1975, 14, 130. l 2 G. C. Lie, E. Clementi and M. Yoshimine, J. Chem. Phys., 1976,64, 2314; G. C. Lie and E. Clementi, J. Chem. Phys., 1975, 62, 2195; H. Popkie, H. Kistenmacher and E. Clementi, J. Chem. Phys., 1973, 59, 1325. l 3 J. A. Barker and R. 0. Watts, Chem. Phys. Letters, 1969,3, 144. l4 A. J. C. Ladd, Mol. Phys., 1977, 33, 1039. l6 J. A. Barker and R. 0. Watts, Mol. Phys., 1973, 26, 789; R. 0. Watts, Mol. Phys., 1974, 28, D. Levesque, G. M. Patey and J. J. Weis, Mol. Phys., 1977, 34, 1077. 1069. S. G. Brush, H. L. Sahlin and E. Teller, J. Chem. Phys., 1966, 45, 2102. G. C. Lie and E. Clementi, J. Chem. Phys., 1976,64,5308. l9 J. S. Rowlinson, Trans. Faraday SOC., 1951, 47, 120. 2 o I. R. McDonald and M. L. Klein, J. Chem. Phys., 1978, 68,4875. 21 S. Swamimathan and D. L. Beveridge, J. Amer. Chem. SOC., 1977,99, 8392. 22 I. R. McDonald and M. L. Klein, J. Chem. Phys., 1976, 64,4790. 23 A. H. Narten and H. A. Levy, J. Chem. Phys., 1971, 55, 2268; A. H. Narten, J. Chem Phys., 24 N. Neto, R. Righini, S. H. Walmsley and S. Califano, Chem. Phys., 1978, 29, 167; R. Righini 25 D. E. O’Reilly, J. Chem. Phys., 1969, 49, 541 6; 1970, 52, 5974. 26 H. Boutin, G. J. Safford and V. Brajovic, J . Chem. Phys., 1969,39, 3135; A. Axmann, W. Biem, P. Borsch, F. HoDfeld and H. Stiller, Disc. Faraday Soc., 1969, 48, 69; J. W. Ring and P. A. Egelstaff, J. Chem. Phys., 1969,51,762; P. S. Goyal, B. A. Dasannacharya, C. L. Thaper and P. K. Iyengar, Phys. Stat. Solidi, 1972,50,701; H. Prask, H. Boutin and S. Yip, J. Chem. Phys., 1968, 48, 3367. 1968, 49, 1692; 1977, 66, 31 17. and M. L. Klein, J. Chern. Phys., 1978, 68, 5553. ’’ M. L. Klein, I. R. McDonald and S. F. O’Shea, J. Chem. Phys., 1978,69,63.

ISSN:0301-7249    

DOI:10.1039/DC9786600048

出版商:RSC    

年代:1978

数据来源: RSC

摘要:

Inclusion of Reaction Fields in Molecular Dynamics : Application to Liquid Water BY WILFRED F. VAN GUNSTEREN, HERMAN J. C. BERENDSEN AND JOHAN A. C. RULLMANN Laboratory of Physical Chemistry, University of Groningen, NijeRborgh 16, 9747 AG Groningen, the Netherlands Received 15th May, 1978 The influence of the inclusion of a reaction field on molecular dynamics simulations of liquid water (ST2 model) has been evaluated, both for a momentary and a delayed reaction field. The influence on radial distribution function and energy is not large, but the spatial dipolar correlation is changed considerably. Total dipole moment fluctuations are much enhanced by the reaction field. If the delay of the reaction field is not taken into account, the reaction field effects are overestimated. Dynamic properties are sensitive to the reaction field : the diffusion constant almost doubles and rotational correlation times decrease.Computer simulations were carried out over a total time of 55 ps with a cut-off radius of 0.578 nm and a time of 0.002 ps. The small cut-off and large time step do not influence the accuracy of statis- tical averages generated by the dynamics run. 1 . INTRODUCTION Molecular dynamics (MD) calculations on polar liquids such as those of Rahman and Stillinger on suffer from an artefact resulting from the use of a cut-off radius R, in the computation of electrostatic interactions. The interaction between polar molecules, and even more so between charged particles, contains long-range components that prescribe the use of very large values of P, or even of infinite-sum- mation methods, such as Ewald S U ~ S ~ - ~ and related or Poisson-Fourier methods.'o-12 While the use of the latter methods for simulations of systems with periodic boundary conditions is mathematically consistent, they produce simulations of a cubic crystal with a large unit cell rather than of a liquid.Thus long-range inter- actions are included for configurations that are essentially non-liquid in character. One cannot expect such long-range terms to contribute in a meaningful way to the simulation of liquid properties. In our opinion MD simulations of liquids should never include direct pair interactions over distances exceeding half the box size in order to prevent the augmentation of artefacts resulting from the use of periodic conditions.How long range the expected electrostatic effects in a liquid really are can best be appreciated if one considers an infinite volume with randomly oriented dipoles. The r.m.s. fluctuation of the field at any given point due to dipoles outside a (large) radius R, is proportional to R,-3/2, while the r.m.s. fluctuation of the potential is pro- portional to &-*. In a system of randomly placed positive and negative charges with equal average densities the convergence is even worse: the r.m.s. fluctuation of the field is proportional to R,- *, while the potential fluctuations diverge. Although the quasi-random components of the long-range interactions are quite large, systematic components ( i e . , components correlated with the configurationW . F .VAN GUNSTEREN, H . J . BERENDSEN AND J . A . C . RULLMANN 59 within the cut-off radius) could seriously influence the structure and dynamics of the simulated ensemble. This systematic component is generally referred to as the reac- tion Jield, which is related to the average dielectric permittivity of the medium outside the cut-off radius. It is the purpose of this paper to investigate the influence of the reaction field on the structure and dynamics of water, using the ST2 model of inter- a ~ t i o n . ~ The reaction field thus has been used as an approximation to the field produced by electrostatic interactions over distances larger than R,. From MD simulations1 it has become apparent that with the use of a cut-off radius strong anticorrelation of dipolar orientations occurs within the cut-off sphere.This anticorrelation severely reduces the dipole moment fluctuation of an ensemble of molecules. This effect is theoretically expected if one considers the relations between dipole moment fluctuations and dielectric constant in a sphere, either isolated in vacuum or embedded in a continuum of the same dielectric ~0nstant.l~ For a sphere of volume V in a continuum the fluctuation of the total dipole moment M is given by* These fluctuations increase roughly proportional to E . In an isolated sphere, relevant for simulations with a cut-off radius, the dipole moment fluctuations are much smaller: ( M 2 ) 9 ( E - 1) -=- eokTV ~ + 2 ’ The fluctuations are now insensitive to the value of E and cannot be used to predict macroscopic dielectric properties.The influence of the reaction field on dipole moment fluctuations has been dis- cussed by Barker and Watts14-17 and It has been shown that incorpora- tion of the reaction field in Monte Carlo (MC) calculations on liquid waterI4-l7 indeed leads to a removal of the dipolar anticorrelation. Monte Carlo calculations, however, answer no questions on dynamic properties. Hitherto they have also dis- regarded the time lag of the reaction field, corresponding to the frequency dispersion of the dielectric constant. In the following we investigate the influence of both a momentary and a delayed reaction field on MD simulations of liquid water, The effects appear to be considerable both on static and dynamic properties (such as the diffusion constant).The delay of the reaction field turns out to be essential. Since reaction field effects are quite slow, simulations over several tens of pico- seconds are necessary. We have therefore decreased the cut-off radius and increased the time step of the simulation to make the simulations feasible. A smaller R, enhances the effects to be studied. From previous evaluations we know that the cut-off radius determines the accuracy of the algorithm to a much larger extent than the time step does, allowing for a substantial increase in time step A t . In section 3 our small R,, large A t run without reaction field is compared with that of Stillinger and Rahman.j The runs with momentary and with delayed reaction field are discussed in section 4 and section 5, respectively.Section 6 contains a general discussion and conclusions. In section 2 the model and computational procedure is described. 2. MODEL A N D COMPUTATIONAL PROCEDURE As a model for liquid water we use the ST2 model of Stillinger and Rahman (SR).3 The MD technique is applied to a set of N = 216 rigid water molecules in a cubical * Rationalized S.I. units are used throughout. E is the relative dielectric constant, E,, is the permit- tivity of vacuum, 8.854 x lo-’’ F m-’.60 REACTION FIELDS I N MOLECULAR DYNAMICS box and subject to periodic boundary conditions. The mass density is fixed at 1 g ~ m - ~ , the cube edge length being 1.862 nm. The easiest way of describing the computational parameters and procedure is to specify the differences between our calculation and that of SR.3 (1) The constraints of the rigid water molecules are satisfied by applying the pro- cedure called SHAKE? They are satisfied with a relative accuracy that has been set equal to lo-’.SR applied the rigid body Newton-Euler equations; our simulations with SHAKE are entirely in Cartesian coordinates. (2) The interaction between two molecules is disregarded when the distance be- tween their centres of mass (SR: oxygen-oxygen distance) is larger than the cut-off radius R,. ( 3 ) In order to be able to describe long time scale events we use a small cut-off radius: R, = 0.578 nm (SR: 0.846 nm). This value is close to the second minimum in the radial distribution function. The influence of the reaction field will show up more clearly when applying a small cut-off radius.(4) The noise in the MD calculation is mainly due to the application of a cut-off in the interaction. Since R, is very small, it makes no sense to choose a short time step A t . Therefore we take At = 0.002 ps, which still provides a stable simulation (SR: A t = 0.000 212 61 ps). (5) The equations of motion are integrated by applying the Verlet algorithmzz instead of a Gear algorithm, since the former performs better for a large time step After 100 steps the average temperature of those 100 steps is calculated. When this value differs more than 5 K from the initial value of 289 K, all velocities are rescaled. Three different MD runs have been performed: (1) Without reaction field, denoted by the symbol NRF. Starting from an equi- librium configuration obtained from Rahman, 10 000 steps have been per- formed, covering 20 ps.Using the same starting configuration 12 500 steps have been performed, covering 25 ps. Starting from the MRF con- figuration obtained after 8000 steps (16 ps), 5000 steps (10 ps) have been per- formed. In all three cases the system needs several picoseconds to “ age ”, since reaction field effects are quite slow. The ageing is monitored by calculating over every 100 steps the mean square total dipole moment of the molecules (M2> and thedipole- dipole correlation function gp(r). The MD statistical averages reported below are calculated over intervals for which these quantities are stable. For the three runs these intervals equal 10 ps (NRF), 10 ps (MRF), 8 ps (DRF).The calculation consumes about 6.7 s CPU time per step on a CDC Cyber 74-26. Autocorrelation functions have been calculated by using fast fourier transform technique^.^^*^^ Velocity autocorrelation functions (fig. 1) have been evaluated from runs of 2500 steps (5 ps) for the NRF and MRF calculations and from a run of 1500 steps (3 ps) for the DRF calculation. In all cases all 648 velocity components have been used every 5th step (0.01 ps). Molecular dipole and quadrupole autocorrelation functions have been derived from 10 ps runs (8 ps for DRF), using averages of the molecular dipole moment components calculated over intervals of 0.2 ps. Relaxa- tion times have been estimated from a logarithmic plot of the autocorrelation func- tions (fig. 5 and 6). (6) The temperature is stabilized during a MD run as follows.(2) With momentary reaction field, denoted by MRF. (3) With delayed reactionfield, denoted by DRF.w. F . VAN GUNSTEREN, H . J . c . BERENDSEN AND J . A . c . RULLMANN 61 h e e. v 0.5 h \ !l c v > 1.c 8 v h > 0.5 0.c - I I I 0.8 t l p s (c) FIG. 1 .-Normalised velocity autocorrelation functions (u(0) . u(t)>/<u’(O)> for the centre of mass motion for the different MD runs: (a) NRF, (6) MRF, (c) DRF. 3. MD WITHOUT REACTlON FIELD (NRF) In this section we compare our NRF run with two runs of Stillinger and Rahman (SR) at nearby temperature^,^ which have also been performed without reaction field. This is done in order to trace the effects of our small R, and large A t on the static and dynamic quantities that we use in sections 4 and 5 in the evaluation of reaction field effects on MD of water.Values for In table 1 energies and pressures are compared for various simulations. TABLE ENERGY AND PRESSURE DATA FOR VARIOUS MD RUNS SR 299 -34.8 (- 35.7) -42.3 (-43.2) -0.6 (- 1 .O) NRF 299 f 6 -34.1 f 0.3 (-37.4) -41.5 5 0.2 (-44.8) -0.3 & 0.4 (-1.8) MRF 308 f 8 -32.4 0.4 (-32.8) -40.1 f 0.2 (-40.5) 1.2 f 0.6 (0.8) DRF 301 f 7 -32.8 & 0.5 (-33.2) -40.3 f 0.4 (-40.8) 0.2 & 0.4 (-0.2) Values in parentheses are corrected for long-range interactions (see text). Errors are r.m.s. fluctuations in the MD runs.62 REACTION FIELDS IN MOLECULAR DYNAMICS the SR reference have been obtained for 299 K by linear interpolation of the data for runs at 283 and 314 K, given in ref.(3). Numbers in parentheses give values corrected for the average interaction beyond A,. The correction3 consists of a Lennard-Jones term in all cases and of an electrostatic term for the SR and NRF data. Our large A t , small R, run does not yield exactly the SR values, as expected. The electrostatic correction terms overestimate the effect of the neglected interactions. The noise generated by the small R, and large A t renders a calculation of quantities that essentially depend on fluctuations in the ensemble,26 such as the specific heat, the thermal pressure coefficient and the adiabatic compressibility, meaningless. TABLE 2.--oXYGEN-OXYGEN RADIAL DISTRIBUTION FUNCTIONS SR 299 0.263 0.319 0.424 0.285 3.09 0.470 1.13 0.353 0.72 NRF 299 0.26 0.32 0.42 0.29 3.03 0.47 1.10 0.35 0.77 MRF 308 0.27 0.33 0.43 0.29 2.95 0.48 1.09 0.36 0.78 DRF 301 0.27 0.32 0.42 0.29 3.03 0.47 1.10 0.35 0.75 ~~ ~ ~~~ The R1 are the zeros, the MI are the maxima and the ml are the minima of the functiongoo(r) - 1.Table 2 shows that the oxygen-oxygen radial distribution function goo(r) is slightly smoother than those calculated with a larger R,. The r-values of the zeros, maxima and minima of goo(r) - 1 are the same in both cases. The velocity autocorrelation function for the centre of mass motion given in fig. 1 looks very much like that of SR, when interpolating the functions for nearby tempera- tures displayed in ref. (3). The self-diffusion constant D can be calculated in two ways, by D, = lim (6t)-I ( [ R ( t ) - R(0)I2) D2 = 1/3 1: ( u ( 0 ) .u ( t ) ) dt t+m and by where R(t) and u ( t ) denote the position and velocity of the molecular centre of mass. The mean square displacement of the molecules is displayed in fig. 2 as a function of time. The D, values in table 3 have been derived from these curves by a least squares TABLE 3 .-SELF-DIFFUSION CONSTANTS, DIELECTRIC QUANTITIES, AND ROTATIONAL CORRELATION TIMES MD run <T)/K D1/10-9 m2 s-I D2/10-9 m2 s-' <M2)/N,u2 a E d P S 4 P S - - - SR 299 3.1 - 0.16 NRF 299 3.1 3.3 3.9 1.5 MRF 308 6.2 5.8 23 f- 8 630 & 200 2.9 1 . 1 DRF 301 4.6 4.0 1.2 rt 0.5 34 & 15 2.9 1.1 0.19 0.07 - Errors are r.m.s. fluctuations in the MD runs. a Values in this column are referred to as " ap- parent Kirkwood factor ". From 8 ps of the equilibrium run.fit. constant. The small R, and large A t appear to have little influence on the resulting diffusionw . F . VAN GUNSTEREN, H . J . c . BERENDSEN AND J . A . c . RULLMANN 63 The static dielectric properties can be evaluated by examining dipole-dipole corre- lation for dipoles at a distance in the interval (r, r + dr): where p denotes the molecular dipole moment. This equation defines the dipole- dipole correlation function gp(r). For evaluation of eqn (3.3) the ST2 dipoles have been put at the centres of mass. Since ref. (3) does not report gP(r), no comparison can be made. ( d o ) m) = 47w2r2gP(r) dr ( 3 - 3) 7 I 5 i 5 10 -15 20 t l p s FIG. 2.-Mean square displacement <[R(t) - R(O)]’)/O.l nm2 of the molecules as a function of time (ps) for the different MD runs: (a) NRF, (6) MRF, (c) DRF.Another measure for the static dielectric structure is the mean square of the total dipole moment N i = l M = 2 pi. (3.4) Values given in table 3 show that we find roughly the same value as SR, both indicating a strong anticorrelation. Dynamic orientational correlations that were evaluated are the rnoZecuZar dipolar autocorrelation function ( d o ) w> (3 * 5)64 REACTION FIELDS IN MOLECULAR DYNAMICS (fig. 5 ) with correlation time z1 (table 3) and the correlation function of the second- order Legendre function P2 of the dipolar orientation, which we shall denote by the molecular quadrupolar autocorrelation function : This function is given in fig. 6; it has a correlation time z2 (table 3). No comparison can be made with SR for the orientational autocorrelation functions, because these have not been evaluated for the ST2 model. For the closely related BNS model' a value of 5.6 ps was found at 307 K, to be compared with our value of 3.9 ps.For z2 the BNS model gave a value of 2.1 ps, compared with our value of 1.5 ps. We do not know if these discrepancies are due to the difference in the models used or to the small R,, large At simulation. Since the orientational correlation functions are not strictly exponential, the values for the relaxation times are not very accurate. They have been derived from a least squares fit of the first 6 (for z,) or 3 (for z,) ps of the autocorrelation functions. Summarizing we may state that the MD statistical averages are not much dis- Only the fluctuations in the average quanti- torted by using a small R, and a large A t .ties have lost their physical meaning. 4. M D WITH MOMENTARY REACTION FIELD (MRF) The momentary reaction field that is experienced by a molecule i is given byl4-I6 where Mi denotes the total dipole moment of a sphere with radius R, surrounding the ith molecule. The energy of dipole and reaction field is equal to -+pi. Ei. (4.2) In the calculation the value27*28 E = 78.4 has been used. In this section we compare the MRF run with the NRF run in order to find the reaction field effects. From table 1 it can be observed that the electric correction term does not adequately simu- late the influence of the reaction field. The pressure appears rather sensitive to inclu- sion of the momentary reaction field.As was also found by Watts,l7 the radial distri- bution function g ( r ) does not change by including the reaction field (table 2). In con- trast to this the ST2 simulation by Ladd,' using infinite summation techniques for the long range interactions, produces a rather different g(r), in which the second neighbour peak is smaller and shifted outwards. The velocity autocorrelation function (fig. 1) is shifted to positive values by the momentary reaction field. This results in a self-diffusion constant D that is about twice as large as without reaction field (table 3). Although the translational tem- perature is slightly larger in the MRF run, the change in D cannot be due to a tem- perature effect alone. As was pointed out in ref.(1 7), the dipole-dipole correlation function changes con- siderably by taking the momentary reaction field into account. The quenching of the total dipole moment is relieved, as can be seen from fig. 3 and table 3. It now appears that a strong positive correlation extends over large distances, causing a build-up of the total dipole moment. The fluctuation of the total dipole moment increases (fig. 4) to values corresponding to extremely large dielectric constants (table 3). The We notice that the temperatures are slightly different.w. F . V A N GUNSTEREN, H . J . c. BERENDSEN AND J . A . c . RULLMANN 65 same behaviour was found in Monte Carlo calculations l7 of the BNS potential, where the system displayed near-ferroelectric properties. Whether the momentary reaction field overemphasises the dipolar correlation ox genuinely expresses a real property of the i ! I! 0.8 I,! 0.6 L h v - 0.4 d ! I -.-._ / . - . / ---‘‘i/ __----_ 0.0 I I \ \.- I - I -1 -- / ’ i. I _- 0.8 0.2 - - 0.2 FIG. 3.-Dipole-dipole correlation functions g d r ) for the different MD runs. - NRF, - - - DRF, _ . _ MRF. model, cannot be decided from these simulations alone. We shall see in section 5 , however, that introduction of a time lag in the reaction field brings the correlation again down to reasonable values. 1 t/ps FIG. 4.-Fluctuation of the total dipole moment M’/Np’ as a function of time (ps) for the different MD runs: (a) NRF, (b) MRF, (c) DRF. The molecular dipolar autocorrelation function (fig. 5) again shows nearly ex- ponential behaviour, but with a significantly shorter correlation time than in the simu- lations without reaction field (table 3).Qualitatively this means that the dynamic66 REACTION FIELDS I N MOLECULAR DYNAMICS counterpart of the dipolar ant icorrelation due to a cut-off is an inhibition and retarda- tion of the molecular reorientation. It is interesting that not only the rotational but also the translational motion is inhibited as judged from the lower value of the self- diffusion constant in the NRF simulation. A. -.-.-. I I I I \.-.- 2 4 6 8 10 tlps FIG. 5.-Normalised molecular dipole autocorrelation functions [eqn (3.5)J for the different NRF, - - - DRF, -. - MRF. MD runs. The molecular quadrupolar autocorrelation function (fig. 6) is changed in a similar way as the dipolar correlation function.The ratio zl/zz, which gives an indication of the type of rotational motion and reaches a value of 3 for a random rotational diffusion, equals 2.6 for both the MRF and NRF case. I 1 I I 1 2 4 6 8 - 0.21 t l p s FIG. 6.-Normalised molecular quadrupolar autocorrelation functions [eqn (3.6)] for the different MD runs. ~ NRF, - - - DRF, - * - MRF.w. F . VAN GUNSTEREN, H . J . c . BERENDSEN AND J . A . c . RULLMANN 67 Summarizing we observe that a momentary reaction field very strongly enhances the dipolar spatial correlation, increases the rate of rotational as well as translational diffusion, but does not change the radial distribution function. The potential energy increases, contrary to what would be expected of an analytical correction term.5. M D WITH DELAYED REACTION FIELD (DRF) The delayed reaction field that is experienced by a molecule i is given by (see ap- pendix) 1 & * - I . {Mi@) + e-Pz' Mi(t - t ' ) df} (5.1) Ei(t) = 4ne,(Em + +)R,3 0 where Here, Td denotes the dielectric relaxation time for which we take the experimental value of 8.20 ps.27 The static dielectric constant E, and the high-frequency value are taken equal to 78.4 and 5.3, r e s p e ~ t i v e l y . ~ ~ ~ ~ ~ The use of experimental values for the dielectric properties is an approximation ; ideally one should adjust the parameters for the continuum in an iterative procedure to the properties of the model. This is clearly an impossible task. Another inconsistency with the model is that the reaction field contains an instantaneous component related to the high-frequency polarisability of the continuum, while the model contains no such polarisability.This inconsist- ency can only be solved by using a polarisable water model. Although a prelimin- ary attempt was made to incorporate polarisability 28 we will not consider the influence of polarisability on the effects of the reaction field in this paper. It is of interest to note that the time constant of the reaction field [P-l in eqn (5.3)] is smaller than the dielectric relaxation time by a factor of 13.6, and is only 0.6 ps. This reduction of relaxation time results from the very non-linear dependence of the reaction field on the dielectric constant. It can be observed from fig. 4 that the total dipole moment fluctuation rapidly decreases when the reaction field delay is switched on.From table 1 it is observed that no appreciable changes occur as compared with the MRF run. The pressure is again reduced. The radial distribution function is not altered (table 2). The influence on the velocity correlation function and diffusion constant is to reduce the effect of the momentary reaction field, the values lying some- where between those of the MRF and those of the NRF runs. From fig. 3 it is seen that the positive dipolar pair correlations produced by the momentary reaction field are inhibited again and a correlation intermediate between NRF and MRF persists. The dipolar fluctuations are very much influenced and reduce to much smaller values. corresponding to an apparent Kirkwood factor of 1.2.The autocorrelation func- tions of the molecular dipolar and quadrupolar orientation (fig, 5 ) remain close to those of the MRF run, and have about the same correlation times. Summarizing we observe that the incorporation of a time lag in the reaction field has an appreciable influence on dipole moment correlations, but the thermodynamic properties remain similar to those of the run with momentary reaction field. The effect of the time lag is to reduce the influence of the momentary reaction field.68 REACTION FIELDS I N MOLECULAR DYNAMICS 6. DISCUSSION AND CONCLUSIONS From this work two conclusions can be drawn: (a) No obvious errors occur in the static and dynamic statistical averages in a MD run when the time step is increased. The conditions seem to be that the algorithm remains below its stability limit and fluctuations of the total energy do not exceed kinetic energy fluctuations.However, no physical quantities can be derived from fluctuations of energy and pressure. (b) The effects of a reaction field (RF) as an approximation to long-range inter- actions are considerable. The structure, as monitored by the radial distribution func- tion, is not changed, but the dipolar correlations are strongly dependent on the RF. Incorporation of the appropriate time delay in the RF is essential for dipolar corre- lations. The isothermal potential energy increases slightly, contrary to analytical corrections that can be made.3 Hence it is not appropriate to apply such apusteriuri corrections. The dynamic properties are influenced : the diffusion constant increases almost two-fold and the rotational correlation time decreases.The time lag of the RF does not have a strong influence on these dynamic properties. These results indicate that the inclusion of delayed reaction fields in the simulation of polar liquids is quite important. The momentary RF overemphasises the dipolar spatial correlations so much that it seems not to be valid to use such instantaneous fields. For Monte Carlo simulations this presents the problem that the expectation value for the delayed reaction field for a given configuration has to be evaluated. This means that in eqn (5.1) Mi(t - t') in the convolution has to be replaced by its expectation value over all configurations with given Mi(t) : (Mi(t - t ' ) ) = Mi(t) exp ( - t ' / T d ) (6.1) where we assume that Mi(t) decays with the dielectric relaxation time zd.The average retarded reaction field now becomes 1 Em-1 (Ei) = - -- This means that the RF due to is fully active and the additional RF due to E, is 93% of the corresponding value for the unretarded case. Thus the momentary RF slightly overemphasises the reaction field effects. The effect of the RF on the dipolar correlations is apparently very sensitive to the magnitude of the RF. Numerically, if we write the RF as Ei = ( ~ ~ E , R , ~ ) - ' L M ~ (6.3) J. = 0.74 if only the high-frequency part of E is taken into account; J. = 0.96 if the averaged delayed field according to eqn (6.2) is considered, and 3, = 0.98 for the mo- mentary value of the RF.Thus the average RF is very similar in the latter two cases, while the apparent Kirkwood factor is very different. It is likely that the behaviour of M(t) at small times contains terms that are not accounted for by the exponential dielectric relaxation, but which have a strong influence on the behaviour of the system. The value for the dielectric constant derived from the DRF case is quite reasonable considering its sensitivity for the RF, the inaccuracy of the value, and the neglect of polarisability in the model. The dipolar correlation can be compared to that derived from purely geometrical considerations for an ice lattice by Rahman and Stillinger30 and by Orban.19 The geometrical restrictions for the proton positions in ice lead to an overall Kirkwood factor of about 2, with a strongly positive contribution from the first neighbours,w.F. VAN GUNSTEREN, H . J . c. BERENDSEN AND J . A . c. RULLMANN 69 a negative contribution from neighbours beyond the first shell up to about 0.65 nm and a positive contribution for larger distances. The same pattern is present in the spatial dipolar correlation (fig. 3), although the strongest negative correlation occurs at shorter range than in ice. This qualitative agreement may indicate that dipolar correlation is to a large extent determined by the statistical distribution of protons over the hydrogen bonds in a largely tetrahedral network. The influence of the R F on dynamical properties like the diffusion constant cannot be traced to an obvious cause. The decrease in the orientational correlation times might indicate a shorter life-time of hydrogen bonds which would also lead to a larger diffusion constant.This point has not been further investigated. The type of rota- tional motion does not seem to depend on the RF: in all cases the ratio of dipolar to quadrupolar relaxation time is 2.6. In terms of an (inadequate) simplified Markovian rotational jump model3' this value indicates a jump angle of <45". Whether a reaction field is included or not does not appreciably influence the radial distribution function goo(r). When long-range interactions in a periodic structure are taken into a c c o ~ n t , ~ the goo(r) changes to the extent that the second neighbour peak is displaced, indicating rather drastic structural changes.We, therefore, would like to stress our remark made in section 1 that the use of infinite summation methods accentuates the artefacts of the periodic boundary conditions. Hence such methods should be avoided if one wishes to simulate the properties of a real polar or ionic liquid rather than those of a mathematical construct. APPENDIX DERIVATION OF DELAYED REACTION FIELD The reaction field E experienced by a molecule in the centre of a sphere with radius R, due to a continuum outside the sphere with dielectric constant E is given by l4-I7 E = FM, (A. 1) where M is the total dipole moment within the sphere and the response F is given by For time-dependent moments M(t), a time-dependent response E(t) results. Since this response is linear, it is determined by the convolution of M with a response func- tion f( z) : E(t) == lorn dzf(z) M(t - z) + g M(t).(A.3) In this notationf(z) contains no &function; the instantaneous reponse g is deter- mined by the high-frequency component of F. We assume that the frequency de- pendence of E is adequately described by &(a) = Ew + ( E , - &W)/(l + iWZd). (A.4) In the experimental relaxation function 28 the term i o T d occurs with an exponent 0.986; this exponent we approximate by 1.70 REACTION FIELDS I N MOLECULAR DYNAMICS From linear response theory it follows that Eqn (A.6) can be rewritten to wherep = ( E , - E ~ ) / ( E ~ - 1) (A.8) and q = (E, - E ~ ) / ( E , + 4) (A.9) which, after evaluation of the integral, yields f(z> = F(co)(P - q)zd-l exp [-(l + q>z/zdl* (A.10) Combination of eqn (A.3), (AS) and (A.lO) immediately results in eqn (5.1). A. Rahman and F. H. Stillinger, J. Chem. Phys., 1971, 55, 3336. ’ F. H. Stillinger and A. Rahman, J. Chem. Phys., 1972, 57, 1281. F. H. Stillinger and A. Rahman, J. Chem. Phys., 1974, 60, 1545. P. P. Ewald, Ann. Phys., 1921, 21, 1087. H. M. Evjen, Phys. Rev., 1932, 39, 675. F. Lantelme, P. Turq, B. Quentrec and J. W. E. Lewis, Mol. Phys., 1974, 28, 1537. V. M. Jansoone, Chem. Phys., 1974,3, 78. E. R. Smith and J. W. Perram, Mol. Phys., 1975, 30, 31. A. J. C. Ladd, Mol. Phys., 1977, 33, 1039. lo R. W. Hockney, Methods Comput. Phys., 1970, 9, 176. l1 R. W. Hockney, S. P. Goel and J. W. Eastwood, Chem. Phys. Letters, 1973, 21, 589. l2 R. W. Hockney, S. P. Goel and J. W. Eastwood, J. Comput. Phys., 1974, 14, 148. l3 H. Frohlich, Theory of Dielectrics (Oxford Univ. Press, London, 2nd edn, 1958). l4 J. A. Barker and R. 0. Watts, Faraday Symp. Chem. Soc., 1972, 6, 161. l5 J. A. Barker and R. 0. Watts, Mol. Phys., 1973, 26, 789. l6 J. A. Barker, Report of CECAM workshop on Molecular Dynamics and Monte Carlo Calcula- tions on Water, ed. H. J. C. Berendsen (Centre Europken de Calcul Atomique et MolCculair, Orsay, 1972), chap. 2, p. 21. l7 R. 0. Watts, Mol. Phys., 1974, 28, 1069. ’* A. Bellemans, Report of CECAM workshop, 1972 [see ref. (16)], chap. 2, p. 6. l9 J. Orban, Report of CECAM workshop, 1972 [see ref. (16)], chap. 2, p. 1 1 . ‘O H. J. C. Berendsen. Report of CECAM workshop, 1972 [see ref. (16)], chap. 2, p. 29. J. P. Ryckaert, G. Ciccotti and H. J. C. Berendsen, J. Comput. Phys., 1977, 23. 327. ’’ L. Verlet, Phys. Rev., 1967, 159, 98. 23 W. F. van Gunsteren and H. J. C. Berendsen, Mol. Phys., 1977,34, 131 1. 24 R. P. Futrelle and D. J. McGinty, Chem. Phys. Letters, 1971, 12, 285. ’’ E. Kerstemont and J. van Craen, J. Comput. Phys., 1976, 22, 451. 26 P. S. Y . Cheung, Mol. Phys., 1977, 33, 519. 27 R. Pottel and U. Kaatze, Ber. Bunsenges. phys. Chem., 1969, 73, 437. 29 H. J. C. Berendsen and G. A. van der Velde, Report of CECAM workshop, 1972 [see ref. (16)], 30 A. Rahman and F. H. Stillinger, J. Chem. Phys., 1972, 57, 4009. 31 J. E. Anderson, Faraday Symp. Chem. Soc., 1972, 6, 82. C. G. Malmberg and A. A. Maryott, J. Res. Natl. Bur. Stand., 1956, 56, 1 . chap. 6, p. 63.

ISSN:0301-7249    

DOI:10.1039/DC9786600058

出版商:RSC    

年代:1978

数据来源: RSC

摘要:

GENERAL DISCUSSION Prof. D. Chandler (IZZinois) said: After reading the first two papers, I feel that perhaps more is being asked of the RISM theory than ought to be expected. It is really a very primitive theory. I developed it to predict what sorts of intermolecular structural effects could be produced by the intramolecular constraints associated with chemical bonding. As a rough estimate for the mechanism of intermolecular site- site correlations, I devised an Ornstein-Zernike-like integral equation which included intramolecular pair distributions in an explicit way. The integral equation pictures atomic pair correlations as being propagated uia simple chains composed of " direct " intermolecular pair correlations and intramolecular pair correlations. If this picture is physically meaningful, the direct correlations should play the role of an effective potential and be short ranged.This idea leads to the RISM closure relation that looks similar to the PY closure for hard sphere molecules. With these comments I want to draw attention to the following points: (1) The comparison of the RISM theory results with those of computer simulations on hard core polyatomics shows that the major qualitative features of the inter- molecular site-site distributions are all correctly predicted by the RISM equation. Since many of these features are very different from those observed in monatomic systems, the agreement between RISM and the computer simulations serves to estab- lish the essential correctness of the physical picture surrounding the theory. (2) To understand the details that are missed by the RISM equation (e.g., the shoulder to the left of the main peak observed in the centre-centre distribution for long diatomics or long linear triatomics), one should consider what physical elements are missing from the basic framework of the equation.I do not think one should artificially graft long ranged tails onto the site-site direct correlation function. I defined this function in the context of a theory which utilizes only pair intramolecular correlations. Corrections to the RISM theory must include higher order intra- molecular distributions. Indeed, interaction site cluster graphs involving intramolecu- lar distributions with three sites give rise to nonanalytic behaviour in the regions where Streett and Tildesley have observed the (fortunately small) deficiencies in the RISM equation. (3) The physical picture described above emphasizes that the RISM equation is not applicable for systems with strong and specific attractive interactions. Thus, I think the presence of hydrogen bonds is the source of the disagreement reported by Sandler between the RISM theory and the experimental results on tertiary butyl alcohol.The disagreement is almost certainly not due to dipole-dipole (or other long ranged electrostatic) interactions since the RISM theory has no difficulty in explaining X-ray and neutron diffraction experiments on acetonitrile.' Dr. M. Rigby and Mr. P. A. Monson (London) (communicated): The analysis of Streett and Tildesley assumes that the resummation of the spherical harmonic expan- sion can accurately represent the pair correlation function for a fixed relative orienta- l C.S. Hsu and D Chandler, Mol. Phys., 1978, 36, 215.72 GENERAL DISCUSSION tion at separations close to discontinuities in the pair potential. In the low density limit, the pair potential and the pair correlation function are related by the equation For hard particle systems, in which the pair energy is either zero or infinite, g(r12, el, 02, $9,,) changes discontinuously at the value of r12 which corresponds to con- tact. This must be true even for orientations for which gooO(rl2) is zero at the contact separation. Dr. D. J. Tildesley and Prof. W. B. Streett, (Zthaca) said : The spherical harmonic expansion is recombined using g(r12, 81, 8 2 , $912) = 471 2 glll2m(r12) Yl,m(m1> Y1,-m(w2>* (1) Lhm The harmonic coefficients are calculated in the sirnulation from Substituting eqn (2) in (1) g(r129 81,829 $912) = 16n2gooo(r,J 2 ( Y;tn(wl) Y ; - m ( d > Ydm,) IJpm Yl~--m(w2) shell r12 to r t a + dr1n (3) Now at closest approach all the indications seem to be that gooO(rl2) is zero.1*2 There- fore for hard diatomics in the cross orientation (6, = 8, = q12 = 4 2 ) and for hard spherocylinders in the crossed and parallel orientations (8, = O2 = 4 2 , v),, = 0), g(rI2, O,, 02, q12) must be zero at contact because of the factor gooO(rl2) = 0 in eqn (3).For these fixed orientations at least, the first peak of the recombined distribution is removed from zero and we speculate that for many other orientations it is as well. Prof.W. R. Smith (Nova Scotia) said: Can the spherical harmonic coefficients be recombined to obtain sections through the full distribution function g(r12, el, 02, p12) with any accuracy? Dr. D. J. Tildesley (Ithaca) said: If the spherical harmonic coefficients are to be recombined to calculate the distribution of a particular orientation, then the accuracy of that recombination will depend on just how many times that orientation, or orienta- tions close to it, were sampled in the simulation where the coefficients were calculated. For example, if the T-orientation is rare at a particular value of rI2, then the recom- bination of the coefficients for the T-orientation at that r12 may give inaccurate or meaningless (negative) value of g(r12, B,, 02, q12).For soft potential models, coeffi- cients based on 25 000-50 000 samplings of the pair distribution are estimated to have an average uncertainty of 0.03 and a maximum of 0.1. The total uncertainty in the recombination then depends on the number of coefficients involved. g(r12, el, 02, p12) is a function which contains an embarrassing wealth of information and although we realise that the most stringent test of a theory is to compare theoretically calculated and computed sections through g(r12, O,, 02, p12), we feel that theories should be cast to calculate particular harmonic coefficients. These are adequate for many purposes and can be obtained accurately from the simulation for the purposes of comparison. Mr. P. A. Monson and Dr.M. Rigby (London) (communicated): Smith has com- mented on the convergence of the spherical harmonic expansion of the pair distribution W. B. Streett and D. J. Tildesley, Proc. Roy. Soc. A , 1976,348, 485. P. A. Monson and M. Rigby, Chem. Phys. Letters, 1978,58, 122.GENERAL DISCUSSION 73 function, especially as revealed by the behaviour of fixed orientation distribution functions close to contact. We have observed behaviour similar to that reported by Streett and Tildesley' in our simulations of hard spherocylinder fluids., For fixed relative orientations, e.g., parallel, crossed, T-shaped, or end-to-end configurations, the distribution functions have maxima at separations which are a little greater than the contact value. The detailed behaviour seen is quite sensitive to the number of terms used in the spherical harmonics summation, and we are now carrying out further simu- lations in order to try to establish directly the experimental behaviour of these func- tions. By analogy with hard sphere results these would be expected to rise smoothly and monotonically to contact, and it seems probable that the qualitatively different results observed by re-summing the spherical harmonic expansion reflect the poor convergence of the expansion in this range.Dr. P. A. Madden (Cambridge) said : The orientational correlation parameter g , is calculated from eqn (5.3) of the paper. I would like to know what the range of the integrand was found to be. My second question could equally well be directed to Chandler, for Hsu, Chandler and Lowden3 proposed that a reasonable and readily calculable representation of the full pair correlation function g could be found from the site-site correlation functions gafi by setting Have Tildesley and Streett, or anyone else, tested this approximation against a direct simulation of g? Dr.D. J. Tildesley (Zthaca) said: Fig. 1 shows the average value of P,(cosy) per molecule as a function of intermolecular separation for a model fluid corresponding to ::::I - 0.5 FIG. 1.-Average value of P2(cosly) per molecule in a particular shell is plotted as a function of r12 (the intermolecular separation). The fluid, model I, corresponds to liquid CS2 at 298 K and 1 atm pressure. W. B. Streett and D. J. Tildesley, Proc. Roy. Soc. A , 1976, 348, 485. P. A. Monson and M.Rigby, Chem. Phys. Letters, 1978, 58, 122. C. S. Hsu, D. Chandler and L. Lowden, Chem. Phys., 1976,14,213.74 GENERAL DISCUSSION liquid CS2. This function, which is a measure of the angular ordering, becomes zero at 254.50,. The range and accuracy of G2 can be determined by truncating the upper integration limit in eqn (5.3) of our paper at various distances and observing the effect on G2. 2.500 0.2007 3.125 0.2002 3.750 0.1641 4.375 0.261 1 5.000 0.3325 5.200 0.3349 Answering Madden’s second query : the site superposition approximation (SSA) of Hsu, Chandler and Lowden enables us to estimate the total distribution function in the liquid [g(r120102~12)] from the site distribution functions gap(raB). That is m m where the products are over all sites in a pair.This approximation has been used to calculate the integrand in eqn (5.3) of our paper using only the site distribution func- tions for the model.’ In calculating this integrand, which can also be obtained exactly from the simulation, the SSA preserves the positions of the peaks and all the significant structure in the function, but produces slightly incorrect peak heights at close approach. For this model SSA enables us to calculate one set of partial structure functions, the gllZsm coefficients, from another set, the gbcc(rap) site distribution functions, with reasonable quantitative accuracy. Prof. D. Chandler (Illinois) said: Hsu, Lowden and I introduced the function G( 1, 2) = C(R1,) gay( Irl(a) - Y 2 (’) I) a,y to provide a qualitative measure of the full two-molecule distribution function., I suspect that there are probably some systems for which g o , 2) = G(L 2) is not nearly as good an approximation as Ladanyi et al.found it to be. correlations contributing to g,. Let The other point I want to make concerns Madden’s question about the range of g , = 1 + ( ( N - l)Pl(U’ U,)), 1 = 1,2, 3 , . . ., where the pointed brackets indicate the equilibrium ensemble average, Pl(x) is the Zth Legendre polynomial of x, and Uj is a unit vector associated with the jth molecuie in an N molecule system. One may relate g , - 1 to the 21th moment of a simple linear combination of gay(r)’s. This fact is quickly verified by examining the small B. M. Ladanyi, T. Keyes, W. B. Streett and D. J. Tildesley, MoZ. Phys., 1979, in press.C. S. Hsu, D. Chandler and L. J. Lowden, Chem. Phys., 1976,14,213.GENERAL DISCUSSION 75 wave vector Taylor series of the fourier transforms of the site-site distributions. Thus, g , and g , are long wavelength properties that are probably not easily interpreted in terms of local structures. To emphasize this point, let me describe a fact that is not well known but may be of interest to many. If one considers a symmetric linear molecule such as CS,, the exact value of g , must differ by some function of the density from the uncorrelated value of unity. However, for this class of molecules, the evaluation of the fourth moment of the appropriate linear combination of gcLy(r) values in the RISM equation approxima- tion gives a g , value of exactly unity. This fact does not mean that the RISM equa- tion is a poor theory for systems composed of symmetric linear molecules.Indeed, Sandler and Tildesley have pointed out that the theory is quite good at describing the local structure of liquid CS,. I think the appropriate conclusion is that local struc- tures of molecular liquids are not usually related in any qualitative way to long wave- length properties. Prof. J. M. Haile (Clemson, S. Carolina) said: In reply to Madden's question con- cerning the long-range nature of the angular correlation parameter G,, we have performed molecular dynamics simulations of a molecular fluid modelled by a Len- nard-Jones atom plus quadrupole-quadrupole interaction. During the simulation the coefficients glllBm(r) for a spherical harmonic expansion of the angular pair correla- tion function g(ro,c~),) were calculated via the method of Streett and TildesleyS2 These coefficients can be used to calculate the static properties of the fluid; thus, G2 is given by: G2 = (C P2(c0s y i j ) ) (1) where y i j is the angle between the symmetry axes of linear molecules i andj, p* = pa3 is the reduced number density, and (..) indicates an equilibrium ensemble average. Fig. 2 shows the integrand of eqn (2); i.e., 2 m = - 2 C g22m(r)r2 plotted against molecular pair separation r for quadrupolar fluids having reduced quadrupolar moments Q* = Q ( E ~ ) - + = 0.5 and 1.0. Also shown, for comparison, is the integrand for the quadrupolar contribution to the internal energy UQQ of the Q* = 1 fluid: UQQINE = 2np*Q* - dm / g,,,(r)C(224; m ~ O ) r - ~ dr (3) 5 m E - 2 where C is a Clebsch-Gordon coefficient in the convention of Rose3 and m_ = -m.Fig. 2 indicates that G2 is systematically nonzero at significantly larger molecular pair separations than most thermodynamic properties, at least for this particular po- tential model. This long range character of G2 and the fact that the integrand is com- posed of both positive and negative parts of approximately equal magnitudes are likely parts of the explanation for the difficulty in obtaining reliable values for G2 from simulation. J. M. Haile, K. E. Gubbins, C. G. Gray and W. B. Streett, Mol. Phys., 1979, in press. W. B. Streett and D. J. Tildesley, Proc. Roy. SOC. A, 1976, 348, 485. M. E. Rose, Elementary Theory of Angular Momentrim (Wiley, New York, 1957).76 GENERAL DISCUSSION -21 I I I 1 I 1 2 3 rlcT FIG.2.-Indication of the long-range nature of the angular correlation parameter G1 from a molecular dynamics simulation of a molecular fluid modelled by a Lennard-Jones atom plus quadrupole- quadrupole potential. Shown is the radial dependence of the integrand for G2 for fluids having quadrupolar moments Q* = 0.5 (a) (pa" = 0.85, kT/& = 1.277) and Q* = 1.0 (b) (po3 = 0.85, kT/& = 1.294) compared with that for the quadrupolar contribution to the internal energy for the fluid having Q* = 1 (c). In response to questions on estimation of the angular pair correlation function g(ro,m2) by means of a superposition of the site-site functions gab(rab), we have tested the approximation using molecular dynamics results.We use the superposition ap- proximation in the form : where gc is the centre's pair correlation function and the angle brackets (. -> indicate an unweighted average over orientations : The right-hand side (r.h.s.) of eqn (4) differs by a constant from that originally sug- gested by Hsu, et al.' With this definition the angle average of the r.h.s. of eqn (4) is g&) which is the result obtained on angle averaging the left-hand side (1.h.s.). C. S. Hsu, D. Chandler and L. J. Lowden, Chern. Phys., 1976, 14, 312.GENERAL DISCUSSION 77 Eqn (4) has been tested using molecular dynamics results for a Lennard-Jones atom plus quadrupole-quadrupole model fluid. During the simulation the site-site and site-centre distribution functions g,a(raa) and the coefficients gzlzo,(r) in a spherical harmonic expansion of g(rc0,wJ were determined.For the purpose of solving the rotational equations of motion, the molecules were considered to be homonuclear diatomics with atom-atom bond length = 0.3292 0. The site-site distribution functions were obtained by sampling distances between " sites " located at the ends of this bond axis. The 1.h.s. of eqn (4) was calculated by recombining the coefficients gZllzm ( r ) for particular pair orientations. The r.h.s. of eqn (4) was calculated directly 0 I .o 1.4 18 2.2 T I C FIG. 3.-Test of the method of estimating g(rwloz) by a superposition of the site-site correlation functions gaB(rtLB) for a Lennard-Jones atom plus quadrupole-quadrupole fluid. Shown areg(rolw2) values for molecular pairs in the " cross " (+ ; O1 = Oz = v, = 90") and " tee " (T; 8, = go", 19~ = 0", v, undefined) orientations. The values were determined by molecular dynamics (MD) and by the superposition approximation (SA) using MD determined gaB(raa).Parameters for the MD calculation were: p 4 = 0.85, kT/e = 1.270, Q* = 1, molecular bond length Z/a = 0.3292, number of molecules N = 256. from the MD determined gab values. Comparison of the 1.h.s. and r.h.s. for a number of pair orientations showed eqn (4) to be in serious error for this model potential. Samples of the results are shown in fig. 3 for molecular pairs in the " tee " and " cross " orientations. Thus, although an approximation of the type given by eqn (4) appears to work for site-site potential models, it does not seem to generalize to other fluids.The last paragraph in section 3 of the paper by Streett and Tildesley is worth emphasizing; namely, that over the last few years a consensus has evolved which states that the short-range, repulsive part of the intermolecular potential largely deter- mines the structure in atomic fluids (argon, for example) but that this is probably an oversimplification if applied to molecular fluids. In studying molecular fluids we must consider not only how the potential varies with intermolecular separation (forces), but also how it changes with relative orientation of the molecules (torques). If it is worth trying to develop similar generalizations for molecular fluids in terms of both forces and torques, perhaps we should seek limiting cases wherein either the forces or the torques are dominant; i.e., study particular molecules in certain regions of the phase diagram where only one of these effects is important.For example, some mole- cules (nitrogen, methane) exhibit a plastic crystalline phase wherein the molecules78 GENERAL DISCUSSION undergo almost wholly rotational motion as opposed to translation. It seems to me that these ideas are worthy of further thought and discussion. Prof. W. A. Steele (Pennsylvania) said: It is interesting to note that the gzzm(r) seems to change from a short-range function for hard-core fluids, as found in several of the simulations of Streett and Tildesley,' (in addition to the work presented here) to a relatively long range function when the model potential includes quad- rupolar terms (see the comment by J.M. Haile in this discussion). This behaviour is quite different from that of the centres correlation [gooo(r)], for dense fluids, which has been shown to be rather insensitive to the long-range forces (angle-dependent or otherwise) by a number of workers. It suggests that theoretical studies of the effects of long-range angle-dependent interaction energies upon molecular correlation func- tions might be well advised to concentrate on those parts of the correlation function tion that have the same angle-dependence as the potential terms. Existing studies of this kind can be found in recent theories of the dielectric constant of polar fluids. Prof. D. Chandler (Illinois) said: As Egelstaff's paper suggests, I think we should focus attention on the site-site distribution functions, gay(r), since it is these correlation functions that we measure experimentally.Prof. J. S. Rowlinson (Oxford) said: The discussion on the advantages of repre- senting molecular structure by site-site distribution functions, gcrB(r), or by the multi- variable function, g(r, mi, w j ) or its expansion into spherical harmonic terms glllSm(r), emphasised both the direct relation ofgab(r) to what is observed in an X-ray or neutron diffraction experiment, and also the convenience of the other functions for pictorial representation. There are, however, experiments which lead directly to certain terms in the spherical harmonic expansions. Thus dielectric and infrared absorption select the component with 2 = 1, and light scattering that with 2 = 2.The advantage of a direct link with experiment is, therefore, not entirely on the side of the site-site functions. Mr. P. A. Monson and Dr. M. Rigby (London) (communicated): We have studied the equation of state and structural properties of fluids of hard prolate spherocylinders of length: breadth (LIB) ratios 2: 1 and 3: 1 using the Monte Carlo method. Equi- molar mixtures of hard spheres and spherocylinders (LIB = 2) have also been studied for the case in which both components have the same breadth (diameter, a). For pure spherocylinders the equation of state was found to agree well with the results of previous worker^.^^^ Results for the mixture are shown in table 1, together with TABLE 1 .-P V/ Nk T FOR THE MIXTURE OF SPHERES AND SPHEROCYLINDERS d Monte Carlo Boublik eqn 0.30 4.10 f 0.04 4.24 0.40 7.31 f 0.06 7.54 0.45 9.87 f 0.10 10.30 Reduced density, d, is the packing fraction.D. J. Tildesley and W. B. Streett, Proc. Roy. SOC. A, 1976, 348, 485. D. W. Rebertus and K. M. Sando, J. Chem. Phys., 1977,67,2585. J . Vieillard-Baron, Mol. Phys., 1974, 28, 809.GENERAL DISCUSSION 79 values calculated from the equation of Boublik.' The agreement is generally good, and the discrepancies may be largely attributed to the contribution from spherocylinder- spherocylinder interactions. We have attempted to use our mixture results to evalu- ate the volume of mixing of spheres and spherocylinders at constant pressure. The excess volume is very close to zero, within the error of the data.Structural information about pure spherocylinder fluids has been obtained by calculating coefficients in the spherical harmonics expansion of the orientation depend- ent pair distribution function.2 Fig. 4 shows gooo(r) for pure spherocylinders at 0.0 I I r l u I I 1 .o 2.0 3.0 FIG. 4.-gOo0(r) for pure spherocylinders as a function of r. (- -) L/B = 2, d* = 0.5; (. - -) LIB = 2, d* = 0.4; (-)LIB = 2, d* = 0.3; (---)LIB = 3, d* = 0.3. several densities and for both shapes. At low density a broad, low first peak is seen, which develops a shoulder at higher densities. For the more elongated molecules, LIB = 3, the main peak is more spread out. The coefficients in the expansion of the angle dependent distribution function were evaluated and were found to resemble qualitatively those reported for hard diatomic molecules.2 For the mixtures, three pair distribution functions were used, representing the three types of pair interactions present.The results showed that the coefficients in the expansion of the spherocylinder-spherocylinder distribution function were almost indistinguishable from those for pure spherocylinders at the same packing fraction. The addition of spheres does not appear to significantly modify the spherocylinder- spherocylinder correlations. The sphere-spherocylinder distribution functions show significant orientation dependence at small separations, but these are beyond 2.5 0. Dr. J. C. Dore and Mr. I. P. Gibson (Canterbury) (communicated): The results shown in fig.4 of Streett's paper for neutron diffraction from liquid CS, were obtained at Harwell (A = 0.84 A) and we have now completed additional measurements (A = 0.7 A) of higher statistical accuracy at the Institut Laue-Langevin. Fig. 5 shows the experimentally measured structure factor with the fitted molecular form-factor. The resulting curves for the intermolecular funclion are shown in fig. 6; the upper curve corresponded to the data presented in the paper by Streett and Tildesley. There is good agreement between both data sets and this should permit detailed comparison with simulations based on modified potential functions. T. Boublik, J . Chem. Phys., 1975, 63, 4084. ' W. B. Streett and D. J. Tildesley, Proc. Roy. SOC. A , 1976, 348, 485.80 GENERAL DISCUSSION 1.00 0.80 0.60 0.40 h 0.20 Y \- 0 Q Y E -2 0.80 0.60 0.40 0.20 0.00 I I I I I I I I I 1 I I I I I I I I 1 I I 2 4 6 8 10 12 14 16 Q /A’’ FIG.5.-Molecular structure factor, S,,,(Q), and the fitted molecular form-factor fl(Q), for carbon disulphide using an incident neutron wavelength of (a) 0.84 and (b) 0.7 A. Dr. J. L. Finney (London) said: McDonald and Klein rightly stress the need to test proposed water potentials against known gas and solid phase properties, They then draw attention to the present inability of any known set of pair potentials to reconcile such experimental data on the gaseous and condensed phases. This failure is to be expected. We have several pieces of experimental and theoretical evidence which lead us to conclude that any pair-additive potential will inherently be unable to account for both gas and condensed phase properties.For example: (1) water molecules in the gas phase are poor electron donors and acceptors,l yet this acid/base property is one of the major attributes of the liquid phase. The isolated water molecule is a weaker Bronsted base than formic acid; thus any poten- tial which attempts to account for those solvent properties of water which are related to this donor/acceptor capability must be able to account for the almost qualitatively different behaviour in the condensed liquid as against the dilute gas; (2) quantum mechanical calculations on single molecules, dimers and small clusters demonstrate the relatively large magnitude of non-pair additive effects in the water-water interac- t i ~ n .~ . ~ For example, Del Bene and Pople2 calculate energies of -6.09, - 14.66 and E. M. Arnett, B. Chawla and N. J. Hornung, J . Soln. Chem., 1977,6, 781. D. Hankins, J. W. Moskowitz and F. H. Stillinger, J. Chem. Phys., 1970,53,4544. * J. Del Bene and J. A. Pople, J. Chem. Phys., 1970,52,4858.GENERAL DISCUSSION 81 0.3 0.2 0.1 0.0 0.4 0.3 0 - 0.2 E 0.0 -0.1 -0.2 -0.3 -0.4 FIG. 6.--lntermolecular contributions, Dm(Q), for liquid carbon disulphide using an incident neutron wavelength of (a) 0.84 and (b) 0.7 A. -23.45 kcal mo1-' for the water dimer, sequential trimer, and open tetramer respectively. Considering the numbers of hydrogen bonds in the three clusters (1, 2 and 3, respectively), it is clear that the potential is highly non-pair additive, and that no effective pair potential could reproduce these energies even approximately; (3) the dipole moment of the isolated water molecule in hexagonal ice I was estimated by Coulson and Eisenberg' to be around 2.6 D, while a recent measurement gives a value between 2.6 and 3.0 D.2 Such a large (40-50%) enhancement of average dipole moment is highly relevant to the dielectric properties of water, and hence to an under- standing of its solvent properties.Yet such enhancements cannot be dealt with within the limitations of a pair-additive formalism. Once the significance of non-pair-additivity is recognised, there are several possible ways of including it in the design of a suitable potential. We have d e ~ e l o p e d ~ - ~ one such approach, which has proved successful in accounting for a variety of both gas and condensed phase properties.In its simplest form, the water molecule is repre- C. A. Coulson and D. Eisenberg, Proc. Roy. Soc. A , 1966, 291,445. E. Whalley, personal communication. P. Barnes, in Report of CECAM Workshop on Molecular Dynamics and Monte Carlo Calcula- tions on Water, ed. H. J. C. Berendsen (Centre Europeen de Calcul Atomique et Moleculaire, Orsay, 1972), p. 77. P. Barnes, in Progress in Liquid Physics, ed. C. A. Croxton (John Wiley, Chichester, 1978). P. Barnes, J. L. Finney, J. D. Nicholas, J. Quinn and A. V. Westerman, Acta Cryst., 1978, A34, S200.82 GENERAL DISCUSSION sented by a spherically-symmetrical repulsive core (presently of LJ 6-1 2 type) together with an electrical multipole expansion, of which the dipole and quadrupole terms are dominant.The non-pair-additivity is handled by allowing the electron distribution so described to be polarisable under the influence of the electric fields of its neighbours. The electrical moments and polarisability are chosen from best known experimental or theoretical values and are not adjustable. The only adjustable parameter allowed is the o of the Lennard-Jones contribution which is fixed with reference to the experi- mental dimer equilibrium separation.’ Although this potential is purely classical, there is good reason from quantum mechanical calculations to expect this to be a remarkably good approximation. A review of such calculations suggests the apparent degree of electron exchange in the water-water hydrogen bond falls as the size of the basis set used in the calculation is increased.Notwithstanding this sug- gestion, the potential could be varied to take account of specific exchange interactions if required. As Hildebrand has stressed in his paper presented at this meeting, if a model dis- agrees with one experimental fact, it is wrong. In this spirit, we have been concerned to test this model potential as fully as possible against known gas and solid phase properties before turning the (expensive) handle on simulation calculations on the liquid phase. So far, we have been able to account for the equilibrium geometry and energy of the water dimer, as well as the bending and stretching force constants. The predicted second virial coefficient is good.For small (3-6 water molecule) clusters, the trends in non-pair-additivity energy contributions predicted by quantum mechani- cal calculation^^*^ are reproduced satisfactorily by this potential. The average dipole moment in hexagonal ice I is predicted to be about 2.8 D, in excellent agreement with the experimental measurement. Calculations on the liquid phase show a similar enhancement. The available evidence demonstrates to our satisfaction that the non-pair-additivity of the water-water interaction is sufficiently large that effective pair potentials cannot be expected to successfully account for both gas and condensed phase properties of water. If we are successfully to simulate water, rather than a hypothetical polar fluid, it would appear necessary to move away from the simpler “ effective” pair potentials to others which can handle the known significant non-pair-additive, or “ cooperative ” effects.This is even more necessary when trying to deal with inter- faces, for example in biological system~.~ These cooperative effects can be handled within our simple classical polarisable electropole formalism, while Berendsen has shown a similar approach is possible using a polarisable four-point-charge model. More work is urgently needed along these or similar lines. Dr. J. C. Dore (Canterbury) said: I would like to comment about the molecular dynamics calculations for water in relation to the structural studies made on heavy water using neutron diffraction techniques. Gibson and I have recently completed a series of measurements at different temperatures to investigate the change in structure for D20 but I wish to compare the experimental results with the simulations at room temperature.Fig. ‘7 shows the intermolecular pair correlation function g(r) obtained from a transform of the neutron data and also the latest Rahman and Stillinger predic- T. R. Dyke, K. M. Mack and J. S. Muenter, J . Chem. Phys., 1977,66,498. J. L. Finney, in preparation. J. L. Finney, in Water: A Comprehensive Treatise, ed. F. Franks (Plenum, New York), vol. 6, in press. H. J. C. Berendsen, in Report of CECAM Workshop on Molecular Dynamics and Monte Carlo CaZcuZations on Water, ed. H. J. C. Berendsen (Centre Europeen de Calcul Atomique et Mole- culaire, Orsay, 1972), p. 22.GENERAL DISCUSSION 83 1.2 1.0 0.8 - L 2 0.6 0.4 0 .2 0 .o 1 2 3 4 5 6 7 r / i FIG.7.-Intermolecular pair correlation function g ( r ) obtained from molecular dynamics at 29 "C (-) and from neutron data at 21 "C (. . .). The g ( r ) function for neutrons is a weighted sum of tions' using the central force model. The g(r) distribution function for neutrons involves a weighted sum of &D(r), go,(r) and goo(r) functions which are shown separ- ately for the MD calculations. It is immediately apparent that the experimental data do not exhibit the sharp features at short distances which are predicted by the model simulation. This is presumably because the neutron data, unlike the X-ray data, are sensitive to orientational correlations between neighbouring molecules. It appears that the potential used in the calculations produces too much tetrahedral ordering in the local environment when compared with the experimental data.I would therefore like to support Finney's comments and add a fourth item to his list of discrepancies. The Rahman-Stillinger model has been of immense value in understanding the proper- ties of water but it is important to recognise its limitations and it may be misleading to put too much reliance on the details of the dynamics until better agreement is obtained with the time-averaged structure. gDD(r) (- * - '>, gOD(r) (- - -) and gOO(r) (- - -). Dr. G. Pilinkis (Budapest) said: The paper of McDonald and Klein classifies the H20 potentials into two classes, those which describe and those which do not describe the second peak in goo.Perhaps a more useful classification could be achieved on the basis of all three pair correlation functions goo, goH, gHH. Last year we determined these three functions experimentally, combining Narten's X-ray and neutron data with our electron diffraction data on heavy water.2 A qualitative comparison between experimental goo, gOD, gDD pair correlation functions and the corresponding theoretical functions calculated by Stillinger and Rahman from ST2 potentials by MD and computed by Lie et al. from CI potentials by MC can be seen in fig. 8. Irrespective of the discrepancies in the 0-D, D-D second peaks, I think there are more serious problems involved in long range 0-D and D-D interactions. Diffrac- J. H. Stillinger and A. Rahman, J . Chem. Phys., 1978, 68, 666.G. Palinkas, E. Kalman and P. Kovacs, Mol. Phys., 1977, 34, 525.a4 GENERAL DISCUSSION ~ ~ ~~~ 2 4 6 8 FIG. 8.-Pair correlation functions (a) goo(r), (b) goD(r) and (c) gDD(r). (-) CI, (- - -) ST2 and ( a . a ) experimental. tion techniques seem to offer more structured water than obtained by computer simulation, which is probably due to long range interactions. I should like to add the following to Dore's comment. If we compare in fig. 8 the first peaks in simulated and separated pair correlation functions, the experimental partials apparently show more ordered structure in the local environment. These later functions reproduce, however, the experimental total neutron g(r) function with- out sharp features at short distances. This proves that the featureless nature of the total neutron g ( r ) is an insufficient criterion for the degree of local order.We studied tetrahedral order in ST2 water by molecular dynamics simulation of aqueous alkali halide solutions with ST2 solvent interacti0ns.l p 2 For the structure determination all pair correlation functions were derived and several orientational probability densities computed from inspection of the frame of particles. The coordinate system for water molecules was chosen with hydrogen atoms in the lower half of the yz plane symmetrically arranged around the - z axis and with the centre of mass at the origin. In this coordinate system the spherical angles of the four tetrahedral directions are 8 = 8T,2, n - 8T,2; q~ = 0, n/2, 3n/2, 2n, respectively. In fig.9 and 10 can be seen G . Palinkas, W. 0. Riede and K. Heinzinger, Z . Naturforsch., 1977, 32A, 1 137. * G. Palinkas and K. Heinzinger, to be published.GENERAL DISCUSSION 85 polar diagrams of statistics on spherical angles 8, q of radius vectors of the four nearest water neighbours of a solvent molecule. These distribution functions are proportional to the probability of finding one of the four nearest neighbours in a given direction. In the figures the full lines correspond to a 0.5 mol kg-l solution. Recalling spherical angles of the tetrahedral directions the figures show clear tetrahedral sym- metry. At the same time, however, the conditional probability of finding all four nearest neighbours coincidentally in tetrahedral directions even at accuracy of 10-1 0" was zero, during the long MD run.FIG. 9.-Polar diagram of Pww(@. n: 2 c 0 3% 2 FIG. 10.-Polar diagram of Pww(q). -86 GENERAL DISCUSSION It may be concluded that neither an exact nor even an approximate tetrahedral structure is attained instantly. The ST2 water shows tetrahedrally symmetrical local order, only in the particle and time average, and this symmetry is the property of only two-particle probability densities. Dr. J. L. Finney (London) said: The 0-D partial correlation function obtained by Dore is consistent with several other pieces of evidence which also suggest that the popular ST2l and Lemberg-Stillinger central force’ potentials overemphasise the tetrahedrality of the water-water interaction, and hence lead to models which are too highly structured.For example : (1) quantum mechanical calculations on the isolated water molecule suggest a very poor separation of the two sets of lone pair electron^,^^^ certainly very much less than is allowed for in classical models. More- over, (2) examination of the geometries of the water molecules in the large number of crystalline hydrates whose structures have been determined at high resolution by neutron diffraction suggests that the structure-determining role of the water lone pairs is much less than that of the hydrogen atoms on the m ~ l e c u l e . ~ ~ ~ If a molecule con- tains a donatable hydrogen, then normally it will participate in a hydrogen bond via that proton. Finally, (3) a detailed examination of the geometries of ordered water molecules found in ex- tensive networks in certain well-ordered protein crystals 4*6 shows extensive 3-co- ordination, consistent again with a lower degree of tetrahedrality in the bonding, and a lack of strong directional control from the poorly-separated lone pair electrons.Moreover, there is somewhat more circumstantial evidence for poor tetrahedrality in the hydrogen-bonding to solvent of the lone pairs of those OH groups exposed on the surfaces of protein^.^ This evidence leads us to question the high degree of tetrahedrality generally accepted to be characteristic of the water-water interaction. In arriving at a highly tetrahedral model, a lot of weight is normally given to the strong tetrahedrality of the coordination in the ices, e.g. by Bernal and Fowler in their classical paper.8 Yet in ice there are other factors to be considered.The water molecules, each potentially capable of participating in four hydrogen bonds-being the proton acceptor in two, proton donor in two-need to organise themselves in a minimum free energy structure with long range order. It seems plausible that the double-donor-double-acceptor capability, under the constraints of long-range order and tetrahedrally-placed protons should adopt the fully tetrahedral structures we observe, in which the lone pairs play an essentially passive role. The directional control is exerted primarily by the protons, with the lone pairs consistent with, but not a significant driving force towards, the resultant tetrahedrally-coordinated crystal structure. The same strong tendency is not borne out for the lone pairs.Prof. D. Chandler (Illinois) said: I have two remarks to make. First, the pheno- mena that Finney has discussed, if true, are described theoretically by the cavity distri- F. H. Stillinger and A. Rahman, J. Chem. Phys., 1972, 57, 1281. A. Rahman, F. H. Stillinger and H. L. Lemberg, J. Chem. Phys., 1975,63, 5223. e.g., see G. H. E. Diercksen, Theor. Chim. Acta, 1971,21, 335. J. L. Finney, in Water: A Comprehensiue Treatise, ed. F. Franks (Plenum New York), vol. 6, in press. I. Olovsson and P-G. Jonsson, in The Hydrogen Bond, ed. P. Schuster, G. Zundel and C. Sandorfy (North-Holland, Amsterdam), 1976, 3, 393. J. L. Finney, Phil. Trans. B, 1977, 278, 3. J. L. Finney, J. Mol. Biol., 1978, 119, 415. J. D. Bernal and R.H. Fowler, J. Chem. Phys., 1933, 1, 515.GENERAL DISCUSSION 87 bution function y*(r) for water. Perhaps we ought to calculate the function from computer simulations. Second, Finney says that real water is less tetrahedral and less structured than com- puter simulated water. Just the opposite conclusion can be developed from the analysis of diffraction data cited by Palinkas. It is my impression that the analysis of neutron diffraction is difficult due to our uncertainty in methods for determining recoil corrections. Perhaps, this uncertainty has something to do with the differences in the conclusions I have cited. Dr. J. C. Dore (Canterbury) said: The inelasticity problems encountered in the analysis of neutron diffraction data arise mainly from the molecular form-factorf,( Q ) or h,(k).The intermolecular contribution is evaluated from the observed diffraction pattern through the relation D,( Q ) = S,( Q ) - f'( Q) and is then transformed to give the weighted pair correlation function d,(r) = 4npr[g(r) - 11. The possible errors arise in fitting the fi( Q ) function to the measured da/dR(O, A) data at high Q-values. In our case, the use of short wavelengths (0.70,0.50 and 0.35 A) may be used to obtain a good fit for.f,(Q) using apparent bond-length parameters which differ from the real values due to molecular recoil effects in the interference function. However, these effects have only a small effect on Dm(Q) since the corrections are much smaller at low Q-values. The lack of structure in g ( r ) at low r-values shown in our data is also apparent in the earlier results of Nartenl although there are some differences in the relative heights of the shoulder and first peak.In order to reproduce the shape given by the molecular dynamics calculations one would expect to see either oscillations in D,( Q) extending to higher Q-values or sharper features in D,( Q) over the main peak. These are almost certainly incompatible with the experimental data but comparisons in Q-space should be made. Palinkas and collab- orators have used their electron data2 with Narten's X-ray and neutron data to obtain the partial functions. In their analysis the g D D ( r ) and & D ( r ) have oscillations extend- ing to large r-values but they have opposite sign and cancel out in the neutron pattern.The electron diffraction experiments are also difficult and it would be useful to see how the errors propagate in the separation of the partial functions. A significant dis- crepancy is present in the comparison of the inter-molecular roD distance for D20 which is found to be surprisngly small, (0.94 8.) for the electron diffraction results but cannot be less than 0.96 A for the neutron data; I hope to discuss these points with Palinkis. A1 ternative methods of determining the partial functions are being investigated and it will probably be necessary to over-determine the solution (i.e., more than three independent measurements) before there is complete confidence in the partial distribution functions. The question of long range structure is a separate question.Dr. G. Pilinkis (Budapest) said: Two questions have been raised: The short range or tetrahedral structure and the long range structure of water. In reply, we have two answers to both of these questions, based on MD, MC computer simulations and diffraction data. As I showed in the example of ST2 water, tetrahedral symmetry appears only in a particle and time average, it is a property characteristic only of two-particle probability densities. Such functions are for example the showed P(O), P(p) functions and gxg(r) pair correlation functions measurable by diffraction techniques. In no case could I think the two answers to the first question are not strictly in contradiction. A. H. Narten, J. Chem. Phys., 1972, 56, 5681. G. Palinkhs, E. Kalman and P. Kovacs, Mol.Phys., 1977, 34, 525.88 GENERAL DISCUSSION we find tetrahedra in ST2 water instantly. It can be stated, however, that ST2 water molecules have a strong tendency to form hydrogen bonds. If we now compare the first peaks in the experimental and simulated pair Correlation functions we may con- clude that they are almost in the same place; the only differences are in the resolution and in the area of the peaks. For example, the 0-0 first coordination number is 4 from experiment and 5.5 from MD. It should be mentioned here, that the coordina- tion number 4 is verified by evidence, established by X-rays. The question of long range structure is more complicated. The differences between simulated and diffraction results are more serious and a new problem arises here, the accuracy of the separated atom-atom correlation functions.I think the third peak about 7 A in the 0-0 pair correlation function is not ques- tionable. It is again a piece of evidence, which has been reproduced from Morgan and Warren (1 938) in all X-ray experiments. Therefore, the correlations between the centres of water molecules are extended at least up to 7 A. Of course, this does not involve the presence of clusters in real water. It is sufficient if there are hydrogen bonded chains from three or four water molecules as was shown in the ST2 film, presented by Berendsen. But if the correlation length for 0-0 interactions is about 7 A and is due to hydrogen bonding, then directing forces in water play an important role in structure formation, and it is not difficult to accept the result that 0-D and D-D correlations extend beyond first neighbour interactions.The experimental pair correlation functions determined are the result of the first attempt to obtain them. The functions determined also reflect the given accuracy of all three types of experiments. Intensity data were measured with an accuracy of one percent. Thus, the Plazcek correction was applied to Narten’s neutron data. Errors due to unknown corrections were eliminated through an iterative procedure requiring the absence of non-physical contributions in the experimental total correlation functions. I think the absence of error trends in all three total structure functions proves the absence of serious systematic errors in the final structure functions.No error propagation analysis was performed because only inaccuracies in intensity measurements were known. The analysis of more serious systematic errors has not been carried out yet for diffrac- tion experiments. Electron diffraction experiments are difficult, but I believe they have now reached the accuracy of the other two experimental techniques. Comparison between ED and ND results for the intramolecular OD distance does not show any significant discrepancy. In the ED experiment this was found to be 0.95 A, so for the time being the inaccuracy in the determination of intramolecular distances is not less than 0.02 A in any diffraction experiment. Dore says that the oscillations in go, and gDD partials cancel out in the total neutron correlation function.This might occur, indeed. Note, however, that separation was carried out in k-space and the separated partial structure functions do not cancel out at low k-range in the total neutron structure function. The non-conditional behaviour of the simultaneous equations was checked through their determinants and by a control separation. From the pair correlation functions determined by Stillinger and Rahman by MD all three MD structure functions (XD, ND, ED) could be derived. Using these as “ measured ” functions with the same simultaneous equations and coefficients, we reconstructed all three correlation functions. In fig. 11 comparison between the original and calculated functions can be seen. The good agreement between these values proves that, except for errors in the experimental structure functions, separation does not involve any major error in the determined pair correlation functions.Finally, I have one more argument concerning the accuracy of determined corre- All three data sets were corrected separately.GENERAL DISCUSSION 89 r/LJ .. :: FIG. 1 1 .-Correlation functions (a) gHH(r), (b) go&) and (c) go&). (-) original and (. - .) calcu- lated functions. lation functions. Namely, these reproduce well all three total experimental correla- tion functions. This has also been verified by the recent results of more precise neutron experiments carried out by Dore and coworkers. Prof. J. L. Rivail (Nancy) said: Berendsen's attempt to include long range forces in molecular dynamic calculations through a reaction field is very interesting and the qualitative improvement reported is encouraging.Nevertheless, it seems rather unlikely that one would get a quantitative fit with experimental data in the case of water because of the oversimplifications of the model used. In particular in the delayed reaction field, Berendsen uses a value of cco = 5.3 which is very controversial and which corresponds, in the liquid, to a dipole moment of 1.85 D with a Kirkwood correlation factor very close to 1. Since the model involves a single Debye process, the corresponding value of should be 4.1.' These two values are quite far from the square of the refractive index in the visible (< 2) because of high frequency relaxa- tion processes. This problem seems closely connected with the existence of hydrogen M.Magat, Disc. Faraday Soc., 1967, 43, 145. J. B. Hasted in Water, a Comprehensive Treatise, ed. F. Franks (Plenum, New York, 1972), p. 255.90 GENERAL DISCUSSION bonds, their dynamics and the strong perturbations they produce in the electronic structure of individual molecules. I would suggest a careful investigation of the role of the reaction field in the case of a simpler liquid (such as HC1 for instance). Prof. H. J. C. Berendsen (Groningen) said: Following discussion remarks made informally by Klein, we would like to refer to Levesque et a2.l who investigated the influence of cut-off radius and system size for dipolar hard-sphere liquids. They in- vestigated in one case the effect of a reaction field and found that the reaction field recovers only a part of the dipolar correlation if compared with theoretical results of linearized hypernetted-chain approximation for an infinite system.A detailed Monte Carlo analysis by Adams and McDonald2 of dipolar correla- tions in polar lattices revealed that dipolar fluctuations are consistently larger in simu- lations using Ewald summations than in simulations using reaction fields. They con- clude that Ewald summations lead to anomalously high dielectric constants. Dr. P. Bordewijk (Leiden) said: I want to put forward a suggestion that may solve the dilemma of either having to choose the reaction field to be proportional to the instantaneous value of the electric moment of the region considered on a molecular basis, or to use a convolution integral like the one in eqn (5.1) of the paper of van Gunsteren et al.Concerning the former choice one would object that it is improbable that the field due to the surrounding molecules adapts itself instantaneously to a change of the moment in the inner sphere. The latter choice, however, implies that the field Ei(t) for given Mi(?) is smaller than according to the static dielectric permittivity, while ( E i ( t ) IMi(t)) should depend only on the potential energy of interaction of the molecules, and not on any kinetic term, so that it should be given by the static permittivity, which also depends only on the potential energy of interaction, whereas the frequency dependent permittivity depends also on the molecular motions. A similar argument in favour of the use of the static dielectric permittivity was given in connection with the calculation of the average moment of a macroscopic sphere for given orientation of one molecule by Ki r kwo ~ d .~ 9 Moreoever, the use of eqn (5.1) has as a consequence that (Mi(t) - Ei(t + z)) is not symmetric with respect to z = 0, although this average is built up from contribu- tions ( p i ( t ) - Tij - p j ( t ) ) that should be symmetric in any event. As an alternative, I suggest one ought to consider, apart from the instantaneous value of the reaction field, i.e., (Ei(t)lMi(t)), which is given by the static permittivity, also the remaining part of Ei at time t. This contribution to the field, although uncorrelated with M,(t), is nevertheless correlated with Mi(t + z) since it influences the evolution of Mi.Writing AE,(t) = Ei(t) - fMi(t) for the latter contribution to the field, one ob- tains : If the delay of the reaction field is accounted for, the right-hand side is non-zero for z > 0, yielding a non-zero correlation between AEi(t) and [Mi(t + z) - Mi(t)], in D. Levesque, G. N. Patey and J. J. Weis, Mol. Phys., 1977, 34, 1077. D. J. Adams and I. R. McDonald, Mol. Phys., 1976,32, 931. J. G. Kirkwood, J. Chem. Phys., 1939, 7 , 911. ' P. Bordewijk, Chem. Phys., 1978, 33, 451.GENERAL DISCUSSION 91 such a way that (M,(t + z) - Ei(t + z)) keeps the value obtained from the static per- mittivity.' Prof. H. J. C. Berendsen (Groningen) said: The question whether time-reversal symmetry should apply to the reaction field is of fundamental importance.There is no doubt that the time correlation function of the dipole pair interaction (,ui(t) Ti,p(t + z)) should be symmetric in z. The reaction field, however, approximates the field XTijp,. In this approximation there is no symmetry between dipoles inside and outside the cut-off radius, and no time reversal symmetry is required for the interaction function (Mi(t)Ei(t + z)}. In fact, if the Onsager reaction field is calcu- lated for a time dependent dipole in a cavity in a medium with frequency-dependent dielectric constant, a delayed reaction field is obtained. In this respect our treatment follows closely that of Linder.2 The question is not whether the reaction field is delayed, but whether the delayed reaction field is a good approximation to the total long-range interaction.The suggestion of Bordewij k that the quasi-random com- ponents of the long-range fields produce a systematic effect because Mi(t) is dependent on previous " random " fields, is worth investigating. Since the delayed reaction field is dissipative, i.e., it removes energy from the system, a fluctuating term is needed to compensate this dissipation in order to generate a stationary process. Whether the combination of delayed reaction field and " random " field is equivalent to the static reaction field is not obvious and remains to be investigated. The argument that (E,(t)M,(t)} should depend only on the potential energy of interaction and not on any kinetic term, although correct in itself, does not prove that our expression (6.2) for the expectation value of the reaction field is incorrect.The dielectric relaxation time does not occur in eqn (6.2). The occurrence of is simply the result of the inconsistent use of For the model E~ = 1 with eqn (6.2) reduces to in a non-polarisable model. 1 &,-1 4mORc3 E, + 2 (El} = - - Mi. This deviates from the Onsager formula but is consistent with the requirement that the expectation value of the reaction field is independent of dynamic quantities. Prof. K. Singer (Egham) said: In a molecular dynamics simulation the periodic boundary conditions alone produce unphysical correlations whether or not the Ewald method is used in the evaluation of the forces. Prof. H. J. C. Berendsen (Groningen) said: This is certainly true; unfortunately there is no alternative without the introduction of other boundary artefacts.But we feel that periodic correlations should not be enhanced by the inclusion of long-range forces from periodic images. Dr. D. J. Adams (Southampton) said: I should also like to make some remarks concerning reaction fields. I am puzzled by the massive difference between the results for the moment fluctuations obtained with the " momentary " reaction field and the delayed reaction field, shown in fig. 4 of the paper by van Gunsteren et al. According to them a change of a factor in the " momentary " reaction field from 0.98 to 0.96 C. J. F. Bottcher, Theory of Electric Polarization, vol. I, Dielectrics in Static Fields (Elsevier, Amsterdam, 2nd edn, 1973), p. 130. B. Linder, in Intermolecular Forces, Adv. Chem.Phys., ed. J. Hirschfelder (Interscience, N.Y., 1967), vol. 12, p. 225.92 GENERAL DISCUSSION should reduce the ensemble average (M2)/Np2 very considerably and I find this difficult to believe. Have they tried this? It could be tried either with their existing molecular dynamics programme or in a Monte Carlo calculation. Bordewijk's suggestion of a stochastic reaction field seems a good idea to me. There is a bias in the delayed reaction field. In an actual polar fluid there will be a strong coupling between a spherical region and its surrounds but this does not take the form of the surrounds always following the sphere moment and lagging behind, sometimes the sphere will lead and sometimes lag. Consider a thought experiment. A spherical region in a polar fluid is observed and the field produced by its surrounds at the sphere centre is monitored as a function of the sphere's momentary dipole moment.In this way an average value of the instantaneous reaction field due to any momentary sphere dipole moment can be found, it will be the equilibrium reaction field of classical electrostatics. This argument was first advanced by Frohlich.' Obviously both the " momentary " and delayed reaction fields are approximations to the true field due to the surrounding fluid but I contend that the " momentary " reac- tion field is the better approximation. Energy will no longer be conserved as there will be dissipation or dielectric loss. This effect is hidden in the present study by frequent velocity scaling to maintain constant energy.The authors have rejected the use of infinite-summation methods. I have made Monte Carlo calculations using both the reaction field and the Ewald-Kornfeld summation for direct comparison with the results of Levesque, Patey and Weis2 for the fluid of dipolar hard spheres. It was found3 that the results for h,(r) and h&) from the reaction field and the Ewald-Kornfeld summation are fairly close to each other and are closer to the results of the LHNC theory than the summation methods evaluated by Levesque et aL2 In particular the reaction field result for ha(r) using 500 dipoles is very close to that found with 256 dipoles using Ladd's method of in- finite ~ummation.~ For the fluid of dipolar hard-spheres at least there is no great difference between reaction field and infinite summation.It has been sug- gested6 that ordinary Monte Carlo is not a very suitable method for studying highly associated liquids. The simulation is liable to become trapped in one region of phase space instead of adequately sampling all important regions. Just possibly this may be the explanation. There is a technical problem also with a delayed reaction field. The goo@) curve for water found by Ladd5 remains a mystery. Prof. K. Singer (Egham) said: The suggestion seems to have been made that the thermodynamic properties of molecules in the condensed phases are rather insensitive to the pair potential. If this is so, it is hard to understand why neither for the simplest type of molecule, N2 (for which both liquid and solid have been extensively investi- gated), nor, to my knowledge, for any other molecule, a pair potential has been found which adequately reproduces the thermodynamic properties of the liquid and the crystalline phases.Prof. J. N. Murrell (Sussex) said: It would seem appropriate at an early point in this discussion to state the present situation regarding intermolecular potential func- tions, particularly as a rather gloomy view has already been presented. H. Frohlich, Theory of Dielectrics (Oxford Univ. Press, London, 2nd edn, 1958). D. J. Adams, Mol. Phys., 1979, to be submitted. A. J. C. Ladd, Mol. Phys., 1978, 36, 463. A. J. C. Ladd, Mol. Phys., 1977, 33, 1039. C. Pangali, M. Rao and B. J. Berne, Chern. Phys. Letters, 1978, 55,413. * D. Levesque, G. N. Patey and J. J. Weis, Mol.Phys., 1977, 34, 1077.GENERAL DISCUSSION 93 Only gas phase experiments or calculations are capable of giving accurate pair potentials. The two experiments which appear to me to have the greatest sensitivity to the potential are low temperature (supersonic nozzle) spectroscopy, which is sensitive to the attractive regions of the potential, and molecular beam inelastic rotational scattering, which is sensitive to the attractive and moderately repulsive parts of the potential. Although neither of these methods can give a direct method of determin- ing the potential from experimental data, the mathematical techniques for predicting experimental results from an assumed potential are now advancing to such a point where trial and error fitting is a feasible approach.For example, the " sudden " approximation for inelastic scattering is computationally economical and appears to be accurate over a wide range of conditions. The non-pair contributions to the potential should be deduced from properties of the solid; primarily structure but perhaps also the phonon distribution. One must start with the pair potentials deduced from gas phase data and deduce three-body terms (in the first place) which are needed to give the correct solid state properties. Except possibly for ionic solids the pair potentials should give a good first approxima- tion to the potential of the solid. Ionic solids may be the exception in that polariza- tion forces, which are not pair additive, may be important. Experiments on the liquid state are not likely to be the source of accurate inter- molecular potentials because the sensitivity of liquid state data to the potential is much less than the best gas and solid data. The ab initio evaluation of accurate intermolecular potentials from quantum mechanics is now feasible for atoms and diatomic molecules. This is not to say that such calculations are trivial or even straightforward. I would place a problem like the NTN2 potential at the limit of present theoretical techniques and computer tech- nology but certainly solvable if the desire for a solution is strong enough. I would be much more sceptical about our ability to produce accurate potentials for more complex systems such as water, although many amongst the computer simulation fraternity appear to be happy with such potentials that have been derived by ab initio met hods. Dr. E. B. Smith (Oxford) said: The dangers inherent in using traditional methods to extract information about intermolecular forces from experimentally measured quantities have been stressed. Often " force fitting " to inadequate analytical func- tions can yield misleading results. However, recently data inversion methods have been devised to unfold the potential energy function, U(r), from thermophysical properties which have a number of distinct advantages over more traditional procedures. Most importantly, they yield only that specific information about U(r) which is contained in the measurements under analysis. Correspondingly, this information is obtained in numerical form, uncon- strained by possibly inadequate analytic functions. Furthermore, unlike more formal inversion procedures, they demand only realistically attainable precision in the input data. When applied to second virial coefficient measurements, they enable the re- pulsive branch of the potential, and the well-width as a function of its depth, to be characterised.l Gas transport property data covering the normally accessible tem- perature range may be inverted to give a major portion of the potential energy func- tion.2 When combined with older inversion techniques such as the Rydberg-Klein-Rees method which may be used to analyse the U.V. spectrum of dimers, the methods al- l D. W. Gough, G. C. Maitland and E. B. Smith, MoZ. Phys., 1972, 24, 151. G. C. Maitland and E. B. Smith, Mol. Phys., 1972, 24, 1185.94 GENERAL DISCUSSION lowed a very direct attack on the intermolecular forces of monatomic species.l They have also been shown to be effective in the study of ion-neutral interactions.2 Attempts are now under way to develop inversion methods for non-spherical potential energy functions. If these are successful, they should place the intermolecular forces of polyatomic species on a more secure basis. For reviews see: E. B. Smith, Physicu, 1974,73,211; G . C . Maitland and E. B. Smith, Seuenth Symposium on Thermophysicul Properties (Amer. SOC. Mech. Eng, New York, 1977). L. A. Viehland, M. M. Harrington and E. A. Mason, Chem. Phys., 1976,17,433.

ISSN:0301-7249    

DOI:10.1039/DC9786600071

出版商:RSC    

年代:1978

数据来源: RSC

摘要:

Molecular Dynamics of Liquid Alkanes BY JEAN-PAUL RYCKAERT~ AND ANDRB BELLEMANS Facultk des Sciences, UniversitC Libre de Bruxelles, Belgium Received 19th May, 1978 The method of molecular dynamics is applied to the simulation of liquid systems of n-alkanes. The model used is a semi-rigid one with fixed C-C bonds and C-C-C angles. In addition to the static and dynamic properties usually deduced for monatomic fluids from such computer experiments, the configurational properties and the internal relaxation behaviour of the alkane chain are also studied. The results of two different simulations of n-butane and of one simulation of n-decane are analysed. The usefulness and the limitations of such computer experiments are discussed briefly. 1. INTRODUCTION Computer simulations of liquid systems at the molecular level are the primary source of non experimental information on their structure and dynamical properties.As far as chain molecules are concerned, Monte Carlo methods were first developed but from such experiments only static properties are accessible, although some dy- namical information can be derived by simulating a brownian motion of the chains. In order to study the dynamical properties of an assembly of chain molecules in a realistic way, it is necessary to use the second traditional method of computer simula- tion, i.e., the method of molecular dynamics which, originally applied to hard spheres,l has been extensively used in the last few years for treating more and more complicated systems, including water.2 In the present paper we report the results of a simulation of short n-alkane liquids which, as far as we know, is the first study of this kind by the method of molecular dynamics. The n-alkane molecules were chosen on account both of their simple structure and their well known thermodynamic properties.In section 2 we describe the model used while section 3 is devoted to a brief de- scription of the various computer experiments: two simulations of liquid n-butane a t different temperatures and one simulation of liquid n-decane. (Results of a prelimin- ary experiment on n-butane with somewhat different interaction parameters have been published previou~ly.)~ Technical details are omitted because the specific methods for simulating flexible chain molecules have been recently described in specialized paper^.^*^ In section 4 we compare the computed thermodynamic properties of n-butane and n-decane with their actual values; section 5 deals with microscopic static properties and section 6 with dynamical properties.2. MODEL OF n-ALKANES The skeleton of a molecule is made of n points coinciding with the C nuclei, all C-C bonds are rigid with a common length I equal to 1.53 A and the angles between adjacent bonds are fixed at 109" 28'. The potential energy associated with the relative rotation of the two parts of a chain adjacent to a C-C bond is a function of the t Present address : Laboratory of Physical Chemistry, The University of Groningen, the Nether- lands.96 MOLECULAR DYNAMICS OF L I Q U I D ALKANES correspodding dihedral angle a ; this function V ( a ) was constructed from the available data on the internal rotation of n-butane: 6 p V(a)/k = [l .116 + 1.462 cos a - 1.578 C062 ct - 0.368 cos3 a + 3.156 c0s4 ct - 3.788 c0s5 a] lo3 K (2.1) where k is the Boltzmann constant. The corresponding extremal energies (expressed in kcal mol-l) are V(0) = 0 (trans configuration) V(&n/3) = 2.95 V(-l2~/3) = 0.70 (gauche configuration) V('Z) = 10.7.All CH2 and CH3 groups are treated as identical united atoms and, as their interaction will be described by means of a Lennard-Jones potential, we call them L-J atoms; e.g., n-butane and n-decane are respectively treated as linear chains of four and ten L-J atoms. Two L-J atoms belonging to different molecules interact by means of the Lennard-Jones 6-12 potential with the following parameters Elk = 72 K, o = 3.923 A (2.2) (These values differ slightly from those used in our previous simulation of n-b~tane.)~ In a similar way L-J atoms belonging to the same molecule and separated by at least three other L-J atoms, interact through this same Lennard-Jones potential.(This obviously does not happen for n-butane). As we do not consider the H atoms explicitly, the number of degrees of freedom of an n-alkane C,H2,+ is equal to n + 3 in our model, corresponding to three translations and n rotations; all of them are classical and for simplicity the mass of CH3 is taken to be equal to the mass of CH2. 3. MOLECULAR DYNAMICS EXPERIMENTS As usual we considered N molecules enclosed in a cube of side L with periodic boundaries.In all experiments we adopted the same cut-off distance rc for the Len- nard-Jones potential: rc = 2.5 O. This distance was always smaller than L/2 in order to prevent a L-J atom from interacting more than once with any other atom and its images. Moreoever, given the largest end-to-end distance r,,, of a molecule, care was taken to fulfill the following relation which guarantees that an atom of a given molecule cannot interact: direclly with any atom pertaining to one of the images of this same molecule. Systems of 64 and 27 molecules were respectively used for simulating n-butane and n-decane. For n-butane the equations of motions in Lagrangian generalized coordin- ates (seven degrees of freedom) were integrated step by step by means of an algorithm due to Gear.* Later it was found that for long chains the Cartesian equations of motion are easier to work with, even though care has to be taken with the constraints arising from fixed bonds and angles at each step of the integration ; it was verified on n-butane that the computing efficiency is about the same as for generalized coordinates.Therefore, we elected to work with Cartesian coordinates for simulating n-decane. A starting configuration is generated by placing the centres of mass of all molecules on a simple cubic lattice, matching the desired density, and by randomly disposing the skeleton of these molecules around that The simulations were initiated as follows.J . P. R Y C K A E R T A N D A . BELLEMANS 97 point (respecting of course the requirement of fixed bonds and angles).All velocities are set equal to zero and the o value of the Lennard-Jones potential is treated as a free parameter. At the beginning o is fixed at such a value that the potential energy of interaction between the molecules is close to zero. As the integration of the equations of motion proceeds, the kinetic energy increases and when it reaches a first maximum, all velocities are reset to zero while a new (larger) value of cr is selected in order to again nullify the potential energy. This process is repeated until 0 reaches its actual value [eqn (2.2)] and the system is then allowed to age to equilibrium. In the next sections we analyse the results obtained from three experiments, two on n-butane (B1 and B2) and one on n-decane. Their main characteristics are listed in table 1.The reduced time-step of integration was equal to 0.002 o(m/&)* TABLE 1 .-CHARACTERISTICS OF THE COMPUTER EXPERIMENTS n-a1 kane thermodynamic number of time step length of method state molecules s simulation n-butane (Bl) liquid at -290 K 3 64 3.9 14 ps generalized coord. n-butane (B2) liquid at -200 K 3 64 3.9 17 ps generalized coord. n-decane liquidat -480K 27 3.9 19 ps Cartesian coord. where rn is the mass of CH,; this corresponds to 3.9 x 1O-l' s in real time. of the total energy was observed with such a time step. No drift 4. THERMODYNAMIC PROPERTIES The two simulations of n-butane were realized respectively near its normal boiling point (273 K) and ai about half way between its melting point (145 K) and its boiling point.The simulation of n-decane corresponds to a temperature somewhat higher than its normal boiling point (447 K). Equilibrium averages of energy, temperature and virial are listed in table 2 for all three experiments. Recall the virial expression for polyatomic molecules A' 1 3(pV - NkT) == (2 Fk Rk) (4.1) where F k and Rk denote the total force acting on and the vector position of the centre of mass of molecule k. Due to the cut-off of the Lennard-Jones interactions, the potential energy and the virial computed at each step of the simulations are over- estimated and a correction was introduced by assuming the density of L-J atoms to be uniform beyond the cut-off. Note also that in table 2 the ratios of the rotational and the translational kinetic energies are close to their theoretical equilibrium values (4/3 for n-butane and 10/3 for n-decane): this can be considered as a good indication that the systems as a whole were at equilibrium.In table 3 we re-express some of the data of table 2 in conventional units by means of the values [eqn (2.2)] of E and cr and we compare them with the actual properties of n-butane and n-decane. The quantities EL and E$ are respectively the molar con- figurational energies of the liquid and of the perfect gas. The experimental values of EL - E$ were obtained from vapour pressure data,9 corrected for gas imperfec- tions.lO*ll The values of EL - E$ corresponding to the simulations are given by U i n t e r + Uintra - Uintra (perfect gas). ( 4 498 MOLECULAR DYNAMICS OF L I Q U I D ALKANES For n-butane we obtained from eqn (2.1) by'means of the Boltzmann formula: [Uintra (perfect gas)]/N~ = 4.21 at 291.5 K = 2.71 at 199.9 K. We note that Uintra is higher for the (simulated) liquid than for the gas, indicating that the probability of gauche conformations is enhanced in the dense phase; this (4.3) TABLE 2.-EQUILIBRIUM AVERAGES OF ENERGY, TEMPERATURE AND VIRIAL (IN REDUCED UNITS) FOR ALL THREE COMPUTER EXPERIMENTS thermodynamic quantities B1 B2 total energy E/NE (constant of motion) intermolecular potential energy Uinter/NE -11.06 intramolecular potential energy Ui"tra/NE 5.14 - 30.37 -32.59* temperature k TIE 4.05 virial ( p V - Nk T)/NE 0.96 ratio of rotational and translational kinetic energies 1.41 -3.48* - 22.95 - 36.47 - 39.41 * 3.80 2.78 4.06 1.35 -1.83* n-decane 13.26 36.03 -68.19 - 74.06* 6.68 -9.83 -21.57" 3.29 ~ * Cut-off correction included.point will be discussed in the next section. gas) was estimated by simulating the gas phase; this led to the approximate value For n-decane the value of Uintra (perfect The experimental values quoted in table 3 for n-butane and n-decane correspond to The values from eqn (2.2) of e and o states situated on the gas-liquid coexistence line. TABLE 3.--COMPARISON OF THE SIMULATED THERMODYNAMIC DATA WITH ACTUAL PROPERTIES OF n-BUTANE AND n-DECANE thermodynamic simulation liquid simulation liquid simulation n-decane at quantities B1 n-butane B2 n-butane of n-decane boiling point temperature/K 291.5 290 199.9 203 481.0 447.3 molar volume/cm3 99.7 99.7 86.7 86.7 236.5 236.2 (EL - G$)/RT -7.82* -8.12 -13.8* -14.3 -10.5* -9.9 * Cut off correction included.were actually adjusted in order to gel. a good fit with experiment B1. (This kind of adjustment is not as simple as for monatomic liquids because E and Q are not the only parameters of the model: there are also the internal potential V(a) and the C-C bond length). The fact that the same values of E and o lead to an equally good fit of the data for experiment B2 is very gratifying. The agreement observed for n-decane is not quite as good but is still reasonable. It could obviously be improved by adopting values of E and Q differing slightly from those of n-butane. A more realistic approach would be to differentiate the end L-J atoms corresponding to CH, from the ones corresponding to CH2.12J .P . RYCKAERT AND A. BELLEMANS 99 5. MICROSCOPIC STATIC PROPERTIES PAIR DISTRIBUTION FUNCTION OF c NUCLEI Fig. 1 and 2 show the pair distribution function of C nuclei (i.e., L-J atoms) as a function of their separation r. Inter- and intra-molecular correlations are mixed. Peaks at 1.53 and 2.49 A are related to rigidly connected pairs of C atoms (first and second neighbours in each molecule) while peaks labelled T and G respectively corre- spond to trans and gauche positions of fourth neighbours inside a given molecule. FIG. 1.-Pair distribution function of C nuclei in n-butane (experiment B1, 291.5 K). FIG. 2.-Pair distribution function of C nuclei in n-decane (481 K). n-decane peaks related to successive TT and GT conformations can also be seen; note the absence of peaks corresponding to two successive gauche conformations : such structures are practically forbidden because they would bring together two L-J atoms at distances much less than o (2.5 A for G+G- and 3.5 A for G+G+ or G-G-).As far as the intermolecular correlations are concerned, it is obvious that the oscilla- tions around g(r) = 1 for r larger than o are much weaker than for monatomic liquids. CONFORMATION OF MOLECULES The distribution of the n-butane molecules as a function of the internal angle cc is shown in fig. 3 for experiment B1, where the peaks observed at cc = 0 and cc = &2n/3 correspond respectively to trans (T) and gauche (G+ and G-) conformations. The magnitudes of the G + and G- peaks are very different, in contrast to the sym- metry of the potential energy V(E).What actually happened is that the initial con- figuration of this experiment was not symmetrical in G+ and G- molecules (excess of1 00 MOLECULAR DYNAMICS OF L I Q U I D ALKANES G+) and this situation was not corrected during the duration of the run (-14ps, real time), in spite of the fact that equipartition was attained between translational and rotational degrees of freedom. Indeed the frequency of transition between T and G molecules (measured by the transfer of a molecule from the central peak to the lateral ones or vice-versa) is rather low: during the entire run, 45 G + T and 43 T --j G transitions were respectively observed for the 64 molecule assembly, i.e., the mean z ; 0. h 0 0 0 0 TT oco -lT FIG.3.-Normalized distribution function of n-butane molecules as a function of the internal angle a (experiment B1, 291.5 K). time interval between two transitions is !z 10 ps. In spite of the fact that a true equilibrium between G+ and G- molecules was not established, it seems reasonable to us to assume that there was equilibrium between T and G molecules on the whole as the observed numbers of G -+ T and T -+ G transitions were practically equal. (Direct transitions between G+ and G- molecules were never observed: they must be x Q 0.5 0 0 0 0 ** O 0 0 0 - T l 2 3 TT 6c 3 0 3 FIG. 4.-Normalized distribution function of n-butane molecules as a function of the internal angle a: 0 , symmetrization of the curve of fig. 3. 0, gaseous phase (291.5 K).J .P . RYCKAERT A N D A . BELLEMANS 101 extremely rare at the temperature of the experiment). The artz$cially symmetrized distribution of molecules as it depends on a is shown in fig. 4, together with the corresponding distribution for the gas phase, computed directly from eqn (2.1) by means of the Boltzmann formula. There appears to be a significant displacement in favour of G molecules when going from the gas to the liquid, in agreement with recent theoretical predictions ; l3 the integration of the peak leads to the following estimations (291.5 K) gas: 66% T, 34% G liquid: 54% T, 46% G (This is in agreement with the preceding observation that Uintra is higher in the liquid than in the gas). The conformation of an n-decane molecule is specified in our model by seven internal angles denoted as al, .. . a7 going from one end of the chain to the other; by symmetry, angles a1 + and a7 - with i = 0, 1, 2 play an equivalent role. As in the case of n-butane, a perfect equilibrium between G+ and G- conformations could not be reached, though the departure was perhaps less since the temperature of the simula- tion was higher (48 1 K). The percentages of T, G+ and G- conformations, estimated as in the case of n-butane, are listed in table 4. An interesting conclusion is that T TABLE 4.-PERCENTAGES OF ttYZnS AND gauche CONFORMATIONS ALONG THE I1-DECANE CHAIN (481 K) G+ T G - Q1, a7 20.4 54.5 25.1 u2, 20.1 63.8 16.1 Q39 Q5 21.6 65.6 12.8 a4 24.1 55.4 20.5 conformations appear less abundant in the centre (a4) and at the end of the chain (a1, a7), presumably on account of distant internal correlations.The conformation of the n-decane molecule can also be characterized globally by its end-to-end distance r and its radius of gyration s around the centre of mass. The equilibrium distribution of r is shown in fig. 5 ; very similar results were obtained for s, from which the following averages were computed: ( r ) = 5.76, ( ( r - cr))')' = 0.71, (r') = 33.64, ((r2 - (r2))')' = 7.73, ( s ) = 2.03 ((s - (s)') = 0.12 (s'} = 4.11 ((s2 - ( s ' ) ) ~ ) = 0.49 (481 K; C-C bond length taken as a unit). ORIENTATJONAL CORRELATIONS BETWEEN MOLECULES Experiments on the depolarisation of diffused light by pure n-alkanes and solutions of n-alkanes have led to the conclusion that to a certain degree there exists a parallel alignment of neighbouring chains at room temperature in the pure liquid phase, when the number of C atoms exceeds six.14 In order to see if such an effect was present in our experiments, we associated a unit vector V' with each molecule, pointing in the longitudinal direction of the chain.For n-butane this vector was oriented along the line joining the mid-points of the two extreme C-C bonds. For n-decane, V,102 MOLECULAR DYNAMICS OF L I Q U I D ALKANES was aligned along the longer principal axis of the inertial ellipsoid associated with the instantaneous conformation of the chain. We next computed the correlation coeffi- cient C = (I V, V,.]) where the average is taken over all pairs of molecules for all configurations of the system. Random orientations of the vectors V, lead to C = 1/2 while strict parallel alignment corresponds to C = 1.Both for n-butane (experience Bl) and n-decane, we obtained values of C which did not deviate significantly from 3 X I cc 7 FIG. 5.-Normalized distribution function f ( r ) of the end-to-end distance of n-decane molecules; r is measured by taking the C-C bond length as unit. Vertical arrows indicate the value of Q and the maximum extension rmax of a molecule; the vertical line corresponds to (Y). 0.5. This was indeed expected for n-butane. The fact that no alignment was ob- served for n-decane molecules either does not really contradict the experimental evi- denceI4 as the temperature of the present simulation (481 K) is much higher. 6. DYNAMICAL PROPERTIES SELF-DIFFUSION The velocity autocorrelation functions of the centre of mass of the molecules are shown on fig.6 , 7 and 8 for experiments B 1, B2 and D, respectively. As already ob- served in our preliminary paper,6 these functions are markedly different from those observed for monatomic liquids like argon or even diatomic and rigid polyatomic ones like nitrogen and water. In all these cases the rapid initial decay of the velocity auto- correlation function is followed by a relatively important negative part, usually inter- preted as a cage eflect: after a time corresponding to the mean collision time, many molecules rebound on their neighbours and reverse their velocity. In the present case, apart from a rather small negative part in experiment B2 (fig. 7), the autocorrelation function remains positive until it enters the region of statistical fluctuations for t larger than We believe that this peculiar behaviour is due to two facts: s.J .P . RYCKAERT A N D A . BELLEMANS 103 1.0 e ” 0 Ob (a) molecular collisions are softened by internal degrees of freedom acting like shock- absorbers and (b) the molecules take advantage of their flexibility to adjust their shape k 0 0 0 O 0 0 0 0 I 0 0 0 0 0 0.5 1.0 1.5 I - ~ ~ - ~ o o o o o o o o o w ~ - - ~ ~ o - 0 4 0 0 0 w o o o FIG. 6.-Velocity autocorrelation function of the centre of mass of n-butane molecules (experiment B1, 291.5 K). FIG. 8.-Velocity autocorrelation function of the centre of mass of n-decane molecules (481 K). to the environment, minimizing the so-called cage effect and easing their diffusion through the surrounding medium.The values of the self-diffusion coefficient D observed for n-alkanes l5 are relatively high if one compares them with the values for spherical molecules under the same conditions (neighbourhood of the normal boiling point). We compare in table 5104 MOLECULAR DYNAMICS OF LIQUID ALKANES the values of D deduced from the computer experiments with those estimated for n-butane and n-decane by means of Houghton’s empirical formula.16 The agree- ment is very satisfactory. TABLE 5.-sELF-DIFFUSION COEFFICENT OF II-BUTANE AND n-DECANE : (a) DEDUCED FROM MEAN-SQUARE DISPLACEMENT OF CENTRE OF MASS,4 (b) OBTAINED BY INTEGRATING VELOCITY AUTOCORRELATION FUNCTION, (c) HOUGHTON’S FORMULA.'^ simulation n-butane (Bl) 291.5 6.1 6.9 6.4 n-butane (B2) 199.9 2.1 2.4 1.9 n- bu tane 48 1 7.5 7.7 6.8 INTERNAL RELAXATION OF THE CHAIN As mentioned in section 5, the internal relaxation of n-butane molecules between T, G+ and G- conformations is a relatively slow process, the mean time interval between transitions being of the order of 10 ps in experiment Bl (291.5 K) which is nearly the same as the total duration of the experiment (14 ps).As the n-decane simulation was performed at a much higher temperature (481 K), the number ofjumps observed between T and G conformations was notably higher, allowing one to dis- cuss the relaxation of the chain in some detail. Table 6 shows the frequencies of TABLE 6.-FREQUENCIES OF INTERNAL TRANSITIONS IN n-DECANE MOLCULES (48 1 K) internal angle average number of average time interval transitions between transisions during lops IPS 5.22 4.86 4.16 5.78 1.92 2.06 2.40 1.73 transitions between G and T conformations or vice-versa, all along the chain.In- tuitively one might expect a larger mobility of the end atoms: this is not confirmed by our results, probably on account of steric effects due to neighbouring molecules and of possible couphgs between internal motions; on the average the mean time interval between two successive G-T transitions is of the order of 2 . 0 ~ ~ . The global relaxation of the n-decane chain may be characterized by the time autocorrelations of the square end-to-end distance and of the square radius of gyra- tion, respectively shown in fig. 9 as pl(t) and p2(t). The corresponding relaxation times are E 3 ps, which implies the occurrence of about 10 internal G-T transitions, a reasonable result.We also studied the angular autocorrelation functions C(t) = (cos [a(?) - a(O)]) and S(t) = (sin [a(t) - a(0)J) for the different angles al, . . . a7. The function C(t) starts from unity and tends ultimately to the asymptotic limit (cos CC)~, This occurred for all angles al, . . . a7 after z 10 ps, i.e., about five G-T transitions. The function S(t) should strictly vanish at all times by symmetryJ . P . RYCKAERT A N D A . BELLEMANS 105 (for an infinite system) and was indeed observed to remain close to zero. As an illus- tration c(t) and s(t) are shown in fig. 10 for the equivalent angles cc2 and FIG. 9.-Normalized autocorrelation functions of the square end-to-end distance square radius of gyration (s2) of n-decane (481 K).<rZ(O)r2(f)) - <r2(0)Y. <s'(0>s2(t)> - <sz(o)>2 ' 9 '1 = +4(0)) - <rZ(O)>z 0, '2 = <s4(0)) - <s2(0)>2 1 (r2) and of the o o o o o o o o o o 0 1 2 3 L 5 6 7 8 9 t / lo& 0. La*," " t o 8 4 8 ' FIG. 10-Autocorrelation functions C ( t ) (full circles) and S(t) (open circles) for the internal angles u2, a6 (n-decane, 481 K). 7. CONCLUSIONS Although it shows that simula- tions of dense systems of chain molecules are feasible with the present generation of computers, it also demonstrates the limitations of this kind of approach. The com- puter time needed is rather prohibitive (more than three hours on an IBM 370/168 for each simulation) so that it seems to be out of the question to construct the whole phase diagram of chain molecules of various length or to investigate the excess thermo- dynamic properties of chain mixtures.To study internal relaxation processes and transport properties seems at present the sensible thing to do : it should permit one to prove or disprove the adequacy of various theoretical models proposed in the litera- ture. Two difficulties are met with longer chain molecules: (a) a broadening of the spectrum of relaxation times, which requires extension of the duration of the simula- tion in real time in order to cope with slow modes and (b) the obvious necessity to en- large the system, which hopelessly increases the number of degrees of freedom to be handled. A possible escape from these difficulties could be found by combining the usual method of molecular dynamics with stoachastic techniques. This investigation is in some sense exploratory. The simulations were made at the CECAM (Facultb des Sciences, Orsay, France) thanks to the hospitality of Prof. Carl Moser, to whom we express our gratitude.106 MOLECULAR DYNAMICS OF LIQUID ALKANES T. Wainwright and B. J. Alder, Nuovo Cimento Suppl., 1958, 9, 116. A. Rahman and F. M. Stillinger, J. Chem. Phys., 1974,60, 1545. J. P. Ryckaert and A. Bellemans, Chem. Phys. Letters, 1975, 30, 123. J. P. Ryckaert, G. Ciccotti and H. J. C. Berendsen, J. Comp. Phys., 1977,23, 327. W. G. Van Gunsteren and H. J. C. Berendsen, Mol. Phys., 1977, 34, 131 1 . R. A. Scott and H. A. Scheraga, J. Chem. Phys., 1966,44, 3054. ' P. B. Woller and E. W. Garbisch Jr, J. Amer. Chem. Soc., 1972, 54, 9310. C. W. Gear, ANL Report No. ANL 7126 (1976). F. D. Rossini, K. S. Pitzer, R. L. Arnett, R. M. Brown and G. C. Pimentel, Selected Values of Physical and Thermodynamic Properties of Hydrocarbons and Related Compounds (American Petroleum Institute, Carnegie Press, 1953). lo J. M. Dymond and E. B. Smith, The Virial Coefficients of Gases (Oxford Clarendon Press, London, 1969). l1 M. L. McGlashan and P. J. B. Potter, Pruc. Roy. SOC. A , 1962,267,478. l2 P. K. Warme and M. A. Scheraga, J. Comp. Phys., 1973,12,49. l3 D. Chandler and L. R. Pratt, J. Chem. Phys., 1976,65,2925. l4 P. Bothorel, C. Such and C. Clement, J. Chim. phys., 1972,69, 1453. D. C. Douglas and D. W. McCall, J. Phys. Chem., 1958,62, 1 102. l6 G. Houghton, J. Chem. Phys., 1964,40, 1628.

ISSN:0301-7249    

DOI:10.1039/DC9786600095

出版商:RSC    

年代:1978

数据来源: RSC

摘要:

Computer Simulation of Gas-Liquid Surfaces : Molecular Fluids BY STEPHEN M. THOMPSON School of Chemical Engineering, Cornell University, Ithaca, New York 14853, U.S.A. Received 5th May, 1978 The molecular dynamics method is applied to the coriiputer simulation of the gas-liquid surface of molecular fluids. Two fluids (nitrogen and chlorine) are modelled with a site-site potential at three temperatures for each system. Surface tensions are in good agreement with experiment and show a significant size effect, while density against orientation profiles show orientational ordering in the surface zone at low temperatures. 1. INTRODUCTION Here we describe the results of molecular dynamics (MD) computer simulations of the gas-liquid surface of two molecular fluids, modelled on nitrogen and chlorine, at three temperatures from the triple point region upwards. Results are presented here for the density against orientation profiles and surface tensions, while results for selected non-equilibrium properties will appear elsewhere at a later date.Computer simulation was chosen as the tool used to obtain these results for two main reasons. First, a number of properties of surfaces are difficult to obtain experimentally, for example, the density against orientation profiles and pair distribution functions. Secondly the theoretical calculation of other surface properties for which defined experi- mental techniques exist, such as the surface tension, is very difficult unless one makes drastic simplying assumptions, as the theory usually involves experimentally inaccessible quantities.All these properties are obtainable in principle by performing computer simulation experiments. So far as we are aware, this is the first work of this type on molecular fluids,l although a body of literature exists on monatomic fluids.2-8 Calculation of the density against orientation profiles enables one to study the preferred orientations (if any) in the surface zone; the effect of molecular shape on surface tension can also be investigated. 2. MOLECULAR MODELS Each fluid was modelled by means of a site-site or atom-atom p~tential.~ Each molecule is assumed to consist of two interaction centres (commonly assumed posi- tionally coincident with the atomic nuclei). The intermolecular potential is then given by 2 2 ~ ( r l 2 ~ 0 1 ~ 0 2 ) = 2 2 uLJS (IY~") - rj(2)l).(2.1) i = l j = 1 Where ri(") is the position vector of site i on molecule k, rI2 is the centre to centre108 COMPUTER SIMULATION OF GAS-LIQUID SURFACES separation, cok = {Ok& are the polar angles specifying the orientation of molecule k, and ~Ljs(r) is a shifted Lennard-Jones (12, 6) interaction lo where uL,(r) is a Lennard-Jones (12,6) interaction uLJ(r) = 4&[(O/r)” - (O/r)6] (2.3) and Y, is the potential cut-off distance. Here E has the dimensions of energy and 0 TABLE 1 .-PAIR POTENTIAL PARAMETERS reference for fluid (&lk)lK alnm Lla parameters NZ 37.3 0 . 3 3 1 0 0.3292 11 Cl* 1 7 3 . 5 0 . 3 3 5 3 0.6080 1 2 has the dimensions of length. The values used are given in table 1. The cut-off distance was given in each case by rc = 2.50 + L.The second term in eqn (2.2) is the usual shifting that results from using a trun- cated p ~ t e n t i a l , ~ ~ ~ while the final term is introduced to ensure that both the potential and its derivatives are continuous at the cut-off. This leads to better dynamics (see below). The results are reduced by means of the parameters E~ and oR, which are defined as the well depth and collision diameter respectively of the spherically averaged pair potential : (2.4) <u(aR, co1c02))wlwz = 0 (2.5) <u(YR, co1c02))wlw2 == -&R (2.6) where rR is the location of the minimum, This enables the results to be directly compared. Values of the parameters These values were obtained by numerical integra- cR, O, and rR are given in table 2.tion. TABLE 2.-sPHERICALLY AVERAGED POTENTIAL PARAMETERS N 2 l3 97.04 0.3754 0 . 4 1 57 0.3292 ClZ 2 3 3 . 3 9 0.4568 0.4984 0.6080 Currently under investigation is the addition of point charges of magnitudes q to the sites and -2q to the centres of mass, where q is chosen to reproduce the known quad- rupole moments of N2 and Cl,.S . M. THOMPSON 109 3. SIMULATION Each system consists of 216 molecules confined to square prisms of dimensions x,, yL (=xL) and z,( > x,) with the usual periodic boundary conditions in the (xy) plane of the surface. The liquid was confined to the centre of the cell with vapour phases at either end. Thus an infinite film is simulated, with two surfaces that are stable without the inclusion of external force^.^-^*^ The cell widths were xL = 2r, (5.65840 and 6.2160 for N2 and C1, respectively) and the cell heights were 190 for both systems at the lower two temperatures and 220 at the highest temperature.The translational and angular velocities of the molecules, ui and mi, respectively, are chosen randomly subject to the constraints that the desired kinetic energy is obtained, this being properly divided between translations and rotations, and that the centre of mass of the system has no resultant translational or angular velocity. At the lowest tem- peratures the molecular centres of mass are assigned coordinates on an f.c.c. lattice structure of the appropriate liquid density, with empty vapour phases. An initial equilibration period of typically lo4 integration steps is rejected, during which the lattice structure melts and the density profile becomes stable.The runs at the higher temperatures are started from the last configuration of the run at the next lowest temperature. A vapour molecule which attempts to leave the cell through one of the end walls is returned to the cell by simply bouncing it off the wall. This procedure is thought to have negligible effect on the surface structure. The position of the centre of mass was monitored during the run, and moved on average a distance of 0.005~ per 10 000 steps. Initially the centre of mass is fixed in the centre (z,,, = zJ2) of the cell. The centre of mass motion was found by solving the equations Fi = M-'Fi (3.1) using a third-order predictor-corrector method. Here ri is the centre of mass vector of molecule i, M is the molecular mass, H i is the force on i due to all the other mole- cules, and the dots indicate differentiation with respect to time.The rotational motion was found by solving the equations hi, = I-'TiP (3.2) using a second order predictor-corrector method, where miP and Ti, are the principal angular velocity of and torque on molecule i respectively, and Zis its moment of inertia. Molecular orientations were expressed by means of quaternion parameter^,'^ giving singularity-free equations of rotational motion with better resultant energy conserva- tion. The equations of motion were integrated with a reduced time step of 0.0015 [ t * = t ( ~ M / a ' > ~ ] . s for N2 and C12, respectively. Total energy was found to oscillate around a fixed value with an amplitude of -0.25% an1 a period of 3000-5000 time steps.The use of a truncated but not shifted pair potential changes this behaviour to a slow monotonic drift. The use of the same algorithm to simulate a bulk phase gave total energy values exact to five significant figures over 2000 steps, thus showing the inhomogeneity of the system to be a cause of numerical instability. The calculations were performed on a PDP 11/70 computer. The computer programs were overlaid and occupied a maximum of 28.5 K 16-bit words of core. Various time-saving schemes were implemented. A run of 50 000 integration steps took 100-170 h, depending on the temperature. The method used to solve eqn (3.2) has been described previou~ly.'~ In real units this is 4.39 x and 3.53 x1 1 0 COMPUTER SIMULATION OF GAS-LIQUID SURFACES 4.PROPERTIES CALCULATED The density against orientation profile p(z, 0) gives the probability density for finding a molecule at height z with the molecular axis at an angle 8 to the z axis.15 It may be expanded as a sum of Legendre polynomials in cos 8: For symmetrical linear molecules pl(z) vanishes when I is odd. The total density profile (centres of mass) is pCOM(z) = [ p(z, 8) sin 8 d8. (4.2) Thus where dl0 is a Kronecker delta. The individual coefficients can be calculated during the simulation from the formulae 3.4) s = l i where the sum over s is over a sequence of Ns time steps, the sum over i is over those molecules with coordinates in the range z - Az < zf < z + Az, and A is the surface area.A value of Az = 0.050 was used in this work. The coefficients I = 0, 2 and 4 were calculated. The surface tension is calculated during the simulation by means of the equation where i and j are molecular subscripts and the r, depend implicitly on i and j . Note that the calculation applies for a truncated Lennard-Jones diatomic fluid, with the exception that the resulting shift in the coexisting densitiesS is not allowed for (see below). If the sum over i a n d j in eqn (4.5) is taken over all molecular pairs, then the result must be divided by two because there are two free surfaces. A correction for the missing tail of the potential was added in a similar manner to that explained previously.8 5. RESULTS The averages obtained during the simulations are given in table 3.The coexisting densities are plotted along with experimental values in fig. 1. Values obtained for a monatomic (1 2, 6) fluid modelling argon are plotted for comparison. Simulation values are low as has been observed previously:s this is due to the shift in the pair potential. While the values for N2 are similar to those for Ar, being only slightly higher, values for C12 are much higher, showing a large size effect. Total density profiles are plotted in fig. 2 and 3 for N2 at 66.5 K and C12 at 171.1 K, respectively. These are averages for the two surfaces in the box in each case; the profiles for the two surfaces were in almost exact agreement with each other. These two profiles are similar in their major details : the " structure " observable in the liquid region changed in form during the simulation and is not considered significant.The estimated error in liquid density in table 3 is based on the size of these oscillations.S. M . THOMPSON 111 Apart from the size of these oscillations being much smaller at higher temperatures, the other profiles are very similar to the ones shown here and are therefore not plotted. Here it is important to note that the real-time duration of the present simulations is less TABLE 3 .-AVERAGES OBTAINED DURING SIMULATIONS Na 42.75 x lo3 66.5 0.887 0.971 1.26 1.20 1.58 1.876 NI 50.0 x 103 77.2 0.803 0.918 0.89 0.94 1.68 2.194 N1 40.0 x lo3 96.6 0.649 0.811 0.46 0.48 3.36 1.755 CI, 50.0 X lo3 171.1 1.335 1.398 2.43 2.54 0.81 1.763 CI, 50.0 x lo3 195.8 1.272 1.350 2.30 2.24 0.99 1.763 C11 50.0 x lo3 218.8 1.211 1.306 2.09 1.95 1.10 1.763 52.0 zt0.035 50.07 h2.2 f0.025 50.04 52.6 zt0.017 50.04 k5.2 h0.066 *0.10 *5.9 f0.036 *0.10 56.9 50.029 f0.05 a Number of time steps used in run proper.Temperature. Reduced liquid density from simulation. * Reduced I Reduced surface tension from experiment. liquid density from experiment. g Reduced surface thickness. * Real time duration of simulation. Reduced surface tension from simulation. 0 a El 0 I . 0.5 0 . 7 '.* kT/ER O.' 1 . o Fig. 1 .-Coexisting liquid densities as a function of reduced temperature (a) Ar - 0 , sirn"; expt; (b) N2 - x , sim"; 0, expt; (c) C12 - x, sim"; 0 expt. than half that of those for Ar in ref. (8) and hence larger statistical error must be ex- pected. A surface thickness may be defined in terms of the gradient of the total density profile at the Gibbs dividing surface.8 Values obtained in this manner are112 COMPUTER SIMULATION OF GAS-LIQUID SURFACES 1 .0 - 0.75 "b" - h, - 5, 0.50- < 0.25 0.0 shown in fig. 4 together with argon values for comparison.* Values for N2 and Ar are very similar, but with an apparently more rapid divergence for the former at - - + 1 higher temperature. Thicknesses for C12 are smaller, which is reflected in the surface tension values (see below). In fig. 5 is plotted the second density harmonic coefficient p2(z) for C12 at 171 K. I /a FIG. 3.-Total density profile for C12 at 171.1 K. The density orientation profile constructed from this coefficient and po(z) is plotted in fig. 6. For N2 there was no evidence of any structure in these curves at any of the temperatures studied although the statistical noise was rather large.This was trueS . M. THOMPSON 4.0 3 . 0 - F -m 2.0 1 . 0 - 0.0 113 - - 0 O x x o I I L 0.7 0.8 0.9 1 .o 8 I * x 0 B FIG. 5.-Second density harmonic coefficient p2(z) for C12 at 171.1 K. The arrow locates the Gibbs dividing surface.114 COMPUTER SIMULATION OF GAS-LIQUID SURFACES 0.0 FIG. I I L 6.-Density against orientation profile po(z) + p2(z)P2(c0s 0) for C12 at 171.1 K. correspond to the asterisks in fig. 5 . to locates the Gibbs surface. The heights P I IS. M . THOMPSON 115 also for the p2(z) coefficients for C12 at 218.8 K and for the p4(z) coefficients for both systems at all temperatures. The p2(z) coefficient for C1, at 195.8 K was similar to that presented here, although of reduced amplitude.It also appeared to be shifted in towards the liquid phase relative to the Gibbs surface. The non-zero values in the bulk liquid region in fig. 5 are statistical noise : here the coefficient differed substanti- ally in both box halves and changed in form as the simulation progressed, while this was not so for the structure in the surface region. A pronounced (height dependent) tendency to adopt preferred orientations is indicated. In the liquid phase at the posi- tive maximum in p2(z), the molecules have a tendency to orient with their axes per- pendicular to the plane of the interface (@ = 0), while at the negative minimum (0.050 on the vapour side of the Gibbs surface) the molecules tend to orient with their axes parallel to the interfacial plane.These results are in qualitative agreement with the predictions of first order per- turbation theory using a spherical reference potential with anisotropic overlap forces.15 A quantiative comparison of simulation and perturbation theory for the potential employed in the simulation is in progress. We find little tendency for the near-spherical molecule N2 to adopt preferred orientations in the surface zone, while this tendency is pronounced for C1, at low temperatures. The surface tensions obtained are plotted in fig. 7 together with those of argon* for comparison. Agreement with experiment is good, although this is probably fortuitous in view of the coexisting densities obtained.Little effect of deviation from sphericity is shown by N,, while again the tendency is pronounced for Cl,. This is intuitively obvious. Probably the inclusion of multipolar forces will significantly alter the surface tension. The author thanks the S.R.C., the National Science Foundation and the Petroleum Research Fund (administered by the American Chemical Society) for support of this work. He also thanks Professors K. E. Gubbins and J. S. Rowlinson and Mr. S. Murad for many discussions. S. M. Thompson and K. E. Gubbins, presented at the Symposium on Computer Simulation of Bulk Matter from a Molecular Perspective, American Chemical Society Meeting, Anaheim, California, March 14, 1978. To appear in the ACS Symposium Series. G. A. Chapela, G. Saville and J. S. Rowlinson, Faraday Disc. Chem. Soc., 1975,59, 22. J. K. Lee, J. A. Barker and G. M. Pound, J . Chem. Phys., 1974,60, 1976. F. F. Abraham, D. E. Schreiber and J. A. Barker, J . Chem. Phys., 1975, 62, 1958. A. C. L. Opitz, Phys. Letters A , 1974, 47, 439. K. S. Liu, J . Chem. Phys., 1974, 60, 4226. M. Rao and D. Levesque, J . Chem. Phys., 1976, 65, 3233. 1133. J. R. Sweet and W. A. Steele, J . Chem. Phys., 1967, 47, 3029. lo W. B. Streett, D. J. Tildesley and G. Saville, as ref. (1). l1 P. S. Y. Cheung and J. G. Powles, Mol. Phys., 1975, 30, 921. l2 K. Singer, A. Taylor and J. V. L. Singer, Mol. Phys., 1977, 33, 1757. l3 S. M. Thompson, D. J. Tildesley and W. B. Streett, Mul. Phys., 1976, 32, 711. l4 D. J. Evans and S. Murad, Mol. Phys., 1977, 34, 327. l5 J. M. Haile, K. E. Gubbins and C. G. Gray, J . Chem. Phys., 1976, 64, 1852. l6 R. C. Reid, J. M. Prausnitz and T. K. Sherwood, The Properties of Gases and Liquids (McGraw- l7 J. J. Jasper, J . Phys. Chem. Ref. Data, 1972, 1, 841. * G. A. Chapela, G. Saville, S. M. Thompson and J. S. Rowlinson, J.C.S. Faraduy ZZ, 1977, 73, Hill, New York, 3rd edn, 1977).

ISSN:0301-7249    

DOI:10.1039/DC9786600107

出版商:RSC    

年代:1978

数据来源: RSC