Odd-even effect in n-alkane systems: A molecular dynamics study

Alternation in various properties of n-alkanes (CnH2n+2) as a function of carbon content (n) is termed odd-even effect. Here, we report a comprehensive molecular dynamics simulation study on n-alkane systems carried out with n ranging between 3 (propane) and 8 (octane), examining the odd-even effect in melting point, density, intramolecular conformational ordering, translational and rotational motion. We observe an odd-even alternation in these properties, but with heptane (n = 7) exhibiting anomalous behavior for all except conformational ordering. Our simulations also show the presence of odd-even behavior in rotational and translational dynamics, below and above the melting point, respectively. The results highlight the role of both molecular shape and the variation in density and their interplay in the origins of the odd-even effect.


Introduction
Several properties of long chain hydrocarbons show an alternating variation with the number of carbon atoms, which is termed the 'odd-even' effect [1][2][3][4][5][6]. In general, oddeven alteration behavior is correlated with the difference in packing density for long ranged periodically packed molecules, wherein the odd numbered alkanes (C2symmetry) do not pack as densely as the even numbered ones (Ci-symmetry) [5]. This effect is correlated with the physical properties of various homologous systems, particularly in n-alkanes (CnH2n+2) and their α-and α, ω-substituents [7 -10]. Also, the effect is more pronounced for shorter chains, with the longer the chain length, the smaller the difference in melting points between consecutive odd and even numbered n-alkanes.
The understanding of odd -even dynamics in alkanes is useful from an industrial application viewpoint. Particularly, in the petroleum industry the transport behavior of alkanes is crucial to understand the chemical processes involved in recovery and refinement [11][12][13][14][15].
Although the odd-even effect is well known and observed in various systems, there is still a lack of understanding of the dynamical basis of odd-even alteration. Many of the experimental studies indicate the requirement for a detailed understanding of the odd-even behavior at molecular level for which molecular dynamics (MD) simulations are well suited [27][28][29][30][31][32][33][34][35][36]. An MD simulation study performed on a 4-n-alkyl-4'-cyanobiphenyl liquid crystal showed a slight odd-even effect in the rotational diffusion coefficient, rotational viscosity coefficient and molecular ordering parameter [27]. Also, an MD study performed on ionic liquid systems is consistent with the behavior observed for the diffusion coefficient and electrical conductivity in ionic liquid systems [30,31]. Another recently reported combined NMR and computational study highlights the importance of dynamic behavior in homologous systems [18]. In the literature, several studies have been reported investigating the correlation between structure and physical properties, such as surface tension, volatility and dissociation energies [6,[32][33][34][35][36].
In this work, our aim is to investigate the odd-even effect on dynamical behavior of molecules in n-alkane systems and further to elucidate the physical properties that drive this effect. We have performed MD simulations on n-alkane systems for n = 3 -8. To understand the dynamical properties of n-alkanes around their melting points, we performed simulations at several temperatures with ± 5 K steps above and below the melting points. Initial estimates of the melting point of each alkane were derived from the experimental work reported in literature. To understand the dynamical properties of n-alkanes around their melting points, we performed simulations at several temperatures with ± 5 K steps above and below the melting points. Initial estimates of the melting point of each alkane were derived from the experimental work reported in literature. With addition of every CH2 group in n-alkane series the intermolecular distance between the end CH3 groups is increased [37]. This in turn influences the molecular structure in the liquid state. Probing several different quantities related to the structure and dynamics pertaining to both orientational as well as positional order, helps us point to those properties that are the most dominant contributors to the origin of the effect.

Simulation Details
We have carried out classical MD simulations using the DL-POLY classic [39,40] on a series of n-alkanes (n = 3 -8), using the TraPPE-UA convention [38]. The n-alkanes (n = 3 -8) selected are listed in Table 1. A united-atom and semi -flexible model was used for all the n-alkane molecules [38], where all the CH3 and CH2 groups in a molecule are treated as a single unit. The TraPPE-UA force field was used, which provides a method for simulating alkanes in a united atom approach. The reason for this is twofold [38]. Firstly, ambiguity related to the hydrogen flexibility in the resulting data can be avoided. Secondly, this also has the advantage of being computationally less expensive in comparison to an all-atom simulation. The initial crystallographic structural information for all the n-alkane molecules was obtained from the Cambridge crystallographic data centre (CCDC) and by the x-ray diffraction work carried out by Boese et al. [5] Each individual simulation cell was built using the VESTA software [41]. The whole system was repeated in such a way that final size of the cell is at least 28 Å along all the directions of the unit cell. Preliminary runs for all the temperatures in the NPT ensemble for 1.2 ns were carried out to estimate the density at the given temperature. Thereafter, these density values were used to carry out the production simulations in the NVT ensemble. The system equilibration was confirmed by ascertaining that the fluctuations in temperature and pressure are below 5 %, and with energies (potential and kinetic) stabilized.
All the NVT simulations were performed in two parts to capture processes occurring at both short (sub picosecond) and long (nanosecond) timescales. At shorter timescale, the simulations were performed for 1.2 ns with trajectories saved for every 0.010 ps, while at longer timescale the simulations were carried out for 10 ns with trajectories written at every 0.5 ps. The initial run time of 0.5 ns and 1 ns for the respective short and long runs is considered as equilibration time. For unwrapping the trajectories the VMD software tool was utilized [42]. Long-range interactions were calculated by a 3D Ewald sum. Three-dimensional periodic boundary conditions were used and following the TraPPE-UA convention 14 Å is taken as cut-off. All simulations employed an integration step of 1 fs. Further simulation details and unit cell dimensions and the related multiplicity values are described in the supplemental and Table S1.

