Thermodynamically effective molecular surfaces for more ecient study of condensed-phase thermodynamics

Evaluation of molecular surfaces plays the key role in a wide range of cutting-edge scientific fields and technologies, due to the well-characterized dependency between molecular surfaces and condensed phase thermodynamics. Numerous methods to evaluate molecular surfaces such as van-der-Waals and solvent accessible surface areas and various parameterizations for each one, have been proposed in the literature and typically yield quite diverse estimations of molecular surfaces. Despite this diversity, numerous successful applications have been reported for each one, which has become possible via ad-hoc modifications and parametrizations employed to accommodate inappropriately defined molecular surfaces. The main aim of the present study is to propose “thermodynamically effective” molecular surface which unlike the conventionally accepted molecular surfaces, can be defined only uniquely, can be measured experimentally for each molecule directly and straightforwardly, is defined based on a well-characterized theoretically described dependency between molecular surfaces and solution thermodynamics, and is highly accurate in evaluating various thermodynamics quantities in solution for a wide temperature range and different types of molecules, without requiring any ad-hoc modification.

the main applications of molecular surface estimation is theoretical evaluation of solution thermodynamics via continuum solvation models 10 which is an extensively applied method in very diverse scientific fields, ranging from catalysis 11,12 , advanced nanomaterials 13 , surface science 14 , or mechanisms of chemical reactions in the condensed phase 15,16 to unraveling the activity mechanism of coronavirus 17 .
Employing molecular surfaces for unravelling scientific challenges in the condensed phase has been an active scientific area with a history of more than a century. One of the earliest attempts in this regard dates back to 1886 and was proposed by Eötvös. He suggested a proportionality between the free energy per unit of interfacial surface, i.e. the surface tension, and the surface area of liquid-phase molecules 18 . This work was indeed one of the earliest examples of experimental evaluation of molecular surfaces, which was achieved by assuming solution-phase molecules as perfect spheres, allowing to evaluate molecular surface area via liquid molar volume. Although the assumption of perfect spheres is the simplest approach to get a rough estimation of molecular surfaces, it satisfactorily holds for mono-atomic molecules. Accordingly, the earliest successful applications of molecular surfaces to study solution thermodynamics both exploit surfaces determined via perfect sphere assumption and are mainly limited to noble gases 19,20 ,21,22 .
In 1964, the van-der-Waals (vdW) surface area concept was proposed in the pioneering work of Bondi 23 , which became the cornerstone of more advanced molecular surfaces such as Solvent Excluded Surfaces (SES) and Solvent Accessible Surface Area (SASA). Since then, a large number of methods and algorithms have been proposed to evaluate molecular surfaces in solution, in particular multiple variants of solvent excluded or solvent accessible surface areas. This wide variety of methods typically yields quite diverse estimations of molecular surfaces, as can be seen in a comparison of vdW and SAS molecular surfaces of ethylene, as depicted in figure 1. For a broader comparison, we provide computed molecular surfaces for 215 molecules via different parameterizations of the vdW method as supplementary material. These data imply that selecting the most appropriate method is not a trivial task.

Figure 1-comparison of vdW (left) and SAS (right) surfaces in ethylene
Surprisingly, despite the diversity of methods and molecular surface approximations they yield, there are numerous examples of reporting successful applications for each one of these methods in studying thermodynamics in solution. This is mainly because the majority of these research works are based on empirically defined relationships between solution thermodynamics and molecular surfaces 24 ,25 .
Inaccuracies due to deviations of the employed molecular surfaces from the actual values are then corrected via ad-hoc modifications and parametrizations, mainly applied to atomic radii. For example, while the Gaussian 03 software package used SES surfaces and UA0 atomic radii as default for computing solvent effects based on the polarizable continuum solvation models, for the latest release of the same software, vdW surfaces and UFF atomic radii are considered as default 26 . Similarly, the most widely applied continuum solvation models exploit their specifically defined molecular surfaces and parameterizations of atomic radii 10 .
One main reason behind this diversity in defining and parameterizing molecular surfaces and requiring such ad-hoc modification is that although the main application of molecular surfaces are commonly for studying solution phase thermodynamics, they are typically parameterized for reproducing other target quantities.
For example, the Bondi parameterization of atomic radii has been done using physical quantities like X-ray diffraction data, gas kinetic collision cross section and liquid density as target quantities 23 while the UFF or UA0 atomic radii are parametrized against bond distances 27 .
It indeed stems from unavailability of a rigorous theoretical method which allows precise and analytical characterization of the relationship between molecular surfaces and solution thermodynamics without adhoc modifications or parameterization.
Obviously, such a theoretical method offers a number of advantages. First and foremost, it allows defining molecular surfaces with physical significance, without requiring ad-hoc modifications and thus uniquely definable. Furthermore, it provides the possibility of evaluating the performance of various methods and parametrizations in reproducing those reference molecular surfaces. Last but not least, it allows for better understanding and treating some of the main challenges in theoretical studies of solvation, such as appropriate treatment of solvent effects in continuum solvation models.
To achieve this goal, in the present study we exploit a theoretically derived relationship describing temperature dependence of vaporization enthalpy to the molecular surfaces in solution which is an extension of a recent study 28 leading to a remarkable improvement of the formerly reported results.
Among all potential thermodynamics quantities of solution which can be analytically related to molecular surfaces for this purpose, the vaporization enthalpy, which is employed in the present study, possesses a number of obvious advantages. The main one is that vaporization enthalpy can be directly determined experimentally, while free energy or entropy can only be determined indirectly and via measuring the temperature dependence of vaporization enthalpy or equilibrium vapor pressures at multiple temperatures, which implies accumulation of errors inherent in both experimental measurements and the subsequent computations. The more convenient experimental procedure for enthalpy measurement has also made accurate benchmark datasets more readily available, which is another advantage of using vaporization enthalpy. Finally, evaluation of molecular surfaces via vaporization enthalpy is not only both more accurate and less challenging but also once it is found, it can be conveniently used to obtain relationships between molecular surfaces and other thermodynamic quantities, via the fundamental thermodynamics relationships, as shown in section 4-3.

