Hindered rotation and bending anharmonicity in aluminum alkyls: implications for methylaluminoxane thermodynamics

Abstract

Accurate thermodynamic calculations for aluminum alkyls require proper treatment of low-frequency vibrations poorly described by the harmonic approximation (HA). Here, we present a systematic investigation of hindered rotation and out-of-plane bending in aluminum trichloride (ATC) and its methyl derivatives, employing advanced computational methods to perform anharmonic entropy corrections, such as torsional eigenvalue summation (TES), the extended two-dimensional torsion method (E2DT), the multi-structural approximation with torsional anharmonicity (MS-T), and Fourier grid Hamiltonian (FGH). Our results reveal distinct structure-dependent behaviors: monomers exhibit near-free methyl rotations where the HA overestimates entropy by 20–30 J K⁻¹ mol⁻¹, while dimers show more hindered rotations adequately described by the HA around room temperature. Out-of-plane bending modes in dimers display increasing energy level spacing that reduces their thermodynamic contribution compared to HA predictions. Simple quasi-harmonic approaches reduce HA inaccuracy for monomers but systematically underestimate the entropy of low-frequency modes in dimers. By applying appropriate methods to ATC and trimethylaluminum (TMA), we achieve excellent agreement with experimental entropy values (within 3.5 J K⁻¹ mol⁻¹). This validation supports extending our approach to methylaluminoxane (MAO), critical for understanding polyolefin catalysis. The dimer species studied here directly relate to MAO edge sites involved in metallocene catalyst activation. Our findings suggest that mode-specific anharmonic corrections are essential for accurate thermodynamic modeling of MAO and similar aluminum-containing systems.

---

Introduction

Aluminum trichloride (ATC) and its methyl derivatives are industrially significant compounds, whose major applications lie in thin film technology and polyolefin catalysis. For the atomic layer deposition of aluminum oxide (Al₂O₃), both ATC^1,2 and dimethylaluminum chloride^3–7 (DMAC) can serve as precursors, although the process based on trimethylaluminum (TMA) and water is considered ideal.^2,8 Additionally, TMA is the key building block of methylaluminoxane (MAO), a staple metallocene cocatalyst in polyolefin catalysis.⁹

A central challenge in the computational modeling of these aluminum systems is their thermodynamic characterization. This challenge is especially relevant to our research interests in MAO, where accurate computation of entropy and free energy has proven critical for understanding its structural preferences. Despite MAO's industrial importance, its precise structure remained elusive for over four decades. Yet, a sheet form, predicted by our electrospray ionization mass spectrometry supported computational studies,^10–13 was recently confirmed by crystallographic data.¹⁴ Crucially, the preference for the sheet structures over previously proposed cages¹⁵ appears to be explained by their higher vibrational entropy.¹⁶

Still, the quantitative accuracy of this observation remains uncertain due to the incapability of the standard computational approach for vibration, the harmonic approximation (HA), to appropriately describe MAO low-frequency vibrations. Beyond structural determination, proper description of these anharmonic vibrations would be essential for understanding the thermodynamic driving forces of MAO reactivity. However, the problem with the HA is not limited to the large MAO structures: previous studies^13,17,18 have shown that the HA fails to capture the thermodynamics of even the simplest MAO building block—trimethylaluminum. More specifically, the HA overestimates the entropy of dimerization of TMA (−238 J mol⁻¹ K⁻¹) by over 30% compared to the experimental value (−181 J mol⁻¹ K⁻¹).

This discrepancy, together with the need to understand the thermodynamic factors governing MAO's structural preferences and reactivity, motivated us to investigate the vibrational entropy of aluminum compounds in more detail, exploring theoretical approaches beyond the HA. For this purpose, ATC and its methyl derivatives provide an ideal starting point as they are the simplest organoaluminum compounds, also exhibiting structural similarity to MAO. The MAO sheet structures contain aluminum sites with bonding environments reminiscent of the bridging motifs found in TMA dimers (see Fig. 1), providing reactive centers¹⁹ for precatalyst activation. The chlorinated species are also relevant for MAO applications, as chloride abstraction from metallocene precatalysts can lead to analogous structural environments.^20,21 Therefore, systematic investigation of HA alternatives for these model systems is essential, beginning with an understanding of the fundamental limitations of harmonic treatments.

Fig. 1 Structural similarity between trimethylaluminum dimer (left) and edge of a methylaluminoxane sheet (right). Al in Pink, O in red, C in grey, H omitted for clarity.

The main problem with the HA is its poor potential energy description of large-amplitude vibrations: torsions and out-of-plane bends. Thermodynamically, this misrepresentation is critical as these vibrations tend to have low frequencies, thus contributing substantially to partition functions. In most cases, this contribution is overestimated by the HA. One solution is to use the quasi-harmonic (qh) approximations which implement an artificial frequency cutoff; vibrations below this threshold are either considered rotations²² or raised²³ to the cutoff value. These approximations, also utilized in our earlier work, can help to prevent overestimating the vibrational contribution by the HA. However, imposing artificial frequency cutoffs lacks physical basis; more accurate results require solving the Schrödinger equation (SE) for torsions and out-of-plane bends with more suitable potential energy description.

Torsional motion can manifest as free rotation, hindered rotation, or vibration, depending on the relation between potential energy barriers and thermal energy. Because the rotational barriers of single bonds are relatively low, torsional motion typically occurs as hindered rotation even at sub-ambient temperatures. Although the first correction models to hindered rotation were developed by Pitzer and Gwinn²⁴ in the 1940s, more practical and accurate approximations^25,26 have surfaced in the past three decades. Originally, Pitzer and Gwinn numerically solved the SE for the hindered rotation of ethane using a simple periodic cosine potential

where V₀ is the barrier height, M is the periodicity of the rotor and ϕ is the torsional coordinate. They also formulated an approximation (PG) that mimics these results by interpolating between harmonic and free rotor (FR) partition functions. Decades later, small adjustments to the PG approximation were made by Truhlar²⁷ and subsequently by McClurg, Flagan, and Goddard^28,29 (MFG) and Ayala and Schlegel³⁰ (AS). For TMA, Ehm et al. tested the MFG method, significantly improving the computational value for the free energy of dissociation.¹⁷

While these approximations work well for symmetric hindered rotors (HRs) with periodic potential, more advanced methods are required for asymmetric HRs. The most accurate results for one-dimensional hindered rotation are obtained via torsional eigenvalue summation (TES).²⁵ Using a higher order Fourier series to describe the potential energy, TES computes the torsional energy levels via diagonalization of Hamiltonian matrix in a plane wave basis.

