Inhomogeneous Monte Carlo simulation of the vapor-liquid equilibrium of benzene between 300 K and 530 K

The inhomogeneous Monte Carlo technique is used in studying the vapor-liquid interface of benzene in a broad range of temperatures using the TraPPE potential field. The obtained values of the VLE parameters are in good agreement with the experimental values as well as with the results from GEMC simulations. In contrast to the GEMC, within one simulation box the inhomogeneous MC technique also yields information on the structural properties of the interphase between the two phases. The values of the vaporization enthalpy and the vapor pressure very well satisfy the Clausius-Clapeyron equation.


Introduction
Benzene belongs to the cyclic hydrocarbons, a class of organic compounds which are frequently studied by molecular simulations due to their great industrial importance.It is also one of the simplest molecular systems in which the hydrophobic interactions play a significant role; many systems of the kind are of use in environmental, clinical and biological applications [1][2][3][4].
Several potential energy models were developed to describe the molecular interactions of hydrocarbons.One of these potential fields for hydrocarbons -TraPPE (Transferable Potentials for Phase Equilibria) -shows good agreement with experimental Vapor-Liquid Equilibrium (VLE) data in a broad range of temperatures [5][6][7].Later on the TraPPE-parameters for alcohols, ethers, glycols, ketones, and aldehydes were determined [8,9].Moreover, a version with explicit consideration of the hydrogen atoms was established [10].
One way of studying the vapor-liquid equilibria of molecular fluids is using the Gibbs Ensemble Monte Carlo technique (GEMC) improved by Siepmann et al. [11] in order to determine the vapor-liquid coexistence curves of several n-alkanes.Another way of studying the vapor-liquid equilibria (VLE) is the direct simulation of two coexisting phases in a box elongated in one direction and containing two phases with different densities.These phases are separated by an interface perpendicular to the elongated edge of the box in order to minimize the surface energy.This technique was used in studying the vapor-liquid interfaces in various systems -Lennard-Jones fluids [12], molten salts [13], water [14][15][16], aqueous solutions [17,18] etc.This method can also be modified in order to study the liquid-liquid equilibria between two immiscible liquids [19,20].The use of boxes with interface in the Gibbs Ensemble Monte Carlo also represents a powerful technique which enables us to study the properties of interfaces between vapor and saturated solution of two immiscible liquids [21][22][23].Molecular Dynamics (MD) simulations in such boxes were also used in studying the free surface of alkane oligomers [24][25][26].
One common problem of all these simulations is the phenomenon that the equilibrium propertiesespecially those of nonpolar systems -significantly depend on the type of the truncation of longrange dispersion interactions and on the applied cut-off radii.
Several approaches were suggested to treat the truncated dispersion interactions during the inhomogenous simulations.Guo et al. [27,28] suggested an approach based on the assumption of the local dependence of thermodynamic properties.Mecke et al. [29] have included the long-range corrections within MD simulation by adding an additive force contribution in the direction perpendicular to the interface.In the case of Lennard-Jones fluids both methods give better results in comparison with simulations without any long-range corrections.Nevertheless the density profiles and the surface tension show at higher temperatures a significant dependency on the cut-off radius.Much better results were obtained by López-Lemus et al. [30,31] using the lattice sum technique; nevertheless this improvement is provided due to a significant increase of the computational expenses.
Recently we have proposed a new method for treating the long-range corrections in systems with planar interface treated by the Inhomogeneous Monte Carlo simulation [32].In the case of Lennard-Jones particles this method works surprisingly well.For cut-off distances of 2.5σ and of 5.0σ this method gives identical density profiles and surface tension values for temperatures up to 95% of the critical temperature.Since the increase of the computational time is only about 20% compared to simulations without any long range corrections.In loc.cit.[33] we applied this recent approach to molecular fluids.
In this work we employ the inhomogeneous simulations with long range corrections for the calculation of VLE at elevated temperatures nearly up to the critical point.As a representative system we take benzene as rather rigid molecule.
In the next section we briefly describe the employed simulation technique, explain the determination of quantities characterizing the phase equilibria, and give a short survey on the technical details of the MC simulations.In the last part of this work we discuss the obtained results.

