Simulation Investigation of Bulk and Surface Properties of Liquid Benzonitrile: Ring Stacking-Assessment and Deconvolution.

The content and the molecular dynamics (MD) simulation analysis here are inspired by our recent ab initio calculation on benzonitrile (BZN), whereas the present results are to expand and develop macroscopic documentation involving data verification. MD simulations of the bulk liquid BZN in the range of 293-323 K unravel the hydrogen bond (-C≡N···H) formation with strength in the order of ortho-H ≫ meta-H ∼> para-H. The possibility for ortho-Hs to get involved in the formation of two bonds simultaneously confirms each having σ- and π-bonding features. Accordingly, we used vast efforts for structural analysis particularly based on the deconvolution of the corresponding complex correlation functions. Specific angle-dependent correlation functions led to the recognition of the molecular stacking with a strict anti-parallel orientation. The in-plane dimer and trimer also take part in the structural recognition. A singularity, found in the trend of the simulated temperature-dependent viscosity and diffusion coefficient of liquid BZN, is centered at about 313 K and quite fascinatingly emulates the reported experiment viscosity. An interplay between a small change in the trend of density and a large change in the corresponding viscosity is a key factor in supporting the singularity. Deconvolution of the simulation results allows attributing the singularity to structural alteration involving H-bonding of different types and extent. Approaching the range of 308-313 K, an alteration between hydrogen bond formation involving mostly ortho-Hs and mixed ortho-Hs + meta-H is possible and supports the singularity.

Introduction

As the simplest aromatic nitrile, benzonitrile (C6H5CN, cyanobenzene) provides a good bridging point between aliphatic nitriles, already extensively researched, and more complex, substituted species that have greater industrial relevance in several sectors of the chemical manufacturing industry such as pharmaceuticals and agrochemicals. Beside a wide distribution of natural products,[1] benzonitrile as a simplest nitrogen-bearing aromatic molecules has been discovered recently in the interstellar medium.[2,3] Benzonitrile and its derivatives are comprising an important class of organic compounds in the chemical industry. This is a stable compound for pyrolysis,[4] an essential process in the synthesis of benzoic acid, benzylamine, benzamide, pesticides, and fatty amines, and dyes could be precursors for the synthesis of benzylamine, which can be used as a pharmaceutical intermediate.[5,6]

Benzonitrile is generally a good solvent for the synthesis of organic, anhydrous inorganic, and organometallic compounds and is generally similar in its behavior as a solvent to acetonitrile (CH3CN).[7] However, benzonitrile differs from acetonitrile in two interesting and significant ways: it has delocalized π-electrons in the phenyl ring constantly interacting with the −C≡N group, and it has a large phenyl group as compared to the methyl group, creating a steric problem.[8] Furthermore, the benzonitrile unlike acetonitrile lacks a hydrogen atom in the alpha position. While acetonitrile and benzonitrile have the same functional group −C≡N, the differences must be sought out in the contribution of the delocalized electrons of the phenyl ring in benzonitrile. More specifically, the cyano group in benzonitriles represents one of the most prevalent functional groups and can feasibly transfer to other functional moieties, such as aldehydes, carboxylic acids, amidines, amides, and amines.[1] The cyano compounds show larger discrepancies relative to the analogous nitro compounds. This is an observable difference between benzonitrile and nitrobenzene even though the dipole moments are similar and the two compounds have equal polarizabilities. This may be explained in terms of a more localized charge in benzonitrile which results in a strongly oriented dipole-moment vector.[9] Benzonitrile is an interesting molecule because of its ability to form a variety of hydrogen-bonded structures, having three different sites through which may form a hydrogen bond: (1) the π-electrons of the CN triple bond, (2) nonbonding electrons of the nitrogen atom of the CN group, and (3) the π-electrons of the phenyl ring.[10] The features of these interactions have been characterized recently[11] by studying the specification of the benzonitrile cluster formation of different sizes by quantum chemical calculations.