With multiple torsional degrees of freedom, obtaining accurate partition functions becomes challenging due to more complicated torsional potential energy surfaces (TPESs) and kinetically coupled rotors. Fernández-Ramos et al. developed the extended two-dimensional torsion method^31,32 (E2DT) to calculate accurate torsional eigenvalues for two-dimensional cases using Fourier potentials, analogous to TES.

Beyond two torsions, solving the SE accurately becomes prohibitive. For these cases, the multi-structural approximation with torsional anharmonicity (MS-T), developed by Truhlar et al.,^33–35 is the most suitable option. The MS-T utilizes the basic interpolating formula of the PG approximations, while also addressing conformations and kinetic energy coupling of rotors. The thermodynamic performance of MS-T was evaluated by Gao et al.,³⁶ who showed MS-T results to align well with experimental values for heat capacity and entropy of hydrocarbons.

For out-of-plane bends, computational approaches are less established than for hindered rotation. Some studies^37–39 have considered these bends in the context of transition state theory, adding quartic corrections to the HA quadratic potential. Solving the SE with more appropriate potential description requires approximations such as the semiclassical Jeffreys–Wentzel–Kramers–Brillouin^40–43 (JWKB) method or the Fourier grid Hamiltonian⁴⁴ (FGH) approach.

In this work, we study the large amplitude motions of ATC and its methyl derivatives, evaluating the thermodynamic effects of different anharmonic corrections. We study torsional vibrations using various methods, including the FR approximation and the HR methods TES, MFG, E2DT, and MS-T. Additionally, we address the low frequency out-of-plane bends of Al₂X₆ structures with more accurate potential descriptions using the FGH method. Our focus is on entropy, the thermodynamic quantity most significantly affected by changes in partition functions, and the property most critical for understanding the behavior of MAO and related systems in catalytic applications. By systematically comparing these correction methods against experimental data where available, we aim to establish best practices for thermodynamic calculations of aluminum-containing systems with significant low-frequency modes. Ultimately, our approach forms the basis for more accurate modeling of MAO and its reactivity, which is required for rational design of improved catalyst systems.

Methodology

General approach to vibrational analysis

In standard quantum chemical calculations, vibrational analysis typically employs the HA, which models each vibrational mode as a quantum harmonic oscillator (QHO) with quadratic potential energy. The vibrational contribution to thermodynamic properties is calculated through partition functions. For a QHO, the vibrational partition function⁴⁵ iswhere h is the Planck constant, ν is vibrational frequency, and β = (kT)⁻¹, where k is the Boltzmann constant and T is temperature. Vibrational entropy⁴⁵ is generally calculated aswhere R is the gas constant and partition function q depends on the approximation used. For torsions, q corresponds to the symmetry-corrected partition function. In this work, q represents the partition function of a single vibrational mode (or of two coupled rotors in the E2DT case), and Q denotes the total vibrational (or rotational) partition function.

While the harmonic potential provides a satisfactory description for high-frequency, small-amplitude vibrations, it misrepresents low-frequency, large-amplitude motions such as torsions and out-of-plane bends, leading to incorrect vibrational energy levels. Additionally, the near-zero frequencies of these vibrations are highly susceptible to numerical noise in quantum chemical calculations: small absolute differences in frequencies between two methods can have large thermodynamic impacts.

The qh approximations address both problems simultaneously by treating vibrations below a frequency cutoff differently from the HA—Grimme's approach interpolates between harmonic and FR entropies, while Truhlar's approach raises low frequencies to the cutoff value. These procedures remove the numerical noise and reduce the overestimation of energy levels by the HA, which is especially the case for very low-frequency torsions.

Still, the accuracy of the qh approximations is limited by their black box nature—obtaining more accurate results requires methods that are based on solutions of the SE with improved potential description. In this work, we evaluate the effect of such methods for large-amplitude vibrations of simple aluminum compounds by comparison to harmonic, quasi-harmonic, and experimental results. The studied compounds, specific methods, and computational details are described in the following sections.

Studied molecules

We investigate ATC and its methyl derivatives in both monomeric and dimeric form. These compounds are characterized by a trigonal planar geometry in the monomer state, with a central aluminum atom bonded to three substituents (Cl or CH₃). The monomers examined are AlCl₃, AlMeCl₂, AlMe₂Cl, and AlMe₃, representing progressive methylation of the parent ATC compound. The corresponding dimers are Al₂Cl₆, Al₂MeCl₅, Al₂Me₂Cl₄, Al₂Me₃Cl₃, Al₂Me₄Cl₂, and Al₂Me₆. The dimerization occurs through bridging substituents, with chlorine bridges being energetically preferred over methyl bridges^46,47 in the mixed alkyl-halide species. In this work, we focus exclusively on chlorine-bridged structures for all mixed species; this approach restricts the number of studied isomers and explains the absence of Al₂Me₅Cl from our series. The only exception to our chlorine-bridge focus is Al₂Me₆, which necessarily contains methyl bridges. The dimer with two methyl groups (Al₂Me₂Cl₄) exists in three distinct isomers shown in Fig. 2.

Fig. 2 Three studied isomers of Al₂Me₂Cl₄. Left: 1,1-Al₂Me₂Cl₄ (C_2v). Middle: 1,2-Al₂Me₂Cl₄ (C_2v). Right: 1,2-Al₂Me₂Cl₄ (C_2h). 1,1 and 1,2 refer to methyl groups that are bonded to the same Al and to different Als, respectively.

This systematic series of compounds allows us to evaluate the thermodynamic impact of increasing methyl rotor count for both monomer and dimer species. The dimers also enable the investigation of low-frequency out-of-plane bending vibrations.

Hindered rotor corrections

Pitzer and Gwinn²⁴ formulated an expression for the HR partition function q^PG by scaling the classical partition function q_class by a correction factor based on the quantum and classical harmonic oscillator (CHO) partition functions (eqn (4)). The QHO partition function was already described in the previous section, while the CHO partition function is simply the inverse of the dimensionless frequency (eqn (5)).

q_CHO = (hνβ)⁻¹ (5)

Including an integral over the torsional coordinate with the potential energy description, the classical partition function^24,30 is computed as

where I is the torsional moment of inertia and ℏ is the reduced Planck constant. For symmetric HRs, this integral simplifies as the moment of inertia becomes constant, resulting inwhere the pre-exponential term equals the FR partition function and J₀(βV₀/2) is a modified zeroth order Bessel function. This formulation provides a smooth transition between the FR behavior at high temperatures or low barriers and QHO behavior at low temperatures or high barriers.

McClurg, Flagan, and Goddard^28,29 modified the PG approximation to correct for the zero-point energy overestimation of the HA:

where the term preceding β in the exponent is a Padé approximant.

