Bridging molecular dynamics and correlated wave-function methods for accurate finite-temperature properties

TL;DR

The paper presents the selPT perturbative approach that combines ab initio molecular dynamics with correlated wave-function methods to accurately predict finite-temperature properties, demonstrated on adsorption enthalpies in zeolites with excellent experimental agreement.

Contribution

The novel selPT method efficiently integrates AIMD with correlated wave-function techniques for precise finite-temperature property calculations.

Findings

Accurate adsorption enthalpies for CH₄ and CO₂ in zeolites at 300 K.

Excellent agreement with experimental data.

Advances toward quantitative AIMD predictions.

Abstract

We introduce the "selPT" perturbative approach, based on ab initio molecular dynamics (AIMD), for computing accurate finite-temperature properties by efficiently using correlated wave-function methods. We demonstrate the power of the method by computing prototypical molecular enthalpies of adsorption in zeolite (CH$_4$ and CO$_2$ on protonated chabazite at 300~K) using the random phase approximation. Results are in excellent agreement with experiment. The improved accuracy provided by selPT represents a crucial step towards the goal of truly quantitative AIMD prediction of experimental observables at finite temperature.

Tables

Table 1. Table 1: Enthalpy of adsorption Δ ads ​ H ​ ( T ) subscript Δ ads 𝐻 𝑇 \Delta_{\text{ads}}H(T) (in kcal/mol) of CH 4 in protonated chabazite at 300 K evaluated at the PBE+D2, PBE, and RPA levels of theory. Results are calculated using full MD trajectories (MD@PBE and MD@PBE+D2) applying the selPT method to a set with N G = 60 subscript 𝑁 𝐺 60 N_{G}=60 configurations (SET1@PBE+D2), N G = 120 subscript 𝑁 𝐺 120 N_{G}=120 (SET2+3@PBE+D2), and N G = 180 subscript 𝑁 𝐺 180 N_{G}=180 (SET1+2+3@PBE+D2). Experimental results at room temperature give Δ ads ​ H ​ ( T ) = − 4.06 subscript Δ ads 𝐻 𝑇 4.06 \Delta_{\text{ads}}H(T)=-4.06  kcal/mol Piccini et al. ( 2015 ); Luo et al. ( 2016 ) .

	MD	MD	SET1	SET2+3	SET1+2+3
	@PBE	@PBE+D2	@PBE+D2	@PBE+D2	@PBE+D2
PBE+D2	$-5.72\pm 0.48$	$-6.09\pm 0.18$	-	-	-
PBE	$-1.13\pm 0.18$	$-1.11\pm 0.25$	$-1.31\pm 0.58$	$-1.84\pm 0.53$	$-1.73\pm 0.43$
RPA	–	–	$-4.20\pm 0.62$	$-4.45\pm 0.48$	$-4.38\pm 0.39$

Equations

$$⟨O{⟩}_H=\underset{\tau \to \infty }{\lim }\frac{1}{\tau }{\int }_0^{\tau }O(\mathbf{q}(\mathbf{t}),\mathbf{p}(\mathbf{t}))dt,$$

$$V(\mathbf{q})=$$

$$H^{\prime }(\mathbf{q},\mathbf{p})=$$

$$⟨V{⟩}_H=$$

$$P_{H^{\prime }}^{V^{\prime }}(E)=$$

$${\Delta }_{\text{ads}}H(T)=⟨V_{S+A}{⟩}_H-\left[⟨V_S{⟩}_H+⟨V_A{⟩}_H\right]-k_{B}T,$$

$$\tilde{P}_{H}^{V}(E)=\frac{1}{M}\sum_{i=1}^{N_{G}}h_{i}\,\exp\bigg{\{}-\frac{\left(E-\epsilon_{i}\right)^{2}}{2w_{i}^{2}}\bigg{\}},$$

$$V({\mathbf{q}}_i)\in [{ϵ}_i-\Delta ,{ϵ}_i+\Delta ],$$

$${{ϵ}_i}^{\prime }=$$

$$\tilde{P}_{H^{\prime}}^{V^{\prime}}(E)=\frac{1}{M^{\prime}}\sum_{i=1}^{N_{G}}h^{\prime}_{i}\,\exp\bigg{\{}-\frac{(E-{{\epsilon}_{i}}^{\prime})^{2}}{2w_{i}^{2}}\bigg{\}}\,,$$

