Dielectric and Transport Properties of Acetonitrile at Varying Temperatures: a Molecular Dynamics Study

Abstract

Use of acetonitrile in electrolytes promotes better operation of supercapacitors. Recent efforts show that electrolytes containing acetonitrile can also function in a wide range of operating temperatures. Therefore, this paper addresses the dielectric relaxation processes, structure and dynamic properties of the bulk acetonitrile at various temperatures. Systems of acetonitrile were modeled using canonical ensemble and simulated by employing Molecular Dynamics method. Results show that interactions among the molecules were correlated within a cut-off radius while parallel and anti-parallel arrangements are observed beyond this radius at relatively high and low temperatures respectively. Furthermore, effects of C-C-N and C-H bending modes were greatly appreciated on the power spectral density of time rate change of dipole-dipole correlations whereas frequency shifts were observed on all modes at the lowest temperature under consideration. Linear variations with temperature were depicted for reorientation times and self-diffusion coefficients. Shear viscosity was also computed with a good accuracy within a certain range of the temperature as well.

Introduction

Acetonitrile (ACN) is an aprotic and polar solvent which is widely used in chemical industry. Supercapacitors which are also known as Electrochemical Double Layer Capacitors (EDLCs) are one of the technological applications in which the ACN is employed as a solvent in electrolytes due to its low viscosity, good conductivity and permittivity despite its very low flash point (6 ℃), high flammability and toxicity.1 Boiling (82 ℃) and melting (−45 ℃) points of ACN have been seemed to be a drawback for a wide range of operational temperatures of supercapacitors. However, there is an increasing scientific interest to enhance the working temperature span. Improvements on the EDLCs’ characteristics can be classified in twofold, namely on electrode materials and on content of electrolytes. System of exohedral nanostructured carbon (nanotubes and onions), activated carbon cloths or graphene nanosheets as an electrode and ACN accompanying ionic liquids allows operation in a wide range of temperatures of EDLCs.2-5 Furthermore, graphene as an electrode material coupled with an electrolyte composed of ACN and an ionic liquid increases electrochemical window that can enhance the energy density and specific power density.6 Use of Room Temperature Ionic Liquids (RTILs) incorporated in molecular solvents, such as ACN, γ-butyrolactone favors thermally stabilized and less toxic electrolytes with a lower cold crystallization temperature, −53 ℃.7,8 Moreover, a promised electrolyte system involving ACN and a low melting co-solvent reduces the operational temperature to −75 ℃ albeit with capacitance reduced to about half of its room temperature value.9 A recent study shows that large residual conductivities in ACN makes low temperature operation of EDLCs possible.10

There are limited number of studies dealing with modeling and simulation of bulk ACN. In an earlier study performed by Bӧhm et al., a six-site model was introduced and it was shown to be reproduced some static and dynamic properties of the real system of ACN.11 Another theoretical accomplishment was based upon RISM calculations with six-site hard sphere model and predicted X-ray structure factors compared to those of experiments very well.12 Jorgensen and Briggs performed Monte Carlo simulations using a three-site model at 25 ℃ and 75 ℃.13 The thermodynamic results of them agreed well with experimental findings. They especially pointed out the importance of electrostatic interactions on the determination of the structure of the liquid ACN. Most of the recent studies about the bulk ACN are directly related to calibration of force field parameters either by three-site or by six-site models. Gee and Van Gunsteren have iterated the model parameters that have originally based upon a model introduced by Edwards et al. until a good match is observed between the predictions of iterated model and experiments.14 Another study by Grabuleda et al. reproduced density, heat of vaporization and isothermal compressibility successfully.15 A flexible six-site model was employed in MD simulations of liquid ACN in that study. A further effort on the calculations of thermodynamic properties constituting with the experiments was reported in the study of Nikitin and Lyubartsev.16 The other calculated quantity, the dielectric constant, was predicted fairly below the experimental one. In a recent study, Albertí et al. proposed a portable model potential to account self-interactions and interactions with ions of the ACN. It anticipates the self-diffusion coefficient and the density well in addition to discussed cation-ACN structure and coordination.17 Edwards et al. developed a three-site model based on the six-site model of Bӧhm et al., and then performed MD simulation to obtain frequency dependent dielectric permittivity, Kirkwood g factor which was close to unity and the other quantities as well.18 Moreover, a study focusing on the calculation of the shear viscosity and the dielectric constant was considered three threesite models; namely the model reported in the study of Edwards et al. and other two variants. The results obtained from the MD simulations of the original three-site model agreed well with those of experiments.19 A further three-site model which was derived by Guàrdia et al. was compared with the available three-site models. In this comprehensive work, it was shown that an ab initio potential improves the thermodynamic and dielectric properties, and discussed the effect of molecular polarizability on the interaction potentials.20