2-Theory
By considering vaporization as a dynamic process at which evaporation and condensation have the same rates and equating the rates of evaporation and condensation described by transition state theory and some manipulations, the ratio of partition functions of the gas and liquid phases is obtained as 28 : where A is Avogadro's constant, CAD is the saturation vapor pressure of the liquid, Δ CG is the energy for moving one molecule from the liquid surface to the gas phase, and I and ℎ are Boltzmann and Planck constants, respectively. Using the statistical thermodynamics relationship between the energy and partition function stated as: and with some algebraic manipulations, it can be shown that the temperature dependence of the vaporization enthalpy follows 28 : where is a constant and Δ SC is the energy required for moving one molecule from the bulk of the liquid to the surface. Evaluation of Δ SC via experimentally measurable quantities can be achieved using the fundamental thermodynamics relationships between energy ( ), Helmholtz free energy ( ) and entropy ( ), which implies 28 : where Δ SC is the free energy change for moving one molecule from the bulk of liquid to the surface.
Another straightforward way to obtain Eq. (4) is using the Gibbs-Helmholtz equation: This then clearly yields Eq. (4) by subtracting the resulting equations for the bulk and surface states.
Exploiting the thermodynamics relationship among Δ SC , surface tension (γ) and the molecular surface area ( C ) which is defined as 28 : Eq. (4) can be rewritten as: Halving the molecular surfaces C in Eq. (7) is considered here because in fact only one half of the molecular surfaces contribute in forming the gas-liquid interface and the other half remains in the bulk of the liquid 28 . the correlation between the surface tension and molar vaporization enthalpy ( QAR ) is obtained as : in which is a constant. Knowing that at the critical temperature both vaporization enthalpy and surface tension approach zero, and due to continuity of the surface tension, the _o _. term also approaches zero, the constant is found as: which by substitution into Eq. (9) finally results in: As discussed earlier, the main advantage of the theoretically derived relationship among vaporization enthalpy, surface tension and molecular surfaces described by Eq. (11)

3-1-Dataset
The theoretically derived methods were benchmarked using thermophysical data of the DIPPR801 database 29 . Screening the initial dataset and selecting only the compounds with maximum uncertainty of 5% for vaporization enthalpy and surface tension resulted in 215 compounds from diverse chemical families, provided as supplementary materials.
For each compound, the experimentally determined data of vaporization enthalpies for 25 temperatures linearly distributed between the melting point and the critical temperature were evaluated using the provided relationships in the DIPPR database. Due to the scarcity of accurate surface tension data at all of the required data points at which vaporization enthalpy data were available, we employed the Guggenheim-Katayama relationship stated as 30 : as a rigorous model for evaluating temperature dependence of surface tension and its temperature derivatives required by Eq. (11). One main advantage of employing the Guggenheim-Katayama relationship is that it perfectly satisfies the boundary conditions, i.e. yielding exactly zero for surface tension and its higher order derivatives with respect to temperature at the critical point, which is commonly violated by other relationships like those proposed in the DIPPR dataset. This is indeed a key requirement, as it was one of the premises of obtaining the constant in Eq. (11). Additionally, compared to various surface tension predictive models such as those used by the DIPPR database and the Eötvös relationship 18 , we found that the Guggenheim-Katayama relationship provides the most accurate evaluation of not only temperature dependence of surface tension but also its derivate for the whole temperature range and hence, the most accurate prediction of vaporization enthalpy via Eq. (11).
To obtain surface tension data at the required temperatures via the Guggenheim-Katayama relationship, for each compound the pre-factor ∘ was calculated by optimization using the available experimental data points of surface tension. The calculated values of ∘ for each compound is reported in the supplementary materials.
The accuracy of the predicted vaporization enthalpy is reported as Average Absolute Deviation (AAD), defined as: Considering that at the critical point, vaporization enthalpy approaches zero and slight deviations in predicted enthalpies results in very large relative errors, AAD provides a more reliable evaluation of the model performances, compared to relative errors.