For systems with multiple torsional degrees of freedom, Truhlar et al.^33–35 developed the MS-T method. The MS-T partition function incorporates rotational and vibrational contributions across multiple conformations:

where Q^rot_j is the rotational partition function of conformation j, U_j is the potential energy of conformation j relative to the lowest energy conformation, Q^QHO_j is the QHO partition function of conformation j, and f_j,η is the interpolating factor including kinetic coupling of torsional motions η by the method of Kilpatrick and Pitzer.⁴⁸

While the interpolating PG approximation and its derivatives offer straightforward approaches, more accurate results for HRs can be obtained by improving the potential energy description. The TES method²⁵ employs a higher-order Fourier series to represent the torsional potential:

where b_k are the Fourier coefficients and n_k are the corresponding integer periodicities. The SE is solved with this potential by diagonalizing the Hamiltonian matrix in a plane wave basis, yielding energy eigenvalues for constructing accurate partition functions bywhere ε_i is the energy of state i.

For systems with two coupled torsional degrees of freedom, the E2DT method³² extends the Fourier potential description into two dimensions using both cosine and sine terms to capture the potential energy coupling between the two torsional coordinates. The two-dimensional SE is then solved with this potential energy function, also incorporating coupled kinetic energy terms.⁴⁸ The resulting torsional eigenvalues are used to construct the two-dimensional non-separable (2D-NS) partition function³¹ that includes the quantum effects of the torsional modes. The E2DT method then uses the ratio of quantum and classical partition functions to correct the extended hindered rotor (EHR) partition function developed by Vansteenkiste et al.⁴⁹

The classical torsional partition function is defined as

where ϕ₁ and ϕ₂ are the torsional angles of the coupled HRs, σ_tor is the torsional symmetry number and D is the Kilpatrick–Pitzer torsional moment of inertia matrix.⁴⁸

Out-of-plane bending corrections

For the out-of-plane bending vibrations of our dimer structures, we utilized the FGH method⁴⁴ developed by Marston and Balint–Kurti. Unlike hindered rotation, which has established correction schemes, approaches for out-of-plane bending are less standardized. The FGH method solves the SE for arbitrary potential using a grid-based approach derived from discrete Fourier transform algorithms. The solved eigenvalues determine the partition function according to eqn (11).

In addition to potential energy, the FGH method requires the kinetic part of the Hamiltonian. For the description of kinetic energy, two options exist: using mass-scaled coordinates with reduced mass or torsional coordinates with moment of inertia. While previous studies on bending vibrations (such as in cyclopentane) have utilized mass-scaled coordinates,³⁸ we chose to use torsional coordinates instead. We justify our choice by the appearance of the bending vibrations of our dimer structures: these vibrations mainly comprise torsional motion along the Al–Cl–Cl–Al dihedral (see Fig. 3). Therefore, compared to mass-scaled coordinates, torsional coordinates provide a more accurate description of kinetic energy.

Fig. 3 Low frequency out-of-plane bending normal mode of Al₂Cl₆. Al in pink, Cl in green.

While using a constant moment of inertia is justified for the hindered rotation of a methyl rotor, these bending vibrations impose a more significant angular dependence on the moment of inertia. For HRs, Ellingson et al.²⁵ presented an averaged approach for the moment of inertia by integrating the Pitzer moment of inertia⁵⁰ over the torsional coordinate. We adjusted their approach to our case by Boltzmann-weighting the Pitzer moment of inertia along the torsional coordinate based on potential energy. This Boltzmann-weighting approach, imposing temperature dependence on the moment of inertia and energy levels, was originally proposed by Reinisch et al.⁵¹

When deciding between rigid and relaxed scans for the dimer TPESs, we opted for the rigid scan approach. While relaxed scans are utilized in HR methods like E2DT, they lead to significant deformations of the rotating groups during bending motions. These deformations would conflict with our chosen kinetic energy approach. Ultimately, we view the choice between relaxed and rigid scans as a trade-off between the accuracy of potential and kinetic energy descriptions.

The bending normal mode of TMA dimer (Al₂Me₆) required special treatment since its bridging structure involves methyl groups rather than chlorine atoms. In this case, the hydrogens of the bridging methyl groups do not rotate with the AlMe₂ groups. To exclude these hydrogens from the torsional moment of inertia calculation as a rotating part, we replaced them with dummy atoms having a mass three times the atomic mass of hydrogen, positioned on the torsional axis at their average distance. This approach provides the correct rotating groups for the torsional moments of inertia while minimizing changes to the overall moments of inertia of the molecule.

Computational details

Geometry optimizations, frequency calculations, and scan calculations of this work were performed at MP2/def2-TZVPP⁵² level of theory using Gaussian 16.⁵³ Tight criteria were used for both optimizations and SCF convergence.

For the last decade, the M06-2X⁵⁴ DFT functional has been the method of choice for aluminum compounds, largely due to the benchmark study of Ehm et al.¹⁷ They found the functional to provide the lowest mean absolute deviation (0.6 kcal mol⁻¹) for the dimerization energy of aluminum alkyls and analogous structures. However, we chose MP2, which also provided sufficient values (1.2 kcal mol⁻¹), based on our recent findings related to MAO structures. We found that M06-2X artificially stabilizes certain structural features in organoaluminum systems,¹⁹ while MP2 energies better align with the reference CCSD(T) results. Therefore, this ab initio approach ensures reliable energetics and direct comparability with our broader research on MAO-related systems.

All optimized structures were confirmed as minima through vibrational frequency analysis. For the hindered rotation TPES scan calculations, we performed relaxed optimizations with a step size of 5 degrees. Rigid TPES scans with the same step size were conducted for the out-of-plane bending vibrations of the Al₂X₆ dimer structures.

We used MSTor 2023 program⁵⁵ to calculate HR corrections for the TES and MS-T methods. To perform the MS-T calculations, we had to modify the covalent radius value of aluminum in the source code; initially the program misunderstood the bonding of our dimer structures. Changes were also made to the TES code, so that the program prints more than one energy eigenvalue (see SI for code modifications). The MFG, E2DT and FGH calculations were performed using Gaussian 16,⁵³ Q2DTor,⁵⁶ and FGH1D⁵⁷ programs, respectively. In the FGH calculations, we used 500 grid points. Moments of inertia for the FR and FGH calculations were determined with MSTor 2023⁵⁵ (see Tables S6 and S7 for the calculated moment of inertia values).

Throughout the thermodynamic calculations of this work, we applied a scaling factor of 0.9749 to the vibrational frequencies. This value was determined for our method and basis set with the FREQ program⁵⁸ by minimizing the root mean square error of zero-point energy in 15 test molecules. Quasi-harmonic entropy corrections were addressed using the Goodvibes script⁵⁹ with a cut-off frequency of 100 cm⁻¹.

