{'fraction_non_alphanumeric': 0.03176958324241763, 'fraction_numerical': 0.011346279729434867, 'mean_word_length': 5.411863458310017, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Partial Hessian vibrational analysis: vapour pressure isotope effects and their relation with non-covalent interactions

In this study equilibrium isotope effects (EIEs) on evaporation to several organic solvents (bromobenzene, dibromomethane, ethanol, methanol, and trichloromethane) in the pure phase are estimated employing Kohn-Sham Density functional theory (KS-DFT) along with full Hessian vibrational analysis (FHVA) and partial Hessian vibrational analysis (PHVA). Both FHVA and PHVA qualitatively predict compounds EIEs. Weakly interacting systems (bromobenzene, dibromomethane, and trichloromethane) display identical EIE values for both, FHVA and PHVA. Whereas for strongly interacting systems (ethanol and methanol), PHVA estimates slightly smaller EIE compared to FHVA. By employing the symmetry-adapted perturbation theory (SAPT) along with independent gradient model (IGM), it was possible to establish how minimal changes in compounds interaction energy affected carbon and bromine EIEs estimation in dibromomethane. Moreover, SAPT and IGM revealed how in ethanol, apparently insignificant intermolecular interactions possessed an enormous impact in carbon position specific isotope effects, which highlight the importance of solvent effects for strongly interacting compounds. The presented results have important implications for correctly comprehending and characterising the role of non-covalent interactions in the prediction of isotope effects, as well as providing a robust test to PHVA for the prediction of EIEs.


Introduction
Equilibrium isotope effects (EIEs) are one of the most successful tools in biochemistry, geochemistry, and chemistry for meticulously studying exchange reactions during processes such as volatilisation, sorption, and evaporation. [1][2][3][4] Given processes are vital to comprehend pollution, biodegradation pathways and remediation of contaminated aquiferous systems like oceans, lakes, and underground reservoirs. Halogenated compounds such as trichloromethane and bromobenzene represent a global concern, given that are not merely perilous pollutants with suspected carcinogen properties, but equally are known for their slow dissolution (per years) in anoxic environments. 4 Technological advances during the last two decades permit now the accurate measuring of position specific isotope effects, hence providing an valuable insight into the nature of bulk isotope effects. However, precise interpretation of experimental results remains a complex task, predominantly owing to physicochemical phenomena that might be easily overlooked such as diffusion or adsorption, but posses a substantial influence over the outcome. 5,6 For instance, it is conceivable that the diffusion effect vanishes from compounds boiling point on, leaving as result mere EIE or so-called vapour pressure isotope effects (VPIE). 2,7 However when direct measurements are performed at standard conditions, to some extend the diffusion effects are present, under those circumstances the term volatilisation isotope effect (VIEs) is employed to emphasise the temperature factor.
Usually, owing to the quantum nature of EIEs, 8 experimental data is adequately depicted via electronic structure calculations. Thus, finding stationary points on a potential-energy surface (PES), the theory of molecular vibrations of molecular systems 9 is employed by the Bigeleisen-Mayer approach for the practical calculation of isotope effects. 8 Predominantly, electron structure calculations based on Kohn-Sham density functional theory (KS DFT) 10 are used. Given choice is motivated mainly by the relative robust performance of density functionals and their affordable computational cost for the full quantum-mechanical treatment of large systems (molecular systems with more than 40 atoms). In practice, the B3LYP functional [11][12][13][14] is likely the most used functional for isotope effects calculation, 9 and although in recent KS DFT benchmarks it does not perform well, 15 its notable ability to predict correctly the isotopic frequency shifts converts it in one of the DFT functionals for excellence. Another exceptional functional is MPW1PW91, 16-18 which has permitted the accurate predictions of hydrogen and carbon VPIEs for fluoroform. 19,20 Important, and a point often overlooked, is the quality of stationary point atomic positions. It must be remembered that the Hessian is the second derivative of the molecular energy with respect to nuclear coordinates (Cartesian coordinates), therefore undergoing an accurate reproduction of them intentionally serves as a pivotal step in order to calculate reliable isotope effects for the precise reason. Moreover, correctly reproducing compounds interatomic interactions conjointly guarantee solvent effects emulation. Goerigk et al. have demonstrated in a recent study that weighted total mean absolute deviation for non-covalent interaction calculated by B3LYP and MPW1PW91 functionals are notably large, especially if dispersion corrections are not included. 15 Therefore, meticulous care during the functional selection is required, given that generally accurate atomic positions implies expensive numerical calculations within the "Jacobs ladder" rungs.
After carefully choosing an adequate DFT functional for the purpose of calculating EIEs, the Hessian matrix must be computed, however, because the Hessian matrix is routinely calculated numerically (finite difference of analytic gradients) 21,22 within KS DFT, even for relatively compact systems its calculation is computationally expensive. 23 In order to surmount given drawback, either software package parallelisation along with more efficient algorithms or rigorous Hessian approximations are been continuously developed. [23][24][25][26][27] The later is appealing, particularly the partial Hessian approach, 28,29 which principle lies on the diagonalisation of a sub-block of the mass-weighted energy second derivative matrix computed for a subgroup of atoms. This pragmatic approach has been successfully employed for characterisation of the infrared spectrum of large transition metals, 30 vibrational analysis of relevant cluster compounds, 31 fingerprint vibrations of adsorbates on surfaces, 28,29,31 and molecular vibrational energy and entropy. 32 Although partial Hessian 28 or adapted versions have been efficiently implemented in the past for estimation of EIEs, 33,34 those studies did not explore at least one of the following issues i.) the quality of intermolecular interactions compared to experimental data or reliable levels of theory, ii.) equivalent configurations of PES by nuclear coordinates, iii.) effect of partial Hessian in the reduced partition function ratio (RPFR) terms, 8 or iv.) practical limitation of partial Hessian approach.
We report on numerical calculation results of EIEs for two organic solvents which mainly interact via hydrogen bonds: ethanol, and methanol, and other four highly polluting halogenated compounds bromobenzene, trichlorofluoromethane, dibromomethane, and trichloromethane.
EIEs results from full Hessian and partial Hessian vibrational analysis are properly compared, as well as their direct effects on the components of the reduced partition function ratio. Moreover, unique insight into intermolecular interactions of certain studied compounds and their relationship with EIEs are carefully analysed and described.

