ACS Publications. Most Trusted. Most Cited. Most Read
My Activity
CONTENT TYPES

Figure 1Loading Img

Chemically Accurate Vibrational Free Energies of Adsorption from Density Functional Theory Molecular Dynamics: Alkanes in Zeolites

Cite this: J. Chem. Theory Comput. 2021, 17, 9, 5849–5862
Publication Date (Web):August 30, 2021
https://doi.org/10.1021/acs.jctc.1c00519

Copyright © 2021 The Authors. Published by American Chemical Society. This publication is licensed under

CC-BY-NC-ND 4.0.
  • Open Access

Article Views

3379

Altmetric

-

Citations

LEARN ABOUT THESE METRICS
PDF (1 MB)
Supporting Info (1)»

Abstract

We present a methodology to compute, at reduced computational cost, Gibbs free energies, enthalpies, and entropies of adsorption from molecular dynamics. We calculate vibrational partition functions from vibrational energies, which we obtain from the vibrational density of states by projection on the normal modes. The use of a set of well-chosen reference structures along the trajectories accounts for the anharmonicities of the modes. For the adsorption of methane, ethane, and propane in the H-CHA zeolite, we limit our treatment to a set of vibrational modes localized at the adsorption site (zeolitic OH group) and the alkane molecule interacting with it. Only two short trajectories (1–20 ps) are required to reach convergence (<1 kJ/mol) for the thermodynamic functions. The mean absolute deviations from the experimentally measured values are 2.6, 2.8, and 4.7 kJ/mol for the Gibbs free energy, the enthalpy, and the entropy term (−TΔS), respectively. In particular, the entropy terms show a major improvement compared to the harmonic approximation and almost reach the accuracy of the previous use of anharmonic frequencies obtained with curvilinear distortions of individual modes. The thermodynamic functions so obtained follow the trend of the experimental values for methane, ethane, and propane, and the Gibbs free energy of adsorption at experimental conditions is correctly predicted to change from positive for methane (5.9 kJ/mol) to negative for ethane (−4.8 kJ/mol) and propane (−7.1 kJ/mol).

This publication is licensed under

CC-BY-NC-ND 4.0.
  • cc licence
  • by licence
  • nc licence
  • nd licence

1. Introduction

ARTICLE SECTIONS
Jump To

The important role of entropy in surface chemistry has long been recognized. (1−3) Prominent examples are the confinement effects in zeolite catalysis such as alkane cracking and dehydrogenation. (3) The quantification of entropy effects on rate and equilibrium constants is crucial for the atomistic understanding of heterogenous catalysis. In turn, this requires the chemically accurate prediction of free energies for different elementary steps such as adsorption, chemical transformation, and desorption.
Statistical thermodynamics is a powerful tool to evaluate the changes in free energies in complex molecular systems from which rate and equilibrium constants are obtained. While the energy barriers and reaction energies for large chemical systems can be predicted ab initio with the desired accuracy and affordable computational cost, (4,5) the prediction of free energies and their variations is still challenging. Substantial progress has been made (3) since early approaches which used the harmonic approximation (6,7) and were often limited to the six rigid body coordinates of the adsorbed molecules relative to the surface. (1,2)
Sampling the potential energy surface (PES) by biased ab initio molecular dynamics (MD), for example, umbrella sampling and thermodynamic integration, or Monte Carlo (MC) simulations requires many steps, for example, long dynamic runs from 50 ps to even 1 ns, (8−15) with millions of energy and force calculations. (9,16−18) This is computationally extremely challenging even if density functional theory (DFT) with the more affordable but less accurate GGA-type functionals (GGA—generalized gradient approximation) is applied in ab initio MD (AIMD) or MC (AIMC) simulations. Metadynamics (19,20) allows one to force the sampling of rare events. However, even for optimal cases simulation times of 50–100 ps are still required for an accurate sampling of a single barrier. (8) Moreover, the computed quantity is the variation of the free energy. It is not straightforward and is extremely expensive to obtain, following this path, absolute free energies or to separate the enthalpic component from the entropic one. As a consequence, AIMD or AIMC simulations for surface problems have been completed only in exceptional cases. (10,18,21−26) While for molecular enthalpies, MD schemes have been proposed to approach chemical accuracy at reasonable computational cost, (27) the entropic contributions are more challenging because they require extensive sampling of the complete phase space with all degrees of freedom.
An alternative approach is to calculate free energies directly from the vibrational partition function of the system. However, standard harmonic calculations are not good enough: it has been shown that anharmonic contributions to vibrational partition functions are crucial for reaching chemical accuracy. (3,21,28) Piccini et al. (29,30) have successfully developed and implemented an approach to predict anharmonic frequencies by finite distortion along the normal modes. (28,29) To minimize the anharmonic couplings between the one-dimensional oscillators, expressing the normal modes in curvilinear coordinates (such as internal coordinates) is mandatory when making finite distortions. (31) The method has been demonstrated to reach chemical accuracy for the predicted free energies, (30,32,33) but it is difficult to automatize and not straightforward to apply. It requires the definition of a set of internal coordinates in a periodic system. (34)
Here, we propose a different approach to get the vibrational partition function and thus the free energies. We use MD to compute the vibrational density of states (VDOS). (3,35−37) In contrast to MD with thermodynamic integration, meta-dynamics, and particle insertion, which require long simulation times to get well-converged total free energies, VDOS based on velocity auto-correlation functions are converged within rather short MD trajectories of roughly 10 ps. (38,39)
The drawback of this method is that the required VDOS for the prediction of vibrational partition functions must contain vibrational fundamental bands only, while the VDOS computed from MD trajectories, by construction, also includes overtones and combination bands. These unwanted contributions must be identified and removed if one wants to make quantitative predictions. One possible strategy is to project the MD trajectories on the vibrational modes and then filter the signal. In the past, for crystalline systems, Zhang et al. (35) and Carreras et al. (36) have made use of a hybrid scheme that combines MD simulations and harmonic phonon calculations. However, the underlying implicit assumption that the anharmonic modes only slightly differ from the harmonic ones is not valid anymore when large amplitude vibrations are involved, for example, when floppy organic molecules or non-rigid solids such as zeolite frameworks are tackled.
Actually, anharmonic effects and large amplitude vibrations have been shown to give non-negligible contributions to the adsorption free energies, (30) and the use of the correct modes is the only way to ensure that the projection allows one to identify the fundamental band without too many spurious contributions, generating uncertainties in the free energy values.
We present here a methodology which is based on harmonic normal modes of a well-chosen reference structure and allows one to predict the evolution of the anharmonic vibrational modes “on-the-fly” along the trajectories. The method has a reduced computational cost compared to recomputing the normal modes at each time step of a simulation and has no issues even when a region of the PES far from one of the minima is explored. For small alkanes in the H-CHA zeolite, we will show that the method yields Gibbs free energies of adsorption in close agreement with the experiment (within 5.6 kJ/mol), whereas the harmonic approximation yields errors up to 15.8 kJ/mol. Approaching chemical accuracy (±4 kJ/mol), our method is a good compromise between accuracy and computational cost.
Here, we study the adsorption of small alkanes (methane, ethane, and propane) on the proton form of zeolite chabazite, H-CHA. Its framework (Figure 1) consists of corner-sharing SiO4 tetrahedra. The acidity of zeolites is due to bridging Si–O(H)–Al hydroxyl groups obtained by the isoelectronic substitution of Al for Si (Si/Al ratio of 11/1) and addition of a proton for charge compensation. The small dimensions of its unit cell make H-CHA suitable as a benchmark for computational methods, (10,18,28−30,40−42) before more extended systems of similar nature and industrial relevance, (43) such as H-MFI, H-ZSM-22, or H-FAU zeolites, will be studied.

Figure 1

Figure 1. Zeolite chabazite, H-CHA with an Si/Al ratio of 11/1, and the proton attached to the O2 position: monoclinic 2 × 1 × 1 supercell (A), adsorption complexes of methane (B), ethane (C), and propane (D). Color code: yellow—Si, red—O, blue—Al, gray—C, and white—H. Panel (E): harmonic VDOSs for the adsorption complex of ethane in H-CHA, modes belonging to adsorbed ethane (top), the O–H of the Si–O(H)–Al bridge (middle), and the framework (bottom), with intensities renormalized for the number of atoms for each subsystem. Panels (A–D) adapted with permission from ref (30).

The monoclinic (2 × 1 × 1) supercell of the bare H-CHA (Figure 1A) has the composition (HAlO2)2(SiO2)22. The proton (one for each unit cell) is attached to the oxygen in position O2 which is shared by two eight-membered rings and one four-membered ring. (30) Together with the bare surface, Figure 1 also shows the adsorption complex formed with methane (1B), ethane (1C), and propane via the primary carbon (1D). The Si–O(H)–Al hydroxyl group points toward the nearest carbon of the adsorbed alkane. In the harmonic VDOS (Figure 1E), the O–H stretching bands of the “free” bridging OH group (3675 cm–1) and of the bridging OH group interacting with ethane (3473 cm–1) are well separated from each other, and also from the CH stretching vibrations of the alkane molecule that are positioned around 3000 cm–1. As expected, the region between 1800 and 2900 cm–1 shows no harmonic vibrations, while the region below 1200 cm–1 is quite crowded. Here, we can find contributions belonging both to the alkane molecules (bending, torsions...), the Si–O(H)–Al bridges (bending, out-of-plane bending, torsion...), and the full set of the Si–O–Si framework vibrations.

2. Methods

ARTICLE SECTIONS
Jump To

2.1. Quasi-Harmonic Approximation

From standard thermodynamics, the adsorption Gibbs free energy, ΔGa, of a system at temperature T can be calculated as
(1)
where ΔHa and ΔSa are the enthalpy and entropy of adsorption. For the vibrational contributions of the gas phase molecule, molec, the bare surface, surf, and the adsorbed molecule on the surface, molec@surf, Gvib is computed directly from the vibrational partition function Zvib of the system; for the rotational and translational contributions of the gas phase molecules see the Supporting Information. In the harmonic approximation, Zvib is the product of the vibrational partition functions of a set {i} of independent vibrational states for which the summation can be performed analytically (44)
(2)
kb is the Boltzmann constant, T is the temperature, ωi is the vibrational frequency of the i-th mode, and ℏ = h/2π with the Planck constant h.
Correspondingly, the vibrational entropies, Svib, and enthalpies, Hvib, are obtained from
(3)
and
(4)
respectively.
The zero-point vibrational energy (ZPE) is given by
(5)
Here, we rely on the quasi-harmonic approximation (QHA) which assumes that eqs 24 are still valid also for a set of independent anharmonic oscillators, where now ωi becomes the anharmonic fundamental frequency of the vibrational mode. Differently from Piccini et al. (30) who solved one-dimensional Schrödinger equations for each vibrational mode, here we obtain the anharmonic frequencies from the VDOS computed by MD simulations. (21,35,36,45−48)
If we consider a system consisting of a set {i} of independent vibrational modes, which form a continuous band of harmonic states, eqs 25 assume the form (49−51)
(2a)
(3a)
(4a)
with
(5a)
where D(ω) is the VDOS.

2.2. VDOS—without and with Projection

For the integration in eqs 2a5a, in the past, (21,37,45−47) use has been made of the VDOS obtained directly from the Fourier transform (FT) of the time autocorrelation function of cartesian velocities
(6)
vα(t) is the velocity of atom α at time t in mass-weighted coordinates and Ek is the average kinetic energy of the system. We will call this method in the following density of state integration (DOS-I).
However, the VDOS required for eqs 2a5a must contain vibrational fundamental bands only, while D(ω), computed according to eq 6 from MD cartesian velocities, by construction includes also overtones and combination bands. The latter arise from the mechanical anharmonicities that are naturally included in the MD trajectories.
As an example, Figure 2 (left panel) shows the VDOS between 2850 and 3000 cm–1 obtained from the cartesian velocity autocorrelation function of a trajectory of gas phase ethane at 300 K: a band can be clearly spotted at 2870 cm–1. At this frequency, no fundamental bands are expected and most probably this contribution arises from a combination band between the CH stretching and the CH3 torsion.

Figure 2

Figure 2. Left panel: VDOS of the ethane molecule in the gas phase at 313 K computed from the autocorrelation function of cartesian velocities (eq 6). Right panel: contribution of the Si–O stretching vibration [bridging Si–O(H)–Al groups] for the bare H-CHA to the total VDOS between 700 and 850 cm–1 computed by projecting the velocities on the harmonic modes without (blue line) or with applying a Lorentzian window filter (red line).

DOS-I implicitly assumes that the artificial contributions of overtones and combination bands to the free energies are negligible. While it is reasonable to assume that these contributions are small compared to the total absolute free energy of the system, it is more difficult to quantify how much they affect free energy differences such as adsorption free energies, especially if the latter are small, a few kJ/mol. If the overtones and combination bands are similar for the bare surface and the surface with the adsorbed molecule, the error will cancel out, but this is not always the case. For this reason, we have decided to follow a different path that allows one to remove these additional contributions.