In the highlight of some recent literature about supercapacitors, one can deduce that the role of ACN in electrolytes is still important and so the operation of EDLCs at extended temperatures is already in progress. Although this brief summary does not give a detailed picture of scientific and technological interests, it shows that analyses of physical behavior of ACN at extended temperatures are required. To the author’s knowledge, there is no any study concentrating on dielectric relaxation phenomena at relatively high and low temperatures. These are the reasons why this study focuses on the dielectric and transport behavior of the ACN at extended temperatures.

The paper is organized as follows. After briefly representing force field parameters in the next section, details of molecular dynamics simulations are presented, and then the transport and dielectric properties of the liquid acetonitrile are given in each separate section followed by conclusions.

 

Modeling and Simulation

A six-site model that consists of three methyl hydrogen atoms (H), a methyl carbon atom (C1), a nitrile carbon atom (C2), and a nitrogen atom (N) is employed for the ACN. Total force acting on an atom is calculated as a sum of inter and intra-molecular forces. Force field parameters of the ACN were adopted from DREIDING/A force field21 which was subjected to some minor changes in this study. The force field was reported to be worked well at various temperatures.39-42 Intra-molecular interactions are the sum of bond stretching (EB), bond-angle bending (EA, three body), dihedral angle torsion (ET, four body) and inversion (EI, four body) interactions. Dihedral angle torsion and inversion terms were not taken into account according to DREIDING/ A. In intra-molecular force field parameters, stretching energy between the C1-H and C1-N were set to 665 and 2500 kcal/ mol Å2 respectively. Furthermore, the angle bending energy for H-C1-H and C1-C2-N are decreased to 64 and 49 kcal/ mol rad2 respectively.

Both pair-wise dispersion and electrostatic interactions were considered in the computation of intermolecular interactions such that

where rαiβj is the distance between the atom i in molecule α and the atom j in molecule β. σ and ε stand for the equilibrium atomic separation and the potential well respectively. q and C are the charge of each atom in the molecule and Coulomb’s constant respectively. Lorentz-Berthelot mixing rule was applied in the evaluation of cross interactions between the atoms. Furthermore, the point charge distribution of the acetonitrile has been modeled by employing optimized molecular charges reported in the study of Cabaleiro-Lago and Rios.22 The charges used in simulations are −0.552, 0.475, 0.190, −0.493 for C1, C2, H and N respectively. The calculated dipole moment, based on these point charges and intra-molecular geometry given in Table 1, equals 3.88 D that is very close to experimental dipole moment, 3.91 D (Ref. 38). In the report of Cabaleiro-Lago et al.,22 optimized point charges have been obtained exclusively from ab initio calculations using the basis set MP2/6- 311+G*. The obtained dipole moment, 3.95 D, agreed well with the mentioned experimental dipole moment. Among the studies reporting these sets of point charges, Grabuleda et al.15 derived a new set in comparison with the pioneering work of Cabaleiro-Lago et al.22 In that study, the point charges, −0.479, 0.481, 0.177 and −0.532 standing for charges of atoms C1, C2, H and N lead to a higher dipole moment, 4.18 D, with respect to experimental dipole moment. Another set of point charges introduced by Nikitin and Lyubartsev16 was based on ab initio calculations using the basis set MP2/6- 311++G(3df,3p). This set of charges gives a dipole moment that is very close to experimental value, and it is very similar to results of Ref. 22. Due to these reasons, this study uses the charge distribution given in the report of Cabaleiro-Lago et al.22