Methods
Molecular dynamic simulations were performed for bromobenzene (BB -C 6 H 5 Br), trichlorofluoromethane (CFC -CFCl 3 ), dibromomethane (DBM -CH 2 Br 2 ), ethanol (ETH -C 2 H 5 OH), methanol (MeOH -CH 3 OH), and trichloromethane (TCM -CHCl 3 ) (table 1) using the opensource software LAMMPS (large-scale atomic/molecular massively parallel simulator). 35 Optimized potentials for liquid simulations (OPLS) force field parameters 36,37 from the LigPar-Gen web server 38 have been used for model simulations. Periodic boundary conditions for a cubic box comprising 90 molecules have been properly minimised and brought to equilibrium at standard conditions. Consequently, for each compound, five clusters containing four molecules were extracted from trajectories of equilibrated MD simulations. Nuclear coordinates of clusters were further optimised using a long-range-corrected (LRC) hybrid functional with non local correlation,  39 Basis set double, and triple Dunning's augmented correlation consistent 40,41 were used. All geometrical optimisations were carefully performed using the SG-2 42 integration grid. Stationary points on the potential surface were confirmed from full Hessian vibrational analysis (FHVA) without imaginary frequencies. During the vibrational analysis the SG-3 grid 42 was employed. All quantum mechanical calculations were carried out in Q-Chem software packages (version 4.4), 43 and were typically restricted to an allowable maximum of 4032 computational hours.
Partial Hessian vibrational analyses 28,31 (PHVA) were performed to clusters for which FHVA was carried out successfull. In PHVA only one molecule was "active", i.e. the second partial derivatives of the potential with respect to displacement of the atoms is calculated, whereas the other three remained "frozen". Given approach was used with each one of the molecules within a cluster.
With the purpose to robustly assess the accuracy of PHVA with respect to FHVA, equilibrium isotope effects (EIE) were properly calculated implementing the Bigeleisen-Mayer approximation from the ratio of the free energies of the isotopologues, 8 thus, using vibrational analysis for Bigeleisen-Mayer approximation gives where h is the Planck's constant, c represents the speed of light in vacuum, k B introduces the Boltzmann constant, T refers to temperature, ν indicates the harmonic frequency, and (s h /s l ) f is the reduced partition function ratio.
It is worth to note that EIE can be expressed in terms of enrichment factor ε as where n represents the number of atoms of the same element in a molecule (?? ) that have been isotopically substituted. EIE and ε have been calculated in the developing version of the open source program isoC. 44 Enrichment factor ε is used along this work for consistency with experimental data.
A decomposition energy analysis is performed in order to obtain an insight into compounds interactions strength and correctly identify patters between FHVA and isotope effects. Clusters were divided in single molecules (four per cluster) and all arranges (six) of possible dimers were promptly formed, consequently, interaction energy E int for dimers were calculated using the symmetry-adapted perturbation theory (SAPT). 45 Within the so-called supermolecular approach that SAPT follows, interaction energy E int of a system typically consisting of two separable subsystems A and B is defined as In parallel, the SAPT interaction energy can be further expanded and decomposed into four physically meaningful components i.e. electrostatic, exchange, induction, and dispersion yielding where E elst is the electrostatic component, E exch represents the exchange term, E ind gives the inductions component, and E disp indicates the dispersion term of interaction energy.
Owing to the relative compact size of the studied molecules, a large basis set (Dunning's correlation-consistent augmented triple) along with accurate higher-order SAPT (SAPT+2) with coupled-cluster doubles (CCD) treatment of dispersion and supermolecular MP2 interaction energy (δMP2) 40 Complementarity to SAPT analysis, non-covalent interactions were visualised and meticulously studied by means of independent gradient model (IGM). 47 Intermolecular interactions described by the density gradient as implemented in the stable release of Multiwfn package software (version 3.6) was used. 48 Table 2 contains average root-mean-square deviation (RMDS) of aligned atomic positions from successfully optimised geometries that also were confirmed to be stationary points without negative frequencies. Although all five BB clusters were geometrically optimised, Table 2: Root-mean-square deviation of atomic positions for geometrically optimised clusters (CFC, DBM, EtOH, MeOH, and TCM) at ωB97X-V/aug-cc-pvdz and ωB97X-V/aug-cc-pvdz levels of theory.

RMSD [ A]
Compound Atoms compared aug-cc-pvdz aug-cc-pvtz stationary points could not be validated from FHVA calculations as those exceeded the computational time limit. It is important to realise that generally speaking, clusters are heterogeneous and possess relatively large RMSD, indicating that VIE and VPIE predictions are obtained to several configurations of the PES. The increase of the basis from aug-cc-pvdz to aug-cc-pvtz, in general only slightly decrease the RMSD.
Two pair of clusters, ones in case of EtOH (EtOH-1 and EtOH-2) and MeOH (MeOH-3 and MeOH-4), resulted in identical atomic positions, which can be interpreted as likely abundant configuration within the PES. Another key point from RMSD is its minimal value for oxygen atoms in EtOH and MeOH, probably owing to their relevance in the formation of hydrogen bonds (HB). It is well recognised that given compounds mainly interact via hydrogen bonding, a weak non-covalent interaction, thus it is reasonably expected that oxygen atomic positions would not variate or at least not as halogen atoms, which noncovalent interactions are significantly weaker.
Average VIEs for studied compounds are given in table 3. Isotope effect results for CFC Table 3: VIEs (ε ) predicted at standard conditions for trichlorofluoromethane, dibromomethane, ethanol, methanol, and trichloromethane, at ωB97X-V/aug-cc-pvdz level of theory.
produce a qualitatively match with respect to experimental data, 3 however it is underestimated. Small mean standard deviation (STD) from position-specific isotope effects (see table S1 of supplementary information) in CFC along with RMSD, indicate that although cluster atomic positions differ appreciably from each other, this variety in the PES does not affect the final prediction. 81 Br isotope effect for DBM is overestimated, while 13 C is rightly predicted as normal, but its magnitude is three timer larger than the one from experimental data. The last from the group of halogenated compounds is TCM, which 13 C isotope effect is only predicted qualitatively, i.e. an inverse isotope effects is predicted, however its value is a factor of magnitude underestimated.
Average position-specific isotope effect for 13 C 1 in EtOH is underestimated, whereas 13 C 2 is not predicted qualitatively either quantitatively. Nevertheless, it is worth mentioning that position-specific isotope effects values for 13 C 2 were in their vast majority close to zero, with only a couple being larger than one. Thus, if 13 C 2 is recalculated with the position-specific values closer to zero, the predicted isotope effect respectively becomes 0.27 . The apparent reason behind such particularly differences in the positions-specific magnitudes is diligently studied and explained in the next section (?? ).
According to our compounds scheme, the second carbon position in EtOH and the carbon position in MeOH are similar, as both atoms are bonded to OH, therefore it is interesting to characterise how different or similar their position-specific isotope effects are. In conjunction with VIEs, VPIE were calculated for compounds with present experimental data. 13 C VPIE for EtOH and MeOH at ωB97X-V/aug-cc-pvdz and ωB97X-V/aug-ccpvtz levels of theory are given in the tables 5 and 6, respectively. Contrary to experimental results, predicted VPIE are notoriously unaffected by the temperature change. In the case of 13 C 2 for EtOH and 13 C in MeOH, predicted values are closer to experimental data, specially 13 C 2 . However, it must be recalled that it is precisely the position-specific isotope effect for carbon atom bonded OH, the one that reflected large STD.     Table 8 contains VIEs predictions for studies systems at ωB97X-V/aug-cc-pvdz level of theory. Predominantly, numerical calculations reveal that unless molecules in the cluster are Table 7: VIEs (ε ) predicted from partial Hessian vibrational analysis at standard conditions for bromobenzene, trichlorofluoromethane, dibromomethane, ethanol, methanol, and trichloromethane, at ωB97X-V/aug-cc-pvdz level of theory. It must be remembered that within the concept of the reduced partition function ratio (RPFR) by the Bigeleisen-Mayer approximation, 8 it is equally possible to describe the isotope effect as the product of three factors, i) the product factor (PF), which contains molecular masses and moments of inertia, ii) the excitation factor (EXP), which introduces a correction for excited molecules, and iii) the zero-point energy contribution (ZPE) factor, which describes the zero-point energy differences for each of the vibrational modes. 9,50 As a rule, PF tends to remain a minor factor, as well as EXP which vanishes at elevated temper-atures, therefore, ZPE represents the factor that typically determines the direction of the isotope effect (normal or inverse).