Full text

An alternative to DFT is to use correlated wave-function methods, which are more accurate and systematically improvable, but involve a higher computational cost. For applications to molecular systems the Møller-Plesset perturbation theory Møller and Plesset (1934) (MP2, MP3 etc.) and coupled-cluster theory Bartlett and Musial (2007) (CC) have long been used Szabo and Ostlund (2012). For condensed phase applications, MP2 has only recently appeared – using plane-wave Marsman et al. (2009); Dixit et al. (2017a), localized Pisani et al. (2008), and hybrid Del Ben et al. (2012) basis sets. Approximations based on coupled-cluster theory have been implemented in the VASP Booth et al. (2013) and PySCF codes McClain et al. (2017); Sun et al. (2018).

A particularly promising correlated approach for condensed phase applications is the random phase approximation (RPA) Bohm and Pines (1953); Gell-Mann and Brueckner (1957); Langreth and Perdew (1975); Dobson and Wang (1999); Furche (2001); Eshuis and Furche (2011); Heßelmann and Görling (2011a); Ren et al. (2012). This approach to the ground state electronic correlation energy can be seen both as a DFT approximation Langreth and Perdew (1975, 1977) or a CC approximation Scuseria et al. (2008). The RPA also has analogies with MP theory Mussard et al. (2016). Unlike CC and MP theories, which are applied exclusively as post Hartree-Fock methods, the RPA is most often evaluated from a DFT starting point.

The success of the RPA is largely related to its accurate description of weak van der Waals (vdW) interactions Dobson and Wang (1999); Furche (2001), but this approximation is sometimes less reliable for short range interactions. Although significantly more expensive than traditional DFT approximations, the RPA can be performed for systems with up to a few hundred electrons and a series of applications to condensed phase systems have been reported in the literature Del Ben et al. (2015a); Marini et al. (2006); Lebègue et al. (2010); Harl and Kresse (2008); Lu et al. (2009); Harl et al. (2010); Kaoui and Rocca (2016); Ren et al. (2009); Schimka et al. (2010); Göltl et al. (2012). When extra precision is required, beyond-RPA methods can be used Furche and van Voorhis (2005); Heßelmann and Görling (2011b); Bates and Furche (2013); Olsen and Thygesen (2012); Ren et al. (2013); Lu (2014); Colonna et al. (2014); Mussard et al. (2016); Dixit et al. (2017b); Hellgren et al. (2018).

Whether via RPA, MP or CC theories, correlated wave-function methods hold the promise of a significantly higher level of accuracy, which is of fundamental importance for design of truly predictive approaches to assist experimental and technological developments. However, they are characterized by an unfavorable scaling, slow basis set convergence Marini et al. (2006); Harl and Kresse (2008); Dixit et al. (2018), and slow cell size convergence Liao and Grüneis (2016), making them expensive.

Despite these difficulties, high-level methodologies are of vital importance in many technologically attractive cases involving e.g. nanoporous materials such as zeolites Speybroeck et al. (2015); Chibani et al. (2016); Bučko et al. (2017); Chibani et al. (2018); Grajciar et al. (2018) or metal-organic frameworks Chibani et al. (2018); Grajciar et al. (2018); Babaei and Wilmer (2016); Adil et al. (2017). These materials are intensively investigated for their ability to selectively adsorb target compounds and this property is used in numerous applications such as depollution Speybroeck et al. (2015); Chibani et al. (2016); Bučko et al. (2017); Chibani et al. (2018); Grajciar et al. (2018); Babaei and Wilmer (2016) or separation of chemicals Adil et al. (2017); Sholl and Lively (2016). Reliable predictions of adsorption properties of these materials can be made only if the electronic structure method used in simulations accurately describes all interactions within the system of interest, and if the thermal effects are properly accounted for by using a suitable statistical mechanics method such as molecular dynamics (MD) or Monte Carlo. Unfortunately, simulations combining high-level methodologies with MD or Monte Carlo approaches, usually requiring ensembles of tens of thousands of configurations, are thus far beyond the reach of the computational power available today. In fact, the forces for correlated methods in the condensed phase have only recently been implemented Del Ben et al. (2015b); Ramberger et al. (2017) and employed in MD simulations Del Ben et al. (2015a); Bokdam et al. (2017).