Potential energy surfaces

The calculated potential energy surfaces are available in numerical form in Tables S3–S5.

Hindered rotation potential energy surfaces

We discovered that the TPESs of aluminum alkyl compounds are notably flat compared to typical organic molecules. This characteristic stems from the longer Al–C bonds compared to C–C bonds in hydrocarbons, providing more space for internal rotation. The simplest case, AlMeCl₂ monomer, exhibits an extremely low barrier to internal rotation of only 19 J mol⁻¹. In contrast, a single methyl rotor in Al₂MeCl₅ dimer possesses a greater barrier height of approximately 3 kJ mol⁻¹ due to more limited space for internal rotation.

To put these barrier heights in context, we can compare them to thermal energy (RT), which equals approximately 2.5 kJ mol⁻¹ at room temperature. When RT ≫ V₀, the motion is effectively free rotation. When RT ∼ V₀, the motion is properly characterized as hindered rotation. Therefore, the torsional vibrations of AlMeCl₂ and Al₂MeCl₅ categorize as free rotation and hindered rotation, respectively.

The introduction of a second methyl rotor impacts the potential energy landscapes by imposing saddle points. In a two-dimensional TPES, the saddle points correspond to the restricting rotational barriers of a single methyl rotor, while the global maxima describe the barriers related to combined rotation of the methyl groups. Fig. 4 illustrates the two-dimensional TPES of the 1,2-Al₂Me₂Cl₄ (C_2h) isomer (see Fig. S1–S3 for other two-dimensional TPESs).

Fig. 4 Two-dimensional torsional potential energy surface of 1,2-Al₂Me₂Cl₄ (C_2h) isomer as a function of methyl torsional angles ϕ₁ and ϕ₂.

Effectively, the relation between saddle points and RT determines the manifestation of torsional motion for these two-dimensional cases. Although the rotational barriers of AlMe₂Cl are sixfold compared to AlMeCl₂, the rotation remains relatively free. The rotations of Al₂Me₂Cl₄ dimers are hindered, as the rotational barriers are comparable to those of Al₂MeCl₅. Combined rotations of the methyl groups encounter barriers approximately double those of the single-rotor case, with the 1,1-isomer showing slightly higher barriers due to steric interactions between adjacent methyl groups.

For structures with more than two methyl rotors (AlMe₃, Al₂Me₃Cl₃, Al₂Me₄Cl₂, and Al₂Me₆), exhaustive mapping of the full TPESs becomes computationally prohibitive. However, barrier heights can be estimated by rotating selected methyl groups to energetically unfavorable positions, with first-order saddle points providing the crucial information. For AlMe₃, all rotational barriers remain below 0.5 kJ mol⁻¹, indicating relatively free rotation of the methyl groups. In the dimer structures, the additional methyl rotors have an insignificant effect on the first order saddle points over the two-dimensional case.

Out-of-plane bending potential energy surfaces

The out-of-plane bends in dimer structures exhibit significantly different potential energy characteristics compared to hindered rotation. Fig. 5 presents the potential energy curves for these bending modes along the torsional coordinates. Unlike the hindered rotation TPESs, these curves show high barriers that completely restrict the bending motion, resulting in well-defined vibrational behavior.

Fig. 5 Torsional potential energy surfaces for out-of-plane bending vibrations in dimer structures. Top: Symmetric potentials. Bottom: Asymmetric potentials.

Depending on the specific dimer structure, the TPES curves exhibit either symmetric or asymmetric profiles. The symmetric curves resemble higher-order monomials rather than simple quadratic potentials, while the asymmetric curves show more complex behavior. The width of these potential wells correlates inversely with the size of the terminal groups—bulkier methyl groups result in narrower potential wells due to steric constraints.

For Al₂Cl₆, shown in Fig. 6, the potential energy increases gradually as the molecule bends away from its equilibrium position. This gradual increase at small displacements is followed by a steeper increase at larger displacements. The HA captures the gradual increase at small displacements but fails to estimate the steepness of the potential at large displacements.

Fig. 6 Potential energy surface of Al₂Cl₆ along the out-of-plane bending torsional coordinate compared to the harmonic approximation (HA).

Energy levels

The calculated eigenvalues with TES, E2DT, and FGH methods are available in numerical form in Tables S9–S11, respectively.

Hindered rotation energy levels

The calculated energy levels for one-dimensional torsional motions reveal significant differences between the monomer (AlMeCl₂) and the dimer (Al₂MeCl₅). For these species, Fig. 7 presents a comparison of the energy levels obtained with three approaches: FR, HA, and TES.

Fig. 7 Energy levels of AlMeCl₂ (left) and Al₂MeCl₅ (right) methyl rotations using the free rotor (FR) approximation, torsional eigenvalue summation (TES), and harmonic approximation (HA). The horizontal energy level lines use solid, dashed, and dash-dot styles to indicate degeneracies of 1, 2, and 3, respectively.

For the monomer, the TES energy levels closely align with the FR energy levels, confirming our characterization of this rotor as essentially free. Both methods share a ground state near zero frequency, followed by diverging energy levels of double degeneracy. Due to the lack of restricting potential energy barriers, the HA fails to estimate the energy levels correctly, predicting unreasonably tight energy level spacing.

In contrast, for the corresponding dimer, the harmonic potential provides a reasonably good fit to the true potential near the minimum, resulting in accurate predictions of the two lowest energy levels. However, as we move to higher energy states, the harmonic energy levels diverge increasingly from the TES states, which adopt the FR behavior.

For structures with two torsional degrees of freedom, interpretation of energy eigenvalues in comparison with the HA becomes challenging. As the HA energy levels originate from two separate one-dimensional potentials, straightforward comparison to the 2D-NS energy levels of the E2DT is practically impossible.

Out-of-plane bending energy levels

The energy levels for out-of-plane bending vibrations display a distinct pattern that differs from both harmonic oscillators and HRs. Fig. 8 presents the seven lowest energy levels for Al₂Cl₆ obtained with the HA and FGH. The FGH energy levels, which account for the true shape of the potential energy curve, exhibit gradually increasing spacing between adjacent levels as energy increases. This behavior, arising from the steepness of the potential energy curve at higher displacements, contrasts with the perfectly even spacing of the HA energy levels.

Fig. 8 Lowest energy levels for Al₂Cl₆ out-of-plane bending calculated using the harmonic approximation (HA) and Fourier grid Hamiltonian (FGH) methods.

This gradually increasing separation between energy levels was also reported in a recent study by Wu et al.³⁸ for the bending vibration of cyclopentane using a quadratic–quartic potential. Although our coordinate system differs (we use torsional angles rather than mass-scaled coordinates), the qualitative behavior is remarkably similar.