For modeling intermolecular force field parameters, the energy parameter of the Nitrogen was adopted from the study of Guàrdia et al. and set to 0.1498 kcal/mol.20 Energy parameter of H and separation distance of N were based on DREIDING and set to 0.0152 kcal/mol and 3.292 Å respectively.21 Separation distance for H was taken account as a free parameter for better agreement between the computed and experimental viscosities. Set of Lennard Jones parameters used in the model are summarized in Table 1.

Table 1.Intra and intermolecular force field parameters used in a six-site representation of the ACN

Simulation cell is subjected to periodic boundary condition in each spatial direction. Simulations were performed in canonical ensemble. 1500 number of molecules matches the experimental density, 777 kg/m3, in a box of 50.887 × 50.887 × 50.887 Å3. All simulations have been carried out on Large Scale Atomic/Molecular Massively Parallel Simulator (LAMMPS).23 After initialization, the Conjugate Gradient Method was employed for minimization with 1 × 105 maximum number of allowable iterations and with convergence criteria of 1 × 10−8 and 1 × 10−10 for energy and force evaluations respectively. Long range electrostatic interactions were computed by the Particle-Particle, Particle-Mesh (PPPM) method24 with a cut-off of 13 Å and 1 × 10−5 relative error in the force evaluations. The Lennard-Jones cut-off was also kept 13 Å. The dielectric constant was set to unity and neighbor lists were updated in every 10 time steps. Canonical ensemble was employed in both equilibration and production phases of all simulations. Time step was 0.5 fs. All systems were warmed up to the desired temperatures; 233, 253, 288, 298, 308, 318 K, from 10 K by employing Noose-Hoover thermostat over 50 ps period of time (with 25 fs coupling time) followed by equilibration processes over 300 ps period of time. In the production stage, simulation times were 2.5 ns. Components of shear stresses per atom were sampled in every time step. Positions and velocities of atoms were dumped in every 10 time steps. Sampling periods of time were 250 and 500 ps for velocities and positions respectively. 1-4 interactions (electrostatic and Lennard-Jones) were also taken account during the simulations. Three sets of simulations for each case were carried out for production of data.

 

Results and Discussions

In this section, the structure, the dielectric and transport properties of the ACN will be exhibited systematically.

Structure. Radial Distribution Functions (RDF) between each pair of atoms excluding H atom were illustrated in Figure 1. An important and close coordination was observed between the atoms C1 and N with amplitude 2.24 at about 3.7 Å. The second coordination shell for this pair interaction was taken place at 7.5 Å. Probability of finding C1 atom in the vicinity of the N atom was due to interactions of positively charged H atoms with the negatively charged N atom. Such coordination is observed in both three and sixsite models reported in other references20,16,13,14 with almost the same amplitude and location. Further examination of Figure 1(a) reveals that the least interaction between C1 and N occurred at radius, 5.4 Å, where other interactions, such as C1-C2 and C1-C1, are still high enough. A monotonic decrease of the RDF between C1-C1 from the first maximum, 1.37, to the other one, 1.20, illustrates a close interaction shells located on 4.7 and 5.95 Å. The second peak observed in this interaction is not significant in three-site models.13,14,20 In contrast to the other interactions, N-N interactions take place relatively are taken place in a long distance compared with those of C1-N. Moreover, N-N interactions form a second shell that is very close to the first one. Locations of the peaks of the first and second coordination are 4.45 and 5.65 Å respectively. This form of N-N interaction is noticeable in six-site models in contrast to three-site ones. However, the second maximum of this interaction has relatively high value in the report of Nikitin and Lyubartsev.16 On the other hand, an important C2-C2 coordination occurs relatively long radius, about 4.95 Å, compared to the other interactions. This distance of interaction with almost the same amplitudes was also reported in other studies. Inclusion of H atoms leads to such a relatively long radius of interaction between the C2 and itself. Furthermore, sum of the bond length between the C1-H and intermolecular equilibrium distance between the H and N almost equals to the location of the first coordination between the C1 and N, 3.7 Å. This is due to interaction of N atom in one molecule with both H and C2 atoms (methyl group) in another molecule, in a sequential arrangement.