2.3. VDOS—Projection Method

The total VDOS of the system D(ω) is equal to the sum of the vibrational densities of states di(ω) of each independent vibration i
(7)
where N-mode = 3N – 6 for a gas phase molecule and N-mode = 3N – 3 for the surface. N is the number of atoms of the system.
The individual VDOS di(ω) can be obtained in turn by the FT of the autocorrelation function of the velocities of the vibrational mode qi
(8)
where εki is the average contribution to the kinetic energy of the i-th mode. Notice that at equipartition εki = kbT/2 for each vibrational mode.
To compute the set of {i(t)} from the fully anharmonic “coupled” trajectories, one possibility is to project the cartesian velocities of the atoms onto the normal modes. (35,36,48) If we call L the eigenvector matrix solution of the harmonic problem
(9)
where vα(t) is the velocity of atom α at time t in mass-weighted coordinates, vαi(t) is the projection of the velocity of α on the normal mode i, and Li is the eigenvector of the normal mode i. In this case
(10)
If the potential energy curve for the i-th vibrational mode is a single-well potential, di(ω) will show a single peak function centered on the anharmonic frequency of the fundamental for the i-th mode. In this case, the anharmonic frequencies can be extracted by fitting di(ω) with, for example, a Lorentzian function. (36) The Hvib, Svib, and Gvib values can then be computed by entering these anharmonic frequencies in eqs 24. If, however, multiple quasi-degenerate minima are present on the potential energy curve as it is the case, for example, for hindered rotations and translations, the VDOS spectra will show broad bands and possibly multiple maxima. Moreover, even in the case of a single-well, the inhomogeneous broadening of the bands can play a role in determining Gvib. To take these effects into account, one possibility is to evaluate the contribution of each vibrational band by eqs 2a5a, integrating this time over the individual VDOS di(ω) instead of the total D(ω)
(2b)
(3b)
(4b)
(5b)
Notice that in eqs 25 and 2b5b, the contributions of the different modes are independent of each other. Moreover, E0vib can be computed separately. For each mode, one can therefore freely choose either to follow the “fitting strategy” to extract the anharmonic frequency for use with eq 2 or to perform the integration on di(ω). As we will discuss in more detail in Section 2.9, this choice can be made depending on the shape of the band and on the accuracy of its prediction.
However, as we have already pointed out, the required di(ω) in eqs 2b5b for Gvib, Svib, and Hvib must contain vibrational fundamental bands only. By construction, the FT of the velocity autocorrelation function from MD trajectories, even after the projection, includes contributions coming from overtones and combination bands. Because combination bands and overtones are most of the time much less intense than the fundamental ones in our classical trajectories, they do not constitute a problem if one decides to extract anharmonic frequencies by fitting di(ω), cf.eq 2. However, they will artificially modify the values of Svib and Hvib if one instead directly integrates di(ω) or D(ω).
Moreover, the harmonic projector, L, does not exactly correspond to the real anharmonic one, especially when cartesian coordinates are used to define the vibrations and large amplitude motions are involved. This can create noise in the projected spectra.
To remove these undesired contributions, in this work we have applied a Lorentzian window to the VDOS di(ω) of each mode. The width of the window is usually chosen in a way to be four times the guessed width for the band (estimated by fitting di(ω) with a Lorentzian function). This way, we are able to both clean di(ω) from the spurious contributions and to preserve the shape of the fundamental bands, that is, to still take the possible inhomogeneous broadening into account.
For example, Figure 2 (right panel) shows the 720–820 cm–1 spectral region of the bare H-CHA surface. The contribution of the SiO stretching band of the Si–O(H)–Al bridge to the total VDOS computed following the above strategy, that is, projection and filter (red line), is compared to the one coming from the projected velocities but without applying any filter (blue line). While the fundamental band of the Si–O stretching, ω0, is clearly recognizable in both cases, a set of additional peaks rise when the Lorentzian window filter is not applied. These additional contributions are due to the couplings of this Si–O stretching with other vibrational modes of the framework, in particular with other Si–O vibrations. For this mode, the difference in the integrated VDOS is around 25%.
We will call this method which is based on projecting on normal modes DOS-P.

2.4. Initial Conditions for the Trajectories