3-2-Computational details
To calculate well-stablished molecular surfaces, the geometry of each molecule was first optimized at the

4-1-Verifying the validity of the theoretically derived relationship
By studying various relationships proposed in the past century to analytically relate solution thermodynamics, surface tension and molecular surfaces, we noticed an obvious inconsistency not only between previously proposed models themselves but also with the theoretically derived relationship proposed in the present study as well. Accordingly, while a large number of studies 19,20 ,21,22,25,32,33 support: where C… †QAD•…7 is the free energy of solvation, is the solvent excluded surface of molecules and ℬ is a constant 34 , many other works employ the very similar relationship to relate vaporization enthalpy, surface tension and molecular surfaces. A well-known example of the latter category is Kabo's method which proposes: and has been extensively applied specially in studying the phase equilibrium in ionic liquids 35  Alongside the mentioned paradoxical deviations between the two models, both of them also show obvious inconsistencies compared to our theoretically derived relationship proposed in Eq. (11) which implies the necessity of a careful and rigorous verification of our model.
To that end, we first evaluated the overall accuracy of vaporization enthalpies predicted via the newly developed relationship. Accordingly, for each compound we optimized the value of the C parameter required by Eq. (11) which yielded lowest error in predicting vaporization enthalpies over the whole temperature range. Via the optimized C parameters, which are in fact our proposed "thermodynamically effective" molecular surfaces, an AAD of 0.188 kcal/mol was obtained for the predicted vaporization enthalpies of the whole dataset. This resulting AAD is within both the chemical accuracy and the reported accuracy of the employed, experimentally determined data.
Interestingly, we noticed that the molecular surfaces evaluated via molar volumes at the melting point based . u 2 terms, which are the most obvious differences between our proposed relationship and the conventionally accepted models, we re-optimized C parameters for the two following relationships: and which are two variants of Eq. (11) obtained by removing the term being studied. Using the re-optimized C parameters, the two abovementioned variants yielded AADs of 2.071 and 0.197 kcal/mol, respectively.
These results clearly show that both terms (

4-2-Evaluation of molecular surfaces estimated via computer algorithms
After verifying the validity and robustness of the theoretically derived relationship in the previous section, this section focuses on studying predictability of vaporization enthalpies obtained via various parameterizations of the vdW and solvent accessible surfaces. To that end, we have studied a total number of 252 differently computed molecular surfaces discussed in section 3-2.
According to the results, while via the solvent accessible surfaces we could not achieve any AAD better than 7.889 kcal/mol, for the vdW surfaces the best results with AAD of 0.568 kcal/mol was obtained for   to molecular surfaces is provided in section 4-3.

4-3-Evaluation of other thermodynamics quantities via thermodynamically effective surfaces
As it was discussed earlier, one of the main advantages of characterizing dependency of solution thermodynamics on molecular surfaces through temperature dependence of vaporization enthalpy is its straightforward transferability to other thermodynamics quantities.
As one of the most important thermodynamics quantities, evaluation of the solvation free energy, which is the primary goal in continuum solvation models 10 , can be achieved via the Gibbs-Helmholtz relationship as follows: The abovementioned dependency between vaporization enthalpy and solvation free energy also allows a quantitative comparison of our theoretically derived relationship with the other conventionally accepted model, which in contrast to our model suggests linear dependency between solvation thermodynamics and the product of molecular surfaces and surface tension via Eq. (14). To that end, we studied the predictability The obtained accuracy of these results shows an AAD of only 0.1215 kcal/mol which is by almost a factor of 2 more accurate than best results obtained via advanced continuum solvation models 10 . .
The saturation vapor pressure of non-ideal gases can also be determined the same way via evaluating the temperature dependence of gas phase molar volumes using appropriate equation of states. Nevertheless, in the present study and only for a proof of concept, we study predictability of saturation vapor pressures of ideal gases (pressures up to 2.5 atm) via Eq. (20). A comparison of predicted and experimentally determined saturation vapor pressures for a number of most widely used solvents is depicted in figure 4. The excellent agreement between the theoretically evaluated and experimental data depicted in figure 4 implies the robustness of our proposed thermodynamically effective surfaces. To sum up, in the context of the present study, we could theoretically derive a relationship which describes dependency between vaporization enthalpy, molecular surfaces and surface tension. We could demonstrate that the proposed dependency between solution thermodynamics and molecular surfaces are remarkably more reliable and rigorous compared to the other models empirically proposed for the same purpose within the last century.
Through our newly derived theoretical approach, we proposed thermodynamically effective surfaces. We demonstrated that the "thermodynamically effective" surfaces are only slightly different than empirically proposed vdW surfaces. However, this slight deviation yields a substantial improvement in predictability of multiple thermodynamics quantities. As a result, we propose the thermodynamically effective surfaces as a reliable alternative for the currently defined molecular surfaces, especially for studying the condensed phase thermodynamics.