Thermodynamics

We have systematically compared the entropy values calculated with different theoretical approaches to understand the impact of anharmonic corrections on thermodynamic properties. Table 1 presents entropy values (in J K⁻¹ mol⁻¹) at standard conditions (298.15 K, 1 bar) for all studied structures using various methods. The harmonic vibrational frequencies are presented in Tables S1 and S2, torsional and out-of-plane bending normal modes being explicitly marked.

Table 1 Gas phase entropy values^a for aluminum compounds calculated using different theoretical methods^b, with experimental values^c where available

Structure	Exp.	HA	qh (Gri.)^d	qh (Tru.)^d	FR	TES	MFG	E2DT	MS-T	FGH
a Entropy values in J K^−1 mol^−1 at 298.15 K and 1 bar.b Scaling factor of 0.9749 was applied for the vibrational frequencies, excluding the TES, 2D-NS, and FGH eigenvalues.c Experimental values derived from literature data as described in the Experimental comparison section.d Quasi-harmonic approximations of Truhlar (Tru.) and Grimme (Gri.).
AlCl3	310.2 ± 5.7	313.6	313.2	313.6	—	—	—	—	—	—
AlMeCl2	—	354.5	344.4	340.1	334.5	334.5	340.5	—	320.1	—
AlMe2Cl	—	368.6	352.0	345.9	346.3	—	355.2	346.3	346.4	—
AlMe3	352.8 ± 0.9	381.9	360.7	354.6	355.0	—	366.7	—	355.1	—
Al2Cl6	471.4 ± 11.6	476.0	458.2	456.2	—	—	—	—	—	472.3
Al2MeCl5	—	501.1	482.9	481.6	501.9	501.3	501.7	—	501.6	496.3
1,1-Al2Me2Cl4 (C2v)	—	511.5	491.4	489.1	512.8	—	512.7	509.1	512.2	502.6
1,2-Al2Me2Cl4 (C2v)	—	507.1	489.3	489.5	508.6	—	508.1	506.2	507.8	502.9
1,2-Al2Me2Cl4 (C2h)	—	512.2	491.7	489.4	513.6	—	513.2	509.7	513.0	503.4
Al2Me3Cl3	—	524.6	506.8	508.1	526.7	—	526.1	—	525.5	520.4
Al2Me4Cl2	—	529.9	509.8	509.9	532.7	—	531.8	—	531.1	520.5
Al2Me6	524.9 ± 1.0	526.0	515.1	520.6	529.6	—	528.0	—	527.0	522.6

---

Monomer results

The anharmonic corrections have a particularly pronounced effect on the entropy of our monomer structures with methyl rotors. In AlMeCl₂, the HA significantly overestimates the entropy (354.5 J K⁻¹ mol⁻¹) compared to more appropriate approaches like the TES method and FR approximation, both of which estimate the entropy at 334.5 J K⁻¹ mol⁻¹. This 20 J K⁻¹ mol⁻¹ difference highlights how poorly the HA captures the behavior of near-free rotors. Interestingly, the estimates of interpolating HR treatments, MFG (340.5 J K⁻¹ mol⁻¹) and MS-T (320.1 J K⁻¹ mol⁻¹), significantly deviate from both the TES value and each other.

For one-dimensional hindered rotation, these two methods should yield matching results. We observed that the entropy underestimation of the MS-T seems to be due to severely overestimated rotational barriers originating from misinterpretation of the Hessian matrix. This issue most likely arises from the extremely low force constant of AlMeCl₂ estimated by MP2; M06-2X, with slightly higher force constant, yields a reasonable value (334.6 J K⁻¹ mol⁻¹) that aligns with the TES value. The rotational barriers estimated by the MFG method were also slightly overestimated, which might contribute to the discrepancy between the MFG and TES methods.

The monomers with two (AlMe₂Cl) and three (AlMe₃) methyl rotors follow similar patterns. For AlMe₂Cl, the HA (368.6 J K⁻¹ mol⁻¹) overestimates the entropy by approximately 22 J K⁻¹ mol⁻¹ compared to our best estimate, the E2DT (346.3 J K⁻¹ mol⁻¹). Remarkably, both the MS-T (346.4 J K⁻¹ mol⁻¹) and FR (346.3 J K⁻¹ mol⁻¹) results are in excellent agreement with the E2DT result. Moreover, the excellent agreement of the MS-T and FR extends to AlMe₃, for which the overestimation of entropy by the HA (26.8 J K⁻¹ mol⁻¹) is the most severe. This remarkable alignment of the FR with the TES, E2DT, and MS-T has two indications: the monomer methyl rotors can be considered free, and more importantly, their coupling effects seem to play an insignificant role at these conditions. Despite this observation, the uncoupled MFG method provides substandard corrections for two and three methyl rotors, as in the one-dimensional case. Therefore, the insufficient results appear to be related to the treatment of low-barrier systems by the MFG method rather than coupling effects between rotors.

The qh approximations also improve upon the HA for monomers, with Truhlar's approach providing values closer to the reference methods than Grimme's approach. This order is somewhat surprising, since Grimme's rotation approach would be more justified, physically. However, we identified that the discrepancy between Grimme's qh approach and the reference methods arises from two factors: estimated moments of inertia deviating from the actual values and missing symmetry corrections. Ultimately, neither qh method fully captures the corrections provided by explicit anharmonic treatments.

Dimer results

The thermodynamic behavior of the dimers notably differs from the monomers. For the dimer with one methyl rotor (Al₂MeCl₅), the HA value (501.1 J K⁻¹ mol⁻¹) is remarkably close to the FR value (501.9 J K⁻¹ mol⁻¹), while the values of the HR methods, TES, MFG, and MS-T, fall between these two extremes. Around the 100 cm⁻¹ frequency range, at room temperature, the FR partition function is nearly twice the HA partition function, but the difference in entropy is small, as it is compensated by the derivative term of eqn (1). Due to the alignment of the HA and FR entropy values, the HR corrections to entropy are rather modest at room temperature, despite significant differences in the partition functions themselves.

For the three isomers of Al₂Me₂Cl₄, the HR and FR deviations from the HA follow a nearly similar pattern. The FR entropy values, on average, lie 1.4 J K⁻¹ mol⁻¹ above the HA values, while the MFG and MS-T values settle between the extremes. Although the E2DT and MS-T values align perfectly in the monomer case, for these dimers the E2DT values are systematically lower than the MS-T values (2.7 J K⁻¹ mol⁻¹ on average).