Lots of (experimental, theoretical, and computational) research works have shown interest in benzonitrile and its derivatives, beside other recent applications such as co-solvent due to having large molecular dipole moment (4.01[12] and 4.18 D[13]) and widespread strong miscibility.[14−16] Because of the three different binding sites, benzonitrile has shown interesting adsorption features on metal and non-metal surfaces. So far, various studies have been performed on the adsorption of benzonitrile on the surface of Ag, Ni, Pd, Pt, and Au metals, and Fe-doped carbon nanostructure surfaces.[17−22] Also, the structure of various benzonitrile clusters was investigated through different theoretical[11] and experimental methods such as laser-induced fluorescence spectra in a free jet.[23] Conclusively, benzonitrile must be characterized for the molecular structure in the bulk media, the possible molecular orientation at the liquid interface and surface energy, and the free energy of solvation with other solute moieties where it is considered as the pure or co-solvent.

Until now, the literature reports indicate that the liquid benzonitrile has not been studied widely for thermophysical and structural properties based on microscopic features tractable by MD simulation. The limited number of cases of considering possible dimer and cluster formation with different possible intermolecular structures and ring stacking has been reported to account for observable properties. However, very recently, a study of the molecular reorientation and local structure in liquid benzonitrile has been done through simulations and analysis of time correlation functions.[24] Other properties such as density, viscosity, and diffusion have been simulated. The single-molecule rotational dynamics also have been simulated and studied, but the notable specific trend of temperature-dependent viscosity, which had been reported earlier experimentally at about 308–313 K, has not been addressed.

The purpose of the present study is to investigate the thermodynamic, structural, and transport properties of BZN liquid in the range of 293–323 K by utilizing classical MD simulation. In particular, the phenyl ring interaction in the stacking mode and the contribution to the liquid state structure, morphology, and thermophysical properties will be studied vastly. Finally, the effort is also concentrated on force field development. To our knowledge, the singularity in the trend of temperature-dependent viscosity observed experimentally for benzonitrile is confirmed here using the simulation-specific analyses for the first time.

Method

To search and investigate the interaction involved and responsible for the BZN liquid structure and thermophysical properties, the simulations are performed at different thermodynamic states (i.e., in the temperature range 293–323 K). Particularly, the stacking phenomenon through the phenyl group (being absent in counterpart acetonitrile solvent) and hydrogen bonding and their role are studied by different methods. These have been the reasons for the initiative of an extensive study in line with our recently conducted studies on clusters of different sizes of benzonitrile molecules in which the effect of different possible hydrogen bond formations were examined on their stability.[11] Additionally, various studies on liquid crystals and stacking,[25−28] and the particular effect of hydrogen bond formation[29] are available in the literature.

Computational Details

The benzonitrile molecular structure was optimized using the Gaussian 09 program.[30] Geometry optimizations of BZN monomer were performed by utilizing B3LYP functional and 6-311++G(3df,3pd) basis set. Vibrational frequencies of the neutral compound were examined to ensure that the BZNʼs structure has a qualified potential energy surface. Using CHELPG procedure,[31] as a Gaussian 09 utility, the partial charges for BZN atoms were determined. Figure displays the BZN molecule and the atoms label used in all computations and simulations. The calculated partial charges are shown in Table S1 (in Supporting Information), and BZN geometrical parameters have been given and verified in the recent study.[11]

Figure 1

Structure of the optimized BZN molecule. The atom labels used throughout are shown. Also, the position of the center of mass (com) and the hypothetical vector from the com to the center of the ring (cor, ×) in the BZN molecule are shown.

All simulations were performed using the DL_POLY program version 2.17[32] to investigate the structural properties of BZN liquid and calculate the bulk thermodynamic (density) and transport properties (diffusion coefficient and viscosity) at ambient pressure (P = 1.01325 × 105 Pa). The non-polarizable optimized potentials for liquid simulation all-atom (OPLS-AA)[33] force fields were used in all simulations. Optimized potential includes intramolecular and intermolecular potential. The interaction potential energy function is given below

Potential energy contribution is accounted as the followingswhere all the terms have their usual meaning. Intermolecular interactions involve atom–atom (12-6) Lennard-Jones (LJ) potential for the van der Waals interactions and Coulombic electrostatic interactions between point charges centered on the atoms

The potential well-depth ε, and the hard-sphere diameter σ, is determined by simple combining and mixing rules, respectively

The partial atomic charges for this solvent were taken from our CHELPG calculations. The ε and σ parameters are given in Table S1 (in the Supporting Information). Equilibrium bond distances and angles were based on optimized geometry. The bond and angle force constants and dihedral parameters were obtained from the literature.[34,35] The force field parameters are presented in Table S2 (Supporting Information).