In this work we present a methodology based on perturbation theory that extends the applicability of correlated methods to compute finite temperature properties at a reasonable computational cost. Our approach starts from MD simulations based on computationally inexpensive DFT functionals. Then, a small number (a few tens) of significant configurations are chosen from the ensemble of structures generated by MD to represent the full probability distribution. Finally, correlated wave-function calculations are performed on these configurations and the results are used to reconstruct the corresponding probability distribution and ensemble averages. Below we will apply this method to the computation of the finite temperature adsorption enthalpies of CH4 and CO2 in the zeolite chabazite, using RPA.

The new approach introduced in this letter will be referred to as perturbation theory on selected configurations (selPT). A detailed discussion of perturbation theory in the context of molecular dynamics and free energy calculations can be found in Ref. Chipot and Pohorille, 2016. The main ideas and working equations will be summarized here.

Within the canonical ensemble, the expectation value of a certain observable $O$ is defined as follows: $\langle O\rangle_{H}=\int O(\mathbf{q},\mathbf{p})\text{exp}\left\{-\beta H(\mathbf{q},\mathbf{p})\right\}d\mathbf{q}\,d\mathbf{p}/Z$, for $Z=\int\text{exp}\left\{-\beta H(\mathbf{q},\mathbf{p})\right\}d\mathbf{q}\,d\mathbf{p}$. Here $\beta=1/k_{B}\,T$ and $\mathbf{q}$ and $\mathbf{p}$ denote nuclear positions and momenta, respectively. In this work the notation $\langle\ldots\rangle_{H}$ represents the canonical ensemble average corresponding to the classical Hamiltonian $H(\mathbf{q},\mathbf{p})=T(\mathbf{p})+V(\mathbf{q})$ with $T(\mathbf{p})$ and $V(\mathbf{q})$ being the nuclear kinetic and potential energies, respectively.

Exploring all $\mathbf{p}$ and $\mathbf{q}$ is generally impossible. An alternative is to assume the ergodic hypothesis. Then, the ensemble average can more usefully be expressed as

[TABLE]

a temporal average over an AIMD trajectory $\dot{\mathbf{q}}=\nabla_{\mathbf{p}}H$, $\dot{\mathbf{p}}=-\nabla_{\mathbf{q}}H$. For infinite time, $\langle O\rangle_{H}$ is independent of trajectory, since the system explores all configurations. However, for finite time, the choice of trajectory is crucial.

Ideally, we could use any method to obtain long AIMD trajectories that ensure a sufficient level of statistical precision. AIMD simulations are time-consuming, however, and thus long trajectories require a numerically efficient model electronic energy Hamiltonian $\hat{H}_{e}$, with potential

[TABLE]

at nuclear configuration $\mathbf{q}$ and repulsive nuclear energy $U_{II}$. Unfortunately, this leaves us with a compromise between accuracy (quality of electronic energies) and efficiency (length of trajectories), for a given precision.

In such cases, perturbation theory can be used to transform from an efficient/inaccurate Hamiltonian $H=T+V$ (e.g. PBE+D2) to an expensive/accurate, Hamiltonian $H^{\prime}$ (e.g. RPA), via the perturbation,

[TABLE]

The trick is to write thermal averages of potentials as,

[TABLE]

where $P_{H}^{V}(E)$ describes the energy $E$ distribution of some potential energy $V$ of interest. The subscript $H$ indicates that it is calculated in the trajectory described by Hamiltonian $H$. Then, one uses the perturbation $\Delta V(\mathbf{q})$ to obtain $\langle V^{\prime}\rangle_{H^{\prime}}=\int EdEP_{H^{\prime}}^{V^{\prime}}(E)$ from

[TABLE]

$P_{H^{\prime}}^{V^{\prime}}$ then gives properties of the trajectory $H^{\prime}$, but is calculated on the trajectory $H$.

This lets us define properties of the system driven by the primed Hamiltonian (e.g. RPA) in terms of trajectories (Eq. 1) on the system described by the original Hamiltonian (e.g. PBE+D2). Typical calculations employing perturbative approaches based on this formula involve large (albeit reduced) numbers of configurations computed in a given MD or Monte Carlo simulation, which are implausible for the expensive methods we seek to use. In this work, we thus introduce an innovation that lets the perturbation be evaluated with very few additional calculations, but without significant loss of precision. We thereby break the compromise between speed and accuracy that hampers the application of AIMD to difficult problems, such as the ones presented below.