The trend between the HA and the rotor approximations stays the same with three (Al₂Me₃Cl₃), four (Al₂Me₄Cl₂), and six (Al₂Me₆) methyl rotors. For practical reasons, we were unable to address the hindered rotation of the two TMA bridging methyl groups due to the limitations of our chosen methods. Both the MFG method and MS-T require torsional axes to be determined between real atoms, while the true torsional axes of the bridging methyl groups point towards the center of the imaginary Al–Al bond. For all these dimer species, the MS-T values slightly exceed the HA values unlike in the monomer cases. For Al₂Me₆, the slight overestimation of FR for each terminal methyl group cumulates to entropy value 2.6 J K⁻¹ mol⁻¹ above the MS-T value. The uncoupled MFG method yields satisfactory results for the dimers as the values fall between the MS-T and FR values.

The out-of-plane bending modes of the dimers represent a distinct contribution to thermodynamic properties. Our FGH calculations consistently yield lower entropy values than the HA, with differences ranging from 3.4 J K⁻¹ mol⁻¹ for Al₂Me₆ to 9.4 J K⁻¹ mol⁻¹ for Al₂Me₄Cl₂. The lowering of entropy is generally caused by the increasing energy level separation of FGH, but the magnitude of the lowering also depends on the location of the FGH ground state relative to the HA counterpart. We observed that the FGH ground states tend to lie above the HA ground states, with the largest difference emerging for Al₂Me₆. This behavior of Al₂Me₆ is expected given that the out-of-plane bending motion in the methyl-bridged species is more complex than in the chlorine-bridged species, where the motion is more localized to the Al–Cl–Cl–Al dihedral. The qh results deviate even more from the FGH results, as the low frequencies are raised to unreasonably high values.

Experimental comparison

To validate our computational approaches, comparison to experimental results is necessary. In our case, experimental entropy data are available for ATC and TMA. For TMA, the gas phase entropy value of the dimer⁶⁰ (524.9 ± 1.0 J K⁻¹ mol⁻¹) can be used together with the entropy of dimerization^61,62 (−180.7 ± 1.4 J K⁻¹ mol⁻¹) to obtain the corresponding value for the monomer. For ATC, we calculated the gas phase entropy values for both forms using solid phase data of the monomer⁶³ (109.3 ± 0.4 J K⁻¹ mol⁻¹), enthalpy of sublimation⁶⁴ (59.9 ± 2.0 kJ mol⁻¹), and entropy of dimerization⁶⁵ (−149.0 ± 2.1 J K⁻¹ mol⁻¹). From the relatively recent work of Brunetti et al.,⁶⁴ we selected the second-law enthalpy of sublimation value that provides better thermodynamic consistency with our computational results. With their third-law enthalpy value, the monomer entropy would lie well below the HA value, despite AlCl₃ lacking problematic low-frequency modes.

As the experimental values are constructed from several parts, propagation of uncertainty is required. For TMA dimer in gas phase, McCullough et al. reported an uncertainty of 1.0 J K⁻¹ mol⁻¹ based on low-temperature calorimetry and vapor pressure studies.⁶⁰ However, their seemingly optimistic estimate includes no individual contributions to the total value. For TMA monomer, we approximated the uncertainty of McCullough et al. as standard uncertainty (one Sigma) and combined it with the standard error (intercept) of Henrickson and Eyman⁶¹ Van’t Hoff plot, assuming the variables to be uncorrelated. In the case of ATC, we utilized the same approach for the three uncertainties. We converted the conservative enthalpy of sublimation uncertainty estimate of Brunetti et al. to entropy by division by sublimation point temperature (350 K).⁶⁴ This uncertainty (5.7 J K⁻¹ mol⁻¹) is significantly higher than the other two, thus practically determining the overall uncertainty for both monomer and dimer.

The calculated experimental gas phase entropy values and uncertainties for monomeric and dimeric forms of ATC and TMA were presented earlier in Table 1 together with our computational values. By using the most suitable approximations for each structure, we obtained values that align well with the experimental values, noting that for ATC we used the second-law enthalpy of sublimation from Brunetti et al.⁶⁴

For ATC monomer, where no significant anharmonic corrections are required due to lack of large amplitude motions, both our HA value of AlCl₃ (313.6 J K⁻¹ mol⁻¹) and the HA value of Konings and Booij⁶⁶ (312.7 J K⁻¹ mol⁻¹), based on experimental frequency and geometry data, are slightly higher than the experimental value (310.2 ± 5.7 J K⁻¹ mol⁻¹). Also, for the ATC dimer, the HA (476.0 J K⁻¹ mol⁻¹) overestimates the entropy compared to the experimental value (471.4 ± 11.6 J K⁻¹ mol⁻¹), while the harmonic value is improved (472.3 J K⁻¹ mol⁻¹) by our anharmonic FGH correction for the low-frequency out-of-plane bending normal mode. However, the HA result of Konings and Booij for the dimer lies slightly lower (470.2 J K⁻¹ mol⁻¹) due to higher frequency of the bending vibration (25 cm⁻¹) compared to our MP2 estimate (15 cm⁻¹).

For TMA, our computational results are in excellent agreement with the experimental values using the MS-T and FGH methods. The HR corrections prove crucial for the monomer, as the HA substantially overestimates the entropy (381.9 J K⁻¹ mol⁻¹) compared to the experimental value (352.8 ± 0.9 J K⁻¹ mol⁻¹). The MS-T treatment reduces the error significantly (355.1 J K⁻¹ mol⁻¹), although a similar result (355.0 J K⁻¹ mol⁻¹) is obtained by the simple FR approximation. For the TMA dimer, the HA value (526.0 J K⁻¹ mol⁻¹) is, surprisingly, in excellent agreement with the experimental value (524.9 ± 1.0 J K⁻¹ mol⁻¹). Still, a result of comparable accuracy (523.6 J K⁻¹ mol⁻¹) is achieved when the FGH correction (−3.4 J K⁻¹ mol⁻¹) for the out-of-plane bend is added to the MS-T result (527.0 J K⁻¹ mol⁻¹). The excellent alignment of the FR and HA results with the experimental values of the TMA monomer and dimer, respectively, indicate that simple methods can yield good numerical results for some systems. However, more robust approaches provide greater confidence and consistency across different molecular systems.

Temperature dependence

Although the standard conditions are in many aspects a great setting for benchmarking, they merely provide a snapshot of the overall thermodynamic behavior. As temperature increases, higher energy states become thermally accessible, significantly contributing to partition functions and thermodynamic properties. For this reason, the use of the HA becomes problematic, as it tends to overestimate the density of energy levels for large amplitude motions. This overestimation then transfers to partition functions and entropy, which is exactly the case for the hindered rotation and out-of-plane-bending of our structures.