Now that we have the full set of equations to compute Hvib and Svib from the MD simulations, next we should clarify at which conditions we run the trajectories. A set of initial positions and velocities for the atom of the system must be chosen. The strategy commonly adopted (21,35,36,45−48) is to prepare the system to be at equipartition at a target temperature and to sample the system at these conditions. This we will do here as well.
However, at 300 K most of the vibrational modes of the kind of systems studied here are in our classical trajectories, much below their quantum vibrational ground state, and sampling the system at equipartition may not fully account for the anharmonic nature of the PES. Therefore, in addition, we use the following alternative protocol. Starting from the equilibrium structure a set of possible initial velocities is obtained by vertically exciting at the same time all harmonic normal modes up to their respective ZPE.
(11)
where vα(0) is the velocity of atom α at time 0 in mass-weighted cartesian coordinates, ωiH is the harmonic frequency of the i-th vibrational mode, and ϕi ∈ {0,π} is a random phase assigned to the mode. Therefore, vα(0) includes contributions from the excitation of all modes to which it contributes.
Changing the set of {ϕi}, an ensemble of replicas for the system is created: each replica will have the same initial positions, the same total energy, but different velocities of the atoms. After running the MD simulations, the atomic velocities vα(t) are projected on the vibrational modes (eq 9) to compute the vibrational mode velocities i(t). Then, for each replica, the velocity autocorrelation function, δi(ti(0) of each mode i, is computed separately. Finally, the i-th velocity autocorrelation function to be used in eq 8 is obtained as the average of the velocity autocorrelation functions computed for each replica.
Not being at equipartition, εki is no longer equal to kbT/2 in eqs 8 and 10. One possibility is to estimate εki from the vibrational temperature of the harmonic mode. Another possible choice is to select εki to ensure that ∫0di(ω) = 1 and, thus, each mode contribute with exactly 1 to the total VDOS, as it should be. This latter strategy allows one to compensate for small errors in the preparation of the system, for example, that the system is excited to the harmonic ZPE and not to the exact anharmonic ZPE.
In the Results and Discussion section, we will compare the values for Hvib and Svib computed with DOS-P according to the two protocols: equipartition at a target temperature (superscript T: TDOS-P) or preparation in the harmonic ZPE state (superscript ZPE: ZPEDOS-P).
Without the projection procedure, there is no way to correct a-posteriori for the equipartition. This is a major disadvantage of the direct use of the cartesian velocities, that is, the DOS-I method. With DOS-I, before starting the production run, a perfect equipartition of the energy between the vibrational modes must be reached. This is the only way to ensure that each mode contributes to the VDOS with exactly 1, as it should. For vibrational spectroscopy, an error of 20–30% in the equipartition is not considered to be a problem because it is still possible to assign the bands. However, for the prediction of free energies, a similar error makes the results meaningless.
For systems as liquid water, to reach equipartition is not much of a problem: the multiple interactions and strong couplings will ensure it in a short time. However, when gas phase molecules are involved, as in the surface–molecule interactions studied here, modes that are completely decoupled from the others exist which means that long equilibration trajectories (from 50 ps to 1 ns) must be used prior to the production run. Notice that the stability of the equipartition must be reached: the thermostat cannot be used in the production run otherwise its coupling with the vibrational modes will be a part of the spectrum which will change the final value for the adsorption free energies. In contrast, with our projection approach, equipartition can be enforced a-posteriori. Therefore, short trajectories on the 1–10 ps time scale can still be used to get accurate free energies.
Finally, as we will see in more detail in the next section, the projection strategy allows one to work at non-equilibrium conditions. This makes it possible to obtain results at ZPE conditions even with “cheap” trajectories obtained with standard DFT-MD using classical nuclei, as done in this paper. Without the projection, to explore the PES at the mode-specific ZPE, the use of the computationally expensive semiclassical MD would be required.

2.5. Vibrational Relaxation Time Scale

Another parameter of the MD simulations that is important to get the correct Hvib and Svib values is the right simulation length/time scale for sampling the PES. The usual perception is that the longer the simulations the higher the accuracy and that the simulation length should be limited only as a compromise between the accuracy and computational cost. Whereas this is still true for our trajectories prepared at equipartition (we have used two simulations of 20 ps each), the ZPE initial condition case is special. The replicas of the system are prepared in a way in which each mode is at its ZPE. Unfortunately, for MD with classical nuclei such as the ones used here, the ZPE is a non-equilibrium condition. The system will spontaneously relax and redistribute the energy between the modes to reach equipartition.
Therefore, the length of the trajectories must be short enough that the system does not depart too much from the ZPE condition. At the same time, the simulations must still be long enough that the sampling of each vibrational mode is meaningful, that is, at least 20 oscillations per mode.
The relaxation can happen on different time scales depending on the couplings between the modes. As an example, Figure 3 shows the frequency shift (compared to the harmonic wavenumber) of the vibrational modes of the bare H-CHA surface at different times after the initial vertical excitation. The shift can be considered a qualitative measure of the level of excitation of each vibrational mode. The figure shows that the high-frequency modes can relax within less than 250 fs, while some lower frequency modes are still fully excited after 2 ps.

Figure 3

Figure 3. Frequency shift (cm–1) compared to the harmonic wavenumber (see figure inset for color coding) for the vibrational modes of the bare H-CHA surface, 0.25, 0.50, 1.0, and 2.0 ps after the initial vertical excitation. The level of excitation is qualitatively measured as a wavenumber shift from the harmonic wavenumbers.

We conclude that, for ZPE initiation, the use of different time scales to compute the velocity autocorrelation function for different modes is the best strategy: the lower the harmonic frequency ωi of the i-th vibrational mode is, the longer are the trajectories used to compute the di(ω). For example, for a mode vibrating at 40 cm–1, we would use a 3 ps trajectory, while for one vibrating at 3000 cm–1 the trajectory is reduced to 100 fs (see Supporting Information for more details, Table S1).
Notice that in any case all the modes are initially excited in all the trajectories. What is changing mode-by-mode is only the portion of the same trajectory used to compute the velocity autocorrelation function.

2.6. Rotations

If we consider again eqs 2b5b, we can realize that an additional problem arises for gas phase molecules where rigid body rotations are strongly coupled to vibrations in the MD trajectories. This can lead to artificial broadening and sometimes even to the splitting of vibrational bands. Simple projection on normal modes does not remove these couplings.
One solution is to clean the trajectory before the projection. If we call ξ(t), the (3N) vector that collects the Cartesian coordinates of the N atoms of the molecule in the center of the mass reference, a rotation-free trajectory can be obtained by
(12)
where R(t) is the rotational matrix that guarantees the best overlap between ξ(t) and the initial geometry ξ(0). R(t) can be obtained by a quaternion fit that minimizes the sum of the squared distances between the mass weight coordinates of corresponding atoms. (52) Such a rotation satisfies the Eckart conditions for small displacements (53) and can still be used for the large ones. (38) Once ξvib(t) is obtained, the velocities can be computed by numerical differentiation. We use a five-point central difference formula (see Section S4 of Supporting Information for details) to guarantee a negligible numerical error on the kinetic energy.

2.7. From a Constant L to L Evolving in Time

Carreras et al. (36) and Zhang et al. (35) assume that the eigenvector matrix L (eq 9) which projects the cartesian coordinates on the vibrational space is constant in time. This can be considered a valid approximation for crystalline solids, such as the ones they were studying, for example, crystalline silicon, (36) MgSiO3 perovskite, (35) and CaSiO3 perovskite. (48) In contrast, we are dealing here with systems that show large amplitude motions and can explore multiple minima on the PES. Therefore, L is not constant in time and a way to compute its changes is needed.
Predicting the time evolution of L by diagonalizing the hessian matrix at each time step is not a viable route. Because for most of the simulation time, the system is out of a minimum of the PES, the solution of the harmonic problem would generate a huge number of negative eigenvalues that one must deal with. Moreover, the computational cost would be enormous.
We propose the following strategy. We assume that L(t), while not constant in time, at instant t can still be properly described by a linear combination of the eigenvector matrices Lref (solution of the harmonic vibrational problem) for a set of appropriate reference structures
(13)
where Pj(t) is the probability that at time t the system is vibrating around the j-th reference structure. The relevant structures can be set up by an educated guess (knowing the structure of the systems) and/or by a preliminary MD run to explore the phase space.
To evaluate Pj(t), we assumed a Gaussian distribution around the reference structures (38,54)
(14)
where mj is a properly defined metric measuring the “distance” between the geometry of the system at instant t and the reference structure j and σc is the width of the Gaussian. The metric mj must be chosen in a way that it is easy and cheap to be evaluated on-the-fly but that it can still capture the fundamental differences between the reference structures. It has been shown (38,54) that the following definition has the required characteristics
(15)
where {γk} can be a set of coordinates of the same type, for example, all interatomic distances, bond angles, torsional angles, or coordination numbers. Notice that in the case of the coordination number, the following continuous definition has been adopted (55)
(16)
where r(t) is the distance between the two atoms at time t and R0 is a cutoff distance for the interactions.
If a set of coordinates of different types is required to univocally define the reference structures (e.g., a combination of coordination numbers and torsional angles), the metric becomes a vector mj = {mjx} and Pj can be generalized to
(17)
with
(18)
where {x} is a subset of coordinates of the same type.
Finally, to describe large amplitude motions such as the CH3 torsions, avoiding the introduction of a huge number of reference structures, a multifragment strategy can be applied. (38) The system is divided in a set of relevant fragments and for each fragment the Lα∈frag matrix is rotated in a way to guarantee the consistency of the reference system at each time step.
(19)
Also in this case, the rotational matrix Rα∈frag can be obtained by minimizing the sum of the squared distances between the mass weight coordinates of the corresponding atoms.

2.8. Divide-and-Conquer Approach for the L Matrix

For bare H-CHA, many vibrational modes of the framework fall in the same frequency region as the torsions, hindered rotations, and translation of the adsorbed alkane molecules (Figure 1E). In the harmonic solution, the degeneracy of the vibrational levels can produce an accidental mixing of the vibrational modes of the adsorbed molecule with the ones of the framework because any linear combination of degenerate modes is still an eigenvector. While this has no impact on the calculation of the harmonic entropy, enthalpy, and free energy, it can bias the quality of the results for the projection procedure.
The anharmonicity of the PES can lift the degeneracy of the two modes that were degenerate on the harmonic potential. In this case, the projection will generate a spectrum of peaks instead of a single well-defined one (see Section S3 of Supporting Information for an example). To bypass the problem, the separation of the modes must be enforced also in the harmonic solution, imposing the localization of the vibrational modes that we use for the projection. To do so, we still use the full Hessian but artificially modify the masses of the parts of the system to shift the vibrational frequency (like in a deuteration). Proceeding part-by-part, we can reshape L in a block matrix in which each block is localized on a particular set of atoms of the system. In all cases, the alkane molecules, both gas-phase and adsorbed, have been modeled with a single block of atoms. The zeolite has been modeled with two separated blocks: the O–H group of the Si–O(H)–Al bridge that interacts with the adsorbed molecule (LOH) and the rest of the silicate framework (Lframe). The H-CHA zeolite presents a huge number of quasi-degenerate harmonic modes below 1000 cm–1, as shown in Figure 1. The correct description of Lframe(t) requires the introduction of multiple blocks and fragments which would increase the computational cost. Because the major contributions to the adsorption free energies are expected to come mainly from the Si–O(H)–Al site onto which the alkane is adsorbed, we decide to choose a simpler path and neglect any direct contribution to the VDOS from the vibrational modes localized on the zeolite but the O–H of the Si–O(H)–Al ones, imposing Lframe = 0.
Notice that the excitation with the full Hessian ensures that even if two coupled vibrations will accidentally decouple by the definition of the blocks in the projection they will still fall on the same (correct) frequency. The filter will then allow to compensate for a possible error in the total intensity, and the correct free energy will be recovered (see Section S7 of Supporting Information, for further discussion on the active space).
To summarize, the vibrational space used for the projection is composed of 6 modes for the bare H-CHA zeolite (the 6 vibrations localized at one of the two O–H), 21 modes for the methane@H-CHA (6 localized on the O–H that interacts with the adsorbed molecule, 15 localized on CH4), 30 for ethane@H-CHA, and 39 for propane@H-CHA. Only this subset of modes is used to compute D(ω) and the vibrational partition function according to eqs 24 or 2b4b, while we neglect the contributions from the modes localized at the rest of the silicate framework.
Note that we reduce the vibrational space in the projection step and when summing up the partition functions but not when we generate the initial conditions for the trajectories (see Section 2.4). The vertical excitations are applied to the original, complete, vibrational space of the system (e.g., 219 modes for the bare H-CHA zeolite). Moreover, we project (on the uncoupled modes) trajectories in which all the atoms are free to vibrate. This is particularly important for including the coupling of low-frequency “adsorbate” modes (OH + alkane) with modes of the silicate framework. (56) For example, it has been shown that soft vibrations can couple in a way that favors the diffusion of the molecule in the zeolite framework. (17,57)

2.9. Summary and Additional Choices

In Section 2.2, we have discussed two possible ways to compute Hvib and Svib from the projected di(ω) as follows: anharmonic frequencies were extracted by fitting di(ω) and using eqs 24, or direct integration of di(ω) (eq 5b). The integration of di(ω) accounts for the inhomogeneous broadening of the bands. However, for the high-frequency modes, the 0.4 fs time step combined with the short length of the trajectories imposes a width of the FT of ∼50 cm–1, artificially inducing a Gaussian shape to di(ω). Therefore, the integration of di(ω) is not meaningful for these cases and it is better to extract the anharmonic fundamental frequencies needed for eqs 25 by fitting. In contrast, for the lower frequency modes, the longer trajectories allow one to reasonably predict the shape of the bands and it is worth to use the full integral of di(ω) in eqs 2b4b to get Hvib and Svib. In this way, the complex shape of the PES, which is due, among other features, to the quasi-degenerate minima of the hindered rotations and translations, can be taken into account.
Consequently, we prefer a “hybrid” strategy with our DOS-P method: if the original harmonic frequency ωiH of a mode i is below 450 cm–1, eqs 2b4b are used, whereas if ωiH is above 450 cm–1, eqs 24 are applied with frequencies extracted by the Lorentzian fitting.
Finally, for both the frequency ranges (above or below 450 cm–1), we prefer to calculate the ZPE contributions employing eq 5 with the frequencies extracted by the Lorentzian fitting:
(20)
Here, Givib is the contribution to the vibrational free energy of the i-th mode and ωianh is the anharmonic frequency obtained by fitting.
The choice of different time scales depending on the frequency range implies that the requested number of replicas is not always the same. The ergodic principle allows one to use a reduced number of replicas for the long trajectories. We choose the number of replicas in such a way that the possible error in the predicted energy due to the missed sampling is less than 1.0 kJ/mol, see the Supporting Information for more details. As expected, the enthalpy and the ZPE are quite fast to converge, while the entropy requires a more extensive sampling.
For trajectories prepared at the target temperature, the statistical uncertainty is slightly higher, due to the different explorations of the “quasi-free” state (the molecule is in the pore but not adsorbed on the acid site) along the trajectories. However, with two trajectories, the statistical uncertainty is below 1 kJ/mol, and we consider this a good compromise between computational cost and accuracy.

2.10. DOS-Projection Method

Our computational protocol for extracting anharmonic vibrational frequencies from MD simulations by projection (DOS-P) for use in partition functions consists of the following steps:
1.

Structure optimization and calculation of the harmonic normal modes for the minimum energy structures

2.

Generation of a set of replicas of the system by vertical excitation of the full set of harmonic modes for the system

2a

At a target temperature (TDOS-P) or

2b

Each up to its own ZPE (ZPEDOS-P)

3.

Running the trajectories

4.

Decision about use of the full vibrational space for the projection or of a subspace only (see Section 2.8). In the case of the sub-space, L was redefined with reference structures for the subspace. Only the modes belonging to the subspace will be used to compute the vibrational partition function.

5.

Prediction of the evolution of the eigenvector matrix along the trajectories by a weighted sum of the L of relevant reference structures (eq 19)

6.

Projection of the velocities on the harmonic modes (eq 9)

7.

Calculation of the FT of the velocity autocorrelation function of the vibrational modes (eq 8)

7a

On the full trajectories (equipartition route, TDOS-P)

7b

On the time scale relevant to each vibration (ZPEDOS-P)

8.

For all modes, use of a Lorentzian to fit di(ω) to extract the maximum position and width of each band.

9.

With these anharmonic frequencies for all the vibrational modes the ZPE was computed according to eq 5.

10.

For modes with an original harmonic frequency

10.1

ωiH > 450 cm–1: the anharmonic frequencies obtained from the fit in step 7 were used to calculate Sivib and Hivib according to eqs 3 and 4.

10.2

ωiH < 450 cm–1: filtering of the VDOS to remove spurious contributions and compute Sivib and Hivib by integrating di(ω) according to eqs 3a and 4a.

11.

Calculation of the thermodynamic functions Svib, Hvib, and Gvib by summing up the contributions according to point 10.1 (ωiH > 450 cm–1) and point 10.2 (ωiH < 450 cm–1).

3. Computational Details

ARTICLE SECTIONS
Jump To

3.1. Systems, Electronic Energies, and MD Simulations

DFT calculations with the VASP code (58) have been performed at periodic boundary conditions at the Γ point only, with a 600 eV plane-wave kinetic energy cutoff. Consistently with previous works, (28,30) the PBE (59) functional, augmented with the Grimme D2 dispersion term, (60) has been chosen.
The H-CHA supercell (Si/Al ratio 11/1, Figure 1) has the dimensions a = 18.90 Å, b = 9.44 Å, c = 9.29 Å, α = 94.0051, β = 94.8903, and γ = 95.3793°. (30) The optimized structures for the adsorption complexes correspond to half coverage (θ = 0.5). This ensures that the lateral interactions between the adsorbed molecules are negligible. For gas phase molecules, a 15 Å × 15 Å × 15 Å cell has been adopted.
For propane, two adsorption configurations are possible via the primary (CH3 group) or the secondary carbon atom (CH2 group). Bučko et al. (10) have shown that at 300 K, the probability for adsorption via the primary carbon is around 7 times higher than via the secondary carbon atom. Piccini et al. (30) found that adsorption via the primary carbon is entropically favored by TΔS = 12 kJ/mol at 313 K. Therefore, in the following, we will consider only the adsorption of propane via primary carbon (shown in Figure 1D), while adsorption via the secondary carbon is further discussed in Section S6 of Supporting Information.
To compute the harmonic normal modes, starting from the normal mode-optimized structures published in ref (30), Hessian matrices in Cartesian coordinates have been obtained by numerical differentiation. An SCF convergence criterion of 10–8 eV/cell has been applied.
The electronic energies ΔEel used in this paper are the ones already published in ref (30), obtained with the hybrid MP2:PDE + D2 approach. They are reported in Supporting Information, Table S3, together with the rotational and translational entropy terms for the gas phase (Table S4).
We use Born–Oppenheimer MD, that is, at each time step the electronic wavefunctions (Kohn–Sham orbitals) have been converged, imposing a threshold for the energy difference between two SCF cycles of 10–6 eV/cell. The classical Newton equations of motion for the nuclei have been integrated through the Verlet algorithm with a time step of 0.4 fs.
All the trajectories have been run in the NVE microcanonical ensemble. The initial positions of the atoms for the dynamics are the ones of the optimized structures discussed above. The initial velocities are obtained by vertical excitation using the protocol discussed in the Methods Section 2.4. Because with standard thermostats reaching equipartition of the energy between the vibrational modes requires long simulation times, we have decided to proceed with the vertical excitation of the modes also in the case of the excitation at a target temperature. In this case, the energy deposited in each mode is kbT, where T is equal to 303 K for methane and 313 K for ethane and propane. Two trajectories of 20 ps each have been run for each system. Notice that in this way also the “quasi-free” state in which the molecule is not directly interacting with the surface but still in the pore is sampled.
The additional reference structures needed for the projection procedure (see Section 2) have been obtained following the normal-mode optimization protocol described in refs (28) and (29).
For all the simulations discussed (prepared at harmonic ZPE or at the target temperature), the correct weight of each vibrational mode to the total VDOS has been enforced a posteriori by imposing εki in a way to ensure ∫0di(ω) = 1.
In the case of gas-phase molecules, the trajectory predicted by the MD contains both rotational, translational, and vibrational contributions. To be able to extract a rotation-free velocity trajectory, we have followed the scheme proposed in Section 2.6, that is, obtaining a rotation-free position trajectory and by the finite difference to recompute also a rotation-free velocity trajectory. The effectiveness of the method has been demonstrated by verifying the complete disappearance of the rotational bands from the VDOS. In the case of the H-CHA and the molecule adsorbed on the H-CHA surface, there is no need to remove the rotational contributions. However, for consistency, we have computed also in this case the instantaneous atomic velocities by the finite difference.

3.2. Reference Structures

For all the systems investigated here, we use as reference structures the normal mode-optimized ones of ref (30) plus any equivalent structure with the permutation of the H atoms that may appear in the dynamic runs due to the rotation of the CH3 group.
To be able to describe the large amplitude torsion of the OH bond around the Si–O bond of the Si–O(H)–Al site, an additional set of reference structures has been generated for the degenerate minimum of the zeolite differing from the first ones only for the H·Si·Al·O improper dihedral angle (11.5° vs −17.1°), that is, the angle between the (Si)–O–H bond and the Si–O(H)–Al plane. Due to the low energy barrier, during the dynamics sometimes the system jumps between the two minima along the path of the H·Si·Al·O out-of-plane bending. (61)
In Table S2 of Supporting Information, we report the number of reference structures and the metric used to describe them.

4. Results and Discussion

ARTICLE SECTIONS
Jump To

4.1. Vertical Excitation at the Target Temperature Versus Harmonic ZPE

As already discussed in Section 2.4, there are two possible choices for the initial conditions of the trajectories: all the modes at the target temperature (TDOS-P) or all the modes excited at their own harmonic ZPE (ZPEDOS-P).
Table 1 compares the results for the two possible initial conditions. As expected, the absolute values for the vibrational entropy, enthalpy, and free energies differ substantially between the two kinds of excitation. When the system is prepared at 300 K (TDOS-P), most of the modes are well below their ZPE, and therefore only a small portion of the phase space is accessible. In particular, the PES near the ZPE level that can show for some modes a strong anharmonic shape is not sampled, shifting the computed frequencies of the fundamental bands. As a consequence, E0vib can be up to 10 kJ/mol higher than in the case of the system initially at the harmonic ZPE level (ZPEDOS-P).
Table 1. Absolute ZPEs, E0vib, Vibrational Enthalpies, Hvib, Vibrational Entropy Terms, −TSvib, and Vibrational Gibbs Free Energies, Gvib, for Alkanes Adsorbed in H-CHA at Experimental Conditions (303 K for Methane and 313 K for Ethane and Propane, 0.1 MPa, θ = 0.5), Calculated for Initially Exciting the System at the Target Temperature (TDOS-P, T = 303 K for Methane and T = 313 K for Ethane and Propane) or for Exciting Each Mode at Its Harmonic ZPE (ZPEDOS-P)a
adsorbate ZPEDOS-PTDOS-PΔ
methaneE0vib150.1155.8–5.7
 Hvib15.615.80.2
 TSvib–37.0–40.23.2
 Gvib128.7131.4–2.7
ethaneE0vib222.6229.2–6.6
 Hvib19.519.10.4
 TSvib–49.8–51.71.9
 Gvib192.3196.6–4.3
propaneE0vib293.3302.8–9.5
 Hvib23.123.00.1
 TSvib–58.3–59.31.0
 Gvib258.1266.5–8.4
a

In the last column, we report the difference Δ between the two. All in kJ/mol.

In contrast, the entropy contribution to the total free energy (−TSvib) is lower for TDOS-P, with a maximum difference of 3.2 kJ/mol for methane. The Hvib values predicted by the two approaches agree within 0.5 kJ/mol. This different sensitivities with respect to the initial conditions of E0vib, Hvib, and −TSvib are related to the fact that the lower frequency modes that dominate Svib and Hvib at 300 K are also in our classical dynamics not so far from their real quantum vibrational ground state levels, whereas the high-frequency ones that carry the major contributions to E0vib are well below this. Svib is more sensitive than Hvib because it is influenced also by the tail of the vibrational partition function.
When free energies for adsorption are considered (Table 2) instead of absolute values, a compensation of errors occurs. For the two kinds of excitation, the differences between the calculated ZPEs of adsorption are between 1.0 and 2.4 kJ/mol. The predicted thermodynamic functions of adsorption have more positive values for ZPE excitation (ZPEDOS-P) than for T excitation (TDOS-P). The ZPE, enthalpy, and entropy terms add up to Gibbs free energies that, therefore, are more positive for ZPEDOS-P than for TDOS-P.
Table 2. Electronic Energy from ref (30), ΔEela, ZPE, ΔE0-viba, Enthalpy, ΔHa, Entropy Term, −TΔSa, and Gibbs Free Energy ΔGa for the Adsorption of Methane (303 K), Ethane (313 K), and Propane (313 K) in H-CHA at Experimental Conditions (0.1 MPa, θ = 0.5), All in kJ/mol, Obtained from MD Runs after Initial Excitation at a Target Temperature, TDOS-P (T = 303, 313, and 313 K for Methane, Ethane, and Propane, Respectively), and for Excitation of Each Mode at Its Harmonic ZPE, ZPEDOS-Pa
  ZPEDOS-PTDOS-PΔ
methane@H-CHAΔEela–25.3–25.3 
 ΔE0-viba2.81.81.0
 ΔHa–20.4–21.81.4
 TΔSa25.323.81.5
 ΔGa4.92.02.9
ethane@H-CHAΔEela–36.2–36.2 
 ΔE0-viba2.50.12.4
 ΔHa–31.1–34.33.2
 TΔSa29.927.22.7
 ΔGa–1.1–7.16.0
propane@H-CHAΔEela–46.7–46.7 
 ΔE0-viba3.51.61.9
 ΔHa–40.5–43.12.6
 TΔSa34.933.51.4
 ΔGa–5.8–9.73.9
a

In the last column, we report the difference between the two.

This behavior can be explained by the portion of the phase space accessible with the two kinds of excitation. Having less energy in the system, T excitation explores less energetic parts of the PES, producing lower ΔHa values. At the same time, the longer trajectories (20 ps) allow the exploration also of the “quasi-free” gas phase state in which the molecule is in the pore but not close to the surface. The contributions of this state, only in a small part explored by ZPE excitation, diminish the entropy loss on adsorption.

4.2. Comparison with Experiments

Figure 4 shows the deviations of the results of our ZPEDOS-P, TDOS-P, and DOS-I simulations from experimental values. (30) As reference, we also include the results of Piccini et al. obtained from harmonic (HARM) and anharmonic (curvilinear distortions, ANH) vibrational partition functions. (30)

Figure 4

Figure 4. Absolute deviation of the computed thermodynamic functions from the experimental values for methane (blue), ethane (red), and propane (green) in H-CHA at experimental conditions (303 K for methane and 313 K for ethane and propane, 0.1 MPa, θ = 0.5). The gray shadows represent the statistical uncertainties (see the Supporting Information for more details). Upper panel: adsorption free energy, ΔGa. Center: adsorption enthalpy, ΔHa. Lower panel: adsorption entropy term, −TΔSa. On the left: the harmonic (HARM) and anharmonic (ANH) predictions of Piccini et al. (30) Center panel: DOS-P derived ΔE0-viba. Right panels: harmonic ΔE0-viba. The average deviation for each method between the three molecules is reported above the columns for the three molecules.

We look first on the columns labeled “ZPE from DOS”. The ΔGa predictions of ZPEDOS-P and TDOS-P both agree within chemical accuracy limits (±4 kJ/mol, red dotted line in Figure 4) with the experimental values and apparently indicate a similarly good, for TDOS-P even slightly better, performance than the anharmonic reference calculations. However, the larger deviations for ΔHa and TΔSa show that this is in part due to an error compensation. For TΔSa, the mean absolute deviations (MAD) of 6.5 and 4.7 kJ/mol for ZPE and T excitations, respectively, are a significant improvement compared to the harmonic approach, but outside the chemical accuracy range and much larger than that of the anharmonic calculations.
Also for ΔHa, DOS-P with ZPE excitation yields an improvement with respect to the harmonic approach but T excitation does not. The origin is not in the thermal part but in the ZPE contribution to ΔHa (Figure 5).

Figure 5

Figure 5. ZPE of adsorption, ΔE0-viba (kJ/mol), for methane, ethane, and propane in H-CHA, calculated from harmonic frequencies (blue dotted line), from anharmonic frequencies derived from curvilinear distortions (30) (red dotted line), and obtained with DOS-P methods by initially exciting the system at 300 K (TDOS-P, red full line) or at each mode’s ZPE (ZPEDOS-P, green full line).

Figure 5 shows that the anharmonic ΔE0-viba reference values (30) decrease from methane to ethane and further to propane. The DOS-P and harmonic results reproduce this trend for methane to ethane but not for ethane to propane for which all three show an increase. The DOS-P results with ZPE excitation agree with the anharmonic results for ethane but show large deviations for methane (−2.7 kJ/mol) and propane (1.6 kJ/mol), whereas the DOS-P with T excitation closely agrees with the anharmonic result for propane but show large deviations for methane (−3.7 kJ/mol) and ethane (−2.7 kJ/mol).
The ZPE is the hardest quantity to predict with the DOS-P method because it is related to the PES around the energy minimum, while MD also explores other accessible states at that energy. For systems prepared by ZPE excitation, the short trajectories partially but not completely avoid the problem (especially when the free energy of adsorption is positive) because the system has less time to sample the configuration space away from the energy minimum structure.
Because on average, the harmonic results show a smaller absolute deviation (1.0 kJ/mol) from the anharmonic predictions (30) than the ZPEDOS-P and TDOS-P methods (1.5 and 2.1 kJ/mol, respectively), we considered also substituting ΔE0-viba in the DOS-P adsorption enthalpies with the harmonic ones. The results are shown in the right panel of Figure 4 labeled “harmonic ZPE”.
Indeed, this substitution increases the accuracy of the predicted ΔHa values: the MAD decrease to 2.1 and 2.8 kJ/mol for ZPE and T excitations, respectively. Because this also diminishes the favorable cancellation of error between ΔHa and TΔSa, the MAD for ΔGa slightly increases to 4.5 and 2.6 kJ/mol for ZPEDOS-P and TDOS-P, respectively.
The small MAD for the entropy term (4.7 kJ/mol)—close to chemical accuracy—and the resulting small MAD for the Gibbs free energy (2.6 kJ/mol)—within chemical accuracy limits—suggests the TDOS-P method with harmonic ZPE as the method of choice. However, ZPEDOS-P has other advantages, see Section 4.4.
For adsorption of methane, ethane, and propane in H-CHA at experimental conditions, Figure 6 shows the final values obtained in this work by DOS projection using both ZPEDOS-P and the TDOS-P excitation and harmonic ΔE0-viba. The results of gravimetric/calorimetric sorption experiments (30) are shown as black dots with vertical bars, indicating the chemical accuracy range (±4 kJ/mol). For comparison, the predictions from harmonic and anharmonic frequencies (curvilinear distortions) of Piccini et al. (30) are also shown (see Sections S11 and table S10 of Supporting Information for the data).

Figure 6

Figure 6. Adsorption Gibbs free energy, ΔGa, enthalpy, ΔHa, and entropy term, −TΔSa, for alkanes in H-CHA at experimental conditions (303 K for methane and 313 K for ethane and propane, 0.1 MPa, θ = 0.5) computed with projection of the VDOS using ZPE excitation (ZPEDOS-P, green full line) or T excitation (TDOS-P, red full line), compared to the harmonic (HARM, blue dotted line), and the anharmonic (ANH, red dotted line) predictions. (30) Experimental values of ref (30) (black dots) are shown with the chemical accuracy range (±4 kJ/mol).

Both DOS-P results reproduce the experimental trend from methane to ethane and propane for TΔSa, ΔHa, and ΔGa. The adsorption enthalpy decreases in this series, whereas the entropy loss (−TΔSa) increases. As a result, for the experimental conditions, ΔGa changes from positive for methane, to slightly negative for ethane, and to strongly negative for propane. Hence, ethane and propane are correctly predicted to form stable adsorbates at 313 K, contrary to the prediction based on the harmonic approximation.
The ΔHa results are within chemical accuracy limits of the experimental values but systematically too negative, whereas the −TΔSa term is systematically predicted too positive, on average 6.5 and 4.7 kJ/mol for ZPE and T excitation, respectively.
Comparison with Piccini et al. (30) shows that the DOS-P results are a major improvement with respect to the use of harmonic frequencies (HARMs) and are getting close to the anharmonic (ANH) results. For ΔGa, the remaining deviations from the ANH results are smaller for the T excitation (3.0 kJ/mol on average) than for the ZPE excitation (5.6 kJ/mol on average).
Compared to the ANH predictions, (30) the DOS-P method underestimates the effect of the anharmonicity on the entropy term, on average by 4.6 and 6.4 kJ/mol for the T and ZPE excitation, respectively. Part of the missing contributions come from the Si–O–Si soft modes of the zeolite framework, which we are not including in the calculation of the thermodynamic functions (see Section 2.8). The better performance of TDOS-P can be due to the larger contributions of the “quasi-free” state of the alkane molecule to the partition function that compensates the missing ones from the zeolite framework.

4.3. Effect of the Projection on the Computed Thermodynamic Functions

A part of the motivation for this study was to find a method that needs less human interference than the calculation of anharmonic vibrational frequencies by the curvilinear distortion of normal modes as suggested by Piccini and Sauer. (28,29) In the previous section of this paper, we have shown that the integration of VDOS is a viable alternative. However, the projection procedure still requires following a complex computational protocol. In Section 2.2, we have discussed possible issues of directly integrating the DOS computed from the MD (DOS-I), without projection. This raises the question if DOS-I is a significant improvement compared to the harmonic approximation and if the additional gain in accuracy with DOS-P is worth the human effort.
First, contrary to the DOS-P, with DOS-I the only available choice is to prepare the system at a target temperature, unless one decides to use the more expensive semiclassical dynamics: eqs 2a5a cannot be used in non-equilibrium conditions. Second, without projection, it is not possible to extract the position of the fundamentals, that is, E0vib can only be computed via eq 5a. Therefore, also in this case, we calculate ΔE0-viba from harmonic frequencies as we do for our DOS-P projection methods in Section 4.2. Third, compared to DOS-P, additional MD trajectories are needed to reduce statistical fluctuations.
We have calculated thermodynamic functions by integrating directly the VDOS computed from the cartesian velocities without any projection or filter scheme (DOS-I). We have used exactly the same two trajectories employed for TDOS-P plus a third one which was needed to reduce the statistical fluctuations for −TΔSa below 10 kJ/mol (see Figure S7 of Supporting Information).
The results are included in Figure 4, whereas the MADs from the experiment are not so different between TDOS-P and TDOS-I, projection on normal modes (TDOS-P) strongly reduces the statistical uncertainty (gray shadows in the figure) compared to the direct VDOS integration (TDOS-I) of the same trajectories. For example, for TΔSa (Figure 4, gray shaded), the statistical uncertainty goes from 10 kJ/mol with TDOS-I to 1 kJ/mol with TDOS-P and to 0.5 kJ/mol with ZPEDOS-P (see Supporting Information, Sections S8 and S9, for more details). Similarly, the statistical uncertainty of ΔHa (when the harmonic ZPE is used) decreases from 5 kJ/mol without projection to less than 0.5 kJ/mol with projection.
Moreover, even if the MAD is not so different between the two schemes, contrary to TDOS-P, TDOS-I is not reproducing the experimental trend with the alkane chain length (see Supporting Information, Section S10, for the full set of data). For example, −TΔSa is predicted to be lower for propane (23.9 kJ/mol) than for ethane (32.3 kJ/mol), whereas the experimental entropy contribution is 27.5 kJ/mol for propane and 23.8 kJ/mol for ethane. (30) One of the consequences is that TDOS-I predicts a positive ΔGa value for ethane (2.5 kJ/mol), while experimentally it is slightly negative (−3.7 kJ/mol). Instead, both TDOS-P and ZPEDOS-P predict the correct sign for ΔGa for all cases considered.
It is difficult to separate the statistical error from the error connected with the contributions of the combination bands and overtones in DOS-I. With long trajectories and long pre-equilibrations, the statistical error can be reduced but the computational cost will be much higher than the one required for the DOS-P method. Keep also in mind that DOS-I includes the DOS of all vibrational modes, whereas DOS-P uses a subspace only. Calculating small differences between the numerous framework, vibrations in DOS-I with limited accuracy may produce larger errors than ignoring the small differences as done in DOS-P.

4.4. Comments on Sampling Statistics

The projection and fitting protocol (DOS-P) proposed in this paper allows to reduce the statistical uncertainty of the computed thermodynamic functions compared to the direct integration of the density of states (TDOS-I), for the same length and number of trajectories, that is, for the same computational cost.
While TDOS-P seems to perform slightly better than ZPEDOS-P in terms of accuracy (see Section 4.2), the latter, with its short trajectories, allows a better control of the sampled phase space and seems to further reduce the statistical uncertainty, for example, below 0.2 kJ/mol in the case of ΔHa. The ZPEDOS-P is particularly suitable for large systems, because it is based on the use of a set of really short simulations (no longer than 3 ps) which can easily run in parallel, reducing the required wall time.
The total simulation times are only 40 and 120 ps for TDOS-P and ZPEDOS-P, respectively, for each system, distributed on multiple short trajectories of maximum 20 ps for the former and 3 ps for the latter. This is much below the usual time scale of 50 ps to 1 ns, (8−14) needed for thermodynamic function calculations (ΔHa,–TΔSa, ΔGa), especially for the entropic contributions, with our statistical precision.
Our data obtained are in line with the results of other MD studies that use much longer time scales for the simulations. For example, Rocca et al. (27) report ΔHa = −25.5 ± 0.75 kJ/mol for CH4 in the H-CHA at 300 K (PBE + D2, 400 eV plane wave cutoff) using simulations of 100 ps, whereas with ZPEDOS-P, using a set of trajectories of maximum 3 ps, we predict ΔHa = −26.9 ± 0.2 kJ/mol. To get a directly comparable ΔHa value, we replaced our hybrid MP2/PBE + D2 electronic energy (−25.3 kJ/mol, Table 2) (30) with the corresponding PBE + D2 value (−34.7 kJ/mol). (30) Our ZPEDOS-P results agree within 1.4 ± 0.95 kJ/mol with that of Rocca et al. (27)

5. Conclusions

ARTICLE SECTIONS
Jump To

The MD simulation times for free energies of adsorption can be significantly reduced without scarifying accuracy when the VDOS of short trajectories are used to calculate vibrational partition functions. The additional advantage compared to methods, such as umbrella sampling, thermodynamic integration, or meta-dynamics, is that enthalpy and entropy terms can be calculated separately with the same accuracy and that absolute values, not only differences, of thermodynamic functions can be obtained.
Direct integration of the VDOS (DOS-I) is a significant improvement compared to the use of the harmonic approximation but suffers from large statistical uncertainties and artificial contributions from combination bands and overtones.
The statistical error can be reduced, and the accuracy improved when the fundamentals are identified by the projection of the VDOS on normal modes (DOS-P). The anharmonic nature of the vibrations along the trajectories is taken into account with well-chosen reference structures at which harmonic normal modes are calculated.
The projection is done in Cartesian coordinates and does not require to work in a space of curvilinear coordinates with all its complexity, as the method proposed by Piccini et al., (30) which solves one-dimensional Schrödinger equations for each mode separately using curvilinear distortions.
The DOS-P projection method allows one to work both at equipartition, that is, at a chosen temperature, or at mode-specific zero-ZPEs without the need of expensive semiclassical dynamics. A set of short trajectories (only two in the case of TDOS-P) is required to reach convergence within 1 kJ/mol: 1–20 ps are sufficient not only for the enthalpy but also for the entropy and free energy instead of the 50 ps to 1 ns usually required by MD.
Our results for the adsorption of methane, ethane, and propane in H-CHA show that the projection can be limited to the subspace of modes localized at the OH group of the adsorption site and the alkane molecule interacting with it. Only these modes are included in the summation for the partition functions. With this approximation, our TDOS-P method yields Gibbs free energies of adsorption in close agreement with the experiment (MAD 2.6 kJ/mol), whereas the harmonic approximation including all vibrations yields a MAD of 12.2 kJ/mol. TDOS-P is almost as accurate as the use of anharmonic frequencies obtained with the curvilinear distortion of (all) individual modes (MAD 2.4 kJ/mol). Reaching or approaching chemical accuracy also for adsorption enthalpies and the entropies, our TDOS-P method is a good compromise between accuracy and computational cost.
The method which we have applied here to the adsorption of short alkanes in zeolite can be easily extended to other low coverage cases by adapting the number of fragments and the active space. Examples are the adsorption of other molecules such as longer alkanes, carbon dioxide, water, and alcohols not only in zeolites and metal–organic frameworks (MOFs) but also on flat (oxide) surfaces.
However, for highly loaded porous materials with weak host–guest interactions and liquid-like behavior of the guest molecules or with strong guest–guest interactions, a (too) large number of reference structures would be needed, and the definition of the fragments would not be straightforward.

Supporting Information

ARTICLE SECTIONS
Jump To

The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jctc.1c00519.

  • Additional equations for the calculation of the free energies, decoupling of the vibrations due to the anharmonicity, additional computational details, data on adsorption of propane via secondary carbon, discussion on the vibrational active space, convergence check, and summary of the computed thermodynamic values (PDF)

Terms & Conditions

Most electronic Supporting Information files are available without a subscription to ACS Web Editions. Such files may be downloaded by article for research use (if there is a public use license linked to the relevant article, that license may permit other uses). Permission may be obtained from ACS for other uses through requests via the RightsLink permission system: http://pubs.acs.org/page/copyright/permissions.html.

Author Information

ARTICLE SECTIONS
Jump To

  • Corresponding Author
    • Daria Ruth Galimberti - Institut für Chemie, Humboldt-Universität, Unter den Linden 6, 10117 Berlin, GermanyInstitute for Molecules and Materials, Radboud University, Heyendaalseweg 135, 6525 AJ Nijmegen, The NetherlandsOrcidhttps://orcid.org/0000-0003-2766-3325 Email: [email protected]
  • Author
  • Notes
    The authors declare no competing financial interest.

Acknowledgments

ARTICLE SECTIONS
Jump To

This work has been supported by the Alexander von Humboldt Foundation with a fellowship to DG and by German Research Foundation (DFG) with a Reinhart Koselleck grant, as well as by the “Fonds der Chemischen Industrie”. The authors gratefully acknowledge the Gauss Center for Supercomputing e.V. (www.gauss-centre.eu) for providing computing time on the GCS Supercomputer SuperMUC-NG at Leibniz Supercomputing Center (www.lrz.de), Project ID pn68qu. Discussions with Dr Marcin Rybicki and Fabian Berger are gratefully acknowledged.

References

ARTICLE SECTIONS
Jump To

This article references 61 other publications.

  1. 1
    Hobza, P.; Sauer, J.; Morgeneyer, C.; Hurych, J.; Zahradnik, R. Bonding ability of surface sites on silica and their effect on hydrogen bonds. A quantum-chemical and statistical thermodynamic treatment. J. Phys. Chem. 1981, 85, 40614067,  DOI: 10.1021/j150626a022
  2. 2
    Sauer, J.; Zahradník, R. Quantum Chemical Studies on Zeolites and Silica. Int. J. Quantum Chem. 1984, 26, 793822,  DOI: 10.1002/qua.560260519
  3. 3
    Collinge, G.; Yuk, S. F.; Nguyen, M.-T.; Lee, M.-S.; Glezakou, V.-A.; Rousseau, R. Effect of Collective Dynamics and Anharmonicity on Entropy in Heterogenous Catalysis: Building the Case for Advanced Molecular Simulations. ACS Catal. 2020, 10, 92369260,  DOI: 10.1021/acscatal.0c01501
  4. 4
    Sauer, J. Ab Initio Calculations for Molecule–Surface Interactions with Chemical Accuracy. Acc. Chem. Res. 2019, 52, 35023510,  DOI: 10.1021/acs.accounts.9b00506
  5. 5
    Brandenburg, J. G.; Zen, A.; Fitzner, M.; Ramberger, B.; Kresse, G.; Tsatsoulis, T.; Grüneis, A.; Michaelides, A.; Alfè, D. Physisorption of Water on Graphene: Subchemical Accuracy from Many-Body Electronic Structure Methods. J. Phys. Chem. Lett. 2019, 10, 358368,  DOI: 10.1021/acs.jpclett.8b03679
  6. 6
    Sauer, J.; Ugliengo, P.; Garrone, E.; Saunders, V. R. Theoretical Study of van der Waals Complexes at Surface Sites in Comparison with the Experiment. Chem. Rev. 1994, 94, 20952160,  DOI: 10.1021/cr00031a014
  7. 7
    Hofmann, A.; Sauer, J. The surface structure of hydroxylated and sulphated zirconia. A periodic density-functional study. J. Phys. Chem. B 2004, 108, 1465214662,  DOI: 10.1021/jp049220f
  8. 8
    Laio, A.; Rodriguez-Fortea, A.; Gervasio, F. L.; Ceccarelli, M.; Parrinello, M. Assessing the Accuracy of Metadynamics. J. Phys. Chem. B 2005, 109, 67146721,  DOI: 10.1021/jp045424k
  9. 9
    Trzesniak, D.; Kunz, A.-P. E.; van Gunsteren, W. F. A Comparison of Methods to Compute the Potential of Mean Force. ChemPhysChem 2007, 8, 162169,  DOI: 10.1002/cphc.200600527
  10. 10
    Bučko, T.; Benco, L.; Hafner, J.; Ángyán, J. Proton exchange of small hydrocarbons over acidic chabazite: Ab initio study of entropic effects. J. Catal. 2007, 250, 171183,  DOI: 10.1016/j.jcat.2007.05.025
  11. 11
    Rey, J.; Gomez, A.; Raybaud, P.; Chizallet, C.; Bučko, T. On the origin of the difference between type A and type B skeletal isomerization of alkenes catalyzed by zeolites: The crucial input of ab initio molecular dynamics. J. Catal. 2019, 373, 361373,  DOI: 10.1016/j.jcat.2019.04.014
  12. 12
    Barducci, A.; Bussi, G.; Parrinello, M. Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method. Phys. Rev. Lett. 2008, 100, 020603,  DOI: 10.1103/PhysRevLett.100.020603
  13. 13
    Muñoz-Santiburcio, D.; Marx, D. Chemistry in nanoconfined water. Chem. Sci. 2017, 8, 34443452,  DOI: 10.1039/c6sc04989c
  14. 14
    Steinmann, C.; Olsson, M. A.; Ryde, U. Relative Ligand-Binding Free Energies Calculated from Multiple Short QM/MM MD Simulations. J. Chem. Theory Comput. 2018, 14, 32283237,  DOI: 10.1021/acs.jctc.8b00081
  15. 15
    Amsler, J.; Plessow, P. N.; Studt, F.; Bučko, T. Anharmonic Correction to Adsorption Free Energy from DFT-Based MD Using Thermodynamic Integration. J. Chem. Theory Comput. 2021, 17, 11551169,  DOI: 10.1021/acs.jctc.0c01022
  16. 16
    Rod, T. H.; Ryde, U. Accurate QM/MM Free Energy Calculations of Enzyme Reactions: Methylation by Catechol O-Methyltransferase. J. Chem. Theory Comput. 2005, 1, 12401251,  DOI: 10.1021/ct0501102
  17. 17
    Smit, B.; Maesen, T. L. M. Molecular Simulations of Zeolites: Adsorption, Diffusion, and Shape Selectivity. Chem. Rev. 2008, 108, 41254184,  DOI: 10.1021/cr8002642
  18. 18
    Bučko, T.; Benco, L.; Hafner, J.; Ángyán, J. G. Monomolecular cracking of propane over acidic chabazite: An ab initio molecular dynamics and transition path sampling study. J. Catal. 2011, 279, 220228
  19. 19
    Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 1256212566,  DOI: 10.1073/pnas.202427399
  20. 20
    Pietrucci, F.; Saitta, A. M. Formamide reaction network in gas phase and solution via a unified theoretical approach: Toward a reconciliation of different prebiotic scenarios. Proc. Natl. Acad. Sci. U.S.A. 2015, 112, 1503015035,  DOI: 10.1073/pnas.1512486112
  21. 21
    Alexopoulos, K.; Lee, M.-S.; Liu, Y.; Zhi, Y.; Liu, Y.; Reyniers, M.-F.; Marin, G. B.; Glezakou, V.-A.; Rousseau, R.; Lercher, J. A. Anharmonicity and Confinement in Zeolites: Structure, Spectroscopy, and Adsorption Free Energy of Ethanol in H-ZSM-5. J. Phys. Chem. C 2016, 120, 71727182,  DOI: 10.1021/acs.jpcc.6b00923
  22. 22
    Cnudde, P.; De Wispelaere, K.; Vanduyfhuys, L.; Demuynck, R.; Van der Mynsbrugge, J.; Waroquier, M.; Van Speybroeck, V. How Chain Length and Branching Influence the Alkene Cracking Reactivity on H-ZSM-5. ACS Catal. 2018, 8, 95799595,  DOI: 10.1021/acscatal.8b01779
  23. 23
    Cnudde, P.; De Wispelaere, K.; Van der Mynsbrugge, J.; Waroquier, M.; Van Speybroeck, V. Effect of temperature and branching on the nature and stability of alkene cracking intermediates in H-ZSM-5. J. Catal. 2017, 345, 5369,  DOI: 10.1016/j.jcat.2016.11.010
  24. 24
    Janda, A.; Vlaisavljevich, B.; Lin, L.-C.; Mallikarjun Sharada, S.; Smit, B.; Head-Gordon, M.; Bell, A. T. Adsorption Thermodynamics and Intrinsic Activation Parameters for Monomolecular Cracking of n-Alkanes on Brønsted Acid Sites in Zeolites. J. Phys. Chem. C 2015, 119, 1042710438,  DOI: 10.1021/acs.jpcc.5b01715
  25. 25
    Janda, A.; Vlaisavljevich, B.; Smit, B.; Lin, L.-C.; Bell, A. T. Effects of Pore and Cage Topology on the Thermodynamics of n-Alkane Adsorption at Brønsted Protons in Zeolites at High Temperature. J. Phys. Chem. C 2017, 121, 16181638,  DOI: 10.1021/acs.jpcc.6b09703
  26. 26
    Fetisov, E. O.; Shah, M. S.; Long, J. R.; Tsapatsis, M.; Siepmann, J. I. First principles Monte Carlo simulations of unary and binary adsorption: CO2, N2, and H2O in Mg-MOF-74. Chem. Commun. 2018, 54, 1081610819,  DOI: 10.1039/c8cc06178e
  27. 27
    Rocca, D.; Dixit, A.; Badawi, M.; Lebègue, S.; Gould, T.; Bučko, T. Bridging molecular dynamics and correlated wave-function methods for accurate finite-temperature properties. Phys. Rev. Mater. 2019, 3, 040801,  DOI: 10.1103/physrevmaterials.3.040801
  28. 28
    Piccini, G.; Sauer, J. Quantum Chemical Free Energies: Structure Optimization and Vibrational Frequencies in Normal Modes. J. Chem. Theory Comput. 2013, 9, 50385045,  DOI: 10.1021/ct4005504
  29. 29
    Piccini, G.; Sauer, J. Effect of Anharmonicity on Adsorption Thermodynamics. J. Chem. Theory Comput. 2014, 10, 24792487,  DOI: 10.1021/ct500291x
  30. 30
    Piccini, G.; Alessio, M.; Sauer, J.; Zhi, Y.; Liu, Y.; Kolvenbach, R.; Jentys, A.; Lercher, J. A. Accurate Adsorption Thermodynamics of Small Alkanes in Zeolites. Ab initio Theory and Experiment for H-Chabazite. J. Phys. Chem. C 2015, 119, 61286137,  DOI: 10.1021/acs.jpcc.5b01739
  31. 31
    Njegic, B.; Gordon, M. S. Exploring the effect of anharmonicity of molecular vibrations on thermodynamic properties. J. Chem. Phys. 2006, 125, 224102,  DOI: 10.1063/1.2395940
  32. 32
    Kundu, A.; Piccini, G.; Sillar, K.; Sauer, J. Ab Initio Prediction of Adsorption Isotherms for Small Molecules in Metal–Organic Frameworks. J. Am. Chem. Soc. 2016, 138, 1404714056,  DOI: 10.1021/jacs.6b08646
  33. 33
    Piccini, G.; Alessio, M.; Sauer, J. Ab Initio Calculation of Rate Constants for Molecule–Surface Reactions with Chemical Accuracy. Angew. Chem., Int. Ed. 2016, 55, 52355237,  DOI: 10.1002/anie.201601534
  34. 34
    Bučko, T.; Hafner, J.; Ángyán, J. G. Geometry optimization of periodic systems using internal coordinates. J. Chem. Phys. 2005, 122, 124508
  35. 35
    Zhang, D.-B.; Sun, T.; Wentzcovitch, R. M. Phonon Quasiparticles and Anharmonic Free Energy in Complex Systems. Phys. Rev. Lett. 2014, 112, 058501,  DOI: 10.1103/PhysRevLett.112.058501
  36. 36
    Carreras, A.; Togo, A.; Tanaka, I. DynaPhoPy: A code for extracting phonon quasiparticles from molecular dynamics simulations. Comput. Phys. Commun. 2017, 221, 221234,  DOI: 10.1016/j.cpc.2017.08.017
  37. 37
    Peters, L. D. M.; Dietschreit, J. C. B.; Kussmann, J.; Ochsenfeld, C. Calculating free energies from the vibrational density of states function: Validation and critical assessment. J. Chem. Phys. 2019, 150, 194111,  DOI: 10.1063/1.5079643
  38. 38
    Galimberti, D. R.; Milani, A.; Tommasini, M.; Castiglioni, C.; Gaigeot, M.-P. Combining Static and Dynamical Approaches for Infrared Spectra Calculations of Gas Phase Molecules and Clusters. J. Chem. Theory Comput. 2017, 13, 38023813,  DOI: 10.1021/acs.jctc.7b00471
  39. 39
    Brehm, M.; Kirchner, B. TRAVIS - A Free Analyzer and Visualizer for Monte Carlo and Molecular Dynamics Trajectories. J. Chem. Inf. Model. 2011, 51, 20072023,  DOI: 10.1021/ci200217w
  40. 40
    Nielsen, M.; Brogaard, R. Y.; Falsig, H.; Beato, P.; Swang, O.; Svelle, S. Kinetics of Zeolite Dealumination: Insights from H-SSZ-13. ACS Catal. 2015, 5, 71317139,  DOI: 10.1021/acscatal.5b01496
  41. 41
    Shi, L.; Yang, J.; Shen, G.; Zhao, Y.; Chen, R.; Shen, M.; Wen, Y.; Shan, B. The influence of adjacent Al atoms on the hydrothermal stability of H-SSZ-13: a first-principles study. Phys. Chem. Chem. Phys. 2020, 22, 29302937,  DOI: 10.1039/c9cp05141d
  42. 42
    Plessow, P. N.; Studt, F. Olefin methylation and cracking reactions in H-SSZ-13 investigated with ab initio and DFT calculations. Catal. Sci. Technol. 2018, 8, 44204429,  DOI: 10.1039/c8cy01194j
  43. 43
    Sarazen, M. L.; Doskocil, E.; Iglesia, E. Effects of Void Environment and Acid Strength on Alkene Oligomerization Selectivity. ACS Catal. 2016, 6, 70597070,  DOI: 10.1021/acscatal.6b02128
  44. 44
    McQuarrie, D. A. Statistical Thermodynamics; Harper & Row, 1973.
  45. 45
    Balog, E.; Becker, T.; Oettl, M.; Lechner, R.; Daniel, R.; Finney, J.; Smith, J. C. Direct Determination of Vibrational Density of States Change on Ligand Binding to a Protein. Phys. Rev. Lett. 2004, 93, 028103,  DOI: 10.1103/PhysRevLett.93.028103
  46. 46
    Lin, S.-T.; Maiti, P. K.; Goddard, W. A. Dynamics and Thermodynamics of Water in PAMAM Dendrimers at Subnanosecond Time Scales. J. Phys. Chem. B 2005, 109, 86638672,  DOI: 10.1021/jp0471958
  47. 47
    Sousa, R. L. d.; Alves, H. W. L. Ab initio calculation of the dynamical properties of PPP and PPV. Braz. J. Phys. 2006, 36, 501504,  DOI: 10.1590/s0103-97332006000300072
  48. 48
    Sun, T.; Zhang, D.-B.; Wentzcovitch, R. M. Dynamic stabilization of cubic CaSiO3 perovskite at high temperatures and pressures from ab initio molecular dynamics. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 89, 094109,  DOI: 10.1103/physrevb.89.094109
  49. 49
    Landau, L. D.; Lifshitz, E. M. Statistical Physics; Elsevier, 2013; Vol. 5.
  50. 50
    Jones, W.; March, N. H. Theoretical Solid State Physics; Courier Corporation, 1985; Vol. 35.
  51. 51
    Berens, P. H.; Mackay, D. H. J.; White, G. M.; Wilson, K. R. Thermodynamics and quantum corrections from molecular dynamics for liquid water. J. Chem. Phys. 1983, 79, 23752389,  DOI: 10.1063/1.446044
  52. 52
    Kearsley, S. K. On the orthogonal transformation used for structural comparisons. Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr. 1989, 45, 208210,  DOI: 10.1107/s0108767388010128
  53. 53
    Kudin, K. N.; Dymarsky, A. Y. Eckart axis conditions and the minimization of the root-mean-square deviation: Two closely related problems. J. Chem. Phys. 2005, 122, 224105,  DOI: 10.1063/1.1929739
  54. 54
    Mathias, G.; Ivanov, S. D.; Witt, A.; Baer, M. D.; Marx, D. Infrared Spectroscopy of Fluxional Molecules from (ab Initio) Molecular Dynamics: Resolving Large-Amplitude Motion, Multiple Conformations, and Permutational Symmetries. J. Chem. Theory Comput. 2012, 8, 224234,  DOI: 10.1021/ct2006665
  55. 55
    Iannuzzi, M.; Laio, A.; Parrinello, M. Efficient Exploration of Reactive Potential Energy Surfaces Using Car-Parrinello Molecular Dynamics. Phys. Rev. Lett. 2003, 90, 238302,  DOI: 10.1103/physrevlett.90.238302
  56. 56
    Cazzaniga, M.; Micciarelli, M.; Moriggi, F.; Mahmoud, A.; Gabas, F.; Ceotto, M. Anharmonic calculations of vibrational spectra for molecular adsorbates: A divide-and-conquer semiclassical molecular dynamics approach. J. Chem. Phys. 2020, 152, 104104,  DOI: 10.1063/1.5142682
  57. 57
    Bates, S. P.; Van Well, W. J. M.; Van Santen, R. A.; Smit, B. Configurational-Bias Monte Carlo (CB-MC) Calculations of n-Alkane Sorption in Zeolites Rho and Fer. Mol. Simul. 1997, 19, 301318,  DOI: 10.1080/08927029708024159
  58. 58
    Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 1116911186,  DOI: 10.1103/physrevb.54.11169
  59. 59
    Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple [Phys. Rev. Lett. 77, 3865 (1996)]. Phys. Rev. Lett. 1997, 78, 1396,  DOI: 10.1103/physrevlett.78.1396
  60. 60
    Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. J. Comput. Chem. 2006, 27, 17871799,  DOI: 10.1002/jcc.20495
  61. 61
    Wilson, E. B.; Decius, J. C.; Cross, P. C. Molecular Vibrations: The Theory of Infrared and Raman Vibrational Spectra; Courier Corporation, 1980.

Cited By

ARTICLE SECTIONS
Jump To

This article is cited by 9 publications.

  1. Kacper Drużbicki, Pablo Gila-Herranz, Pelayo Marin-Villa, Mattia Gaboardi, Jeff Armstrong, Felix Fernandez-Alonso. Cation Dynamics as Structure Explorer in Hybrid Perovskites─The Case of MAPbI3. Crystal Growth & Design 2024, 24 (1) , 391-404. https://doi.org/10.1021/acs.cgd.3c01112
  2. Veronique Van Speybroeck, Massimo Bocus, Pieter Cnudde, Louis Vanduyfhuys. Operando Modeling of Zeolite-Catalyzed Reactions Using First-Principles Molecular Dynamics Simulations. ACS Catalysis 2023, 13 (17) , 11455-11493. https://doi.org/10.1021/acscatal.3c01945
  3. Céline Chizallet, Christophe Bouchy, Kim Larmier, Gerhard Pirngruber. Molecular Views on Mechanisms of Brønsted Acid-Catalyzed Reactions in Zeolites. Chemical Reviews 2023, 123 (9) , 6107-6196. https://doi.org/10.1021/acs.chemrev.2c00896
  4. Fabian Berger, Marcin Rybicki, Joachim Sauer. Molecular Dynamics with Chemical Accuracy─Alkane Adsorption in Acidic Zeolites. ACS Catalysis 2023, 13 (3) , 2011-2024. https://doi.org/10.1021/acscatal.2c05493
  5. Massimo Bocus, Veronique Van Speybroeck. Insights into the Mechanism and Reactivity of Zeolite-Catalyzed Alkylphenol Dealkylation. ACS Catalysis 2022, 12 (22) , 14227-14242. https://doi.org/10.1021/acscatal.2c03844
  6. Tom Braeckevelt, Ruben Goeminne, Sander Vandenhaute, Sander Borgmans, Toon Verstraelen, Julian A. Steele, Maarten B. J. Roeffaers, Johan Hofkens, Sven M. J. Rogge, Veronique Van Speybroeck. Accurately Determining the Phase Transition Temperature of CsPbI3 via Random-Phase Approximation Calculations and Phase-Transferable Machine Learning Potentials. Chemistry of Materials 2022, 34 (19) , 8561-8576. https://doi.org/10.1021/acs.chemmater.2c01508
  7. Daria Ruth Galimberti. Vibrational Circular Dichroism from DFT Molecular Dynamics: The AWV Method. Journal of Chemical Theory and Computation 2022, 18 (10) , 6217-6230. https://doi.org/10.1021/acs.jctc.2c00736
  8. Kristof De Wispelaere, Philipp N. Plessow, Felix Studt. Toward Computing Accurate Free Energies in Heterogeneous Catalysis: a Case Study for Adsorbed Isobutene in H-ZSM-5. ACS Physical Chemistry Au 2022, 2 (5) , 399-406. https://doi.org/10.1021/acsphyschemau.2c00020
  9. Marcin Rybicki, Joachim Sauer. Rigid Body Approximation for the Anharmonic Description of Molecule–Surface Vibrations. Journal of Chemical Theory and Computation 2022, 18 (9) , 5618-5635. https://doi.org/10.1021/acs.jctc.2c00597
  • Abstract

    Figure 1

    Figure 1. Zeolite chabazite, H-CHA with an Si/Al ratio of 11/1, and the proton attached to the O2 position: monoclinic 2 × 1 × 1 supercell (A), adsorption complexes of methane (B), ethane (C), and propane (D). Color code: yellow—Si, red—O, blue—Al, gray—C, and white—H. Panel (E): harmonic VDOSs for the adsorption complex of ethane in H-CHA, modes belonging to adsorbed ethane (top), the O–H of the Si–O(H)–Al bridge (middle), and the framework (bottom), with intensities renormalized for the number of atoms for each subsystem. Panels (A–D) adapted with permission from ref (30).

    Figure 2

    Figure 2. Left panel: VDOS of the ethane molecule in the gas phase at 313 K computed from the autocorrelation function of cartesian velocities (eq 6). Right panel: contribution of the Si–O stretching vibration [bridging Si–O(H)–Al groups] for the bare H-CHA to the total VDOS between 700 and 850 cm–1 computed by projecting the velocities on the harmonic modes without (blue line) or with applying a Lorentzian window filter (red line).

    Figure 3

    Figure 3. Frequency shift (cm–1) compared to the harmonic wavenumber (see figure inset for color coding) for the vibrational modes of the bare H-CHA surface, 0.25, 0.50, 1.0, and 2.0 ps after the initial vertical excitation. The level of excitation is qualitatively measured as a wavenumber shift from the harmonic wavenumbers.

    Figure 4

    Figure 4. Absolute deviation of the computed thermodynamic functions from the experimental values for methane (blue), ethane (red), and propane (green) in H-CHA at experimental conditions (303 K for methane and 313 K for ethane and propane, 0.1 MPa, θ = 0.5). The gray shadows represent the statistical uncertainties (see the Supporting Information for more details). Upper panel: adsorption free energy, ΔGa. Center: adsorption enthalpy, ΔHa. Lower panel: adsorption entropy term, −TΔSa. On the left: the harmonic (HARM) and anharmonic (ANH) predictions of Piccini et al. (30) Center panel: DOS-P derived ΔE0-viba. Right panels: harmonic ΔE0-viba. The average deviation for each method between the three molecules is reported above the columns for the three molecules.

    Figure 5

    Figure 5. ZPE of adsorption, ΔE0-viba (kJ/mol), for methane, ethane, and propane in H-CHA, calculated from harmonic frequencies (blue dotted line), from anharmonic frequencies derived from curvilinear distortions (30) (red dotted line), and obtained with DOS-P methods by initially exciting the system at 300 K (TDOS-P, red full line) or at each mode’s ZPE (ZPEDOS-P, green full line).

    Figure 6

    Figure 6. Adsorption Gibbs free energy, ΔGa, enthalpy, ΔHa, and entropy term, −TΔSa, for alkanes in H-CHA at experimental conditions (303 K for methane and 313 K for ethane and propane, 0.1 MPa, θ = 0.5) computed with projection of the VDOS using ZPE excitation (ZPEDOS-P, green full line) or T excitation (TDOS-P, red full line), compared to the harmonic (HARM, blue dotted line), and the anharmonic (ANH, red dotted line) predictions. (30) Experimental values of ref (30) (black dots) are shown with the chemical accuracy range (±4 kJ/mol).

  • References

    ARTICLE SECTIONS
    Jump To

    This article references 61 other publications.

    1. 1
      Hobza, P.; Sauer, J.; Morgeneyer, C.; Hurych, J.; Zahradnik, R. Bonding ability of surface sites on silica and their effect on hydrogen bonds. A quantum-chemical and statistical thermodynamic treatment. J. Phys. Chem. 1981, 85, 40614067,  DOI: 10.1021/j150626a022
    2. 2
      Sauer, J.; Zahradník, R. Quantum Chemical Studies on Zeolites and Silica. Int. J. Quantum Chem. 1984, 26, 793822,  DOI: 10.1002/qua.560260519
    3. 3
      Collinge, G.; Yuk, S. F.; Nguyen, M.-T.; Lee, M.-S.; Glezakou, V.-A.; Rousseau, R. Effect of Collective Dynamics and Anharmonicity on Entropy in Heterogenous Catalysis: Building the Case for Advanced Molecular Simulations. ACS Catal. 2020, 10, 92369260,  DOI: 10.1021/acscatal.0c01501
    4. 4
      Sauer, J. Ab Initio Calculations for Molecule–Surface Interactions with Chemical Accuracy. Acc. Chem. Res. 2019, 52, 35023510,  DOI: 10.1021/acs.accounts.9b00506
    5. 5
      Brandenburg, J. G.; Zen, A.; Fitzner, M.; Ramberger, B.; Kresse, G.; Tsatsoulis, T.; Grüneis, A.; Michaelides, A.; Alfè, D. Physisorption of Water on Graphene: Subchemical Accuracy from Many-Body Electronic Structure Methods. J. Phys. Chem. Lett. 2019, 10, 358368,  DOI: 10.1021/acs.jpclett.8b03679
    6. 6
      Sauer, J.; Ugliengo, P.; Garrone, E.; Saunders, V. R. Theoretical Study of van der Waals Complexes at Surface Sites in Comparison with the Experiment. Chem. Rev. 1994, 94, 20952160,  DOI: 10.1021/cr00031a014
    7. 7
      Hofmann, A.; Sauer, J. The surface structure of hydroxylated and sulphated zirconia. A periodic density-functional study. J. Phys. Chem. B 2004, 108, 1465214662,  DOI: 10.1021/jp049220f
    8. 8
      Laio, A.; Rodriguez-Fortea, A.; Gervasio, F. L.; Ceccarelli, M.; Parrinello, M. Assessing the Accuracy of Metadynamics. J. Phys. Chem. B 2005, 109, 67146721,  DOI: 10.1021/jp045424k
    9. 9
      Trzesniak, D.; Kunz, A.-P. E.; van Gunsteren, W. F. A Comparison of Methods to Compute the Potential of Mean Force. ChemPhysChem 2007, 8, 162169,  DOI: 10.1002/cphc.200600527
    10. 10
      Bučko, T.; Benco, L.; Hafner, J.; Ángyán, J. Proton exchange of small hydrocarbons over acidic chabazite: Ab initio study of entropic effects. J. Catal. 2007, 250, 171183,  DOI: 10.1016/j.jcat.2007.05.025
    11. 11
      Rey, J.; Gomez, A.; Raybaud, P.; Chizallet, C.; Bučko, T. On the origin of the difference between type A and type B skeletal isomerization of alkenes catalyzed by zeolites: The crucial input of ab initio molecular dynamics. J. Catal. 2019, 373, 361373,  DOI: 10.1016/j.jcat.2019.04.014
    12. 12
      Barducci, A.; Bussi, G.; Parrinello, M. Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method. Phys. Rev. Lett. 2008, 100, 020603,  DOI: 10.1103/PhysRevLett.100.020603
    13. 13
      Muñoz-Santiburcio, D.; Marx, D. Chemistry in nanoconfined water. Chem. Sci. 2017, 8, 34443452,  DOI: 10.1039/c6sc04989c
    14. 14
      Steinmann, C.; Olsson, M. A.; Ryde, U. Relative Ligand-Binding Free Energies Calculated from Multiple Short QM/MM MD Simulations. J. Chem. Theory Comput. 2018, 14, 32283237,  DOI: 10.1021/acs.jctc.8b00081
    15. 15
      Amsler, J.; Plessow, P. N.; Studt, F.; Bučko, T. Anharmonic Correction to Adsorption Free Energy from DFT-Based MD Using Thermodynamic Integration. J. Chem. Theory Comput. 2021, 17, 11551169,  DOI: 10.1021/acs.jctc.0c01022
    16. 16
      Rod, T. H.; Ryde, U. Accurate QM/MM Free Energy Calculations of Enzyme Reactions: Methylation by Catechol O-Methyltransferase. J. Chem. Theory Comput. 2005, 1, 12401251,  DOI: 10.1021/ct0501102
    17. 17
      Smit, B.; Maesen, T. L. M. Molecular Simulations of Zeolites: Adsorption, Diffusion, and Shape Selectivity. Chem. Rev. 2008, 108, 41254184,  DOI: 10.1021/cr8002642
    18. 18
      Bučko, T.; Benco, L.; Hafner, J.; Ángyán, J. G. Monomolecular cracking of propane over acidic chabazite: An ab initio molecular dynamics and transition path sampling study. J. Catal. 2011, 279, 220228
    19. 19
      Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 1256212566,  DOI: 10.1073/pnas.202427399
    20. 20
      Pietrucci, F.; Saitta, A. M. Formamide reaction network in gas phase and solution via a unified theoretical approach: Toward a reconciliation of different prebiotic scenarios. Proc. Natl. Acad. Sci. U.S.A. 2015, 112, 1503015035,  DOI: 10.1073/pnas.1512486112
    21. 21
      Alexopoulos, K.; Lee, M.-S.; Liu, Y.; Zhi, Y.; Liu, Y.; Reyniers, M.-F.; Marin, G. B.; Glezakou, V.-A.; Rousseau, R.; Lercher, J. A. Anharmonicity and Confinement in Zeolites: Structure, Spectroscopy, and Adsorption Free Energy of Ethanol in H-ZSM-5. J. Phys. Chem. C 2016, 120, 71727182,  DOI: 10.1021/acs.jpcc.6b00923
    22. 22
      Cnudde, P.; De Wispelaere, K.; Vanduyfhuys, L.; Demuynck, R.; Van der Mynsbrugge, J.; Waroquier, M.; Van Speybroeck, V. How Chain Length and Branching Influence the Alkene Cracking Reactivity on H-ZSM-5. ACS Catal. 2018, 8, 95799595,  DOI: 10.1021/acscatal.8b01779
    23. 23
      Cnudde, P.; De Wispelaere, K.; Van der Mynsbrugge, J.; Waroquier, M.; Van Speybroeck, V. Effect of temperature and branching on the nature and stability of alkene cracking intermediates in H-ZSM-5. J. Catal. 2017, 345, 5369,  DOI: 10.1016/j.jcat.2016.11.010
    24. 24
      Janda, A.; Vlaisavljevich, B.; Lin, L.-C.; Mallikarjun Sharada, S.; Smit, B.; Head-Gordon, M.; Bell, A. T. Adsorption Thermodynamics and Intrinsic Activation Parameters for Monomolecular Cracking of n-Alkanes on Brønsted Acid Sites in Zeolites. J. Phys. Chem. C 2015, 119, 1042710438,  DOI: 10.1021/acs.jpcc.5b01715
    25. 25
      Janda, A.; Vlaisavljevich, B.; Smit, B.; Lin, L.-C.; Bell, A. T. Effects of Pore and Cage Topology on the Thermodynamics of n-Alkane Adsorption at Brønsted Protons in Zeolites at High Temperature. J. Phys. Chem. C 2017, 121, 16181638,  DOI: 10.1021/acs.jpcc.6b09703
    26. 26
      Fetisov, E. O.; Shah, M. S.; Long, J. R.; Tsapatsis, M.; Siepmann, J. I. First principles Monte Carlo simulations of unary and binary adsorption: CO2, N2, and H2O in Mg-MOF-74. Chem. Commun. 2018, 54, 1081610819,  DOI: 10.1039/c8cc06178e
    27. 27
      Rocca, D.; Dixit, A.; Badawi, M.; Lebègue, S.; Gould, T.; Bučko, T. Bridging molecular dynamics and correlated wave-function methods for accurate finite-temperature properties. Phys. Rev. Mater. 2019, 3, 040801,  DOI: 10.1103/physrevmaterials.3.040801
    28. 28
      Piccini, G.; Sauer, J. Quantum Chemical Free Energies: Structure Optimization and Vibrational Frequencies in Normal Modes. J. Chem. Theory Comput. 2013, 9, 50385045,  DOI: 10.1021/ct4005504
    29. 29
      Piccini, G.; Sauer, J. Effect of Anharmonicity on Adsorption Thermodynamics. J. Chem. Theory Comput. 2014, 10, 24792487,  DOI: 10.1021/ct500291x
    30. 30
      Piccini, G.; Alessio, M.; Sauer, J.; Zhi, Y.; Liu, Y.; Kolvenbach, R.; Jentys, A.; Lercher, J. A. Accurate Adsorption Thermodynamics of Small Alkanes in Zeolites. Ab initio Theory and Experiment for H-Chabazite. J. Phys. Chem. C 2015, 119, 61286137,  DOI: 10.1021/acs.jpcc.5b01739
    31. 31
      Njegic, B.; Gordon, M. S. Exploring the effect of anharmonicity of molecular vibrations on thermodynamic properties. J. Chem. Phys. 2006, 125, 224102,  DOI: 10.1063/1.2395940
    32. 32
      Kundu, A.; Piccini, G.; Sillar, K.; Sauer, J. Ab Initio Prediction of Adsorption Isotherms for Small Molecules in Metal–Organic Frameworks. J. Am. Chem. Soc. 2016, 138, 1404714056,  DOI: 10.1021/jacs.6b08646
    33. 33
      Piccini, G.; Alessio, M.; Sauer, J. Ab Initio Calculation of Rate Constants for Molecule–Surface Reactions with Chemical Accuracy. Angew. Chem., Int. Ed. 2016, 55, 52355237,  DOI: 10.1002/anie.201601534
    34. 34
      Bučko, T.; Hafner, J.; Ángyán, J. G. Geometry optimization of periodic systems using internal coordinates. J. Chem. Phys. 2005, 122, 124508
    35. 35
      Zhang, D.-B.; Sun, T.; Wentzcovitch, R. M. Phonon Quasiparticles and Anharmonic Free Energy in Complex Systems. Phys. Rev. Lett. 2014, 112, 058501,  DOI: 10.1103/PhysRevLett.112.058501
    36. 36
      Carreras, A.; Togo, A.; Tanaka, I. DynaPhoPy: A code for extracting phonon quasiparticles from molecular dynamics simulations. Comput. Phys. Commun. 2017, 221, 221234,  DOI: 10.1016/j.cpc.2017.08.017
    37. 37
      Peters, L. D. M.; Dietschreit, J. C. B.; Kussmann, J.; Ochsenfeld, C. Calculating free energies from the vibrational density of states function: Validation and critical assessment. J. Chem. Phys. 2019, 150, 194111,  DOI: 10.1063/1.5079643
    38. 38
      Galimberti, D. R.; Milani, A.; Tommasini, M.; Castiglioni, C.; Gaigeot, M.-P. Combining Static and Dynamical Approaches for Infrared Spectra Calculations of Gas Phase Molecules and Clusters. J. Chem. Theory Comput. 2017, 13, 38023813,  DOI: 10.1021/acs.jctc.7b00471
    39. 39
      Brehm, M.; Kirchner, B. TRAVIS - A Free Analyzer and Visualizer for Monte Carlo and Molecular Dynamics Trajectories. J. Chem. Inf. Model. 2011, 51, 20072023,  DOI: 10.1021/ci200217w
    40. 40
      Nielsen, M.; Brogaard, R. Y.; Falsig, H.; Beato, P.; Swang, O.; Svelle, S. Kinetics of Zeolite Dealumination: Insights from H-SSZ-13. ACS Catal. 2015, 5, 71317139,  DOI: 10.1021/acscatal.5b01496
    41. 41
      Shi, L.; Yang, J.; Shen, G.; Zhao, Y.; Chen, R.; Shen, M.; Wen, Y.; Shan, B. The influence of adjacent Al atoms on the hydrothermal stability of H-SSZ-13: a first-principles study. Phys. Chem. Chem. Phys. 2020, 22, 29302937,  DOI: 10.1039/c9cp05141d
    42. 42
      Plessow, P. N.; Studt, F. Olefin methylation and cracking reactions in H-SSZ-13 investigated with ab initio and DFT calculations. Catal. Sci. Technol. 2018, 8, 44204429,  DOI: 10.1039/c8cy01194j
    43. 43
      Sarazen, M. L.; Doskocil, E.; Iglesia, E. Effects of Void Environment and Acid Strength on Alkene Oligomerization Selectivity. ACS Catal. 2016, 6, 70597070,  DOI: 10.1021/acscatal.6b02128
    44. 44
      McQuarrie, D. A. Statistical Thermodynamics; Harper & Row, 1973.
    45. 45
      Balog, E.; Becker, T.; Oettl, M.; Lechner, R.; Daniel, R.; Finney, J.; Smith, J. C. Direct Determination of Vibrational Density of States Change on Ligand Binding to a Protein. Phys. Rev. Lett. 2004, 93, 028103,  DOI: 10.1103/PhysRevLett.93.028103
    46. 46
      Lin, S.-T.; Maiti, P. K.; Goddard, W. A. Dynamics and Thermodynamics of Water in PAMAM Dendrimers at Subnanosecond Time Scales. J. Phys. Chem. B 2005, 109, 86638672,  DOI: 10.1021/jp0471958
    47. 47
      Sousa, R. L. d.; Alves, H. W. L. Ab initio calculation of the dynamical properties of PPP and PPV. Braz. J. Phys. 2006, 36, 501504,  DOI: 10.1590/s0103-97332006000300072
    48. 48
      Sun, T.; Zhang, D.-B.; Wentzcovitch, R. M. Dynamic stabilization of cubic CaSiO3 perovskite at high temperatures and pressures from ab initio molecular dynamics. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 89, 094109,  DOI: 10.1103/physrevb.89.094109
    49. 49
      Landau, L. D.; Lifshitz, E. M. Statistical Physics; Elsevier, 2013; Vol. 5.
    50. 50
      Jones, W.; March, N. H. Theoretical Solid State Physics; Courier Corporation, 1985; Vol. 35.
    51. 51
      Berens, P. H.; Mackay, D. H. J.; White, G. M.; Wilson, K. R. Thermodynamics and quantum corrections from molecular dynamics for liquid water. J. Chem. Phys. 1983, 79, 23752389,  DOI: 10.1063/1.446044
    52. 52
      Kearsley, S. K. On the orthogonal transformation used for structural comparisons. Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr. 1989, 45, 208210,  DOI: 10.1107/s0108767388010128
    53. 53
      Kudin, K. N.; Dymarsky, A. Y. Eckart axis conditions and the minimization of the root-mean-square deviation: Two closely related problems. J. Chem. Phys. 2005, 122, 224105,  DOI: 10.1063/1.1929739
    54. 54
      Mathias, G.; Ivanov, S. D.; Witt, A.; Baer, M. D.; Marx, D. Infrared Spectroscopy of Fluxional Molecules from (ab Initio) Molecular Dynamics: Resolving Large-Amplitude Motion, Multiple Conformations, and Permutational Symmetries. J. Chem. Theory Comput. 2012, 8, 224234,  DOI: 10.1021/ct2006665
    55. 55
      Iannuzzi, M.; Laio, A.; Parrinello, M. Efficient Exploration of Reactive Potential Energy Surfaces Using Car-Parrinello Molecular Dynamics. Phys. Rev. Lett. 2003, 90, 238302,  DOI: 10.1103/physrevlett.90.238302
    56. 56
      Cazzaniga, M.; Micciarelli, M.; Moriggi, F.; Mahmoud, A.; Gabas, F.; Ceotto, M. Anharmonic calculations of vibrational spectra for molecular adsorbates: A divide-and-conquer semiclassical molecular dynamics approach. J. Chem. Phys. 2020, 152, 104104,  DOI: 10.1063/1.5142682
    57. 57
      Bates, S. P.; Van Well, W. J. M.; Van Santen, R. A.; Smit, B. Configurational-Bias Monte Carlo (CB-MC) Calculations of n-Alkane Sorption in Zeolites Rho and Fer. Mol. Simul. 1997, 19, 301318,  DOI: 10.1080/08927029708024159
    58. 58
      Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 1116911186,  DOI: 10.1103/physrevb.54.11169
    59. 59
      Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple [Phys. Rev. Lett. 77, 3865 (1996)]. Phys. Rev. Lett. 1997, 78, 1396,  DOI: 10.1103/physrevlett.78.1396
    60. 60
      Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. J. Comput. Chem. 2006, 27, 17871799,  DOI: 10.1002/jcc.20495
    61. 61
      Wilson, E. B.; Decius, J. C.; Cross, P. C. Molecular Vibrations: The Theory of Infrared and Raman Vibrational Spectra; Courier Corporation, 1980.
  • Supporting Information

    Supporting Information

    ARTICLE SECTIONS
    Jump To

    The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jctc.1c00519.

    • Additional equations for the calculation of the free energies, decoupling of the vibrations due to the anharmonicity, additional computational details, data on adsorption of propane via secondary carbon, discussion on the vibrational active space, convergence check, and summary of the computed thermodynamic values (PDF)


    Terms & Conditions

    Most electronic Supporting Information files are available without a subscription to ACS Web Editions. Such files may be downloaded by article for research use (if there is a public use license linked to the relevant article, that license may permit other uses). Permission may be obtained from ACS for other uses through requests via the RightsLink permission system: http://pubs.acs.org/page/copyright/permissions.html.

Pair your accounts.

Export articles to Mendeley

Get article recommendations from ACS based on references in your Mendeley library.

Pair your accounts.

Export articles to Mendeley

Get article recommendations from ACS based on references in your Mendeley library.

You’ve supercharged your research process with ACS and Mendeley!

STEP 1:
Click to create an ACS ID

Please note: If you switch to a different device, you may be asked to login again with only your ACS ID.

Please note: If you switch to a different device, you may be asked to login again with only your ACS ID.

Please note: If you switch to a different device, you may be asked to login again with only your ACS ID.

MENDELEY PAIRING EXPIRED
Your Mendeley pairing has expired. Please reconnect