Figure 1.RDFs of heavy atoms of the ACN at 298 K. Panel (a) shows C1-C1, C1-C2 and C1-N interactions. In panel (b) C2-C, C2-N and N-N interactions were illustrated.

The vibrational spectroscopy of the ACN is another point discussed in this study. There are four non-degenerate vibrations of a1 symmetry and four degenerate vibrations of esymmetry.25 In a1 symmetric group, a symmetric C-H stretch, a C≡N stretch, a symmetric C-H bend and a C-C stretch are designated as ν1, ν2, ν3 and ν4 respectively. Four doubly degenerate e-symmetry vibrations, ν5, ν6, ν7 and ν8 stand for a C-H stretch, a C-H bend, a methyl rock and a C−C≡N bending respectively. The calculated spectra is illustrated in Figure 2 and compared with the experimental one in Table 2. The maximum and minimum relative differences, 0.05016 and 0.00266, between the calculated and experimental frequencies are observed for C-C stretching, ν4, and C-N stretching, ν2 modes respectively. The calculated mode frequencies agree well with the experimental one reported in the study of Deak et al.25

Table 2.aFor temperature 233, 298 and 318 K. bFor 233 K. cFor 298 and 318 K

Figure 2.Vibration spectra of acetonitrile at 233, 298 and 318 K. No significant change is observed with respect to the temperature.

Transport Properties. Due to significance of the transport properties at extended temperature, both shear viscosity and self-diffusion coefficient was investigated.

The shear viscosity was obtained from Green-Kubo relation given as

where τab corresponds to molecular stress tensor and calculated as the following,

In Eq. (3), rbij = rbi – rbj , (bth component of the vectorial distance between atom i and j) Faij = Fai – Faj (ath component of the sum of forces acting on atom i and atom j) and a, b = 1, 2, 3 stand for components of Cartesian coordinate system. Autocorrelation function of the stress tensor was computed precisely by considering a suggestion of Nevins and Spera.26 Time interval between each start of successive time windows was chosen sufficiently small, and set to 7-10 fs, to assure the accuracy of computations. Correlations were computed over 2.5 ns, and the correlation time was kept 10 ps because of rapid decay of it. For integration of correlation functions, the composite Simpson rule is employed. Time variation of the shear viscosity in each temperature is shown in Figure 3 and the averages are tabulated in Table 3.

Figure 3.Variation of shear viscosity with correlation time. The correlation function of molecular stress tensor rapidly becomes asymptotes to the mean value it.

Relative differences between the theoretical and experimental predictions for each temperature from 233 K to 318 K are 0.5006, 0.3310, 0.1210, 0.0029, 0.0283 and 0.1513. The predictions of this study are still acceptable for given temperatures excluding those at 233, 253 and 318 K. At 298 and 308 K, the calculated shear viscosity agrees well with experimental values. The shear viscosity at 298 K predicted by other theoretical efforts was 2.9 × 10−4 (Ref. 14), 3.5 × 10−4 (Ref. 19) in shear viscosity unit.

Furthermore, self-diffusion coefficients for each temperature are calculated from integration of velocity auto correlation function,