Evolution of structure with temperature
The odd-even behavior is depicted in thermodynamic properties, such as melting point and density, of n-alkanes with increasing chain length (n). To obtain the melting point of each system from the simulations at each temperature, we calculated the orientational order parameter as a function of time using equation 1, where, θ(t) is the angle made by the molecular axis with the reference direction shown in the inset to Figure 1(a) and〈 〉denotes an ensemble average. The temperature at which the order parameter values change from -0.5 or 0.2, i. e. from anti-parallel or partial alignment with respect to the reference direction, to 0, i. e. with completely random alignment, was denoted as the melting point. As an example the calculated order parameter for pentane for different temperature is shown in Figure 1(a). The calculated orientational order parameter as a function of time at various temperatures for all the alkanes is shown in the supplemental in Figure S1. The orientational structure was also probed by calculating the distributions of molecules with their molecular axis aligned at an angle (θ) with respect to the Cartesian directions, as depicted in the inset to Figure 1(a). This distribution, henceforth referred to as the orientational distribution function (ODF), along with the isotropic distribution expected in the absence of orientational order is shown in the supplemental section in Figure S2 for all the molecules. The transition from an ordered state represented by sharp peaks in the ODF (such as for pentane at 120 K in 1(b)), to a disordered state marked by an ODF that closely follows the isotropic distribution, i.e. forms a perfectly random distribution can be seen. This transition also marks the melting point (e.g., inset to 1(b) at 125 K).
The density values were extracted from the simulations performed in NPT ensemble, as described in the simulation details section. In Figure 1(c) the melting point and density as a function of n are shown. For lower n-alkanes (n = 3 -6) clear odd -even behavior in melting point and density is evidenced, while at higher values of n = 7 and 8, this effect is less prominent. Weakening of odd -even behavior with increasing value of n has been reported for other systems, like diamides [9]. However, the deviation for n = 7 from the expected odd-even behavior, as seen in Figure 1(c), may also be influenced by to the parameterization of force field for n = 3, 4, 5 and 8 alkanes, as explained in the discussion section.

Positional order
To further examine structural information as a function of temperature, radial distribution functions (RDF) for CH3 -CH3 pairs belonging to different molecules are plotted in Figure 2. The temperature dependence of the RDF for CH3 -CH2 and CH2-CH2 pairs for all the molecules are shown in the supplemental Figures S3 and S4, respectively. For all the alkanes, the transition from an ordered state represented by the presence of sharp peaks in RDF to a disordered state shown by shallower and broader peaks can be seen. In the ordered state, i. e. below their respective melting points for all the molecules, the slight broadening observed in the sharp RDF peaks can be attributed to the presence of rotational disorder, as well as angular and dihedral vibrations in the system. Increasing the temperature displays a trend towards the disordered state.
In Figure 2(c) the temperature dependent RDF for pentane CH3-CH3 is shown. Interestingly, unlike the shorter alkane systems, the second nearest neighbor peak at 5.14 Å in RDF is considerably sharper than the first nearest neighbor peak. We note the presence of a prominent layering structure in n > 4 systems, as illustrated in the respective insets to Figure 2. It is possible that intra-layer molecules might be rotating in synchronization and not contribute to the peak broadening in RDF, while this may not be the case for inter-layer molecules. In contrast, for n < 4 systems molecules are relatively more aligned along the b-axis and small enough that the layers are comparable with the inter-layer gap. Similar to pentane, a sharp second nearest-neighbor peak in the RDF for hexane ( Figure 2(d)), ascribed to the rotational motion of intra-layer molecules, is observed at 6.47 Å. For hexane a shift in second nearest neighbor peak position and together with a reduction in intensity is observed. This can be correlated with increase in the intralayer molecular spacing with the addition of CH3. Therefore, as a function of increasing n a continuous shift in the second nearest neighbor peak position with concomitant reduction in intensity is observed (Figure 2(e) for n = 7 and 2(f) for n = 8).
To understand the conformational order in alkanes, we calculated the distribution of CH3 -CH3 intramolecular distances. This CH3 -CH3 distance can also be viewed as an effective length of the molecules. In the insets to Figure 3(a) and 3(b) CH3 -CH3 intramolecular distances at (i) in the crystalline state and (ii) at respective melting points is shown, respectively.
A bimodal distribution in the distances is observed for all the alkanes, except for propane (n =3). To probe the odd -even behavior in the bond length distribution below and above the respective melting points we calculated the full width at half maximum (FWHM) of the sharper peaks and distance between the two peaks, as shown in Figure  3(a) and 3(b), respectively. A clear odd -even behavior with no significant variation with temperature, is observed. Interestingly, here heptane (n = 7) appears to follow the odd -even trend, unlike its behavior observed in other calculated properties. This may indicate that translational and rotational dynamic behavior is strongly influenced by intermolecular interactions, and hence density as compared to the molecular shape.