For monomer hindered rotation and dimer out-of-plane bending, HA entropy overestimation already applies at low temperatures and further increases with rising temperature. However, for dimer hindered rotation (Al₂MeCl₅), the behavior is somewhat different, as seen in Fig. 9. The FR partition function and entropy lie above the HA and TES values at low temperature. As temperature rises, the TES entropy approaches the FR entropy, while the overestimation of the HA becomes more significant above 400 K. The reasonable agreement between the HA and TES below this temperature proves useful for many applications.

Fig. 9 Temperature dependence on the Al₂MeCl₅ partition functions and entropy values obtained with the harmonic approximation (HA), free rotor (FR) approximation, and torsional eigenvalue summation (TES). The values correspond to the contribution of the torsional mode.

For two-dimensional hindered rotation of the dimers, the thermodynamic behavior of the HR methods (E2DT and MS-T) resembles that of TES for Al₂MeCl₅. We noted earlier that the dimer entropy values obtained by the MS-T lie slightly above the E2DT values at standard conditions. Here, we observed that this trend applies throughout the temperature range. This consistency between the E2DT and MS-T indicates that the MS-T is a reliable alternative to the E2DT and the best available option moving forward to multi-dimensional hindered rotation.

Conclusions

In this work, we systematically investigated anharmonic corrections for low-frequency vibrations in aluminum trichloride (ATC) and its methyl derivatives, with particular focus on entropy. By applying state-of-the-art methods for hindered rotation (TES, E2DT, MS-T) and out-of-plane bending (FGH), we obtained results that are in excellent agreement with the experimental values and identified crucial structure–property relationships that impact the computational accuracy of these compounds. These findings also provide valuable insights into the thermodynamic characterization of methylaluminoxane (MAO).

The most striking observation is the distinct behavior of low-frequency torsional modes in methyl containing monomers versus dimers. In monomers, the barriers restricting methyl rotations are extremely low, leading to practically free rotor behavior poorly captured by the harmonic approximation (HA). This mischaracterization leads to substantial entropy overestimation by the HA, which is only partly corrected by the quasi-harmonic (qh) approximations. In contrast, the same methyl groups in dimers experience greater steric constraints, resulting in higher rotational barriers that align better with the HA, particularly at standard conditions.

For out-of-plane bending modes in dimers, our improved potential energy descriptions, featuring a steeper rise at larger displacements, result in a characteristic pattern of increasing energy level spacing inadequately captured by both the HA and qh approximations. While the HA overestimates the energy level density and thus the entropy for these modes, the qh approximations underestimate the entropy by raising the frequencies to artificially high values.

Our observations on the torsional and bending vibrations prove critical for the large number of complex MAO low-frequency vibrations, which exhibit both motional motifs. These vibrations, especially those of larger MAO sheets, comprise bending and twisting of the sheet frame, accompanied by more localized motions such as bending of the edge sites reminiscent of trimethylaluminum (TMA) dimers and methyl rotations. For the vibrations without methyl rotations, which also possess the lowest frequencies, we believe that the problem of the harmonic potential is the same as for the out-of-plane bends of the present work: it increases too slowly at larger displacements. While the entropy deviations of the harmonic and quasi-harmonic approximations are relatively minor for a single mode, they become cumulatively significant when considering all such modes. We believe that improved results could be obtained by a simple but well-parametrized mathematical model, which would mimic the diverging energy levels of the out-of-plane bends observed in the present work. This approach would also correct inaccuracies in harmonic ground states, which can cause substantial energy state overcounting.

Looking forward, we plan to extend our anharmonic analysis to small-scale MAO structures, formulating a new mode specific approach to correct the large-displacement behavior of the HA for low frequency modes with predominant bending character. For the MAO–catalyst interaction, we consider our vibrational characterization of the mixed alkyl-halide species crucial, as chloride abstraction from the precatalyst can create structurally analogous local environments. In this context, finding a thermodynamically sound solvation model for aluminum alkyls is critical, given the solution-phase environment of the catalytic application.

By establishing a more accurate framework for thermodynamic calculations of aluminum systems, this work provides a foundation for improved computational modeling of MAO–catalyst interactions, addressing a long-standing challenge in computational chemistry of polyolefin catalysis.

Author contributions

P. H.: conceptualization, formal analysis, investigation, visualization, writing – original draft. M. L.: conceptualization, funding acquisition, project administration, supervision, writing – review & editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the SI. See DOI: https://doi.org/10.1039/d5cp02487k

Acknowledgements

P. H. acknowledges the support of the Research Council of Finland, decision 357509. M. L. acknowledges the Academy of Finland Flagship Programme, Photonics Research and Innovation (PREIN), decision 320166. The authors acknowledge CSC – IT Center for Science, Finland, and Finnish Grid and Cloud Infrastructure (urn:nbn:fi:research-infras-2016072533) for computational resources. The authors thank Dr Scott Collins (University of Eastern Finland) for helpful discussions regarding the experimental data.

References