Integrations were performed again by the Composite Simpson rule. Self-diffusion coefficients obtained from simulations vary linearly with respect to the temperature as illustrated in Figure 4. Calculated self-diffusion coefficient at 298 K, 3.98 × 10−9 m2/s is very close the experimental finding, 4.04 × 10−9 m2/s. Self-diffusion coefficient for each temperature is also tabulated in Table 3.

Figure 4.Linear variation of the self-diffusion coefficient with respect to the given temperatures.

Table 3.Calculated shear viscosity and self-diffusion coefficient for varying temperatures

Many other studies predicting the self-diffusion coefficient from simulations of either three or six site models are all below the experimental value. These theoretical estimations span from 2.6 × 10−9 (Ref. 20) to 3.55 × 10−9 m2/s (Ref. 19).

Dielectric Properties. In this section, dielectric properties of the ACN will be discussed in detail. For a nonconducting, isotropic and polarizable system, the frequency dependent dielectric constant is calculated as a function of dielectric susceptibility and given as

where ε(ω) is a complex function, ε(ω) = ε′(ω) − iε″(ω) and its real and imaginary parts correspond to the dielectric permittivity and the dielectric loss respectively. The dielectric susceptibility is obtained from the equation,

Here, β = 1/(kBT) , and M and kB stand for the total dipole moment per unit volume of the system and Boltzmann constant respectively. In the classical limit, the tilde denoting Kubo transform is dropped and thus Eq. (5) is rewritten as

If one considers Hamilton’s equations are invariant to time origins of a correlation function, then Eq. (6) is represented as a time rate change of a dipole-dipole correlation function and reorganized as . Thus; the expression for the complex dielectric permittivity is read as

In this study, two different ways on the calculation of frequency dependent dielectric permittivity are followed. In the first approximation, systems were treated as dilute ones. Therefore, the cross terms are negligible in this case, and the correlation function converts to a correlation function of individual molecules that is given as;

where N and μ stand for number of molecules in the system and dipole moment of a molecule respectively. For further details, readers are referred to McQuarrie.31

Calculated dielectric permittivity for each temperature based on correlated system of dipoles is shown in Figure 5. At first glance, similar variations of the dielectric permittivity and the loss are observed overall. However, some differences between those at different temperatures arise at some frequencies. The first observation was that the dielectric loss increases with decrease in temperatures. The most prominent increase is observed at about 380 cm−1 (see inset a in Figure 5) where a sharp decrease of the dielectric permittivity takes place. At this frequency, peak values of dielectric losses, 15.6, 12.5 and 11.89, vary linearly with the temperature as well.

Second, a decrease in all mode frequencies appeared in Figure 5 were observed compared to vibration spectra of the ACN. The magnitude of such decrease, frequency shift, at high temperatures is not so significant with respect to those of colder ACN. The first and easily observable change was on the C-C stretching mode frequency that decreased to 716 cm−1. Another obvious change was at the frequency of ν6 mode. Although this frequency is at about 1465 cm−1 for higher temperatures, it is calculated as 1263 cm−1 at 233 K. Such a frequency shift tends to increase with increase in the frequency of the modes.

Figure 5.Frequency dependent dielectric constant and loss at 233, 298 and 318 K. Inset a shows peaks of each the dielectric loss variations. In inset b, a contradictory change of both the dielectric permittivity and the loss with respect to varying temperature.

Figure 6.The dielectric constant and loss as function of frequency. Calculated quantities are based on the correlated system of dipoles approximation.

In the second approach, the system is assumed to be a system that composes of molecules interacting with each other, and these interactions are correlated. In this case, if one calculates the total dipole moment of the system as , then the cross terms of the form 〈μα(t) ⋅ μβ(t)〉are taken into account in the calculation of the dipole-dipole autocorrelation function that is given as

The calculated frequency dependent dielectric permittivity based on correlated interactions of dipoles for given temperatures is shown in Figure 6.