To understand the challenges involved in the direct application of perturbation theory, we now present a specific numerical example that will also be used later to discuss our new methodological developments. Let us consider the finite temperature adsorption of methane on protonated chabazite (later, results will also be presented for the adsorption of carbon dioxide). The model used for this system, which involves 200 valence electrons per unit cell, is shown in Fig. 1 and further discussed in Supplemental Material. Using MD in the NVT ensemble, the adsorption enthalpy can be computed as

[TABLE]

where subscripts $S$ and $A$ denote the clean substrate (protonated chabazite), and the adsorbate (methane), respectively, while $S+A$ indicates the system of substrate interacting with adsorbate. In order to evaluate the ensemble averages in Eq. 6, MD simulations at 300 K have been performed with the VASP code Kresse and Hafner (1993) using two different approximations, namely PBE and PBE+D2 (a vdW corrected version of the PBE functional Grimme (2006)). Computational details are provided in Supplemental Material. These two MD simulations will be denoted as MD@PBE and MD@PBE+D2, respectively.

In Tab. 1 the PBE (PBE+D2) enthalpy of adsorption corresponding to MD@PBE (MD@PBE+D2) is an “exact” estimate, namely an estimate where the averaged potential and the Hamiltonian used for the sampling are at the same level of theory; the PBE (PBE+D2) value evaluated from MD@PBE+D2 (MD@PBE) is based on the perturbative approach in Eq. 5. Both PBE+D2 and PBE values evaluated perturbatively from MD@PBE and MD@PBE+D2, respectively, are close to their “exact” values, demonstrating the reliability of our approach. The PBE+D2 perturbative estimate is slightly worse, showing the importance of starting from an approach that correctly samples the key adsorption region of phase-space.

From a practical point of view, this brute force application of perturbation theory involves the reevaluation of the energy for thousands of configurations (specifically, the reevaluation of the electronic part of Eq. 2). If, starting from MD@PBE or MD@PBE+D2, we applied this perturbative procedure to compute the RPA enthalpy of adsorption, the corresponding computational time would amount to several millions of CPU hours in our supercomputer center. While technically feasible, such high demands are of limited practical interest.

We now proceed to show that it is possible to reduce the number of times the electronic part of Eq. 2 is recomputed to a very small fraction of the total number of configurations. This is achieved by finding a semi-analytic model of $P_{V}^{H}(E)$, and using it to build an approximate model of $P_{V^{\prime}}^{H^{\prime}}(E)$ by assuming that the distributions generated by $H$ and $H^{\prime}$ (e.g., $H_{\text{PBE+D2}}$ and $H_{\text{RPA}}$, as used here) are similar. We employ the following procedure to approximately reconstruct the correct high-level energy distribution $P_{H^{\prime}}(V^{\prime})$ and corresponding average values:

Step 1:  MD driven by the Hamiltonian $H(\mathbf{q},\mathbf{p})$ is performed to generate the $P_{H}^{V}(E)$ distribution (eq. 4).

Step 2:  $P_{H}^{V}(E)$ is approximated by a set of $N_{G}$ Gaussian functions with $3N_{G}$ fitting parameters $\left\{h_{i},w_{i},\epsilon_{i}\right\}$

[TABLE]

with normalization constant $M=\sum_{i=1}^{N_{G}}\sqrt{2\pi}h_{i}w_{i}$. Parameters are selected to represent key statistical features of the distribution. The tilde introduced in Eq. 7 emphasizes the approximate character of the probability distribution. It is crucial here to forbid any single or a limited number of Gaussians to dominate; for this reason the contribution of each single Gaussian is limited to a certain threshold. We note that this representation of potential energy density is related to the kernel density estimation method of probability density functions Parzen (1962). Technical details on the fitting with Gaussian functions are further discussed in Supplemental Material.

Step 3:  $N_{G}$ configurations with positions $\mathbf{q}_{i}$ are selected from the ensemble of step 1 to satisfy

[TABLE]