For MD simulation, a geometry-optimized BZN structure was enfolded to make a simulation box [50 × 50 × 50 (Å3)] containing 343 BZN molecules. The equation of motion was solved using the Verlet–Leapfrog integration algorithm. Periodic boundary conditions are applied in three dimensions. By using Ewald summation with the precision of 1 × 10–5, coulombic pair interaction in the form of non-bounded LJ interaction was calculated. The potential cut-off distance of 24 Å was set appropriately for the simulation box involved.

The simulation was initially performed in the NVE ensemble with a short time step (0.0001 ps) for 1 ns for BZN molecules to be randomly distributed without overlapping. Then, after the heating (up to 423 K) and cooling process, the simulation was continued for 5 more ns with 0.001 ps time step in NPT ensemble (for simulation of density) under the Hoover thermostat and Andersen barostat algorithm. Then, the simulation was continued, to calculate transport properties, for 10 more ns in the NVT ensemble. The simulation was started from the upper temperature of the temperature range of interest (293–323 K) in 5 K intervals. To simulate the equilibrium and transport property, the output of one NPT ensemble was used as the input of the initial ensemble at the next lower temperature successively. Self-diffusion coefficient and viscosity were determined using mean square displacement (MSD) (initially calculated from the simulated trajectories of the NVT ensemble). The specification of the simulation boxes is shown in Table S3 (Supporting Information).

Also, to investigate the BZN molecular surface properties, 540 BZN molecules were simulated in a cubic box of the sizes 50 × 50 × 50 (Å3). After equilibration, the ensemble was extended in the z-direction to 150 Å to create a liquid slab with a final thickness of about 50 Å. Similarly, MD simulation was performed in the NVE ensemble for 1 ns before performing the heating–cooling cycle. After cooling the ensemble down to 298 K, MD simulation was continued in the NVT ensemble for 20 more ns.

Results and Discussion

The liquid bulk properties including structure, thermodynamics, and transport as well as the surface properties have a fundamental relationship with the molecular structure. We present the result of computational studies that incorporate the stacking phenomena of the BZN phenyl ring.

Bulk Liquid Benzonitrile

Density, Surface Tension, and Force Field Assessment

The liquid density (equivalently the hard-sphere diameter) and intermolecular forces are the most relevant quantity in the thermodynamics of fluids from which the basic parameters of the equation of state can be derived. Fortunately, liquid density can be measured accurately by experimental methods easily. Because the liquid density also depends on the intermolecular forces, this single parameter is usually utilized to test the validity of the force field applied in the simulation. Table contains the simulated density of liquid benzonitrile in the range of 293–323 K at one atm pressure. The densities obtained are in good agreement (with a maximum deviation of 0.21%) compared to the experimental densities.[36,37] The result of the simulation can be tested visually and verified numerically with regard to literature experimental data (see Table and Figure ). The low deviation between simulation and the reported experimental densities warrants the simulation procedure and validates the OPLS-AA force field applied.

Table 1

Comparison of Simulated Densities (at 1 atm) with Experimental Values at Different Temperatures[36,37]a

	ρ, (% Dev.)b (g/cm3)
		experimental
T (K)	simulation [this work]	(36)	(37)
283.15			1.0140
288.15			1.0096
293.15	1.0041		1.0051, (0.09)
298.15	0.9990	1.0009, (0.19)	1.0007, (0.16)
303.15	0.9951	0.9964, (0.13)	0.9963, (0.12)
308.15	0.9898	0.9919, (0.21)	0.9918, (0.19)
313.15	0.9866	0.9872, (0.06)	0.9874, (0.08)
318.15	0.9827		0.9830, (0.03)
323.15	0.9783

The deviation with respect to two sets of experimental densities is also shown.

% Dev = 100 × (|ρexp – ρsim|)/ρexp.

Figure 2

Comparison of the simulated temperature-dependent density of the BZN with literature simulation and experimental. Lines through the point are used to visualize trends more clearly. No trendline is shown for data points of ref (36) to better visualize its slight fluctuation at 313 K from linear behavior while extensively overlapping with ref (37). The plot for ref (24) included is based on data we extracted (ref (24), Figure ).