The variation and shape of the dielectric permittivity curves are very similar to those of the former approximation but the magnitude of the static dielectric constants is smaller. All static dielectric constants obtained from both assumptions, namely correlated and uncorrelated systems of dipoles, are tabulated and compared with experimental data in Table 4. The results obtained from Eqns. (7) and (8) are higher than the experimental data whereas those of Eqns. (7) and (9) are substantially smaller. The minimum relative differences between the current prediction and experimental one were 0.027, 0.057, 0.040 and 0.051 at 288, 298, 308 and 318 K respectively.

Table 4.Dielectric constants obtained from various approaches. The results are also compared with experimental data available in the literature

Figure 7.Time variation of the static dielectric constant.

The other approach is fluctuation formula which evaluates the static dielectric constant as frequency goes to zero. In this case, it can be shown that, in conjunction with Eqns. (7) and (9),

where M is the total dipole moment of the system per unit volume. Since the dielectric constant converges the solution in a long period of time, the calculations were conducted over 400 ps that seems to be enough for ensemble averaging. Time variation of the static dielectric constant is illustrated in Figure 7. The converged values of the dielectric constant are also summarized in Table 4. As shown in Table 4, all results obtained from the fluctuation formula, Eqn. 10, were smaller than the experimental data but agree well with the results of correlated system of dipoles. Dielectric constants calculated by use of several theoretical approximations and reported in various studies were given in Table 4. Gee and Van Gunsteren14 have been successfully anticipated the value of dielectric constant as 35.3 that is very close to the experimental one.

Power spectral density of the time rate change of the dipole-dipole correlation function given in Eqn. 8 is also analyzed to clear out the effect of mode frequencies on the relaxation phenomena. Power spectral densities of the systems at given temperatures were drawn in Figure 8. ν8, C–C≡N bending mode frequencies did not change with temperature at all as illustrated in Figure 8. Although low mode frequencies at each temperature are almost the same, some differences among them arise for frequencies higher than 500 cm−1. All mode frequencies at the lowest temperature under consideration are lower than corresponding mode frequencies at higher temperatures. Some frequency shifts were observed for smaller values of temperatures. The difference between the ν4 modes at 233 K and 298-318 K is at about 111 cm−1. The shift of mode frequencies also increases with temperature. Further shift, 202 cm−1, takes place at C-H stretching, ν6 mode. Such a slowed down relaxation modes have not been observed in the vibration spectra of ACN given in Figure 2.

Figure 8.Power spectral density of time rate change of dipoledipole correlation based on the Eqn. (8). Effect of modes are also illustrated accordingly.

Analyses of relaxation and reorientation times were the other subject of the current study. The molecular reorientation dynamics is analyzed through time correlation function of Legendre polynomials of the first and second kind that are given as where Pl (l = 1, 2) stands for the first and second type Legendre polynomial, and u(t) is the unit vector of the molecular dipole moment defined as u(t) = μ(t)/||μ(t)||2. Correlation time was 20 ps and all correlation functions were computed over 500 ps. Time correlation functions of molecular orientations were shown in Figures 9 and 10.

Figure 9.First-rank orientation time correlation functions Cl(t) for ACN at several temperatures. The inset shows a stretched exponential fit to the computed data at 233 K.

Figure 10.Second-rank orientation time correlation functions C2(t) for ACN at several temperatures. A stretched exponential function is fitted to the computed reorientation correlation function at 233 K.

Both reorientation times deviate from exponential decay pattern and can be described by Kohlrausch-Williams-Watts (KWW) law also known as the stretched exponential law which is read as where 0 ＜ γ ≤ 1 . Mean value of reorientation time is computed from

Analysis shows that the relaxation behavior of ACN does not obey Debye theory. The relaxation best described by a bi-exponential curve such that,

and a mean relaxation time is obtained from

