Near Critical Coexistence for an AUA Model of Thiophenes

Near critical coexistence has been determined by means of parallel tempering coupled with grand canonical Monte Carlo simulations which were later recombined by using histogram reweighting techniques. The data collected during the simulations is not only useful to determine accurately the critical point but also to provide estimates for the coexistence density jump between the phases in equilibrium. A recently introduced algorithm by Kim [Kim Y.C. (2005) Phys. Rev. 71, 051501; Kim Y.C. (2005) Comput. Phys. Commun. 169, 295], based on the scaling of the positions of the different minima found for the Binder parameter, has been applied to the case of a realistic model of thiophene consisting of different Lennard Jones sites. Contrary to the case of the Hard Core Square Well (HCSW) and Restricted Primitive Model (RPM) systems, significant corrections to scaling are found in this case. By readapting the algorithm we are able to calculate the coexistence in the critical region.


INTRODUCTION
Molecular simulation is rapidly becoming a very important tool for the prediction of thermodynamical properties of a wide variety of compounds and mixtures.In particular, thiophene compounds are nowadays of great interest due to the increasing requirements for sulphur removal from fuel being dibenzothiophene (DBT) and its derivatives [3] the most difficult compounds to remove from the streams.Hence, there is a strong need for thermodynamic data for these compounds for which only limited vapour-equilibrium data can be found in the literature.For example, experimental data of Vapour-Liquid Equilibrium (VLE) for mixtures containing thiophene + alkanes + CO 2 have been published recently [4].In addition, mixtures of ethanol+thiophene have also been considered [5].A much older study on systems consisting of aromatic hydrocarbons and thiophene in polar solvents also exists [6,7].The appearance of powerful algorithms like the Gibbs ensemble, introduced in 1987 by Panagiotopoulos [8] allowed the direct determination of phase coexistence without considering the simulation of the interface.
Nevertheless, studying the critical region of the different systems still constitutes a challenge due to the fact that the correlation length diverges when the critical point is approached.Having a good prediction for the critical point is crucial when equations of state or group contribution methods are applied to predict phase equilibria.For realistic potentials, the law of rectilinear diameters combined with the correct scaling relationship has commonly been used to extrapolate phase coexistence data obtained from molecular simulation to estimate critical points [9].Understanding the near-critical region behaviour can be particularly useful for the case of binary mixtures where the application of the law of rectilinear diameters is limited.Furthermore, the introduction of the Histogram Reweighting (HR) techniques [10,11] which allow the estimation of a joint Probability Distribution Function (PDF) has lead to several well-based studies and determinations of the critical point [12][13][14].The use of HR techniques allows the calculation of the PDF at any combination of intensive properties, a fact which becomes very interesting when these techniques are combined with Finite Size Scaling (FSS) studies.
FSS techniques allow the study of the influence of the size of the simulation box on the location of the critical point.For instance, in a recent work [15,16] the use of the fourth order cumulant calculation, also known as the Binder parameter, was combined with FSS studies.Although the use of these methodologies constitutes a significant progress in the location of critical points, less attention has been paid to the region close to the critical point.To give insight into this problem, recently Kim [1] has proposed a new algorithm able to yield precise estimates of the near critical coexisting liquid and gas densities ρ + and ρ − which he applied to the Restricted Primitive Model (RPM) and Hard Core Square Well (HCSW) systems.Due to the previously mentioned correlation length divergence appearing in the near-critical region, these values cannot be reliably estimated by using the conventional criteria of equal-area construction for the density equilibrium distribution function P(ρ; T ).For large enough system sizes and temperatures sufficiently far away from the critical temperature, a double peaked PDF is found providing accurate estimates for the liquid and vapour coexisting densities.On the other hand, when the critical point is approached, these two peaks become indistinguishable making the estimates for ρ + and ρ − unreliable.At this point, a study of the Binder parameter Q L reveals important information which is useful to determine accurate values for the coexisting vapour and liquid densities.In this work, we apply the biased algorithm introduced by Kim [1] for the case of the HCSW and the RPM, to the realistic intermolecular potential developed for thiophenes [17].By using this algorithm, we are able to determine the coexistence density jump Δρ ∞ ≡ ρ + − ρ − while we adopt the hypothesis of the law of rectilinear diameters to calculate ρ diam ≡ (ρ + + ρ − )/2 where it is important noting that both values of Δρ ∞ and ρ diam need to be determined in order to calculate the upper part of the phase diagram.In the next section we give details about the simulations and the model used in this work, while in Section 2 we provide the reader with a theoretical background about the Binder parameter, Q L .