It can be seen that our simulated densities are closer to the experimental data reported in refs (36) and (37) than those simulated in ref (24). Both the trend of our data and those simulated in ref (24) show fluctuations and deviations from the linear behavior, though ours is appearing with enhancements. The reported densities of ref (36) are in a limited range (298–313 K) and it shows some fluctuation (at 313 K), while the experimental data of ref (37) are linear in an extended range of 283–318 K (see Figure ). This fluctuation enhancement can be attributed to the long cut-off distance of 24 Å (we used in the DL_Poly environment) compared to the 12 Å (ref (24) used in the Gromacs environment). Therefore, such fluctuation enhancement could be attributed to the finite size and the simulation conditions (this work and ref (24)) while is likely observable in experimental data (ref (36)). Whether these fluctuations occurring around 308–313 K has something to do with the fluctuation of the viscosity (as we will discuss in subsequent sections) is subjected to the more detailed simulation and experimental measurements.

MD simulation could be used as a prevailing tool in the prediction of surface properties of organic liquids. The surface tension, γ, is proportional to the integral of the difference between the normal components of the pressure tensor: γ = −b(P + P + P)/4, where b is the length of the simulation box along the box principle axis (the z-axis) perpendicular to the liquid/vapor interface, and Ps are the principal components of the pressure tensor. The result of surface tension simulations (for 20 ns in NVT ensemble) of liquid BZN at 298 K was found to be 38.35 mN/m, being in good agreement with the experimental data (38.33 mN/m)[38] within 0.1 mN/m. Such a high level of agreement ensures the accuracy of the refined force field parameters.

Liquid Structure

The investigation of structural properties of liquid benzonitrile involves the simulation of the atom–atom radial distribution function (RDF). We proceed under the condition that the particular hydrogen bonding structures have been already established and confirmed based on the cluster formation of different sizes,[11] as shown in Figure a. The g(r)’s for the −CN group and different H atom sites at equilibrium are shown in Figure b. The first tall and narrow peak indicates that the correlation of the N atom with different H’s (ortho, meta, and para) are rather strong and of low dynamics at short ranges. The first peaks of N···H1(H2) (interacting at ∼2.95 Å) indicate (both equally) strong correlations with low dynamics in the short-range (compared to other cases in Figure b) and persist well in the long-range (second peak). The fact that a conformer comprising, for instance, one benzonitrile dimer forming two strong ortho-H bonds simultaneously [like −C≡N···H1(H2)] can thoughtfully account for this.

Figure 3

(a) Hydrogen bonding is identified in the dimer and trimer BZN clusters. (b) Simulated RDF of N and different H atomic sites for liquid BZN at 298 K.

The first peak is generally well-shaped in a similar way for the three different H atoms, while peak height is slightly different in the order of meta-H ≈> para-H ≈> ortho-H having a width in the order meta-H ∼> para-H ≫ ortho-H. These peaks located in the 2.8–2.9 Å domain could be attributed to the short-range electrostatic interactions enhanced due to hydrogen bondings of moderate strength, being in good agreement with the recent ab initio studies performed by considering benzonitrile clusters up to tetramers.[11]

Owing to such a behavior, we also investigated the temperature dependence of g(r)’s in the range of 293–323 K for different H-bondings and focused on peak heights, as plotted in Figure . As can be seen, the peak height is not considerably different among different −C≡N···H cases (being different approximately by about 0.05–0.07). Therefore, the lower dynamics of −C≡N···H1(H2) brought up a conclusive assessment of the possible H-bonding’s strength and its influential value on liquid structure in the order of ortho-H ≫ meta-H ∼> para-H. The subsequent peaks (particularly, the second peaks) are an indication of rather good coordination of dipoles at long ranges, whereas −C≡N···H1(H2) shows a second peak appreciably stronger comparatively. The fluctuation seen in Figure implies fluctuations in the simulated densities (Figure ), being arisen from the superpositions of the effect of different H-bond types for the temperature-dependent density. For further investigations, other analyses such as dimer existence autocorrelation functions (DACFs) and the evaluated structural condition were calculated by the TRAVIS package.[39]

Figure 4

Temperature dependence of the peak heights of the RDFs. ortho-H[=H1(H2)], meta-H[=H3(H4)], and para-H[=H5].