Mean reorientation times computed from Eqn. (11) and relaxation (based on Eqn. 8) times are listed in Tables 5. The predicted reorientation time at 298 K, 〈τ1〉, is below experimental values, 3.28 ps and 3.68 ps, as reported in Ref. 35 and 36 respectively. The relative difference between 〈τ2〉 values of the current study and experimental one30 is −0.225.

An experimental relaxation time at 298 K has been reported33 as 3.4 ps that is higher than the value tabulated in Table 5. Variation of reorientation times are also drawn in Figure 11. A linear correlation was observed between the temperature and reorientation times.

Table 5.Calculated reorientation and relaxation times

Figure 11.Linear variations of reorientation times with respect to the temperature.

Simulation results for the reorientation time of first and second ranks have been reported in both Ref. 20 and 14. While the former study predicts 〈τ1〉 and 〈τ2〉 as 3.8 and 1.6 ps the predictions of latter one were 1.93 and 0.80 ps that are considerably small compared to the experimental value. Among the theoretical studies, Guardia et al.20 successfully predicted the relaxation time as 3.3 ps. The other values were 4.3 (Ref. 18) and 2.2 ps (Ref. 14).

Although relaxation behavior of the ACN did not decay exponentially, Debye analysis is still significant in terms of examining frequency response of the dielectric relaxation. Debye theory of dielectric relaxation suggests an exponential decay of dipole-dipole autocorrelation function by assuming that ε(∞) goes to unity, and it can be written as

Thus, calculated frequency dependent relaxation times, based on the Debye theory, at 233, 298 and 318 K were shown in Figure 12.

As seen from the Figure 12, degenerate e-symmetric C–C≡N bending mode leads a substantial decrease in the relaxation time. Furthermore, a significant decrease was observed at e-symmetric C-H angle bending mode at about 1345 cm−1.

Orientation of ACN molecules within a sphere gives valuable information about the intermolecular interactions. For this purpose, The Kirkwood correlation factor was calculated as a function of varying radii of sphere,

where z is the number of neighbor molecules enclosing a reference molecule at the center of sphere. If there is a difference between orientations of the reference molecule and neighboring molecule, the correlation factor should not equal to zero. It is clear that, g greater than unity refers to molecules directing themselves parallel to the reference molecule, and vice versa. In this respect, the effect of some temperatures on the Kirkwood correlation factor was shown in Figure 13.

Figure 12.Frequency dependent Debye relaxation times. Effects of vibration modes on the relaxation phenomena are clearly visible in this graph.

Figure 13.Kirkwood correlation factor as a function of the radius of spherical shells, g(r).

As seen from Figure 13, there is an anti-parallel organization just in the vicinity of the reference molecule starting from r = 3.2 Å. Preferred orientation of the molecules turns into parallel arrangement for radius greater than r = 5.2 Å until the role of the arrangement changes in a certain radius due to sinusoidal behavior of g(r). The same continuous change of the g(r) with similar amplitudes is also reported by Guàrdia et al.20 After certain radii, the correlation factors start to fluctuate and they cannot be represented as continuous functions, and then they converge certain values, and become invariant to the radius of spherical shells.