Translational motion
To understand the translational diffusive behavior of the alkanes, the time dependent mean square displacement (MSD) of the centres of mass of the molecules was calculated. Figure 4(a) depicts the temperature dependent MSD for pentane at both short and long (inset to Figure 4(a)) simulated time scales. The time dependence of the MSD for n = 3 -8 as a function of temperature at shorter and longer (inset) timescales is shown in the Figure S5 of supplemental. At shorter timescales, two different time regimes in MSD are observed. Initially, in the first regime, below 0.3 ps, sharp increase in MSD as a function of time is observed. In this regime MSD can be characterized by MSD ∝ t 2 , corresponding to ballistic motion of the molecules. Here molecules move freely without any interaction with the neighbouring molecules and therefore, no significant influence of temperature on MSD in the probed time scale is observed. In the second regime, above the 1 ns time scale, when the system is below its respective melting point a plateau-like trend with no significant variation in MSD with time is observed. This indicates highly restricted motion under the influence of neighboring molecules, also referred to as a 'caging' region [43]. We observe that at these temperatures, molecules are unable to escape the 'cage' formed by the surrounding molecules even at longer times (up to 10 ns). This implies an absence of the long-range diffusion typical of liquids. At higher temperatures, at and above the respective melting points of the molecules, in the second regime (above 1 ns) not only does the MSD exhibits nearly linear time dependence (MSD ∝ t 1 ), but a strong influence of temperature is also seen. The MSD displaying linear behavior as a function of time indicates pure diffusive motion, while sub-or super-diffusive motion would arise from departure from linearity.
where D is the diffusion coefficient and nd (= 3) represents the dimension. The diffusion coefficient as a function of n is shown in Figure 4(c). As a function of n, D exhibits an odd -even behavior both at and above Tm, for Tm + 5 and Tm + 10.

Rotational motion
The temperature variation of orientational dynamics can be investigated by studying the evolution of a unit vector u (shown in the inset to Figure 5(a)) along the position vector of a pseudo-atom comprising the molecule with respect to its centre of mass. The orientational correlation functions (OCF) are calculated using equation 3.

( ) = 〈 ( + ) • ( )〉 (3)
A representative figure depicting the typical OCF behavior with time for pentane at several temperatures is shown in Figure 5(a) for the short and longer run simulations (inset to Figure 5(a)). The time dependent OCF behavior is shown for all the alkanes in Figure S6 of supplemental.
Below their respective melting points, all alkanes show non-zero residual correlations at longer times, indicating a severely constrained rotation. These residual correlations suggest that the corresponding molecules are unable to span the entire orientational space and are restricted to undergo libration-like reorientational motion, a characteristic of crystalline structure. This behavior is in agreement with the high degree of anisotropy observed in the ODF at low temperatures ( Figure S2 of supplemental). Above the respective melting points, most of the correlations decay to zero within 100 ps. This suggests that at these temperatures the molecules span the entire orientational space via rotational motion, a behavior akin to liquids. It is important to note that for pentane and heptane at 125 K and 190 K respectively, we observe an OCF trend of an intermediate state. This behavior is consistent with the RDF and MSD results described earlier.
Further, the shape of the OCF facilitates identification of two time regimes, one below ∼ 1 ps, where the decay of OCF has a relatively steep slope, and another above ∼ 1 ps where the slope has flattened. Such distinct regimes in the OCF have been identified earlier with a 'free-rotor' like motion (first, sub-picosecond regime) and a slower rotational motion where inter-molecular interactions start influencing the rotations (second, longer time regime) [44]. The time scales (τ) are obtained by fitting the OCF in the regime above ∼ 1 ps with an exponential function. The odd -even behavior is most pronounced at temperatures below Tm, particularly at Tm-5.
We obtained the correlation time scales (τ) by fitting the OCF in the regime above ∼ 1 ps with an exponential function. The corresponding variation of τ as a function of n in CnH2n+2 above and below melting temperature (Tm) is depicted in Figure 5(b). Interestingly, the time constant τ displays odd-even behavior which is distinct at temperatures below Tm, particularly at Tm-5. At temperatures Tm and above the odd -even behavior appears less pronounced. This behavior is also consistent with the diffusion coefficient obtained from MSD.