1. T. S. Suntola, A. J. Pakkala and S. G. Lindfors, US Pat., 4389973A, 1983.
2. R. L. Puurunen, J. Appl. Phys., 2005, 97, 121301.
3. K. Kukli, M. Ritala, M. Leskelä and J. Jokinen, J. Vac. Sci. Technol., A, 1997, 15, 2214–2218.
4. R. Matero, A. Rahtu and M. Ritala, Langmuir, 2005, 21, 3498–3502.
5. S. Li, Y. Bao, M. Laitinen, T. Sajavaara, M. Putkonen and H. Savin, Phys. Status Solidi A, 2015, 212, 1795–1799.
6. Y. Bao, M. Laitinen, T. Sajavaara and H. Savin, Adv. Electron. Mater., 2017, 3, 1600491.
7. H. Saare, G. Dianat and G. N. Parsons, J. Phys. Chem. C, 2022, 126, 7036–7046.
8. S. M. George, Chem. Rev., 2010, 110, 111–131.
9. H. S. Zijlstra and S. Harder, Eur. J. Inorg. Chem., 2015, 19–43.
10. A. Joshi, H. S. Zijlstra, E. Liles, C. Concepcion, M. Linnolahti and J. S. Mcindoe, Chem. Sci., 2021, 12, 546–551.
11. A. Joshi, S. Collins, M. Linnolahti, H. S. Zijlstra, E. Liles and J. S. McIndoe, Chem. – Eur. J., 2021, 27, 8753–8763.
12. S. Collins, A. Joshi and M. Linnolahti, Chem. – Eur. J., 2021, 27, 15460–15471.
13. S. Collins, G. Hasan, A. Joshi, J. S. McIndoe and M. Linnolahti, ChemPhysChem, 2021, 22, 1326–1335.
14. L. Luo, J. M. Younker and A. V. Zabula, Science, 2024, 384, 1424–1428.
15. M. Linnolahti and S. Collins, ChemPhysChem, 2017, 18, 3369–3374.
16. S. Collins and M. Linnolahti, ChemPhysChem, 2023, 24, e202300342.
17. C. Ehm, G. Antinucci, P. H. M. Budzelaar and V. Busico, J. Organomet. Chem., 2014, 772, 161–171.
18. M. Linnolahti and S. Collins, Dalton Trans., 2022, 51, 11152–11162.
19. X. Mao, M. Bharti, S. Collins and M. Linnolahti, Chem. – Eur. J., 2025, 31, e202404642.
20. M. Ystenes, J. L. Eilertsen, J. Liu, M. Ott, E. Rytter and J. A. Støvneng, J. Polym. Sci., Part A:Polym. Chem., 2000, 38, 3106–3127.
21. M. S. Kuklin, J. T. Hirvi, M. Bochmann and M. Linnolahti, Organometallics, 2015, 34, 3586–3597.
22. S. Grimme, Chem. – Eur. J., 2012, 18, 9955–9964.
23. R. F. Ribeiro, A. V. Marenich, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. B, 2011, 115, 14556–14562.
24. K. S. Pitzer and W. D. Gwinn, J. Chem. Phys., 1942, 10, 428–440.
25. B. A. Ellingson, V. A. Lynch, S. L. Mielke and D. G. Truhlar, J. Chem. Phys., 2006, 125, 084305.
26. E. Dzib and G. Merino, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2022, 12, e1583.
27. D. G. Truhlar, J. Comput. Chem., 1991, 12, 266–270.
28. R. B. McClurg, R. C. Flagan and W. A. Goddard, J. Chem. Phys., 1997, 106, 6675–6680.
29. R. B. McClurg, J. Chem. Phys., 1999, 111, 7163–7164.
30. P. Y. Ayala and H. B. Schlegel, J. Chem. Phys., 1998, 108, 2314–2325.
31. A. Fernández-Ramos, J. Chem. Phys., 2013, 138, 134112.
32. L. Simón-Carballido, J. L. Bao, T. V. Alves, R. Meana-Pañeda, D. G. Truhlar and A. Fernández-Ramos, J. Chem. Theory Comput., 2017, 13, 3478–3492.
33. J. Zheng, T. Yu, E. Papajak, I. M. Alecu, S. L. Mielke and D. G. Truhlar, Phys. Chem. Chem. Phys., 2011, 13, 10885–10907.
34. J. Zheng and D. G. Truhlar, J. Chem. Theory Comput., 2013, 9, 1356–1367.
35. W. Chen, P. Zhang, D. G. Truhlar, J. Zheng and X. Xu, J. Chem. Theory Comput., 2022, 18, 7671–7682.
36. Y. Gao, T. He, X. Li and X. You, Phys. Chem. Chem. Phys., 2019, 21, 1928–1936.
37. B. C. Garrett and D. G. Truhlar, J. Phys. Chem., 1979, 83, 1915–1924.
38. J. Wu, L. G. Gao, W. Ren and D. G. Truhlar, Chem. Sci., 2020, 11, 2511–2523.
39. W. Chen, X. Guo, L. Chen, R. Zhang, Y. Li, H. Feng, X. Xu and X. Zhang, Phys. Chem. Chem. Phys., 2021, 23, 7333–7342.
40. H. Jeffreys, Proc. London Math. Soc., 1925, s2–23, 428–436.
41. G. Wentzel, Z. Phys., 1926, 38, 518–529.
42. H. A. Kramers, Z. Phys., 1926, 39, 828–840.
43. L. Brillouin, C. R. Hebd. Seances Acad. Sci., 1926, 183, 24–26.
44. C. C. Marston and G. G. Balint-Kurti, J. Chem. Phys., 1989, 91, 3571–3576.
45. D. A. McQuarrie, Statistical Thermodynamics, Harper & Row, New York, NY, 1973.
46. A. Tarazona, E. Koglin, F. Buda, B. B. Coussens, J. Renkema, S. van Heel and R. J. Meier, J. Phys. Chem. B, 1997, 101, 4370–4378.
47. B. B. Champagne, D. H. Mosley, J. G. Fripiat, J.-M. André, A. Bernard, S. Bettonville, P. François and A. Momtaz, THEOCHEM, 1998, 454, 149–159.
48. J. E. Kilpatrick and K. S. Pitzer, J. Chem. Phys., 1949, 17, 1064–1075.
49. P. Vansteenkiste, D. Van Neck, V. Van Speybroeck and M. Waroquier, J. Chem. Phys., 2006, 124, 044314.
50. K. S. Pitzer, J. Chem. Phys., 1946, 14, 239–243.
51. G. Reinisch, J.-M. Leyssale and G. L. Vignoles, J. Chem. Phys., 2010, 133, 154112.
52. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
53. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16 (Revision A.03), Gaussian, Inc., Wallingford CT, 2016.
54. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 120, 215–241.
55. J. Zheng, W. Chen, S. L. Mielke, J. L. Bao, K. L. Clarkson, R. Meana-Pañeda, X. Xu and D. G. Truhlar, MSTor 2023, Tsinghua University, Beijing and University of Minnesota, Minneapolis, 2023.
56. D. Ferro-Costas, M. N. D. S. Cordeiro, D. G. Truhlar and A. Fernández-Ramos, Comput. Phys. Commun., 2018, 232, 190–205.
57. R. D. Johnson III, FGH1D Program for one-dimensional solution of the Schrödinger equation (Version 1.01), Gaithersburg, MD, 1999.
58. H. S. Yu, L. J. Fiedler, I. M. Alecu, S. Kanchanakungwankul and D. G. Truhlar, FREQ – version 2, University of Minnesota, Minneapolis, 2021.
59. G. Luchini, J. V. Alegre-Requena, I. Funes-Ardoiz and R. S. Paton, F1000Research, 2020, 9, 291.
60. J. P. Mccullough, J. F. Messerly, R. T. Moore and S. S. Todd, J. Phys. Chem., 1963, 67, 677–679.
61. C. H. Henrickson and D. P. Eyman, Inorg. Chem., 1967, 6, 1461–1465.
62. M. B. Smith, J. Organomet. Chem., 1972, 46, 31–49.
63. D. R. Stull and H. Prophet, JANAF thermochemical tables, National Bureau of Standards, Washington, DC, 2nd edn, 1971.
64. B. Brunetti, V. Piacente and P. Scardala, J. Chem. Eng. Data, 2009, 54, 940–944.
65. G. E. Vrieland and D. R. Stull, J. Chem. Eng. Data, 1967, 12, 532–535.
66. R. J. M. Konings and A. S. Booij, J. Chem. Thermodyn., 1992, 24, 1181–1188.

---