Converged values of the correlation factor for finite systems considered in this study were 1.408, 1.47 and −2.475 at 318, 253 and 233 K respectively. The predicted correlation factors for given finite systems at 253 K and 318 K are very close the value, 1.48, at 298 K as predicted by Edwards et al.18 On the other hand, such a behavior was not observed for systems at 233 K. For the limiting values of the radius of spherical shells in this system, molecules tend to direct themselves anti-parallel to the reference molecule. The mechanism lying behind the anti-parallel arrangement just above liquid-β solid transition temperature is a decrease in the volume of ACN as indicated by pressure drop in the system. A key role of the volume contraction as a result of phase transitions in molecular arrangement is a change of dominant molecular forces affecting orientations and librational motions.37 Magnitudes of repulsive and attraction forces directly affect the nature of the repulsive and attractive interactions. Another factor affecting the nature of the intermolecular behavior is the electrostatic interactions. To elucidate the reason of the anti-parallel arrangement of the molecular dipoles, one must be consider both electrostatic and Van der Waals interactions. Repulsive pair interactions between N atoms try to direct one N atom apart from the other one in order to construct a parallel arrangement. On the other hand, it is very clear that attractive forces between the N-N atoms must be the main reason associated with the anti-parallel organization of molecules. The attractive force did not cause considerable change in the frequency of the C≡N bond stretching mode as shown in Figure 2. Although, such change is not apparent in the vibration spectra, it was pretty apparent in both complex relative permittivity and also power spectral density graphs illustrated Figures 5, 6 and 8 respectively. In other words, increasing value of intermolecular attractions between N atoms strongly drives the mechanism of anti-parallel organization of molecules and so result in a different fashion of relaxation process just above the liquid-β solid transition temperature compared to higher temperatures.

 

Conclusions

Recent advancements in building supercapacitors make use of acetonitrile not only at room temperature but also extended temperatures. Its low viscosity contributes to better operation of room temperature ionic liquids at low temperatures. In this study, structure, dielectric and transport properties of acetonitrile were highlighted using molecular dynamics method at several temperatures. The number of the particles and the size of control volumes were kept constant during the simulations.

The length of the first coordination between the methyl group carbon and nitrogen was the first important finding of the study. This length is almost equal to sum of the bond length between the C1 and H and intermolecular equilibrium distance between the H and N. This is due to the tendency of molecules directing themselves in a sequential order as expected. The vibrational spectra of acetonitrile are also examined with a good accuracy compared to the experimental findings. In some modes, a frequency shift is observed at the lowest temperature under the consideration. The most probable reason for these changes is to be volume contraction of the system just above liquid−solid phase transition temperature. Such frequency shifts is due to the change in the magnitudes of repulsive and attractive forces resulting from volume contraction.

Transport and dynamic properties of acetonitrile were also successfully anticipated. The predicted viscosity at moderate temperatures were really good agreement with the experimental data. Furthermore, self-diffusion coefficient varies linearly with the temperature. The calculated self-diffusion coefficient at 298 K is very close to experimental datum available in the literature.

The other focus of the study was dielectric relaxation phenomena of the acetonitrile. For this purpose, mainly two ways of calculating complex dielectric permittivity were followed. In the first one, molecules are assumed to be independent from each other and interactions among them are uncorrelated. The predicted dielectric constant by this approach agrees well with the experimental data. On the other hand, the dielectric constant calculated from the correlated interactions approximation is below the experimental findings but very close the other prediction of this study resulting from fluctuation formula. In contrast to results obtained from the correlated molecules assumption, computed dielectric constants by uncorrelated molecules assumption agree well with the experimental data. In order to highlight the effect of mode frequencies on the dielectric relaxation phenomena, power spectral density of the time rate change of dipole-dipole autocorrelation functions was also analyzed. Effects of modes on the relaxation were identified at some frequencies that change with low values of the temperature. The frequencies of all modes at 298 and 318 are almost the same. However, at 233 K, the mode frequencies shift considerably as the frequency of the mode increases. The magnitude of the phase shift was not the same as that of vibrational spectra as discussed.

Another main finding of this study is arrangement of the molecules by probing the temperature. For this purpose, the Kirkwood correlation factor was computed at 233, 253 and 318K, and found to be good agreement with the other studies. It should be noted that at higher temperature, acetonitrile molecules prefer to direct themselves in parallel fashion. As the temperature approaches to phase transition temperature, such preferred direction turns out to be antiparallel arrangement. One possible reason for the change of the directions of the dipoles seems to be attractive forces coming to effect due to contraction in the volume at about phase transition temperatures. The nitrogen atom is prone to be the strongest attraction among the other atom types due to its relatively high potential energy well. The key role of nitrogen atom manages the direction of the dipoles.