MODEL AND SIMULATION DETAILS
To describe the dispersion interactions, thiophene is represented by a set of interacting Lennard-Jones sites for each CH or S group.The interactions between two different united atoms, i and j, from different molecules is calculated according to the 6-12 Lennard Jones potential: To calculate the parameters between unlike united atoms, we use the Lorentz-Berthelot combining rules: (2) To model the aromatic rings we have taken the AUA 4 CH intermolecular potential for polyaromatic compounds by Contreras et al. [18,19], and the S group for thiophenes by Pérez et al. [17].The parameters for the different potentials, which were obtained from fitting to selected equilibrium properties of aromatic and thiophenic compounds, are shown in Table 1.In Table 2 we show the bond angles and distances used for the different bonds.These parameters have been taken from the experimental geometry of thiophene [20] as well as from ab initio calculations at the Density Functional Theory level.The molecule of thiophene has been assumed to be rigid, as for the AUA 4 model for aromatics while no electrostatic charges have been included due to the moderate dipole moment of the thiophene molecule μ = 0.54 D [21].
Regarding the simulations, three different system sizes L = 20, 25, 28 Å have been investigated by means of five grand canonical simulations at different temperatures coupled with parallel tempering exchanges.The cutoff length was set always to half of the simulation box.Periodic boundary conditions were applied.The length of each simulation was at least 10 8 MC steps previously equilibrated for at least 2 × 10 6 MC steps.The frequencies for the different MC moves were set to 0.2 in the case of translations and rotations and 0.6 for insertions or deletions of the particles.

ANALYSIS OF THE BINDER PARAMETER (Q L )
As has been previously shown, the Q L parameter provides useful information when analysing the critical region.As is shown in Figure 1, the study of the Q L parameter along the equilibrium line defined as: where m = ρ − ρ is the order parameter and where it is important noting that in this case U L = Q −1 L , provides a direct route to estimating the critical point through the intersections of the curves plotted for different system sizes L. The U L parameter provides a dimensionless measure of the shape of the order parameter distribution function and is expected to have a universal value U L = 1.6035 at the critical point [22] for the Ising universality class.This means that the infinite system size critical point can be identified as the point where U L becomes system-size independent, being this point where cumulants from different sizes intersect.
Although this estimate is in principle a direct estimate of the infinite system size critical point, there are corrections to scaling which cause the intersections of different system sizes, and hence this estimate of the infinite system size critical point, to be also sensitive to the system size, albeit less so than an estimate taken from just one system size.These corrections to scaling can be used to scale various intersections from different system sizes with system size.
On the other hand, the Q L parameter investigated at constant temperature as a function of ρ, also embodies important information.As is shown in Figure 2, when the Q L parameter is plotted at isothermal conditions, two minima can be appreciated.These two minima give the finite size values of the coexisting vapour and liquid densities ρ ± L (T ).Below the critical temperature, and for system sizes going to infinity, where the finite-size effects are negligible, the values of ρ ± L (T ) for the minima approach the true ρ ± values while their respective heights Q − min and Q + min vanish.Under these conditions, the values obtained for ρ ± are fully consistent with the use of the equal-weight construction criteria [23].Close to the critical temperature, a finite-size scaling approach for Q L , the theoretical basis of which is briefly described here, is proposed by Kim.For further details, the reader is referred to [1] and [2].The following finite-size scaling expression is proposed for Q L : being Q the scaling function, Δρ = ρ − ρ c , t a proportionality constant, and β and ν critical exponents.According to this expression, the average of the minima heights and the average of the normalized density deviations follow respectively the scaling expressions: being M and N the appropriate universal scaling functions.Since both averages depend on the same scaling variable tL 1/ν , Δy min can be expressed as a universal function R(Q min ) of the average of the minima heights.As can be seen from Equation 7, the knowledge of the universal scaling function R(Q min ) leads to the determination of Δρ ∞ , thus the calculation of the critical exponent β as well as the determination of the critical temperature which can be identified as the point where Δρ ∞ vanishes.In order to determine R(Q min ), an algorithm has recently been proposed [2] based on a recursive numerical construction of the function taking as a starting point the exact expansion obtained for the case of the double-peaked Gaussian distribution of the density PDF found well below the critical temperature.The use of this algorithm requires generating very accurate data for a variety of system sizes, at least three, in the critical region.
To overcome this limitation, a second algorithm was proposed by assuming the universality class to which the systems belong to.This assumption is well justified for a wide variety of systems, such as the HCSW fluid, the Lennard-Jones (LJ) model and even the RPM system where the universality class is already known.This second algorithm uses as input the R(Q min ) previously determined by using the unbiased algorithm.This allows Δρ ∞ to be determined by considering only a unique system size thanks to the a priori knowledge of the universal scaling function R(Q min ), that can be presented now as a term which corrects the finitesize jump in density (ρ + min − ρ − min ) transforming it into the infinite system size difference between coexisting densities.
To determine the value of R(Q min ), it has been fitted [2] to the following expression defined in terms of the auxiliary variable q ≡ Q min ln(4/eQ min ): where a 2 = 1.829, a 3 = 1.955, b 2 = 2.340, b 3 = −1.388and q r = q/q c being q c the value at which R(Q min ) vanishes and consequently indicating the critical point.In the absence of corrections to scaling, the value of q c is expected to be universal.For the case of the HCSW and the RPM it was found that q c 0.286.