Simulation details
In the recently suggested method of treating the long-range interactions in inhomogeneous simulations the simulation box consists of three compartments separated by interfaces perpendicular to the z-coordinate.The outer two compartments are treated as identical phases.In order to evaluate the density and the energy profiles, the box is subdivided into strips ∆z of equal thickness.Then, the basic idea consists in the separation of the interaction energy U i of molecule i into the contribution of the interactions with the molecules within the cutoff sphere and the long range part depending on the z-coordinate and on the distance dependent density of all components in the system where u ij is the pair potential between molecules i and j and U lrc i (z) represents the long-range correction term.In the case of a one-component system formed by molecules containing N s sites the general relation [33] is reduced to where the first summation runs over all strips into which the box is divided and the second over all sites of molecule a and the last one over all different sites.ρ β (z k ) is the local density of site β in a strip of thickness ∆z centred at the position z k .In the case of the interaction of Lennard-Jones type the function w αβ (ξ) has the following form Clearly, for cut-off radii larger than 2σ, the long range correction of the repulsion interactions (first term in the brackets) is negligible compared to that of the dispersion interactions.When the long range corrections of the energy are involved in the MC algorithm, the behavior of the system rapidly tends to reach the behavior of the system with non-truncated interactions.The critical temperature of Lennard-Jones fluid with R c = 2.5σ without long-range corrections is T * = 1.186 [12], while the inclusion of long-range corrections enables us to perform the inhomogeneous simulations of Lennard-Jones fluid at T * = 1.25 (which represents 95% of the critical temperature of the fluid interacting via the full Lennard-Jones potential).The values of the coexisting densities and of the surface tension are in agreement with results obtained by other different simulation techniques [32].The equilibrium and interfacial properties of fluids can thus be obtained using the systems of considerably smaller size than it was so far.
The methods for determining the equilibrium properties from inhomogeneous simulations are straightforward.The densities of coexisting phases can be directly obtained from the density profiles which represent the basic output of the inhomogeneous simulations.This is demonstrated in figure 1a for the results of TraPPE benzene at T = 423.15K.The densities of both phases are the averages of the respective local densities.A snapshot of the corresponding simulation box is included in this figure under the density profile.The vapor-liquid equilibrium at a given temperature in a one-component system is also characterized by the pressure.In the case of inhomogeneous simulations the pressure acting on the simulation box is different in the z-direction perpendicular to the interface (normal pressure) and in the lateral direction (transversal pressure); this asymmetry is caused by the non-zero value of the surface tension.The pressure acting on the walls of the system (simulation box) is related to the diagonal components of the virial tensor Π where the term N k B T represents the ideal gas contribution while the short-range part of the non-ideal contribution is given by the sum of products of the components of the force acting on particle i, F µ i , multiplied by the corresponding component of its position vector, r i ; µ = x, y, z.The expressions for the long-range corrections of the components of the virial tensor Π lrc µµ (z i ) can be expressed in a form similar to those for the energy (equation 2).We do not repeat these formulas explicitly; interested readers can find them in [32] or [33].Since the normal component of the pressure is constant in each part of the system, which is the condition for the mechanical stability of the interface, it can be used in evaluating the vapor pressure, P • ≡ P zz .With respect to large fluctuations of the forces within the liquid phase it is more convenient to determine the vapor pressure from the profiles of the transversal or normal pressure along the z-axis.The normal pressure profile in the z-direction can be constructed by the contributions of all interactions of molecules i and j placed in two different strips and all the strips between them.Similarly the transversal (also called tangential) pressure profile is calculated as average from the interactions acting in x-and y-directions (see e.g.[12] for simple fluids or [25] for the generalization for molecular liquids).The most direct way of including the long-range corrections is to evaluate them from the corrections of the components of the virial tensor and add them to the strip in which the given molecule is placed; in a one-component system we can write for the transversal pressure profile where ρ(z) is the density profile, A is the area of the interface, x ij , y ij and z ij are the components of a vector pointing from the center of molecule i to j, f ij is the force acting between these two particles and θ(x) is the Heaviside step function.The normal pressure profile, P N (z) is given by a formula of the same form but with factor z 2 ij instead of (x 2 ij + y 2 ij )/2 and Π lrc zz (z) instead of Π lrc xx (z) and Π lrc yy (z).The vapor pressure is then represented by the average of these pressure profiles within the vapor phase where the scatter is much smaller than in the liquid phase (see figure 1b).
Another quantity related to the vapor-liquid equilibrium is the enthalpy of vaporization.Since the interaction energy of molecules is calculated in every MC step it does not cost any additive work to also construct profiles of the interaction energy as a function of the z-coordinate.From such a profile the excess internal energy of molecules in both phases can be determined similarly to the determination of the coexisting densities (see figure 1c).The vaporization enthalpy can then be obtained as where the subscript m denotes the molar quantities (H: excess enthalpy, U : excess internal energy and V : volume) while the superscripts v and l identify the vapor and the liquid phase.The molar volumes are obtained from the density profiles.At low temperatures it is possible to replace V v m by the expression of the ideal gas, because it is impossible to estimate the gas-phase density from the simulation.
The surface tension and the thickness of the interface can be calculated in the same way as in [33].
The simulation conditions are almost identical with those used in our previous study [33].The system contains 343 molecules of benzene.These are represented by planar rigid hexagons with the vertices representing the united atoms CH; the distance between these united atoms is r CC = 1.4 Å and the sites of two different molecules interact through the Lennard-Jones potential with parameters σ CC = 3.695 Å and CC /k B = 50.5K.The initial configurations were obtained by NPT equilibration of a simple cubic lattice containing 343 molecules.The recommended value of the cut-off radius for the TraPPE potential, used at all simulations, is R c = 14.0 Å.This demands a box length of at least 28 − 30 Å in the x− and y− directions.A cubic box of this length contains about 160 benzene molecules in the liquid phase at T = 298.15K.It is necessary to increase this number in the simulation with respect to the free space for the gas phase and the molecules in the interfacial region.Therefore we used the twofold number of particles.Because this number depends on the temperature it has to be verified in each simulation.
During the equilibration the lateral dimensions of the simulation box were kept constant, 30×30 Å, while the z-dimension was allowed to vary.After this phase the z-dimension was increased to 180 Å and an equilibration phase of 250000 Monte Carlo loops (250000 × 343 generated configurations) was run.Three production runs of the same length were then performed.The maximal translation and rotation parameters were adjusted in order to get an acceptance ratio between 30 and 50%.The simulations at temperatures T 498.15K were performed with 512 molecules in a box with dimensions 30 × 30 × 240 Å.After every 10th MC cycle the density profiles for sites and molecular centers along the z-axis were recorded in strips of thickness 0.25 Å after every 10 MC cycles.With the same frequency, the components of the virial tensor and profiles of normal and transversal pressure were also measured.
The molecular truncation [33] of the intermolecular energy was applied.This means that the truncation was based on the center-center separation of the molecules rather than on the separation of the sites.