where we use $\Delta=0.02$ kcal/mol. Since several time-uncorrelated configurations can satisfy Eq. 8, the specific configuration $\mathbf{q}_{i}$ is randomly chosen with uniform probability. This is a reasonable procedure since configurations with the same energy ($\epsilon_{i}$) have equal probability to be sampled. For the selected configurations $\mathbf{q}_{i}$ a different (higher) level of theory is used to evaluate the energies $V^{\prime}(\mathbf{q}_{i})$ to obtain $\Delta V_{i}\equiv\Delta V(\mathbf{q}_{i})=V^{\prime}(\mathbf{q}_{i})-V(\mathbf{q}_{i})$.

Step 4:  Next, we calculate

[TABLE]

and use them to obtain:

[TABLE]

where $M^{\prime}=\sum_{i=1}^{N_{G}}\sqrt{2\pi}h^{\prime}_{i}w_{i}$. This equation follows from Eqs. 5 and 7, with reweighting factor $\exp\{-\beta\Delta V\}$.

Step 5:  Finally, we determine the estimate of $\langle{V^{\prime}}\rangle_{H^{\prime}}$ by combining eq. 10 with eq. 4. The accuracy of this ensemble average can be evaluated by computing the standard error according to Eq. S17 (see Supplemental Material). If a higher level of accuracy is required the original set of selected configurations can be extended to decrease the statistical error of the final result, as done later.

We label the five step procedure described above by the acronym “selPT”, which stands for perturbation theory on selected configurations. In order for this approach to be useful, statistical convergence should be achieved for $N_{G}\ll N_{\text{MD}}$, where $N_{\text{MD}}$ denotes the total number of MD steps. This is not an unreasonable expectation, as the configurations produced within a certain MD simulation are not completely independent due to the finite correlation time. Below, we show that even $N_{G}=60$ configurations can be sufficient for a reasonable estimate and $N_{G}=180$ can provide highly accurate results – both are vast improvements on the $N_{\text{MD}}\approx 190,000$ runs required to recompute a full MD trajectory. Such a dramatic reduction is possible because $\tilde{P}_{H^{\prime}}^{V^{\prime}}(E)$ is a smooth, representative function that captures the important properties of the true distribution $P_{H^{\prime}}^{V^{\prime}}(E)$.

To establish the accuracy and the efficiency of the selPT procedure we come back to the example of CH4 adsorbed on chabazite. Before computing RPA adsorption energies, the approach is validated using PBE and PBE+D2: First, the “exact” distribution $P_{H_{\text{PBE+D2}}}^{V_{\text{PBE+D2}}}$ generated by the previously described MD@PBE+D2 simulation is fitted by 60 Gaussian-shaped functions (eq. 7); second, 60 configurations $\mathbf{q}_{i}$ are chosen in order to satisfy the condition in Eq. 8; finally, Eq. 10 is used to evaluate $\tilde{P}_{H_{\text{PBE}}}^{V_{\text{PBE}}}$ and the corresponding ensemble averages. The procedure is repeated three times to generate three uncorrelated sets containing 60 configurations each, that will be called SET1@PBE+D2, SET2@PBE+D2, and SET3@PBE+D2. This is helpful to study the behavior of the selPT method as a function of the number of configurations. Indeed, the adsorption enthalpies have been computed using the three sets with 60 configurations, from the three sets with 120 configurations obtained by merging pairs of the original sets (SET1+2@PBE+D2, SET2+3@PBE+D2, and SET1+3@PBE+D2), and with a set that includes all the configurations (SET1+2+3@PBE+D2). Since the enthalpy of adsorption Eq. 6 requires ensemble averages for the three different systems (S, A, S+A), the selPT procedure is applied independently on each one of them.

Fig. 1 demonstrates the fitting quality. The middle panel demonstrates the ability of selPT to reproduce key features of the sampling data, and to perturb them using a more accurate method, by comparing $P_{H}^{V}$ from the full MD@PBE+D2 run, with its parametrized distribution $\tilde{P}_{H}^{V}$ and RPA counterpart $\tilde{P}_{H^{\prime}}^{V^{\prime}}$ approximated from Eq. 10 using SET2@PBE+D2. The right panel shows a phase space diagram for the full MD run, and the 60 samples selected by selPT, reported in terms of variables representing the speed of CH4 and its distance from the surface of the chabazite. Note that selPT samples the probability distribution accurately, by design; as shown in the right panel of Fig. 1 a representative part of phase space is also covered for the system under consideration.

