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.