Emily
Gross‡
a,
Mark D.
Driver‡
b,
Areesha
Saif
b,
Oliver N.
Evans
b and
Christopher A.
Hunter
*b
aDepartment of Chemistry, University of Regensburg, Universitaetsstrasse 31, 93053 Regensburg, Germany
bYusuf Hamied Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge CB2 1EW, UK. E-mail: herchelsmith.orgchem@ch.cam.ac.uk
First published on 9th April 2025
The surface site interaction model for liquids at equilibrium (SSIMPLE) is a method for calculating thermodynamic properties in a fluid phase based on the use of surface site interaction points (SSIP) to represent all of the non-covalent interactions that molecules make with the environment. Interactions between the SSIPs of two different molecules are governed by a non-polar term and a polar term. Here the formulation originally made for room temperature liquids is generalized to any temperature. We show that the non-polar interaction term is temperature independent while the polar interaction term depends on temperature. This formulation was used to develop a description of the temperature dependence of fluid phase density in terms of an expansion energy, which is based on net intermolecular SSIP interactions. The method is shown to accurately model the temperature dependence of experimentally measured association constants for the formation of 1:
1 H-bonded complexes in carbon tetrachloride. The atomic interaction point (AIP) version of the SSIP descripiton of 171 different compounds was used in SSIMPLE to calculate room temperature liquid densities that are in good agreement with experimental data. Since non-covalent interactions in the vapour phase can be treated in the same way as liquid phase interactions, SSIMPLE can also be used to calcuate vapour–liquid equilibria (VLE). Experimental VLE data for 196 binary mixtures of 30 miscible compounds was collected, and SSIMPLE was shown to reproduce the experimental behaviour well.
Empirical approaches are either based on the creation of equations of state (EoS) or on quantitative structure property relationships (QSPR). EoS approaches have been in development for nearly two centuries, since Clapeyron first proposed the ideal gas law,7 and the development of EoS for the description of real gases has been an area of active research since the first equation by van der Waals.8–11 QSPR approaches use the correlation of a series of molecular descriptors to physical properties. The UNIFAC approach is a group contribution method, which was first applied to the generation of relationships to predict VLE and liquid–liquid equilibria (LLE) in 1977.12 Iterative improvements in the descriptions of the fragments have led to increased accuracy in the predictions.13–18
First principles methods based on quantum mechanics have also been used to develop models for predicting VLE. The sigma-profile description of molecules generated by COSMO19 has been used in the creation of semi-empirical relationships to predict VLE and phase diagrams. COSMO-RS has been used to calculate activity coefficients,2,20 and free energies in pure and mixed solvents as a function of temperature.6,21–28 Molecular simulation of liquid–vapour systems can be used to predict phase properties through propagation of the components through phase space. Monte Carlo (MC),29–32 molecular dynamics (MD),33–36 and dissipative particle dynamics (DPD)37–40 have all been used to predict phase compositions. Statistical associating fluid theory (SAFT) is a methodology for including the effects of association into a given theory and develop EoS to describe VLE and LLE.41–43 SAFT extends the thermodynamic perturbation theory developed by Wertheim44,45 to incorporate the treatment of mixtures.31,39–41 The SAFT approach can be used in MC simulations39 and MD simulations31,32 and has been applied to alkanes and polymer mixtures.46–53
The surface site interaction model for liquids at equilibrium (SSIMPLE) provides solvation energies based on the description of molecules as a set of surface site interaction points (SSIP) calculated using DFT.54 The dimensions of an SSIP are defined based on the 0.002 e bohr−3 electron density isosurface of a water molecule, assuming that water should be represented by one SSIP for each of the four H-bonding sites. Each SSIP is assigned an interaction parameter, εi, which is equivalent to the experimentally measured H-bond donor parameter (α) for sites with positive electrostatic potential or the H-bond acceptor parameter (−β) for sites with negative electrostatic potential.55 The H-bond scales were defined such that the product -αβ (or εiεj) corresponds to a free energy change for the interaction of two SSIPs at room temperature in units of kJ mol−1. Empirical relationships between experimentally determined H-bond parameters and the calculated molecular electrostatic potential surface (MEPS) can be used to obtain an SSIP description of any molecule as described previously.56
To describe a liquid phase in SSIMPLE, interactions between SSIPs are treated in a pairwise manner. The association constant for the interaction between SSIP i and SSIP j, Kij, is given by eqn (1).
![]() | (1) |
The interaction energy is made up of a polar term, Epolarij, which is equal to the product εiεj, and a non-polar term, EvdWij, which is the energy of the van der Waals interaction between two SSIPs, previously determined to be −5.6 kJ mol−1 based on VLE data for non-polar molecules.54 van der Waals interactions between non-polar molecules are a linear function of surface area,57 so by choosing a description that gives all SSIPs the same surface area footprint, a constant value can be used for EvdWij. In the SSIP approach, the difference between different atom types is described by differences in the number of SSIPs rather than a variable van der Waals interaction parameter.54 The liquid phase is treated as a Boltzmann ensemble of all possible SSIP interactions, and the matrix of Kij values calculated using eqn (1) can be used to determine the speciation of intermolecular contacts for any mixture of solvents and solutes. The free energy of solvation of SSIP i, ΔGS(i), calculated using eqn (2) corresponds to the free energy of transfer of SSIP i from a dilute vapour where there are no SSIP interactions.
![]() | (2) |
The first term in eqn (2) describes the interactions made by SSIP i with the other SSIPs in the liquid. The second term in eqn (2) corrects for the increased probability of interaction between SSIPs when they are confined to a condensed phase.54 Summation of ΔGS(i) over all SSIPs in a molecule yields the free energy of solvation of the molecule, ΔGS. This approach has been used to successfully predict phase transfer free energies for a wide range of different solutes between two liquid phases.58
SSIMPLE requires the concentrations of all species to be known and is only applicable at room temperature, because the interaction parameters, εi, were empirically parameterised using measurements of the free energy changes for formation of 1:
1 complexes at 298 K. Here, we describe a generalisation of the SSIMPLE method to include the temperature dependence of SSIP interactions and vapour–liquid equilibria.
![]() | ||
Fig. 1 Temperature dependence of ΔH* for argon (black), methane (red) and water (blue). ΔT is the increase in temperature relative to the melting temperature. |
Fig. 2 provides an interpretation of these observations. The bound state is represented by two different species. Molecular rotation within the bound state can cause the loss of the directional polar interactions formed in the H-bonded state, whilst maintaining van der Waals interactions in the non-H-bonded state. Molecular translation breaks all interactions, leading to the free state. Non-polar contacts are temperature-independent, because these interactions are not affected by temperature-induced rotational motion, but the contribution from polar interactions depends on the relative populations of the non-H-bonded and H-bonded states, which are sensitive to temperature. If the relative populations of the two bound states in Fig. 2 are determined by a Boltzmann distribution, the temperature-dependence of the polar interaction term in eqn (1), Epolarij, is given by eqn (3).
![]() | (3) |
Since the value of Epolarij(T) is known at 298 K for any pairwise SSIP interaction, eqn (3) can be solved to obtain the value of Epolarij(0). In practice, a lookup table of Epolarij(0) values is more straightforward to implement.
To assess the validity of this temperature-dependent formulation of the pairwise SSIP interaction energy, we compare calculated equilibrium constants with experimental data measured at different temperatures for the formation of 1:
1 complexes that make a single H-bond, i.e. a single SSIP contact.61–94 The ratios of the association constants experimentally measured in carbon tetrachloride at 293 K and 323 K were obtained for complexes of different solutes with known experimental values of the H-bond parameters, α and β. The H-bond parameters were used as the SSIP values, εi and εj, for the two solutes, and eqn (3) was used in SSIMPLE to calculate the equilibrium constants for the solute-solute interactions, as described previously.95 The calculated and experimental results are compared in Fig. 3. Although there is considerable scatter, the correlation in Fig. 3 validates the proposed formulation for the temperature dependence of SSIP interactions.
![]() | ||
Fig. 3 Comparison of experimental measurements of the ratio of association constants at two different temperatures with the corresponding values calculated with SSIMPLE and eqn (3) for 169 1![]() ![]() |
![]() | (4) |
![]() | (5) |
The zero point number density, ρ0, corresponds to the most concentrated state, which in practice is never reached because the liquid freezes. The lowest possible number density, ρc, occurs at the critical point, when vapour and liquid number densities are both equal (eqn (6)).57
![]() | (6) |
Thermal expansion results in the incorporation of additional void volume, such that the total number density of the liquid and vapour phases decreases with increasing temperature. The expansion energy which governs this process, Eexp, is the average energy required to break interactions between molecules. The expansion energy can be decomposed into contributions from non-polar van der Waals interactions, EvdWexp, and polar interactions, Epolarexp, in the same way as the pairwise SSIP interactions EvdWij and Epolarij that define Kij in eqn (1).
We can use non-polar molecules to examine the relationship between the expansion energy and the total van der Waals interaction energy available to a molecule in the condensed phase. At the critical temperature, Tc, the liquid and vapour number densities are both equal to ½ ρc. Substitution of eqn (5) and (6) into eqn (4) provides a relationship between the expansion energy and the critical temperature (eqn (7)).
![]() | (7) |
Values of Eexp were obtained for 190 alkanes by using the experimental values of Tc and calculated values of VvdW in eqn (7) (see ESI† for details). Fig. 4 shows the relationship between Eexp and the number of SSIPs per molecule, N, which was calculated from the surface area of the 0.002 e bohr−3 electron density isosurface. Although there is some variation, eqn (8), plotted in black in Fig. 4, provides a reasonable description of the data.
![]() | (8) |
![]() | ||
Fig. 4 Relationship between the expansion energy Eexp and the number of SSIPs N for 190 alkanes. The black line is eqn (8). |
If all SSIPs were paired, then there would be N/2 interactions, and the total van der Waals interaction energy would be ½ NEvdWij. The energy barrier to expansion for alkanes is half of this total, but reduced by an additional factor of (1 + N/12), such that there is an upper limit of −3EvdWij on the value of Eexp for large N (the asymptote in Fig. 4).
In order to determine the contribution that polar interactions make to the expansion energy, the experimental temperature dependence of the total number density of water in the liquid and vapour phases was examined. An apparent value of the expansion energy, Eappexp, can be computed using the experimental values of ρl and ρv in eqn (4) (eqn (9)).
![]() | (9) |
Fig. 5 shows that Eappexp changes significantly with temperature for water, i.e. water does not obey the law of rectilinear diameter. The value of Eappexp increases dramatically at low temperatures, where H-bonded states are more highly populated. The more H-bonding interactions are made, the greater the energy barrier to pulling molecules apart in the expansion process.
![]() | ||
Fig. 5 Relationship between apparent expansion energy Eappexp calculated using eqn (9) and temperature for water. |
It is possible to calculate the extent of H-bond formation in the liquid or vapour phase of water as a function of temperature using SSIMPLE. However, the equilibrium constant for pairwise interactions between SSIPs involves contributions from both van der Waals interactions and polar interactions (Eqn (1)). In order to determine the population of the bound state that is H-bonded in Fig. 2, the total population of the bound state calculated using SSIMPLE must be compared with the population of the bound state that is non-H-bonded, which can be calculated using just the van der Waals term in eqn (1) (i.e. by setting Epolarij = 0). For SSIPs that only make van der Waals interactions, the probability that a H-bonded or a non-H-bonded state is populated would be equal because the energy difference is zero. The extent to which polar interactions perturb these populations can therefore be determined using eqn (10) to calculate the fraction of the bound state that is H-bonded in the liquid and vapour phases, ϕlb and ϕvb.
![]() | (10) |
![]() | (11) |
Since SSIPs that make only van der Waals interactions can be considered equivalent, the speciation of SSIP contacts for a phase containing only non-polar molecules is described by eqn (12), which can be solved to obtain an analytical expression for ψvdWf (eqn (13)).
[i] = [if] + 2KvdW[if]2 | (12) |
![]() | (13) |
It is possible to use the temperature dependent formulation of Kij in eqn (1) in conjunction with the experimental number densities of the liquid and vapour phases of water to calculate the speciation in both phases as a function of temperature using SSIMPLE. The fraction of bound states that are H-bonded in the liquid and vapour phases, ϕlb and ϕvb, from eqn (10) can then be used to obtain the total fraction of bound states that are H-bonded for the vapour–liquid system, ϕb, using eqn (14).
![]() | (14) |
Fig. 6 shows that the value of Eappexp calculated from the experimental data for water is closely related to the fraction of the bound state that is H-bonded, ϕb, calculated using SSIMPLE. The line of best fit through the data obtained above 313 K has a correlation coefficient of R2 = 0.9996. The four lower temperature datapoints in Fig. 6 deviate from this line, presumably because there is a change in the structure of the liquid as it approaches the freezing point due to the geometric constraints imposed by highly H-bonded networks i.e. there is a decrease in density associated with formation of ice-like structures.
The line of best fit in Fig. 6 has a slope of 5.3 kJ mol−1 and an intercept of 2.1 kJ mol−1. These values can be compared with the value of EvdWexp calculated using eqn (8) for a molecule with four SSIPs that makes only van der Waals interactions (4.2 kJ mol−1) and with the total polar interaction energy for a water molecule at the zero point (twice Epolarij(0) for a water–water H-bond, which is −25.6 kJ mol−1). Some of the values of Eappexp are two to three times larger than these energies. In order to obtain an expression for Eappexp that is consistent for both polar and non-polar molecules, it is required that if polar interactions were to make no contribution to the expansion energy of water, then ϕb would be 0.5 and Eappexp would be equal to EvdWexp, i.e. 4.2 kJ mol−1. Noting that the intercept in Fig. 6 is exactly half of 4.2 kJ mol−1, these conditions are satisfied by eqn (15).
![]() | (15) |
The value of Epolarexp can therefore be determined from the slope of the line of best fit in Fig. 6, as 6.4 kJ mol−1, which is a quarter of the total polar interaction energy for water at the zero point, i.e. the polar interaction energy per SSIP at the zero point. These results suggest that a suitable formulation of the expansion energy involves the sum of the van der Waals term given by eqn (8) and the polar interaction energy per SSIP at the zero point calculated using eqn (16).
![]() | (16) |
In order to calculate Epolarexp, we therefore require a method to determine which SSIPs interact at the zero point. The hierarchical SSIP pairing strategy described previously for cocrystal prediction is used.98 In this approach, all SSIP interactions are assumed to be independent and free to find the best partner, such that there are no steric constraints on contacts, no packing effects, and no cooperativity between sites that are close in space on the surface of a molecule. This means that the best H-bond donor pairs with the best H-bond acceptor, and the second best H-bond donor with the next best H-bond acceptor etc. Any remaining SSIPs of the same sign that have not been paired are assumed to not form polar interactions, and therefore do not contribute.
Eqn (15) shows that the contribution of polar interactions to the expansion energy is weighted by the fraction of bound states that are H-bonded, i.e the fraction of polar interactions that are actually made at any given temperature determines how much polar interaction is available to oppose expansion. The denominator in eqn (15) can be rationalised by reformulating eqn (4) that was used for non-polar molecules as eqn (17) for polar molecules.
![]() | (17) |
Eexp = EvdWexp + ϕbEpolarexp | (18) |
The 2(1 − ϕb) term in eqn (17) is a factor that scales the input of thermal energy that drives expansion. For molecules that make only van der Waals interactions, this term is equal to one, so eqn (17) becomes identical to eqn (4). For molecules that make polar interactions, this term is less than one. When polar interactions are present, some of the thermal energy goes into rotational motion that breaks directional H-bonding interactions, and the rest goes into translational motion that leads to expansion (see Fig. 2). In effect, the polar interactions must be broken before expansion can take place, so the amount of thermal energy that goes into expansion is proportional to the population of bound states that are non-H-bonded.
The total number of SSIPs obtained by this procedure, NAIP, may differ from the total number of SSIPs calculated using the total surface area of the molecule. Since the number of SSIPs is important for an accurate description of the non-polar van der Waals contribution to the solvation energy, a correction is applied to solvation energies calculated using the AIP description in SSIMPLE (Eqn (19)).
![]() | (19) |
For common solvents that are liquids at room temperature and pressure, the vapour pressure is sufficiently low that we can assume ρv = 0 and ρl = ρ. An iterative algorithm was used to find the number density of the liquid that satisfies eqn (17). The number density of the liquid was initially set at ρc calculated using eqn (6). The SSIMPLE algorithm was then used to calculate the speciation of SSIP contacts and hence the fraction of bound states that are H-bonded, ϕb. Eqn (17) was then used to calculate a new number density. The average of the new and initial number densities was used in the SSIMPLE algorithm for the next iteration of this process, gradually concentrating the liquid until the number density converged.
Fig. 7 compares the experimentally measured densities of 171 liquids at 298 K (see ESI† S1 for details) with the calculated values. The agreement is generally excellent, and the RMSE of 0.87 M is largely down to the two outliers indicated in Fig. 7, hydrogen peroxide and hydrazine. For both of these liquids, the calculated density is significantly lower than the experimental value. This discrepancy could be due to inaccuracies in the calculated SSIP values, which are not polar enough, or due to cooperative multipoint interactions between molecules which are not considered in the SSIMPLE approach, which treats all SSIPs as independent.
![]() | ||
Fig. 7 Relationship between the experimentally measured densities at 298.15 K, ρexp, and the values calculated using eqn (17), ρcalc, for 171 liquids (RMSE = 0.87 M). The two outliers are hydrazine and hydrogen peroxide. |
![]() | (20) |
The value of C0 can be estimated from experimental data. For some liquids, e.g. water and toluene, C0 was found to be approximately zero at room temperature,97 and if we assume that this value does not change significantly with temperature, eqn (20) can be used to estimate the distribution of molecules between the liquid and vapour phases for the pure fluids at the saturated vapour pressure. An iterative algorithm was used to solve eqn (20) for the density of the vapour phase in equilibrium with the liquid phase. The number density of the vapour phase was initially set to zero, and the number density of the liquid phase was initially set to ρ, which was calculated as described above. The speciation of SSIP contacts in the liquid was used to calculate ΔGlS, and the number density of the vapour phase was calculated using eqn (20) starting with ΔGvS = 0. The new number density of the vapour phase was used to obtain a new number density for the liquid phase using eqn (4). SSIMPLE was then used to calculate the speciation of SSIP contacts in both phases, and the solvation energies were used in eqn (20) to obtain the number density of the vapour phase for the next iteration of this process, gradually concentrating the vapour and diluting the liquid until the number densities converged.
Fig. 8 and 9 show the results for water and toluene calculated for the entire liquid temperature range. The calculated mean density of the two phases, ½ ρ, is in excellent agreement with the experimental data for all temperatures. The temperature dependence of the number densities of the liquid and vapour phases are also reproduced well by the calculation. The only significant deviations between calculation and experiment occur for temperatures that are close to the critical point, where the liquid structure begins to break down.
![]() | (21) |
![]() | (22) |
Assuming that the C0 values are independent of the composition of the phases, these parameters can be used in SSIMPLE to calculate number densities in the two components of the vapour phase by solving eqn (20) for the binary mixture using the iterative procedure described above for VLE. The SSIMPLE vapour density is calculated at the saturated vapour pressure when the liquid and vapour pressure are in equilibrium, but the experimental values were measured at pressures of 101325–101
330 Pa. When the saturated vapour pressure is low, intermolecular interactions in the gas phase are negligible, so although the density of the vapour is reduced at lower pressures, the composition is unaffected. VLE data for binary mixtures are presented in the form of X–Y plots, where X(A) is the mole fraction of component A in the liquid phase, and Y(A) is the mole fraction of component A in the gas phase. Fig. 10 illustrates typical results for mixtures of two compounds that contain the same functional groups. These mixtures show close to ideal behaviour (X = Y), which is well-reproduced by the SSIMPLE calculation (see ESI† S8 for more examples). Fig. 11 shows examples of binary mixtures of dissimilar compounds. In some cases, there are significant deviations from ideality, which are again well-described by the calculated values.
The agreement between the SSIMPLE calculation and the experimental data was evaluated using the RMSE of the vapour phase composition (Y), and the results are tabulated in Table 1. The average RMSE between calculation and experiment across the entire dataset is 7.6%. The examples in Fig. 10 and 11, which have RMSE values that range from 1% to 9%, illustrate the quality of the predictions obtained using of SSIMPLE. Fig. 12 shows the X–Y plots for the three binary mixtures with largest values of RMSE (21%). The results in Table 1 are organised so that compounds with similar functional groups can be compared, but it is difficult to discern any clear patterns. The values of RMSE appear to be slightly higher for the chlorinated solvents. Two of the 21% RMSE examples in Fig. 12 involve chloroform, but there are other chloroform mixtures with low values of RMSE: acetic acid 6% and dichloromethane 4%.
Here the SSIMPLE approach is generalized to a temperature-dependent form. The polar interactions are temperature-dependent and can be described using a Boltzmann distribution between two different bound states for the interaction of two SSIPs. In one bound state, polar directional interactions are made, and in the other bound state, rotational motion breaks these interactions, leaving only the non-polar interaction. The approach provides an accurate description of the effect of temperature on experimentally measured association constants for the formation of 1:
1 H-bonded complexes in carbon tetrachloride. This formulation has been extended to describe the temperature dependence of the density of the fluid phase in terms of an expansion energy, which is made up of a constant non-polar contribution and a temperature-dependent polar contribution. The method shows very good agreement with the experimentally determined densities of 171 different liquids at room temperature.
Since non-covalent interactions in the vapour phase can be treated in the same way as in the liquid phase, the speciation of SSIP interactions can be used to calculate solvation free energies in both phases. However, treatment of vapour–liquid equilibria requires an additional compound-dependent free energy that is currently only available from experimental measurements. Using these empirical parameters, SSIMPLE was used to make accurate predictions of vapour–liquid equilibria for 196 binary mixtures at different temperatures. The temperature-dependent version of SSIMPLE described here thus represents a promising general method for predicting the equilibrium and solvation properties of liquids.
The generality of the approach stems from the fact that the calculations are based on the calculated ab initio molecular electrostatic potential surfaces of the molecules together with a relatively small number of empirically established parameters:
1. The van der Waals interaction energy between two SSIPs is 0.3 kJ mol−1 Å−2, based on experimental VLE data for noble gases, alkanes and perfluorocarbons.54
2. The relationship between the polarity of an SSIP and the calculated MEPS, based on experimental data on H-bonded complexes in carbon tetrachloride.53,99
3. The packing coefficient, r, and the zero point void, v, are 90% and 5.0 Å3, based on experimental data on noble gases.57
4. The relationship between the number of SSIPs in a molecule and the non-polar component of the fluid phase expansion energy, based on experimental data on alkanes described here.
The resulting speciation of SSIP interactions in a phase follows directly from these parameters. An additional empirical correction is required to describe vapour–liquid equilibria, but this free energy term can be estimated using the Antoine coefficients for pure compounds. Predictions for liquid and vapour densities, solvation energies, partition coefficients, and solute association constants can thus be made for a wide range of compounds.
Footnotes |
† Electronic supplementary information (ESI) available: Details on the convergence criteria and SSIP descriptions used for calculations, data on pure liquid concentrations, thermal expansion information, association constant data, VLE data for single component systems, details on the conversion of AIPs to SSIPs, information on the experimental VLE data, VLE X-Y plots. See DOI: https://doi.org/10.1039/d5cp00635j |
‡ These authors contributed equally. |
This journal is © the Owner Societies 2025 |