The row corresponding to PBE in Table 1 shows the practical reliability of selPT in computing enthalpies of adsorption $\Delta_{\text{ads}}H(T)$. Specifically, the similarity of the energies across all calculations highlights the success of the approach. If our trajectories sampled an overly-limited set of configurations (whether through a too-short trajectory or too-inaccurate functional), or if our selections sets were too poor, the numbers would differ significantly across columns. The consistency of the values justifies the length and reasonableness of our initial AIMD calculations, and the accuracy of the selection process.

The bottom row of Table 1 shows RPA enthalpies of adsorption, which would be basically impossible to obtain using a full AIMD treatment due to the inherent cost of RPA (full results considering all the possible combinations of sets are provided in Tables S9 and S10 of Supplemental Material). Numerical results are based on the RPA implementation in the VASP code Harl and Kresse (2008) and computational parameters are provided in the Supplemental Material. The statistically most accurate selPT estimate of the RPA enthalpy of adsorption is obtained from the largest set SET1+2+3@PBE+D2 ($N_{G}=180$) and provides a value of $-$4.38 $\pm$ 0.39 kcal/mol in excellent agreement with the experimental value of $-4.06$ kcal/mol Piccini et al. (2015); Luo et al. (2016). From Table 1 (and Tables S9-S10) it can be noticed that the errors present a clear trend to decrease as a function of the number of configurations included in the selPT set. In Table S9 it is shown that RPA results from as few as 60 samples are all within 0.92 kcal/mol of one another, and within 0.83 kcal/mol of experiment. This is already a considerable improvement with respect to the PBE and PBE+D2 estimates, which present deviations of at least 2 kcal/mol well beyond the chemical accuracy.

Key for this high level of accuracy is that the selPT calculations are based on a reasonable PBE+D2 trajectory which, unlike PBE, samples phase space similarly to RPA. This ensures that the perturbation from PBE+D2 to RPA is correctly biased to the weak vdW interaction between the zeolite and molecule. This intuitive point is further discussed in a quantitative way in Supplemental Material where RPA results from a PBE starting point are presented (Secs. SIII and SIV) and it is shown how PBE leads to significantly different sampling of the configurational space with respect to PBE+D2 and RPA (Sec. SV).

To further establish the accuracy and reliability of the selPT approach, a second system is discussed, CO2 adsorbed in protonated chabazite. Computational details are analogous to those of CH4 in chabazite. The adsorption enthalpy of CO2 is sizeably larger than that of CH4, with an experimental value of -8.41 kcal/mol at room temperature Pham et al. (2012). While the PBE+D2 result deviates from this value by 1.3 kcal/mol, our most accurate estimate from the largest set SET1+2+3@PBE+D2 provides an highly accurate enthalpy of adsorption of $-$8.74 $\pm$ 0.39 kcal/mol. Detailed results for all the test sets are presented in Table S9 ($N_{G}=60$) and S10 ($N_{G}=120$) in Supplemental Material. The errors are analogous to those for adsorbed CH4; by considering that the adsorption enthalpy of CO2 is about double than that of CH4, the relative error in this case decreases much more rapidly as a function of the set size and $N_{G}=60$ provides already an accuracy within $\pm 10\%$.

From a practical point of view, the selPT approach has been instrumental to achieve the computation of RPA enthalpies of adsorption. By decreasing the number of configurations to 60, the RPA results are obtained in a CPU time that is only about 3-4 times the CPU time involved in the full MD@PBE calculation, or around ten times for the more accurate simulation with 180 configurations.

In conclusion, we introduced the selPT method to compute ensemble averages of energies and, more generally, energy distributions for correlated wave-function methods. As proof of principle, this methodology was applied to compute the finite-temperature adsorption energy of CH4 and CO2 in protonated chabazite at the RPA level of theory, where it successfully reproduced known experimental values within chemical accuracy. selPT is not limited to the RPA but could be interfaced with other sophisticated approximations (e.g. MP2, double-hybrids or coupled-cluster theory), as long as computations on few tens of configurations can be afforded.

Beyond providing an accurate and efficient methodology to address challenging problems in physics and material science, selPT could be also used to develop finite-temperature test sets from high level theories, to test and develop new DFT functionals and semi-classical models.