A set of molecular models based on quantum mechanical ab initio calculations and thermodynamic data
Pfaffenwaldring 9, 70550 Stuttgart, Germany
A parameterization strategy for molecular models on the basis of force fields is proposed, which allows a rapid development of models for small molecules by using results from quantum mechanical (QM) ab initio calculations and thermodynamic data. The geometry of the molecular models is specified according to the atom positions determined by QM energy minimization. The electrostatic interactions are modeled by reducing the electron density distribution to point dipoles and point quadrupoles located in the center of mass of the molecules. Dispersive and repulsive interactions are described by Lennard-Jones sites, for which the parameters are iteratively optimized to experimental vapor-liquid equilibrium (VLE) data, i.e. vapor pressure, saturated liquid density, and enthalpy of vaporization of the considered substance. The proposed modeling strategy was applied to a sample set of ten molecules from different substance classes. New molecular models are presented for iso-butane, cyclohexane, formaldehyde, dimethyl ether, sulfur dioxide, dimethyl sulfide, thiophene, hydrogen cyanide, acetonitrile, and nitromethane. Most of the models are able to describe the experimental VLE data with deviations of a few percent.
Keywords: Molecular modeling; modeling strategy; vapor-liquid equilibrium; critical properties; iso-butane; cyclohexane; formaldehyde; dimethyl ether; sulfur dioxide; dimethyl sulfide; thiophene; hydrogen cyanide; acetonitrile; nitromethane
Molecular modeling and simulation is a progressive approach for describing and predicting thermophysical properties of both pure substances and mixtures of technical interest. Many authors have shown the excellent capabilities of this molecular approach for different applications [1, 2]. Unfortunately, the more widespread use of molecular methods for engineering applications is still restricted by the scarcity of suitable molecular models.
For many substances, transferable molecular models have been developed, i.e. force fields for classes of compounds like alkanes or alcohols. Thereby, it is assumed that the parameters for functional groups are valid for different molecular species. The main disadvantage of transferable models is that their parameters were adjusted for the whole substance class and may not be optimal for a specific substance. Furthermore, they do not cover substances outside the modeled class. An overview and assessment of some commonly used transferable potentials can be found in .
Molecular models that were optimized for specific substances are available only for selected compounds and, particularly older ones, sometimes do not show the desired quality. Therefore, from an engineering point of view, a method for the rapid development of new molecular models for specific substances is of great interest. Moreover, the molecular models should be accurate, simple, and computationally efficient. In the present paper, a systematic strategy is proposed to develop such models.
The present work is restricted to rigid, non-polarizable models for comparatively small molecules. The models have state-independent parameters throughout. For computational efficiency, the united-atom approach is used, i.e. hydrogen atoms bonded to carbon are not modeled explicitly. The proposed strategy uses information determined by quantum mechanical (QM) ab initio calculations to include physically sound molecular properties and to reduce the number of adjustable parameters. A remaining subset of model parameters – typically two to four – is subsequently optimized by adjustment to experimental data on vapor-liquid equilibria (VLE) of the pure substances. The aim is to achieve deviations to experimental values for the vapor pressure, saturated liquid density, and enthalpy of vaporization in the range from triple point to critical point of below 5, 1, and 5%, respectively.
Such accurate models are known to have an excellent extrapolative and predictive power. This was recently shown, e.g., by the prediction of 17 different thermophysical properties for ethylene oxide, covering phase equilibria, thermal, caloric, transport properties, and surface tension , or for ammonia including structural quantities .
Molecular models that were developed on the basis of QM calculations stand between strictly ab initio models and fully empirical models. The present strategy is based on the idea to include substantial ab initio information without giving up the freedom to reasonably optimize the model to important macroscopic thermodynamic properties. Thus, for the modeling process some experimental data is needed for optimization. All three chosen properties, mentioned above, have the advantage to be well available for numerous engineering fluids and to represent dominant features of the fluid state.
The parameters of a molecular model can be separated into three groups. Firstly, geometric parameters specify the locations of the different interaction sites of the molecular model. Secondly, electrostatic parameters define the interactions of static polarities of the single molecules. And finally, dispersive and repulsive parameters determine the attraction by London forces and the repulsion by overlaps of the electronic orbitals. Here, the Lennard-Jones 12-6 (LJ) potential [6, 7] was used to assure straightforward compatibility with the overwhelming majority of the molecular models in the literature.
To describe the intermolecular interactions, a varying number of LJ sites and superimposed ideal point dipoles and/or ideal linear point quadrupoles were used. Point dipoles or quadrupoles were employed for the description of the electrostatic interactions to reduce the computational effort significantly. A point dipole may, e.g. when a simulation program does not support this interaction site type, be approximated by two point charges separated by a distance . Limited to small , one is free to choose this distance as long as holds. Analogously, a point quadrupole can be approximated by three collinear point charges , , and separated by each, where .
The total intermolecular interaction energy writes as
where , , are the distance, the LJ energy parameter, and the LJ size parameter, respectively, for the pair-wise interaction between LJ site on molecule and LJ site on molecule . is the permittivity of vacuum, whereas and denote the dipole moment and the quadrupole moment of the electrostatic interaction site on molecule , and so forth. are expressions for the dependency of the electrostatic interactions on the orientations and of the molecules and , cf. [8, 9]. Finally, the summation limits , , and denote the number of molecules, the number of LJ sites, and the number of electrostatic sites, respectively.
2 Molecular properties from QM
In a recent publication, Sandler et al.  gives a brief overview on the use of QM for the calculation of thermophysical properties. By numerically solving Schrödinger’s equation, it is nowadays possible to calculate different molecular properties for technically relevant components in a quite standardized way. Many different QM codes are available for this task. For license reasons, the open source code GAMESS(US)  was used in the present work.
All geometric data of the molecular models, i.e. bond lengths, angles, and dihedrals, were directly taken from QM calculations. Therefore, a geometry optimization, i.e. an energy minimization, was initially performed using GAMESS(US) . The Hartree-Fock level of theory was applied with a relatively small (6-31G) basis set. Alternatively, density functional theory (DFT) methods, e.g. BLY3P, can be used, as they are known to give reasonable results for the molecular structure .
The resulting configuration of the atoms was taken without subsequent modification to specify the position of the LJ sites in space, except for the hydrogen atoms. As the united atom approach was used to obtain computationally efficient molecular models, the hydrogen atoms were modeled together with the carbon atom they are bonded to. For the methylene (CH) and methyl (CH) united atom sites, the LJ potential was located at the geometric mean of the nuclei, while the methine (CH) united atom site was located at 0.4 of the distance between carbon and hydrogen atom, cf. Figure 1. These empirical offsets are in good agreement with the results of Ungerer et al. , which were found by optimization of transferable molecular models for n-alkanes.
Intermolecular electrostatic interactions mainly occur due to static polarities of single molecules that can well be obtained by QM. Here, the Møller-Plesset 2 level was used that considers electron correlation in combination with the polarizable 6-311G(d,p) basis set.
The purpose of the present work is the development of effective pair potentials with a state-independent set of model parameters. Obviously, the electrostatic interactions are stronger in the liquid state than in the gaseous state due to the higher density. Furthermore, the mutual polarization raises their magnitude in the liquid. Thus, for the calculation of the electrostatic moments by QM a liquid-like state should be considered. This was done here by placing the molecule within a dielectric continuum and assigning the experimental dielectric constant of the liquid to it, as in the COSMO method .
From the resulting electron density distribution for the small symmetric molecules regarded here, ideal point dipoles and ideal linear point quadrupoles were estimated by simple integration over the orbitals. Magnitudes and orientations of these electrostatic interaction sites were used in the molecular models without any modification.
For other, more complex molecules, more sophisticated methods like CHELP , CHELPG , or the distributed multipole analysis  are available in the literature. These methods adjust a set of partial charges or higher order electrostatic sites to the electrostatic potential around the molecule calculated by QM. Although they are able to reflect the electrostatic interactions with higher accuracy, they are not considered here. They always yield a larger number of interaction sites if they are not co-located with other sites and would thus lead to computationally more expensive molecular models.
2.3 Dispersion and Repulsion
It would be highly desirable to also calculate the dispersive and repulsive interactions using ab initio methods as well. This approach was followed by different authors in the past, e.g. for neon [20, 21, 22, 23], argon [21, 23, 24], krypton , nitrogen , carbon dioxide , hydrogen chloride , acetonitrile , methanol , acetylene , and methanethiol . However, from an engineering point of view, this leads to difficulties.
For an estimation of dispersive and repulsive interactions at least two molecules must be taken into account. To properly scan the energy hyper surface, many QM calculations at different distances and orientations of the molecules have to be performed. As the dispersive, and partly also the repulsive, interactions are just a very small fraction of the total energy calculated by QM, highly accurate methods like coupled cluster (CC) with large basis sets or even extrapolations to the basis set limit must be used for this task .
Due to the fact that this is computationally too expensive for engineering purposes, we used the parameters for the dispersive and repulsive interactions for an initial model from similar sites of other molecular models. Some of these parameters were subsequently fitted in the optimization process to yield the correct VLE behavior of the modeled substance.
3 Optimization to VLE data
The optimization was performed using a Newton scheme as proposed by Stoll . The applied method has many similarities with the one published by Ungerer et al.  and later on modified by Bourasseau et al. . It relies on a least-square minimization of a weighted fitness function that quantifies the deviations of simulation results from a given molecular model compared to experimental data. The weighted fitness function writes as
wherein the -dimensional vector represents the set of model parameters to be optimized. The deviations of results from simulation to experimental data are weighted with the expected simulation uncertainties . Equation (4) allows simultaneous adjustment of the model parameters to different thermophysical properties (saturated liquid densities , vapor pressures , and enthalpies of vaporization at various temperatures in the present work).
The unknown functional dependence of the property on the model parameters is approximated by a first order Taylor series developed in the vicinity of the initial parameter set
Therein, the partial derivatives of with respect to each model parameter , i.e. the sensitivities, are calculated from difference quotients
Assuming a sound choice of the model parameter variations , i.e. small enough to ensure linearity and large enough to yield differences in the simulation results significantly above the statistical uncertainties, this method allows a step-wise optimization of the molecular model by minimization of the fitness function . Experience shows that an optimized set of model parameters was usually found within a few iterative steps when starting from a reasonable initial model.
Correlations for vapor pressure, saturated liquid density, and enthalpy of vaporization, taken from the DIPPR database , were used as ”experimental data” for model adjustment and evaluation. This was done even in cases where the correlation is based on no or only few true experimental data points, as the correlations were regarded best practice. The comparison between simulation results and experiment was done by applying fits to the simulation data according to Lotfi et al. . The relative deviation between fit and correlation was calculated in steps of 1 K from 55 to 97% of the critical temperature and is denoted by ”mean unsigned error” in the following.
Vapor-liquid equilibrium simulations were performed using the Grand Equilibrium method by Vrabec et al. , technical simulation details are given in the appendix.
4 Molecular Models
The selected ten molecules belong to different substance classes to show the wide applicability of the proposed strategy. In the present work, we restricted ourselves to small molecules, where the internal degrees of freedom may be neglected. Thus, the molecular models are rigid, using the most stable configuration determined by QM.
The optimized parameter sets of the new molecular models are summarized in Table 1. Table 2 compares the critical properties from simulation to experimental data. The critical properties from simulation were obtained through fits to VLE simulation results as suggested by Lotfi et al. . The estimated uncertainties of critical temperature, critical density, and critical pressure from simulation are 0.5, 2, and 2%, respectively. A very good agreement between simulation and experiment was reached, being predominantly within the combined error bars.
In the following sections, substance specific details are briefly discussed. Furthermore, references to alternative models from the literature are given and the simulation results from the present models are compared to simulation data from the literature where available. Numerical VLE simulation results are given as Supplemtary Material.
4.1 Iso-Butane and Cyclohexane
The branched alkane iso-butane and the cyclic alkane cyclohexane show only very weak static polarities. Here, the main contributions to the intermolecular interaction are dispersion and repulsion. The electrostatic interactions have only a minor influence but should not be neglected completely.
In the literature different molecular models for iso-butane can be found, which are mostly based on force fields for branched alkanes. The well-known OPLS force field by Jorgensen et al. is available in two versions for iso-butane, one using the united-atom approach  and one using an all-atom description . Both OPLS force fields were optimized to liquid density and enthalpy of vaporization at 293 K only. The model of Poncela et al.  was adjusted to yield correct second virial coefficients. For a better applicability in a wider range of states, recent developments were optimized to experimental VLE data. Examples from this group are the force fields presented by Nath and de Pablo , Martin and Siepmann , Bourasseau et al. , or Chang and Sandler .
For the present model, four LJ sites, one for each methyl group and one for the methine group, were used to describe dispersion and repulsion of iso-butane. The polarity was modeled by a single (weak) dipole (0.1347 D) located in the center of mass. Orientation and magnitude of the dipole were taken from the QM calculation. For the initial model, the LJ parameters were taken from Ungerer et al. . It was sufficient to adjust a single parameter, the offset distance of the methine group. It was optimized to 0.4 of the carbon-hydrogen distance, cf. Figure 1 and held constant subsequently. The parameters of the present model are given in Table 1.
Figures 2 to 4 show saturated densities, vapor pressure, and enthalpy of vaporization, respectively, from the present iso-butane model in comparison to experimental data . Figure 5 shows a deviation plot between simulation and experimental data. In the deviation plot also simulation results from Martin and Siepmann, using their TraPPE force field,  and from Nath and de Pablo  are included. A very good agreement was obtained for the present model yielding mean unsigned errors in vapor pressure, saturated liquid density, and enthalpy of vaporization of 4.2, 0.6, and 1.8%, respectively, in the temperature range from 55 to 97% of the critical temperature, which is about 225 to 395 K. For vapor pressure, the present model yields significantly better results than the TraPPE force field, while no simulation data is available from Nath and de Pablo for this property. For saturated liquid density, all three models yield comparable results within 1% deviation to experimental data. No comparison between the models was possible for enthalpy of vaporization due to missing numerical data in [48, 49].
For cyclohexane different molecular models are available in the literature [52, 53, 54, 55] which all account for the internal degrees of freedom. Nevertheless, the present cyclohexane model was assumed to be rigid and in its most stable configuration, i.e. the saddle shape. The molecular model consists of six LJ sites, one for each methylene group. The static polarity was modeled by a single quadrupole parameterized according to QM. The two LJ parameters for methylene were optimized to experimental VLE data. The parameters of the present model are given in Table 1.
Figures 2 to 4 again show present VLE simulation results in comparison to experimental data for cyclohexane. Figure 6 shows the deviation plot including simulation results from Bourasseau et al. . The mean unsigned errors in vapor pressure, saturated liquid density, and enthalpy of vaporization for the present model are 0.9, 0.5, and 5.6%, respectively. Simulation results for the enthalpy of vaporization show systematic relative deviations towards the critical point, while the absolute deviation is below 2.5 kJ/mol.
Compared to the model of Bourasseau et al. improvements in the description of the saturated liquid density were achieved. The simulation results for vapor pressure agree within their assumed simulation uncertainties which were not reported by Bourasseau et al. An interesting point is that both molecular models yield the same deviations from experimental data on the enthalpy of vaporization, while the DIPPR database reports true experimental data up to . For all other molecular models for cyclohexane mentioned above, no numerical VLE simulation data is reported in the literature. Thus, no comparison is made here.
4.2 Formaldehyde and Dimethyl Ether
The present molecular model for formaldehyde consists of two LJ sites, one for the oxygen atom and one for the methylene group, as well as one dipole. The dipole is located in the center of mass and its moment was specified according to QM results, cf. Table 1. All four LJ parameters were adjusted to experimental VLE data and are given in Table 1.
Figures 7 to 8 compare simulation results and experimental VLE data for formaldehyde, Figure 9 shows the relative deviations. Mean unsigned errors in vapor pressure, saturated liquid density, and enthalpy of vaporization are 4.3, 0.9, and 8.4%, respectively.
Note that, in contrast to iso-butane or cyclohexane, the available experimental data base is very weak here. In fact, for vapor pressure only a single data set from the year 1935  is available. For saturated liquid density and enthalpy of vaporization, respectively, a single data point at 254 K is “accepted” by the DIPPR database . Thus, no further optimization of the molecular model was attempted although the desired quality seems not to be fully reached. Hermida-Raón and Ríos  published a molecular model based on QM calculations of formaldehyde dimers and trimers. They applied their model to liquid phase simulations but report no results on VLE properties. Thus, no comparison is presented here.
Dimethyl ether was modeled with three LJ sites in the present work, one for the oxygen atom and one for each methyl group. A dipole was located in the center of mass and oriented along the symmetry axis of the molecule, where the dipole moment was again taken from QM calculation. For an initial model, the same LJ parameters of the methyl groups were taken as for the iso-butane model, i.e. those by Ungerer et al. . An adjustment of the two LJ parameters of the oxygen site was sufficient to reach the desired quality. Alternative models for dimethyl ether are given in [58, 59, 60, 61].
The simulation results for the present model of dimethyl ether in comparison to the experimental VLE data are shown in Figures 2 to 4. Figure 10 shows the relative deviations between simulation and experiment also including simulation results from Stubbs et al.  and very recent results from Ketko and Potoff . For dimethyl ether a good experimental data base is available for optimization of the molecular model. The simulation data of the present model is in very good agreement with the correlations of experimental data. The mean unsigned errors in vapor pressure, saturated liquid density, and enthalpy of vaporization are 2.6, 0.4, and 1.0%, respectively.
For vapor pressure, both molecular models specifically adjusted to dimethyl ether, i.e. the present model and the model by Ketko and Potoff, yield better results than the transferable molecular model by Stubbs et al., while for saturated density all models perform similarly. For enthalpy of vaporization, the present model and the model by Ketko and Potoff also yield comparable results. Note that the parameters for the electrostatic interactions of the model by Ketko and Potoff were adjusted to experimental VLE data as well.
It can be summarized that the model by Ketko and Potoff  and the present dimethyl ether model are of similar quality and outperform the transferable model by Stubbs et al. . While Ketko and Potoff adjusted four LJ parameters and the point charge magnitudes for their electrostatic interactions, following the proposed modeling strategy an optimization of only two Lennard-Jones parameters was sufficient to reach the same quality.
4.3 Sulfur Dioxide, Dimethyl Sulfide, and Thiophene
For sulfur dioxide, a molecular model was published by Sokolić et al. [62, 63] which was optimized to total energy and pressure in the liquid state. It was recently reviewed by Ribeiro . Alternatively, the commercial force field COMPASS  reports parameters for sulfur dioxide.
For the present molecular model, the intermolecular interactions of sulfur dioxide were modeled with three LJ sites, i.e. one per atom, plus one dipole and one quadrupole. The electrostatic sites are located in the center of mass and parameterized according to the results of QM calculation. The four parameters of the LJ sites, i.e. , , , and , were adjusted to experimental VLE data. All parameters of the molecular model are given in Table 1.
Present simulation results for sulfur dioxide are compared to experimental VLE data in Figures 11 to 13. Figure 9 shows the relative deviations between simulation and experiment for the present model, while no numerical VLE simulation data for sulfur dioxide from other authors was available to us. Mean unsigned errors in vapor pressure, saturated liquid density, and enthalpy of vaporization are 4.0, 0.9, and 1.6%, respectively. The good experimental data base, with more than 60 individual experimental VLE data points, allows a thorough optimization of the molecular model to the desired quality.
Literature models for dimethyl sulfide are given in [66, 67, 68]. The present dimethyl sulfide model consists of three LJ sites, one for the sulfur atom and one for each methyl group. The electrostatic interactions are modeled by one dipole and two quadrupoles oriented perpendicularly to each other. This description of the electrostatics was chosen, as QM yields a charge distribution, which can not properly be described with a lower number of electrostatic sites. The LJ parameters of the methyl groups were assumed to be the same as for iso-butane and dimethyl oxide, while the parameters of the sulfur group were adjusted to experimental VLE data.
Figures 11 to 13 show the present simulation results for dimethyl sulfide in comparison to experimental VLE data. Figure 14 shows the relative deviations between simulation and experiment for the present model and for the model of Lubna et al. . Note that simulation results on enthalpy of vaporization were not included in . For the present model of dimethyl sulfide, mean unsigned errors in vapor pressure, saturated liquid density, and enthalpy of vaporization are 4.0, 0.7, and 3.8%, respectively. Particularly the vapor pressure is better described than by the model of Lubna et al.
Also thiophene was described in the present work by a rigid model in its most stable conformation, as for cyclohexane. Five LJ sites, one for each of the four methylene groups and one for the sulfur atom, as well as one dipole and one quadrupole were used. The electrostatic parameters were directly passed on from QM. A total of four LJ parameters, i.e. , , , and , were adjusted to experimental VLE data. The optimized parameters are given in Table 1. Alternative molecular models for thiophene can be found in the literature [68, 69, 70]
Figures 11 to 13 compare simulation results to experimental VLE data for thiophene, while the relative deviations are shown in Figure 15. In the deviation plot also simulation results from Lubna et al. , Juárez-Guerra et al. , and Pérez-Pellitero et al.  are included. Mean unsigned errors of the present model in vapor pressure, saturated liquid density, and enthalpy of vaporization are 3.8, 1.2, and 3.2%, respectively.
The thiophene model by Juárez-Guerra et al.  shows significant deviations in both vapor pressure and saturated liquid density while no simulation results for the enthalpy of vaporization were given by the authors. The anisotropic united atoms (AUA) potential by Pérez-Pellitero et al.  overpredicts the vapor pressure over the complete temperature range by up to 20%. Simulation results for saturated liquid density and enthalpy of vaporization are in very good agreement with the DIPPR correlation for low temperatures but give noticeably higher values than the DIPPR correlation for temperatures above . The TraPPE force field by Lubna et al.  yields vapor pressure results that agree very well with the DIPPR correlation and the simulation results within their scatter. Saturated densities are higher than the correlation towards the critical point for all four molecular models. For enthalpy of vaporization no numerical data were given by Lubna et al.
It should be noted that for thiophene experimental vapor pressure data is available for temperatures up to around , experimental saturated liquid densities and enthalpies of vaporization are only available up to approximately . Thus, an assessment of the different molecular models regarding density and enthalpy of vaporization above on the basis of the DIPPR correlations is questionable.
4.4 Hydrogen Cyanide, Acetonitrile, and Nitromethane
Hydrogen cyanide was modeled in the present work with two LJ sites, one for the methine group and one for the nitrogen atom. Electrostatic interactions were modeled by one dipole and one quadrupole oriented along the symmetry axis, where the parameters were passed on from QM. All four LJ parameters were adjusted to experimental VLE data. The optimized parameters can be found in Table 1. For hydrogen cyanide no other molecular models were found in the literature.
Figures 11 to 13 compare simulation results and experimental VLE data for hydrogen cyanide. Figure 9 shows the relative deviations. Unfortunately, the DIPPR database  contains no true experimental data on the enthalpy of vaporization for hydrogen cyanide. Consequently, in the optimization process, only minor attention was paid to the enthalpy of vaporization. Mean unsigned errors in vapor pressure, saturated liquid density, and enthalpy of vaporization are nominally 7.2, 1.0, and 12.2%, respectively.
Several molecular models for acetonitrile are available in the literature. Jorgensen and Briggs , Price et al. , Guárdia et al. , and Nikitin and Lyubartsev  present models that were optimized to the liquid density and enthalpy of vaporization at 293 K. Hloucha and Deiters  give a polarizable molecular model for simulations in the liquid state while Hloucha et al.  published a model that is based on ab initio calculations. Finally, Wick et al.  proposed an extension of their TraPPE force field that was optimized to VLE data to cover acetonitrile.
In the present work, acetonitrile was modeled using three LJ sites, one for the methyl group, one for the central carbon atom, and one for the nitrogen atom. The electrostatic interactions were modeled by one dipole and one quadrupole, located in the center of mass and parameterized strictly according to QM results. The parameters of the LJ sites of the methyl group and the nitrogen atom were adjusted to experimental VLE data. The LJ parameters of the central carbon atom were taken from unpublished work on carbon dioxide and excluded from optimization, only a very weak sensitivity of the VLE simulation results on these parameters was found. All parameters of the molecular model are given in Table 1.
Simulation results for acetonitrile are compared to experimental VLE data in Figures 7 to 8. The relative deviations between simulation and experiment are shown in Figure 16 next to results of the OPLS-UA force field by Jorgensen and Briggs  and the TraPPE force field by Wick et al., which were reported for both models in . Despite the good experimental data base, the desired quality was not achieved by the present optimization. Only a fair description of the experimental VLE was reached. Mean unsigned errors in vapor pressure, saturated liquid density, and enthalpy of vaporization are 19.7, 0.9, and 5.4%, respectively.
The large relative deviations in vapor pressure result from systematic underestimations at low temperatures. In this region also difficulties were encountered in the simulative calculation of the chemical potential in the liquid phase to determine the phase equilibrium. This was even the case when the more sophisticated gradual insertion method  was used, as described in the appendix in greater detail.
The OPLS-UA force field by Jorgensen and Briggs  significantly overestimates the vapor pressure of acetonitrile and shows deviations in the saturated liquid density up to 4%, cf. Figure 16. Simulation results for vapor pressure obtained with the TraPPE force field by Wick et al.  show a very good agreement with the DIPPR correlation, while the saturated density is slightly overpredicted towards the critical point. No comparison of the literature models for acetonitrile is possible regarding the enthalpy of vaporization due to the lack of numerical simulation results.
The present molecular model for nitromethane consists of four LJ sites, one for the methyl group, one for the nitrogen atom, and one for the two oxygen atoms each. The electrostatic interactions were modeled by one dipole and one quadrupole oriented along the symmetry axis, where the parameters were specified according to QM results. All six LJ parameters were adjusted to experimental VLE data. Alternative models are available in the literature [75, 76, 77, 78].
Figures 11 to 13 compare the present simulation results for nitromethane with experimental VLE data. Figure 17 shows the relative deviations obtained with the present model, the OPLS-AA force field by Price et al. , and the TraPPE force field by Wick et al. . For the present model of nitromethane, mean unsigned errors in vapor pressure, saturated liquid density, and enthalpy of vaporization are 18.7, 0.2, and 7.0%, respectively. Again, a systematic underprediction of the vapor pressure at low temperatures was found, leading to the high relative deviations, as for acetonitrile.
The OPLS-AA force field overpredicts the vapor pressure by about 80% and underpredicts the saturated liquid density by up to 6%. The TraPPE force field of Wick et al. that was optimized to experimental VLE data, yields very good results in vapor pressure for low temperatures while the result for the highest simulated temperature deviates by +34% from the DIPPR correlation. Regarding saturated liquid density, strong scatter and large statistical uncertainties of the simulation results are observed for both literature models.
A strategy was proposed for the rapid development of molecular models for engineering applications. The strategy relies on results from QM ab initio calculations to include physically sound molecular properties and to reduce the number of adjustable parameters. Dispersive and repulsive interactions were modeled by LJ sites. Thus, the LJ interaction sites were located according to atom positions obtained by QM energy minimization on Hartree-Fock level. For the parameterization of the electrostatic interactions, QM calculations were performed using the Møller-Plesset 2 level of theory and the COSMO method. The resulting electron density distribution was reduced to ideal point dipoles and ideal linear point quadrupoles located in the center of mass. The moments and orientations of the dipoles and quadrupoles were passed on to the molecular models without any modification.
A united-atom approach was used for methine, methylene, and methyl groups to reduce the total number of interaction sites. The parameters of the LJ sites were initially taken from similar sites of other molecular models and were subsequently optimized to reproduce experimental VLE data, i.e. vapor pressure, saturated liquid density, and enthalpy of vaporization. It was aimed to achieve deviations between simulation and experiment of below 5, 1, and 5% in vapor pressure, saturated liquid density, and enthalpy of vaporization, respectively.
The new modeling strategy was successfully applied to ten molecules from different substances classes, i.e. iso-butane, cyclohexane, formaldehyde, dimethyl ether, sulfur dioxide, dimethyl sulfide, thiophene, hydrogen cyanide, acetonitrile, and nitromethane. Simulation results for the different substances agree well with correlations of experimental data taken from the DIPPR database , with noticeable deviations in the vapor pressure at low temperatures for acetonitrile and nitromethane.
For the two elongated molecules acetonitrile and nitromethane, the reduction of the electrostatic interactions to sites located in the center of mass seems to be an over-simplification, as the optimization of the molecular models was not fully successful. Also an adjustment of the electrostatic parameters, while keeping the ratio of the polar moments constant (not reported here in detail), did not yield significant improvements in the quality of the molecular models. A further study, e.g. using two or more spatially distributed electrostatic sites, is beyond the scope of this paper.
For all other molecules a significant improvement compared to available molecular models from the literature was achieved. Through an optimization to substance specific experimental VLE data, a very good description of the phase equilibria was obtained. Furthermore, in case of dimethyl ether it was shown that with the proposed modeling strategy it was possible to reach the same model quality by optimization of just two model parameters compared to five optimized parameters for the model by Ketko and Potoff .
The authors gratefully acknowledge financial support by Deutsche Forschungsgemeinschaft, Schwerpunktprogramm 1155 ”Molecular Modeling and Simulation in Process Engineering”. The simulations were performed the national super computer NEC SX-8 at the High Performance Computing Center Stuttgart (HLRS) under the grant MMHBF and on the HP XC6000 super computer at the Steinbuch Centre for Computing, Karlsruhe under the grant MMSTP.
The authors want to thank Inga Utz for her modeling work on the sulfur and nitrogen containing compounds.
The Grand Equilibrium method  was used to calculate VLE data at seven to thirteen temperatures from 50 to 97% of the critical temperature during the optimization process.For the liquid, molecular dynamics simulations were performed in the isobaric-isothermal () ensemble using isokinetic velocity scaling  and Anderson’s barostat . There, the number of molecules was throughout and the time step was to fs depending on the molecular weight and the magnitude of the intermolecular interactions. The initial configuration was a face centered cubic lattice, the fluid was equilibrated over time steps with the first time steps in the canonical () ensemble. The production run time span was to time steps with a membrane mass of kg/m. Widom’s insertion method  was used to calculate the chemical potential by inserting up to test molecules every production time step.
In cases where Widom’s insertion method yielded large statistical uncertainties for the chemical potential, i.e. at high densities for strongly interacting molecules, Monte Carlo simulations were performed in the ensemble for the liquid. Then, the chemical potential was calculated by the gradual insertion method [82, 79]. The number of molecules was . Starting from a face centered cubic lattice, Monte Carlo cycles were performed for equilibration and for production, each cycle containing translation moves, rotation moves, and volume move. Every cycles, fluctuating state change moves, fluctuating particle translation/rotation moves, and biased particle translation/rotation moves were performed, to determine the chemical potential. These computationally demanding simulations yield the chemical potential in dense and strong interacting liquids with high accuracy, leading to reasonable uncertainties in the VLE.
For the corresponding vapor, Monte Carlo simulations in the pseudo- ensemble were performed. The simulation volume was adjusted to lead to an average number of molecules in the vapor phase. After initial Monte Carlo cycles, starting from a face centered cubic lattice, equilibration cycles in the pseudo- ensemble were performed. The length of the production run was cycles. One cycle is defined here to be a number of attempts to displace and rotate molecules equal to the actual number of molecules plus three insertion and three deletion attempts.
The cut-off radius was set to Å throughout and a center of mass cut-off scheme was employed. Lennard-Jones long-range interactions beyond the cut-off radius were corrected employing angle averaging as proposed by Lustig . Electrostatic interactions were approximated by a resulting molecular dipole and corrected using the reaction field method . Statistical uncertainties in the simulated values were estimated by a block averaging method .
-  Ungerer, P.; Nieto-Draghi, C.; Rousseau, B.; Ahunbay, G.; Lachet, V. J. Mol. Liq. 2007, 134, 71–89.
-  Ungerer, P.; Nieto-Draghi, C.; Lachet, V.; Wender, A.; Di Lella, A.; Boutin, A.; Rousseau, B.; Fuchs, A.H. Mol. Sim., 2007, 33, 15–30.
-  Martin, M.G. Fluid Phase Equilib., 2006, 248, 50–55.
-  Eckl, B.; Vrabec, J.; Hasse, H. Fluid Phase Equilib. 2008, to appear.
-  Eckl, B.; Vrabec, J.; Hasse, H. Mol. Phys. 2008, accepted.
-  Jones, J.E. Proc. Roy. Soc. 1924, 106A, 441–462.
-  Jones, J.E. Proc. Roy. Soc. 1924, 106A, 463–477.
-  Allen, M.P.; Tildesley, D.J. Computer simulations of liquids; Clarendon Press: Oxford, 1987.
-  Gray, C.G.; Gubbins, K.E. Theory of molecular fluids. 1. Fundamentals; Clarendon Press: Oxford, 1984.
-  Lorentz, H.A. Ann. d. Phys. 1881, 12, 127–136.
-  Berthelot, D. Compt. Rend. Ac. Sc. 1889, 126, 1703–1706.
-  Sandler, S.I.; Castier, M. Pure Appl. Chem. 2007, 79, 1345–1359.
-  Schmidt, M.W.; Baldridge, K.K.; Boatz, J.A.; Elbert, S.T.; Gordon, M.S.; Jensen, J.H.; Koseki, S.; Matsunaga, N.; Nguyen, K.A.; Shujun, S.; Windus, T.L.; Dupuis, M.; Montgomery, A.M. J. Comput. Chem. 1993, 14, 1347–1363.
-  Leach, A.R. Molecular Modeling. Principles and Application; 2nd ed.; Prentice-Hall: Englewood Cliffs, 2001.
-  Ungerer, P.; Beauvais, C.; Delhommelle, J.; Boutin, A.; Rousseau, B.; Fuchs, A.H. J. Chem. Phys. 2000, 112, 5499–5510.
-  Klamt, A. J. Phys. Chem. 1995, 99, 2224–2235.
-  Chirlian, L.E.; Francl, M.M. J. Comput. Chem. 1987, 8, 894–905.
-  Breneman, C.M.; Wiberg, K.B. J. Comput. Chem. 1990, 11, 361–373.
-  Stone, A.J.; Alderton, M. Mol. Phys. 2002, 100, 221–233.
-  Eggenberger, R.; Gerber, S.; Huber, H.; Welker, M. Mol. Phys. 1994, 82, 689–699.
-  Vogt, P.S.; Liapine, R.; Kirchner, B.; Dyson, A.J.; Huber, H.; Marcelli, G.; Sadus, R.J. Phys. Chem. Chem. Phys. 2001, 3, 1297–1302.
-  Garrison, S.L.; Sandler, S.I. J. Chem. Phys. 2002, 117, 10571–10580.
-  Nasrabad, A.E.; Laghaei, R.; Deiters, U.K. J. Chem. Phys. 2004, 121, 6423–6434.
-  Ermakova, E.; Solca, J.; Huber, H.; Welker, M. J. Chem. Phys. 1995, 102, 4942–4951.
-  Nasraba, A.E.; Deiters, U.K J. Chem. Phys. 2003, 119, 947–952.
-  Leonhar, K.; Deiters, U.K. Mol. Phys. 2002, 100, 2571–2585.
-  Welker, M.; Steinebrunner, G.; Solca, J.; Huber, H. Chem. Phys. 1996, 213, 253–261.
-  Naicker, P.K.; Sum, A.K.; Sandler, S.I. J. Chem. Phys. 2003, 118, 4086–4093.
-  Hloucha, M.; Sum, A.K.; Sandler, S.I. J. Chem. Phys. 2000, 113, 5401–5406.
-  Garrison, S.L.; Sandler, S.I. J. Phys. Chem. 2004, 108, 18972–18979.
-  Garrison, S.L.; Sandler, S.I. J. Chem. Phys. 2005, 123, 054506.
-  Stoll, J. Molecular Models for the Prediction of Thermophysical Properties of Pure Fluids and Mixtures; Fortschritt-Berichte VDI, Reihe 3, 836; VDI-Verlag: Düsseldorf, 2005.
-  Ungerer, P.; Boutin, A.; Fuchs, A.H. Mol. Phys. 1999, 97, 523–539.
-  Bourasseau, E.; Haboudou, M.; Boutin, A.; Fuchs, A.H.; Ungerer, P. J. Chem. Phys. 2003, 118, 3020–3034.
-  Rowley, R.L.; Wilding, W.V.; Oscarson, J.L.; Yang, Y.; Zundel, N.A.; Daubert, T.E.; Danner, R.P. DIPPR® Data Compilation of Pure Compound Properties; Design Institute for Physical Properties, AIChE: New York, 2006.
-  Lotfi, A.; Vrabec, J.; Fischer, J. Mol. Phys. 1992, 76, 1319–1333.
-  Vrabec, J.; Hasse, H. Mol. Phys. 2002, 100, 3375–3383.
-  Daubert, T.E. J. Eng. Chem. Data 1996, 41, 365–372.
-  Reid, R.C.; Prausnitz, J.M.; Sherwood, T.K. The Properties of Gases and Liquids; 3rd ed.; McGraw-Hill: New York, 1977.
-  Kudchadker, A.P.; Ambrose, A.; Tsonopoulos, C. J. Chem. Eng. Data 2001, 46, 457–479.
-  Mathews, J.F. Chem. Rev. 1972, 72, 71–100.
-  Tsonopoulos, C.; Ambrose, D. J. Chem. Eng. Data 2001, 46, 480–485.
-  Marsh, K.N.; Young, C.L.; Morton, D.W.; Ambrose, D.; Tsonopoulos, C. J. Chem. Eng. Data 2006, 51, 305–314.
-  Kudchadker, A.P.; Alani, G.H.; Zwolinski, B.J. Chem. Rev. 1968, 68, 659–735.
-  Jorgensen, W.L.; Madura, J.D.; Swenson, C.J. J. Am. Chem. Soc. 1984, 106, 6638–6646.
-  Jorgensen, W.L.; Maxwell, D.S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225–11236.
-  Poncela, A.; Rubio, A.M.; Freire, J.J. Mol. Phys. 1997, 91, 189–201.
-  Nath, S.K.; de Pablo, J.J. Mol. Phys. 2000, 98, 231–238.
-  Martin, M.G.; Siepmann, J.I. J. Phys. Chem. B 1999, 103, 4508–4517.
-  Bourasseau, E.; Ungerer, P.; Boutin, A.; Fuchs, A.H. Mol. Sim. 2002, 28, 317–336.
-  Chang, J.; Sandler, S.I. J. Chem. Phys. 2004, 121, 7474–7483.
-  Errington, J.R.; Panagiotopoulos, A.Z. J. Chem. Phys. 1999, 111, 9731–9738.
-  Neubauer, B.; Boutin, A.; Tavitan, B.; Fuchs, A.H. Mol. Phys. 1999, 97, 769–776.
-  Faller, R.; Schmitz, H.; Biermann, O.; Müller-Plathe, F. J. Comput. Chem. 1999, 20, 1009–1017.
-  Bourasseau, E.; Ungerer, P.; Boutin, A. J. Phys. Chem. B 2002, 106, 5483–5491.
-  Spence, R.; Wild, W. J. Chem. Soc. 1935, 1, 506–509.
-  Hermida-Raón, J.M.; Ríos, M.A. J. Phys. Chem. A 1998, 102, 10818–10827.
-  Jorgensen, W.L. J. Am. Chem. Soc. 1981, 103, 335–340.
-  Lin, B.; Halley, J.W. J. Phys. Chem. 1995, 99, 16474–16478.
-  Stubbs, J.M.; Potoff, J.J.; Siepmann, J.I. J. Phys. Chem. B 2004, 108, 17596–17605.
-  Ketko, M.B.H.; Potoff, J.J. Mol. Sim. 2007, 33, 769–776.
-  Sokolić, F.; Guissani, Y.; Guillot, B. Mol. Phys. 1985, 56, 239–253.
-  Sokolić, F.; Guissani, Y.; Guillot, B. J. Phys. Chem. 1985, 89, 3023–3026.
-  Ribeiro, M.C.C. J. Phys. Chem. B 2006, 110, 8789–8797.
-  Yang, J.; Ren, Y.; Tian, A. J. Phys. Chem. B 2000, 104, 4951–4957.
-  Jorgensen, W.L. J. Phys. Chem. 1986, 90, 6379–6388.
-  Delhommelle, J.; Tschirwitz, C.; Ungerer, P.; Granucci, G.; Millié, P.; Pattou, D.; Fuchs, A.H. J. Phys. Chem. B 2000, 104, 4745–4753.
-  Lubna, N.; Kamath, G.; Potoff, J.J.; Rai, N.; Siepmann, J.I. J. Phys. Chem. B 2005, 109, 24100–24107.
-  Juárez-Guerra, F.; Rivera, J.L.; Zúñiga-Moreno, A.; Galicia-Luna, L.A.; Rico, J.L.; Lara, J. Sep. Sci. Technol. 2006, 41, 261–281.
-  Pérez-Pellitero, J.; Ungerer, P.; Mackie, A.D. J. Phys. Chem. B 2007, 111, 4460–4466.
-  Jorgensen, W.L.; Briggs, J.M. Mol. Phys. 1988, 63, 547–558.
-  Guàrdia, E.; Pinzón, R.; Casulleras, J.; Orozco, M.; Luque, F.J. Mol. Sim. 2001, 26, 287–306.
-  Nikitin, A.M.; Lyubartsev, A.P. J. Comput. Chem. 2007, 28, 2020–2026.
-  Hloucha, M.; Deiters, U.K. Mol. Phys. 1997, 90, 593–597.
-  Alper, H.E.; Abu-Awwad, F.; Politzer, P. J. Phys. Chem. B 1999, 103, 9738–9742.
-  Sorescu, D.C.; Rice, B.M.; Thompson, D.L. J. Phys. Chem. A 2001, 105, 9336–9346.
-  Price, M.L.P.; Ostrovsky, D.; Jorgensen, W.L. J. Comput. Chem. 2001, 22, 1340–1352.
-  Wick, C.D.; Stubbs, J.M.; Rai, N.; Siepmann, J.I. J. Phys. Chem. B 2005, 109, 18974–18982.
-  Vrabec, J.; Kettler, M.; Hasse, H. Chem. Phys. Lett. 2002, 356, 431–436.
-  Andersen, H.C. J. Chem. Phys. 1980, 72, 2384–2393.
-  Widom, B. J. Chem. Phys. 1963, 39, 2808–2812.
-  Nezbeda, I.; Kolafa, J. Mol. Sim. 1991, 5, 391–403.
-  Lustig, R. Mol. Phys. 1988, 65, 175–179.
-  Flyvbjerg, H.; Petersen, H.G. J. Chem. Phys. 1989, 91, 461–466.
|continued on next page|
|continued from previous page|
|continued on next page|
|continued from previous page|
List of Figures
- 1 Geometry of the united-atom sites: (a) The methine (CH) site is located at 0.4 of the carbon-hydrogen distance, (b) the methylene (CH) site is located at the geometric mean, and (c) the methyl (CH) site is located at the geometric mean.
- 2 Saturated densities. Present simulation data: iso-butane, cyclohexane, dimethyl ether. — correlations of experimental data , + critical points derived from simulation, experimental critical points.
- 3 Vapor pressure. Present simulation data: iso-butane, cyclohexane, dimethyl ether. — correlations of experimental data , + critical points derived from simulation, experimental critical points.
- 4 Enthalpy of vaporization. Present simulation data: iso-butane, cyclohexane, dimethyl ether. — correlations of experimental data .
- 5 Relative deviations of vapor-liquid equilibrium properties between simulation data and correlations of experimental data  () for iso-butane: present model, Martin and Siepmann , Nath and de Pablo . Top: vapor pressure, center: saturated liquid density, bottom: enthalpy of vaporization.
- 6 Relative deviations of vapor-liquid equilibrium properties between simulation data and correlations of experimental data  () for cyclohexane: present model, Bourasseau et al. . Top: vapor pressure, center: saturated liquid density, bottom: enthalpy of vaporization.
- 7 Saturated densities. Present simulation data: formaldehyde, hydrogen cyanide, acetonitrile, nitromethane. — correlations of experimental data , + critical points derived from simulation, experimental critical points.
- 8 Enthalpy of vaporization. Present simulation data: formaldehyde, hydrogen cyanide, acetonitrile, nitromethane. — correlations of experimental data .
- 9 Relative deviations of vapor-liquid equilibrium properties between simulation data and correlations of experimental data  () for present models: formaldehyde, hydrogen cyanide, sulfur dioxide. Top: vapor pressure, center: saturated liquid density, bottom: enthalpy of vaporization.
- 10 Relative deviations of vapor-liquid equilibrium properties between simulation data and correlations of experimental data  () for dimethyl ether: present model, Stubbs et al. , Ketko and Potoff . Top: vapor pressure, center: saturated liquid density, bottom: enthalpy of vaporization.
- 11 Saturated densities. Present simulation data: sulfur dioxide, dimethyl sulfide, thiophene. — correlations of experimental data , + critical points derived from simulation, experimental critical points, cf. Table 2.
- 12 Vapor pressure. Present simulation data: sulfur dioxide, dimethyl sulfide, thiophene. — correlations of experimental data , + critical points derived from simulation, experimental critical points.
- 13 Enthalpy of vaporization. Present simulation data: sulfur dioxide, dimethyl sulfide, thiophene. — correlations of experimental data .
- 14 Relative deviations of vapor-liquid equilibrium properties between simulation data and correlations of experimental data  () for dimethyl sulfide: present model, Lubna et al. . Top: vapor pressure, center: saturated liquid density, bottom: enthalpy of vaporization.
- 15 Relative deviations of vapor-liquid equilibrium properties between simulation data and correlations of experimental data  () for thiophene: present model, Lubna et al. , Juárez-Guerra et al. , Pérez-Pellitero et al. . Top: vapor pressure, center: saturated liquid density, bottom: enthalpy of vaporization.
- 16 Relative deviations of vapor-liquid equilibrium properties between simulation data and correlations of experimental data  () for acetonitrile: present model, Wick et al. , Jorgensen and Briggs [71, 78]. Top: vapor pressure, center: saturated liquid density, bottom: enthalpy of vaporization.
- 17 Relative deviations of vapor-liquid equilibrium properties between simulation data and correlations of experimental data  () for nitromethane: present model, Wick et al. , Price et al. . Top: vapor pressure, center: saturated liquid density, bottom: enthalpy of vaporization.
|continued on next page|
|continued from previous page|
|continued on next page|
|continued from previous page|