NEAR-CRITICAL COEXISTENCE APPROACH AND CONCLUSIONS
Since the main objective of this work is to determine the near-critical phase coexistence with a reasonable investment of computational time, in this section we show how to adapt the algorithm proposed by Kim and Fisher to the case of our realistic intermolecular potentials where the computational time is increased and furthermore we find corrections to scaling which are not found in the case of the HCSW fluid or the RPM model.Although the biased algorithm is able to provide the user with accurate estimates of Δρ ∞ , in order to close the phase diagram, the value of ρ diam is also necessary.A precise and complex algorithm has already been developed [1] based also on a numerical recursive construction taking as initial point the double-peaked Gaussian distribution found below the critical point, however in this case there is no associated universal function that can be developed and each system can have a different form.By applying this scaling algorithm one can notice that the applicability of the law of rectilinear diameters extends up to the region very close to the critical point.Thus, in order to close the phase diagram, we propose the combination of the assumptions of the law of rectilinear diameters with the predictions for Δρ ∞ obtained from the biased algorithm.To construct the upper part of the phase diagram, two points are defined to calculate the values of ρ diam (T ).First of all, the critical temperature and density are calculated from a FSS analysis of the intersections of the Binder parameter calculated using different system sizes with its universal value [24].Then, a second point for ρ diam (T ) can be fixed at a temperature sufficiently below the critical temperature such that the FSS effects are negligible.Then, by applying the biased algorithm mentioned in Section 2, we can calculate the scaled values of Δρ ∞ which allow us to calculate finally ρ ± L (T ).Due to the observed corrections to scaling appearing in the case of the Lennard-Jones fluid, we have observed in the limit of vanishing Δρ ∞ a value of q ≈ 0.31, slightly different from the expected universal value q ≈ 0.286 calculated by Kim.This can be appreciated in Figure 3 where, as in the case of U L in Figure 1, the values of q L have been found for different temperatures and system sizes.Although in this case the limited quality of the simulation data prevent an accurate estimate of the intersections being made, however, we can conclude that for the case of the Lennard-Jones thiophene system, the intersections occur at a significantly higher value than the one calculated by Kim.To determine from our simulation data, the universal value of q, and according to the already observed behaviour of the intersections of the fourth order cumulant curves with its universal value U L = 1.6035 [22], we have plotted the intersections of the q L curves with the apparent critical temperatures calculated from the U L analysis [24].These values, indicated in Figure 3 by means of squares, seem to take place at an approximate value q L ≈ 0.294 close to the value calculated by Kim.The discrepancy between these values can probably be explained taking into account the error that could be committed when applying the full numerical method as well as the error in the calculation of the q L curves of this work.After this analysis, for the case of the Lennard-Jones potentials used in our model, the behaviour observed suggests that we need to readapt the biased algorithm simply by replacing in the calculation of q r the critical universal value q c = 0.294 with the approximate value at which the intersections take place q ≈ 0.31.Finally, in Figure 4   minima.We can conclude that the use of this methodology allows a well-founded calculation of the upper part of the phase diagram without the necessity of investing additional computational time.

3 q
, we show how the phase diagram of thiophene is closed by applying the combination of HR techniques with a FSS study of the Q L L intersections at different system sizes.Dotted, dashed and solid line represent L = 20, 25, 28 Å systems respectively.

Figure 4 Thiophene
Figure 4Thiophene phase diagram: coexistence densities close to the critical point.