The autocorrelation function is useful in the calculation of the hydrogen bond population based on the structure and geometry of the liquid at a certain time and the dimer existence lifetime population then follows. We calculated the lifetime that two specific atoms involved in hydrogen bonding by integrating the DACFs at 298 K, as shown in Figure S1 (in the Supporting Information).[39] The dimer existence due to the possible hydrogen bond formation by the three kinds of Hs can be assessed, as seen intuitively are decaying at different rates. Assuming the dimer existence decay exponentially over time, the calculated lifetimes are ∼9.9, 2.5, and 0.23 ps for the correlation between N and ortho-, meta-, and para-H atoms, respectively. Therefore, the lifetime for the ortho-H···N interaction is well above others, confirming the importance of ortho-H in the intermolecular interaction. Also, the percentage of the simulation time that two atoms experience a strong (−C≡N···H) correlation has been calculated by the evaluation of structural conditions at all temperatures (Figure S2 in the Supporting Information). These calculations were based on the evaluation of the H-bonding lifetime when any N···H has a distance within in 0–3.0 Å interval optionally, though, conclusively representing the lifetime at the most probable distance of ∼3.0 Å, which can be deduced from Figure . This structural condition has been evaluated by considering the time that each kind of H-bonding is persisted longer. By combing the results in Figure for bond length with those of Figure S2 for the lifetime, an important conclusion brought about by the MD simulation of bulk BZN is an assessment of the (−C≡N···H) interaction. Hence, the relative structural condition that determines the liquid BZN with an overwhelming structure is in the order of ortho-H > meta-H ≫ para-H.

Therefore, considering the results in Figures and 4, it can be concluded that the most strong interaction between benzonitrile molecules occurs through ortho-H···N with low dynamic while lasting longer. Earlier ab initio calculations have shown that the H-bonding in BZN through ortho-H···N is well-contributed specifically by both σ- and π-bonding.[11] Therefore, the possible H-bonding types that have been already demonstrated, for instance, based on the formation of clusters of different sizes (as shown in Figure a),[11] also confirm the above trend in the correlation strength.

Structural Properties: Combined Distribution Function

Further structural analysis was elaborated by utilizing the combined distribution function. Such a goal was achieved by mapping combined radial and angular distribution functions (ADFs) using the TRAVIS package.[39] Therefore, a counter of ADFs versus RDF made it possible to map and understand by what angle (θ) does the −CN group approaches (r) the H atom (Figure a) while is keeping an appreciable correlation. The results illustrated in Figure b is representing a well-defined and precise ortho-H···N correlation, with a maximum between θ = 120 and 135°. For the meta and para cases, the characteristic simulated angles θ ∼ 180° (Figure c,d) are revealing and confirming the earlier finding[11] that the para-H···N and the meta-H···N interactions occur through the σ-bond and π-bond formation, respectively.

Figure 5

(a) Representation of angle θ and distance r, and (b–d) the combined distribution function for hydrogen bonding between N and ortho-, meta-, and para-H atoms.

Transport Properties

In addition to calculating the equilibrium thermodynamic properties (leading to force field development) and the liquid structural feature, the dynamics of the BZN molecules in the liquid state can be supplemented by investigating the simulation of transport properties. Therefore, the MSD and the self-diffusion coefficient of benzonitrile liquid were calculated and shown in Figure a,b, respectively. The MSD was calculated from the knowledge of coordinates of the particles during the simulation in the NVT ensemblewhere r(t) represents the com location of the particle i at time t, and N is the number of particles.

Figure 6

(a) Simulated MSD of the com of benzonitrile molecules in liquid state at different temperatures. The thick lines are composed of a copious number of overlapping data points. Dotted lines are the linear fits. (b) Temperature dependence of the calculated self-diffusion coefficient of benzonitrile liquid. The high and low values of uncertainties are 0.04 and 0.01(error bars), respectively.

Self-diffusion coefficient (D) represents the microscopic liquid dynamics and were obtained from the linear regime of MSD curves at long simulation time using the Einstein relation