Results
The values of the densities and the excess energies for the coexisting phases are listed in table 1.These values were determined in the way shown in figures 1a and 1c.The error bars were estimated as three times the standard deviation of a set of 5 block averages of thickness 10 Å which were chosen uniformly within the bulk phases.Figure 2 shows the comparison of our results with the GEMC data obtained by Wick et al. [7] and with the experimental data [34].Since the TraPPE potential field was parameterized to reproduce densities of coexisting phases in a wide temperature range, the perfect agreement with the experimental data is not surprising.The coincidence with the GEMC results shows that the inhomogeneous simulations can provide equivalent description of the phase equilibria if the finite size effects are correctly involved.
The values of the vapor pressure are summarized in the first two columns in table 2. In the second column the average values of the normal pressure profile within the vapour phase are shown.The third column contains data which were determined from the zz-components of the virial tensor.Evidently, if the vapor pressure is estimated by the average value of the transversal or normal pressure profile in the vapor phase, the statistical uncertainties are one order of magnitude smaller compared to the second method.The dependence of the vapor pressure on the temperature is plotted in figure 3. Again the results of this work (circles with error bars) are compared with GEMC results (diamonds) [7] and with the experimental data (solid line) [34].A slight overestimating of the vapor pressure is one shortcoming of the TraPPE potential field which was already recognized in the paper by Wick et al. [7].The knowledge of the equilibrium pressure enables us to calculate the heat of vaporization.Its values are listed in the fourth column in table 2 and the comparison with the experimental data is shown in figure 4. One can see that the employed potential model for benzene underestimates the values of the heat of vaporization; the relative difference is about 10 % which is similar to the extent in which the vapor pressure is overestimated.The main reason for this insufficiency of the TraPPE model may be the formation of stable dimers of benzene molecules [35] which cannot be described accurately by pairwise additive interaction site potential models.The enthalpy of vaporization is related to the temperature dependence of the vapor pressure through the Clausius-Clapeyron equation With respect to non-ideal behaviour of the vapor phase in the investigated temperature range, the volume change of vaporization cannot be approximated by the molar volume of ideal vapor.The simulated values of the vapor pressure were fitted to the Antoine equation ln where P • 0 = 1 bar and the fitted constants are A = (10.2± 0.5) and B = (−3450 ± 250) K; the differential quotient on the left hand side of the Clausius-Clapeyron equation can be then expressed as ∂ ln(P The consistency between the vaporization enthalpy and the vapor pressure according to equations ( 7) and ( 9) is shown in figure 5  For completeness, the dependence of the surface tension of benzene on the temperature is shown in figure 6.The solid line represents the values obtained by an extrapolation method suggested by Somayajulu [36] and the circles are the results of our simulations.We note that the surface tension values evaluated according to the method of Somayajulu are in agreement with the experimental data up to the normal boiling point [37].

Conclusion
Inhomogeneous Monte Carlo simulations were performed in order to investigate the vaporliquid equilibrium properties of benzene using the TraPPE potential field.The perfect agreement of the results of this simulation technique with those obtained by the GEMC as well as with experimental data is a profound affirmation of the reliability of this technique for the calculation of phase equilibria.The enthalpy of vaporization and the vapor pressure very well satisfy the Clausius-Clapeyron equation.

Figure 1 .
Figure 1.Profiles of thermodynamic properties of TraPPE benzene at 423.15 K along the zcoordinate of the simulation box (vap ≡ v, liq ≡ l).a) Upper part: Density profile and snapshot of the coexisting liquid and vapor phases.b) Middle part: Profile of normal and transversal components of pressure.The normal pressure profile is represented by the solid line while the dashed one is the transversal component.c) Lower part: Profile of the energy felt by one molecule.U v and U l are explained in the text.

Figure 5 .
Figure 5. Representation of the Clausius-Clapeyron equation (CC).The figure shows the right hand sides of equations (7) and (8) as functions of temperature.The uncertainty range of the expression ∆Hm/(P T ∆Vm) is shown by the vertical error bars.The dashed lines indicate the uncertainty range of −B/T 2 .
where the values of the quotient −B/T 2 and ∆H m /(P • T ∆V m ) are plotted against T .The dashed lines indicate the uncertainty range of −B/T 2 while those of the second expression are shown by the vertical error bars.

Table 2 .
Vapor pressure, enthalpy of vaporization and surface tension for TraPPE benzene; results of inhomogeneous simulations.