Effects from Hessian matrix approach
Given that only one molecule in the cluster was active during PHVA, the PF in all performed calculations is close to one. The effect of PHVA in the remaining terms (EXP and ZPE) was studied by calculating those directly within isoC. 44 For direct reference TCM and CFC from FHVA and PHVA were used, as precisely named systems did not present differences in the averaged VIE. There were not found relevant patterns for ZPE from PHVA.
The only observed characteristic was that this term seemed to fluctuate in a convenient way to properly compensate the fact that EXC is always larger than its respective term from FHVA, thus ZPE values from PHVA were shifted. This observation is clearly evident when statistical parameters such as STD, skewness and kurtosis were compared between PHVA and FHVA, indicating close or identical values.

Role of non-covalent interactions in the prediction of isotope effects
Total interaction energies via SAPT along with visualisation of non-covalent interactions from IGM analyses, was used in order characterise the magnitude of interactions in the first cluster of each compound optimised at the ωB97X-V/aug-cc-pvdz level of theory (see fig. S1 in supplementary information).
As each cluster consists precisely of four molecules, a total of six unique dimer arrangement can be formed, thus for instance to ethanol cluster the following dimers were formed:  Results demonstrate that CFC possesses the weakest interactions among studied compounds. In increasing order, CFC is followed by TCM, DBM, BB, MeOH, and EtOH. Although from SAPT calculations BB does not seem to typically reflect relatively strong interaction energies, the orange bar stands out over the rest. A comprehensive analysis of given result indicates it corresponds to a π-π stacking (parallel displaced), which is promptly confirmed after properly comparing dispersion and electrostatics components from the interaction energy decomposition (see fig. S3 in supplementary information). Consequently, interaction energy components (electrostatics, exchange, induction, and dis- Table 8: Average interaction energies from SAPT analysis for geometrically optimised clusters (BB, CFC, DBM, EtOH, MeOH, and TCM) at ωB97X-V/augcc-pvdz and ωB97X-V/aug-cc-pvdz levels of theory. interactions, as the key reason for the immoderate change of isotope effects prediction.
Using IGM for EtOH cluster, we detected that apparently a moderate van der Waals force (red circle in fig. 2) is responsible for the large STD during the calculation of 13 C 2 VIE. Given intermolecular interaction was only observed over/around C 2 atomic positions that, correspondingly, after being isotopically substituted to 13 C, exhibited relative large magnitudes (¿ 1.0 ). Because EtOH molecules interact massively via hydrogen bonds, it is unreasonable to dismiss the fact that the intensity or lack of hydrogen bonds might be equally contributing to the overall STD. Which is precisely the case of MeOH cluster.

Conclusions
We have calculated average volatilisation-pressure isotope effects for trichlorofluoromethane, dibromomethane, ethanol, methanol, and trichloromethane from full Hessian vibrational analysis at the ωB97X-V/aug-cc-pvdz level of theory. Generally, isotope effects were predicted in qualitative agreement with experimental results. Isotope effects calculations with a larger basis set (aug-cc-pvtz) exhibit smaller standard deviation errors, and in most of the cases their magnitudes did not change appreciably compared to the smaller basis set (augcc-pvdz) results. Vapour-pressure isotope effects were in better agreement with experimental data.
As for equilibrium isotope effects calculation, partial Hessian vibrational analysis 28,31 provides a robust, computational inexpensive and easy approach for their accurate calculation.
We observed that unless intermolecular interactions of the studied compound are carried out via via hydrogen bonds, accuracy of predicted isotope effects is uncompromised and equals isotope effect magnitudes from full Hessian vibrational analyses. It is substantial to note that using partial Hessian vibrational analysis, average computational time was reduced by a 75% in comparison to full Hessian vibrational analysis.
Furthermore, we noticed that the product factor and the excitation factor of the reduced partition function ratio by the Bigeleisen-Mayer approximation, 8   General theory and an application to the cyclobutene-butadiene isomerization using  Table S1: VPIEs (ε ) predicted on evaporation using aug-cc-pvdz basis set for trichlorofluoromethane, dibromomethane, ethanol, methanol, and trichloromethane.

Compound
isotope ε STD Exp.   Figure S5: Interaction energy for studied compounds of cluster one optimised at ωB97x-V/aug-cc-pvtz. Each colour bar indicates one dimer combination, within a four molecule cluster. Considering position-specific isotope effects smaller than one. 3 Considering position-specific isotope effects larger than one.

S7
ωB97x-V/aug-cc-pvdz ωB97x-V/aug-cc-pvtz Figure S6: Visualisation of DBM intermolecular interactions by IGM analysis to optimised atomic positions at the ωB97x-V/aug-cc-pvdz (left) and ωB97x-V/aug-cc-pvdz (rigth) level of theory. BGR colour space describes intermolecular interactions, where blue highlights strong attraction, green are van der Waals interactions, and red are strong repulsive.
ωB97x-V/aug-cc-pvdz ωB97x-V/aug-cc-pvtz Figure S7: Visualisation of MeOH intermolecular interactions by IGM analysis to optimised atomic positions at the ωB97x-V/aug-cc-pvdz (left) and ωB97x-V/aug-cc-pvdz (rigth) level of theory. BGR colour space describes intermolecular interactions, where blue highlights strong attraction, green are van der Waals interactions, and red are strong repulsive.