The self-diffusion coefficients for the benzonitrile liquid were calculated (by performing 10 ns simulation) in the temperature range of 293 to 323 K as plotted in Figure b. As can be seen, self-diffusion for BZN liquid increases (almost quadratically) with increasing temperature while showing a fluctuation at 313 K. The statistical uncertainties for self-diffusion were determined by repeating simulation in ensembles with particles having different initial coordinates. Uncertainties reported for self-diffusion were calculated based on three sets of such data collection. In addition, we calculated the viscosity of benzonitrile liquid (Figure ). For characterization and investigation of the microscopic dynamics of the liquid bulk, the viscosities of BZN were estimated using the Stokes–Einstein equationwhere kBT is the thermal energy at temperature T and α is the molecular diameter. By selecting c = 6 for benzonitrile, the results are shown in Figure . To determine a mean value for the length parameter α, we employed the results of quantum mechanical calculation [b3lyp/6-311++g(3df,3pd)] using the Gaussian 09 program[30] and determined the volume of the benzonitrile molecule. The simulated viscosities are plotted in Figure and compared with the experimental data. Similar to all liquids, the simulated BZN viscosity decreases monotonically as the temperature increases, though a fluctuation is seen centered at about 313 K (consistent with the profile of the diffusion coefficient (Figure b). Quite interestingly, a fluctuation observed in the experimental viscosity of benzonitrile,[23] in the same temperature range as demonstrated in Figure , is consistent with the present simulation results. Therefore, the simulation we performed with precaution in a wide range of temperatures in the reduced steps of 5 K could unravel the singularity at 313 K prominent. Also, in a more recent experimental study,[40] despite the small temperature range and the limited number of data points made available, a fluctuation is likely at 308 K.

Figure 7

Calculated viscosity of benzonitrile liquid at different T values. Lines are trend lines.

In addition, this singularity (albeit weak) has been indicative not only in experimental studies but also in a reported simulation study, though not addressed and left unattended.[24] The limitation in data acquisition and presentation (i.e., the long step size of 10 K) looks to be a sort of obstacle for the singularity not to be seen strictly. To remove such a hindrance, we coined more data points by averaging viscosity values (Figure , ref (24)) and diffusion values (Figure S8, ref (24)) between two of their consecutive measurements. Therefore, by such extension of data points, the singularity in the trend of both viscosity and diffusion can be observed clearly, as shown in Figure S3 (Supporting Information),[24] which are confirming our findings in the present work accordingly.

Benzonitrile liquid properties have been an effort and studied widely as a commitment in literature;[24] nonetheless, the simulation for understanding viscosity and diffusion coefficient in detail requires choosing a wide range of temperatures in small steps. Especially, when transport properties are simulated through the time autocorrelation function, sampling a large number of simulated trajectories is essential for the best result to match with the experimental data.[41] On the other hand, the cut-off distance must be selected long enough to account for the LJ interaction perfectly and, hence, to enable distinguishing a real singular trend from the fluctuations in the calculations.

Therefore, these fluctuations are part of the viscosity profile, confirming the validity and accuracy of the simulated viscosity in detail. Under the condition that no interpretation has been given to the observed experimental viscosity to date, we got the impression that the nature of these fluctuations has roots in the factors that can be clarified by the simulation results.

Our experience in this group from the simulation of liquid BZN (by Gromacs 4.5.5[42]) with a large size NVT ensemble (up to 2575 molecules, cut-off distance of 14 Å, OPLS-AA force field[33]) was appealing. Although these simulations were done for slabs of BZN with a large size (liquid/vapor) surface area (80 × 80 Å2), the trend of temperature-dependent viscosity confirms the same fluctuation in the range of 308–313 K.

Specific Molecular Structures in Liquid BZN

To date, several studies have been performed on π–π interactions in liquid benzene and some of its derivatives, such as benzonitrile.[43] In these compounds, due to the dispersion attraction between the positively charged rings and the π electrons, the two rings may combine in a face-to-face sandwich configuration. Certainly, factors such as electron-withdrawing and electron-donating of the benzene ring having substituted with functional groups would affect this behavior. For example, in benzonitrile, the electron-withdrawing substituent −C≡N group helps to enhance the configuration with π-stacking interaction as a result of electron withdrawal from the conjugated π bonds and thus reducing the electrostatic repulsion. For this reason, a 100 ps simulation was performed (Hyperchem 7)[44] on benzonitrile clusters of various sizes at ambient temperature. As shown in Figure , the final structures led to clusters, all interacting in a ring stacking mode.

Figure 8

Appearance of the stacked clusters of different sizes of benzonitrile produced at equilibrium from the initial configuration, (a) dimer-1, (b) dimer-2, (c) trimer-1, (d) trimer-2, (e) tetramer-1, and (f) tetramer-2.

For verifying the behavior of BZN molecules in the bulk environment including possible stacking, the RDFs between the com of molecules were analyzed by utilizing a 10 ns simulation at different temperatures. Although the stacking assessment may be done by using RDF between different parts of the molecules while are correlated in the stacking mode, we found that the assessment based on the coms’ (imaginary sites) correlation produces results for comprehending best in our presentations. This choice indeed puts all things together for visual inspection of the basic facts, relevantly and correctly. In another sense, this produces a module for angle-dependent pair correlation function in a relevant way (see Figure for the position of the com).

Accordingly, rather well-defined RDFs are produced with a sharp peak (shown in Figure a for different temperatures) that implies a good correlation between the centers of mass of the BZN molecules. Notably, the RDFs of Figure are in complete agreement with other simulation reported.[24]

Figure 9

RDFs for com–com of BZN molecules in liquid bulk at different temperatures. The g(r)’s values are extensively overlaid at all different temperatures considered.

A shoulder inappropriately masked underneath the main peak contains interesting information due to its distinct position at about 3.5–4.0 Å. The distance of the closest approach of the masked peak (at about 3.32 Å, Figure ) is notably close to the spacing between two graphene sheets.[45]

To study the structure and relative orientation of molecules in liquid BZN (especially in the two ranges centered at ∼4 and 6 Å), the analysis of the combined distribution function would be appealing. Accordingly, different RDFs are combined and their most probable correlations, where they match and indicate simultaneous occurrences, are examined. As shown in Figure a, both RDFs between the com and the cor of BZN molecules (see Figure ) show the most probable peaks that intuitively match at about ∼4 Å. Therefore, stacking of the two BZN may occur, while the two ring plans are grafted together (say in the antiparallel mode) through the correlation of (the two imaginary points) com with cor at about 4 Å. Therefore, such a plot could be utilized to deconvolute a shoulder peak as an existing peak.

Figure 10

Combined radial–RDF (at 298 K) for (a) cor–cor, (b) N···H1, (c) N···H3, and (d) N···H5 RDFs vs com–com RDF.

Figure b shows the combined distribution function of the RDFs between the com of BZN molecules and the RDFs of ortho-H···N hydrogen bonding. In this form, the most probable state matches with the distance of about 6.23 Å for the com–com of BZN, and the ortho-H···N makes hydrogen bonding at the distance of 2.8 Å. That is, considering the long com–com distance, BZN molecules can also form hydrogen bonds through ortho-H without conforming to the structures that lead to the stacking mode.

Further abstract correlations are combined in Figure c,d, making the use of possible hydrogen bond formation involving meta-H and ortho-H. As shown in Figure c,d, under the condition that the com of two BZN molecules are drifted away from the effective possible stacking mode, strict hydrogen bond formation (e.g., the meta-H···N and para-H···N) is possible with high probability. However, in considering both meta- and para-H involvements, a less likely domain is noted to emerge at about 4 to 5 Å. This can be attributed to the formation of intermediate structures with the two BZN almost flat and able to form hydrogen bonding, while stacked in parallel (antiparallel) structures.

Further analysis has included studying the hypothetical vector defined from the com to the cor laying in the BZN molecule plane (Figure ). Then, the angular distribution among different molecules in the ensemble was calculated to estimate the extent and role of the relative molecular orientation of the possible clusters. Combining this vector orientational distribution with the RDF produces maps suitable for deconvoluting the complex RDF of Figure , which would be otherwise very difficult (if not impossible) by the conventional deconvolution mapping popularly done for the analytical curves.

The combined distribution function shown in Figure (for 298 K) indicates that the correlation between the com of two adjacent rings in the stacking mode at about 4 Å (in a dimer, trimer, ...) is strongest where the most probable angle between their in-plane vectors (Figure ) is centered strictly at 180° and 0°. However, the relative orientation of two vectors at 180° (anti-parallel) are much more probable than at 0° (parallel). The concluded structures are schematically demonstrated in Figure , which are in good agreement with the quantum mechanical studies.[46] More precisely, the simulation results are distance-wise in good agreement with the 3.2–3.4 Å intermolecular plane distance reported by quantum mechanical calculation.[46]

Figure 11

Combined angle distribution function for the (com–com) in-plane-vector vs distance corresponding to RDF (of Figure ) at 298 K. The molecular relative structures that are implied consistent with the mapped angle distribution are also illustrated.

Figure 12

Schematic representation of the most probable states for BZN with molecules oriented in the plane-stacking mode of (a) anti-parallel (180°) and (b) parallel (0°). The coms are marked by a pink sphere. (c) Temperature dependence of the peak heights of the RDFs (Figure ).

According to Figures and 12, geometrical states with the planes of the BZN molecule stacked together are characterized here by the angle between the in-plane-vectors, having occurrence probability in the order of 180° ≫ 0° (Figure ). In short, based on the deconvolution scheme we followed, the two adjacent BZN molecules could correlate and form a stable dimer (alone or participating in trimer, tetramer, and...) effectively by stacking in an antiparallel orientational conformation.

For BZN molecules having coms at ∼6 Å apart (Figure ), the most probable angles formed between the vector (shown in Figure ) are in the order of 180° > 90° > 0°. At this distance, as shown in Figure b–d, BZN molecules give rise to non-stacking structures and form hydrogen bonding with a coplanar configuration. Such structures, along with their in-plane vectors, are given as an example in Figure . Therefore, where the two coms are ∼6.23 Å apart, the most probable angle formed between the two in-plane vectors is 180°. This matches with the possibility of hydrogen bond formation of ortho-H···N between two BZN molecules placed coplanar. Also at this distance, the probability of two in-plane-vectors at 90° and 0° angles indicate the formation of a hydrogen bond of meta-H···N and para-H···N, respectively (see Figure ).

Due to the effective repulsion between cyano groups at short ranges (∼4 Å in stacking mode), the benzonitrile molecules in the predicted most probable orientation of 180° do not overlap exactly and show slight shifting, as shown in (Figure a). From RDFs for −C≡N···C(1–7) in the anti-parallel orientation (Figure S4 in the Supporting Information), the correlation of C2 and C3 (as well as C4 and C5) atoms with the N atom is unraveled by the characteristic peak maximum at 3.6–3.7 Å and the distance of the closest approach of ∼3.0 Å. More importantly, the rising edges of correlation functions (at short range) for N atom with C atoms of the ring (i.e., −C≡N···C(2–6)) have all coincided with the same slopes and distances of closest approaches. This feature boosts the conjecture that all the N atoms approaching ring C atoms are positionally and orientationally the same and happens mostly when they are in the stacking mode with the molecular planes taking a parallel orientation. Such interactions can arise practically through a simultaneous bilateral interaction, involving the N atom of either BZN molecules dimerizing in the most probable antiparallel 180° configuration.

To try reasoning viscosity acting singular at 308–313 K, we did other analyses (similar to those shown in Figures c and 11 for 298 K) at the whole range of temperatures. These results are given in Figures S5 and S6. As can be seen in these two figures, the amplitude of the probabilities decreases somewhat at about 308 K and increases back again at 313 K. Correspondingly, the peak heights we extracted from Figure (the correlation function between the BZN com with the first peak at ∼6 Å representing the conformations for intermolecular interaction through hydrogen bonding) show fluctuation (Figure c) similar to the temperature-dependent viscosity in the range 303–313 K. Therefore, in this temperature range, the system is more prone to hydrogen bond formation. Note that similar changes in the trend of correlation strengths between individual N···H atoms (especially meta and para-H) occur in the temperature range 308–313 K as demonstrated in Figure . Conclusively, the singular behavior of viscosity in the temperature range of 303–313 K is due to the change in the strength of different N···H hydrogen bonds, especially at 313 K temperature, and hence the resulting structural alterations. Conclusively, the singular behavior of viscosity in the temperature range of 303–313 K occurs due to the change in the strength of different N···H hydrogen bonds, leading to structural alterations, especially at 313 K.

Conclusions

MD simulations remarkably have enabled us to determine the details of bulk liquid benzonitrile structural properties at equilibrium and utilize them to enlighten microscopic specific phenomena such as temperature-dependent (293–323) density, viscosity, and diffusion. After elucidating cluster formation in liquid BZN recently,[11] the deconvolution and extraction of the orientational stacking mode were performed by appropriate implementation of the angle-dependent pair correlation function. Accordingly, two BZN molecules prefer interacting in the stacking mode while their principle molecular axis is angled, centered at about 180°. Simulation at different temperatures (within 293–323 K) unravels, consistent with experimental data, a specific singularity in the trend of the temperature-dependent-viscosity centered at about 313 K. According to the further simulation analysis, switching between the configuration involving ortho-H···N interaction (between dimer molecules) and the configuration involving simultaneous ortho-H···N and meta-H···N interactions (between trimer molecules) are responsible for the occurrence of the singular behavior of viscosity.