Discussion
In a study reported by Gautam et al. [45] fully atomistic (explicit hydrogen model) simulations were performed for propane. When comparing the MSD between two simulation models, not much difference has been observed, implying that flexibility in C -H bonds and H -C -H angles may not have a significant influence on the oddeven behavior.
We compare the melting points of the alkanes obtained from our simulations with experimental data by Boese et al. [5] and Yang et al. [37]. Overall, good agreement between the experimental and simulation results is obtained. Both the melting point and densities follow similar trends. The melting points of even numbered alkanes are higher than the odd ones while, again as expected, the respective densities for even numbered alkanes are lower in comparison to odd alkanes. The translational behavior shown in Figure 4 is in good agreement with the quasi-elastic neutron scattering experimental study reported by Yang et al. [37]. The odd even effect in translational motion obtained from the MSD (Figure 4) is retained at higher temperatures at and above (Tm). Also noteworthy is the existence of an odd-even alternation of MSD in the ballistic regime (t=0.1 ps). In this regime, the molecules are under free motion and are not affected by the surrounding molecules. An existence of an odd-even effect even for this ballistic motion implies existence of this effect even independent of the nearest neighbor influences and hence the origin of the effect may be the molecular shape itself as suggested by [9]. Similar to this, the existence of odd -even alternation in the molecular size attributes shown in Figure 3 also suggests the importance of molecular shape/conformation in the origin of odd -even effect. This is consistent with [5,9,10] that suggest the molecular geometry -trapezoid vs parallelogram as the origin of this effect.
For the rotational dynamics obtained from the OCF ( Figure 5) the odd even behavior is only found at temperatures below Tm. Therefore, the simulations exhibit the oddeven effect in translational motion even in liquid state, while not in the rotational motion. This implies that odd-even behavior is more prominent at longer time scales.
In further comparison of the structural parameters such as the RDF and ODF, no signature of structural ordering (long or short ranged) is observed above Tm. However, the presence of odd-even behavior in translational properties even above Tm indicates that conformational order may not be the only driving force. The intermolecular interactions also contribute to the effect via density. To conclude, both the molecular conformation / geometry, and the intermolecular interactions are important and their interplay gives rise to the odd -even effect.

Conclusions
To obtain a molecular level understanding of the origin of odd-even alteration in the properties of n-alkanes, we carried out MD simulations of alkanes with a united atom model. We investigated the structure and dynamics of both positional and orientational order in six n-alkanes with n = 3 (propane) to 8 (octane). For each alkane, MD simulations were carried out at several temperatures below and above the melting point, thereby detailing the evolution of the properties and hence their odd-even alternation in temperature across the melting point. Detailed calculations of several properties from the same simulation data helps us see correlations between the oddeven alterations in different properties and thus helps understand the origin of the odd-even effect. It can be concluded from the results reported here that both the molecular shape and density of alkanes play a significant role in the origin of the oddeven effect.

Structural Details
The orientation order parameter, s(t), is a calculation characterizing the alignment of a system, suitable for intrinsically anisotropic materials, such as liquid crystals. S is given by the equation below, where (t) denote the angle made by the molecular axis with respect to the reference direction and  The radial distribution function (RDF) describes the probability distribution of molecules at a given distance, hence providing the local structural information. The temperature dependence of RDF for CH3 -CH2 and CH2 -CH2 (CH3 -CH3 is shown in the main manuscript) pairs for all molecules are shown in Figure S3 and S4, respectively. Similar for all alkanes, we observe a temperature dependent transition from an ordered crystalline state represented peaks. The slight broadening of RDF peaks in the ordered state, i. e. below their respective melting points for all the molecules, can be attributed to the presence of rotational disorder, angular and dihederal vibrations in the system.

Dynamical Details
The mean squared displacement (MSD) is a measure of the deviation of a molecule position with respect to its starting position over time. Therefore, the time dependent MSD from the center of mass of the molecules can be used to understand the translational diffusive behavior of alkanes. The overall MSD for n = 3 -8 as a function of temperature at shorter and longer (inset) timescales is shown in the Figure S5. To investigate rotational motion, the orientational correlation function (OCF) as a function of temperature was used. For all the alkane the time dependent first order (l = 1) OCF behavior is shown in Figure S6.