Applicability of Tail Corrections in the Molecular Simulations of Porous Materials.

Molecular simulations with periodic boundary conditions require the definition of a certain cutoff radius, rc, beyond which pairwise dispersion interactions are neglected. For the simulation of homogeneous phases the use of tail corrections is well-established, which can remedy this truncation of the potential. These corrections are built under the assumption that beyond rc the radial distribution function, g(r), is equal to one. In this work we shed some light on the discussion of whether tail corrections should be used in the modeling of heterogeneous systems. We show that for the adsorption of gases in a diverse set of nanoporous crystalline materials (zeolites, covalent organic frameworks, and metal-organic frameworks), tail corrections are a convenient choice to make the adsorption results less sensitive to the details of the truncation.

Introduction

Microporous materials are characterized by nanosized pores and large internal surface areas. These features have been exploited to separate[1] and store gas molecules[2] because of the attractive interactions between the molecules and the atoms of the frameworks. Further gas-related applications include catalysis[3,4] and sensing.[5]

Metal–organic frameworks (MOFs), covalent organic frameworks (COFs),[6] and zeolites[7] are three classes of microporous materials that also have the quality of being crystalline, i.e., their unit cells repeat periodically, and the bulk materials can be typically represented by just a few hundreds of atoms. Tens of thousands of these microporous crystals have been reported in the literature,[8−10] and many more have been generated in silico as potentially synthesizable structures.[11−13] Finding the best material for any application related to the adsorption of some specific molecules, or understanding their property–performance relations, appears to be an increasingly vast combinatorial problem.[14,15] Molecular simulations are a valuable tool to screen these materials, but they rely on accurate modeling of the gas–framework interaction and the methods of statistical thermodynamics to predict the thermodynamics properties that are used to identify promising materials among numerous candidates.[16]

It is common in the modeling of gas adsorption in crystalline microporous materials to keep the geometry of the framework fixed[17,18] and to model the adsorbate–gas and gas–gas interactions by dispersion interactions (e.g., Lennard-Jones or Buckingham potential) combined with Coulombic interactions between partial charges.[19] When a reliable protocol is found for the modeling, it is important to provide a detailed description of this protocol,[20] ideally in a findable, accessible, interoperable, and reusable (FAIR) way,[21,22] including the parameters used for the dispersion interactions and their mixing rules and the charge-derivation method; this is referred to as the “force field” used in the model. Other important details should be included to fully describe the force field. These are the cutoff distance for the dispersion interactions, which is used to avoid computing the interaction between atoms that are so far apart that their interactions can be ignored; the truncation method (the choice of setting to zero the potential beyond the cutoff, shifting the entire potential to have a zero value at the cutoff distance, or the usage of a switching potential); and finally, whether tail corrections are applied or not. On this last point we focus our attention, to highlight the influence of tail corrections on the results and convergence behavior (with respect to the cutoff radius) in typical gas adsorption calculations. In this way, we address a question that is still debated in the literature: should tail correction be used for simulations of microporous crystals?[23−25]

In periodic simulations, one would like to avoid a particle interacting twice with the same portion of the space, and therefore, one typically imposes that the size of the simulation box is bigger than two times the cutoff radius.[26] Consequently, in the case of microporous materials, one has to replicate the number of unit cells until this criterion is met, with a significant increase of the computational time for each replication due to the larger number of particles and therefore pairwise distances to compute. Thus, it is advisable to choose an adequate cutoff radius, which is a compromise between the desire to avoid large simulation boxes while preventing the neglect of significant interactions. For a long time, a cutoff radius not greater than 6–8 Å (i.e., around 2.5 times the typical Lennard-Jones σ) was normally employed,[27] whereas in most of the recent simulations common cutoffs are in the range of 12–18 Å,[28] thanks to faster computers.

For homogeneous systems, like liquids, it is usually assumed that one can correct for contributions above the cutoff using tail corrections, which are based on the assumption that the radial distribution function g(r) ≈ 1 for r ≥ rc.[26] This means that one assumes that beyond the cutoff, particles are randomly displaced, and one can analytically compute their overall dispersion interactions with the bead. For example, the correction that applies to the popular Lennard-Jones 12–6 potential energy is[29]and similar expressions exist for the pressure and chemical potential. Even though there are special long-range correction schemes for surfaces[30−32] and other geometries,[33] tail corrections of the general form of eq are still the most widely used ones.

It has been shown that these tail corrections can have a pronounced effect on the computed thermodynamics properties, such as the critical temperature of homogeneous particles.[27,34] For heterogeneous systems, such as microporous crystals, however, the radial distribution function does not necessarily converge to one at rc, and therefore, the validity of the tail-correction approximation is matter of debate: Siperstein found that tail corrections might be applicable,[23] whereas other authors[24,25] state that they should not be used in microporous frameworks. This controversy motivated this work. To address the question of whether the use of tail corrections for adsorption in porous media is advisable, we consider the simple case of methane adsorption in microporous materials, where methane is modeled as a single bead particle interacting with a 12–6 Lennard-Jones potential, where Coulomb interactions can be neglected because they are small for methane adsorption. Here, we systematically study for a large set of structures (including MOFs, COFs, and zeolites) the effect of these tail corrections on commonly computed properties (e.g., the Henry coefficient, gas uptake at different pressures, and deliverable capacity).

Methods

Structure Selection

The simulations shown in this study were performed on a selection of 40 microporous crystals: 10 COFs from the revised CoRE-COF database, which Ongari et al. released as the CURATED-COFs database;[10,35] 10 zeolites from the IZA database;[8] 10 MOFs from the CoRE-MOF database;[9] and 10 other common MOFs, chosen because of their popularity in computational studies.

We characterized all materials in these databases considering the following properties as descriptors of our pore–property space: accessible surface area, density, largest free sphere, number of channels, number of pockets, pocket surface area, and largest included sphere. These properties were calculated using the Zeo++ code,[36] using a probe diameter of 1.86 Å and atomic diameters corresponding to the Lennard-Jones σ parameter of the force field. To rationally choose a diverse set of porous structures for our investigation, the first 30 frameworks were selected from their databases based on k-nearest neighbor clustering (with k = 10) in pore-property space, as implemented in our structure-comp python package.[37] After knn clustering in property space, the samples closest to the k centroids are selected for the study, which ensures that we do not miss any cluster that may exist in the hyperspace of geometrical properties.

Monte Carlo Simulations

We used the AiiDA package[38] to design a reproducible workflow for the Monte Carlo (MC) simulations, which were performed using the RASPA[39] code. In our workflow, we performed all the simulations at 298 K, considering a single-bead methane particle (with methane–methane interactions described by the TraPPE united-atom model).[40]

Simulations, with and without the tail corrections, were run at different cutoffs, in the range of 6–30 Å. All relative errors shown in the following are computed relative to the values at a cutoff radius of 24 Å. The absolute values of all properties are shown in the Supporting Information and indicate that most simulations are well-converged for cutoff radii above 20 Å. We performed 300 000 Widom insertions[41] to calculate the Henry coefficients KH at infinite dilution. Grand-canonical Monte Carlo (GCMC) simulations were used to evaluate the uptake at 5.8, 35, and 65 bar. We used the Peng–Robinson equation of state to convert the pressure to the corresponding fugacity.[42] We used 15 000 cycles for equilibration and 15 000 for production. We found the number of MC cycles sufficient to give small error bars compared to the error due to the truncation of the potential (compare the absolute values in the Supporting Information). NVT ensemble simulations with a single methane bead were employed to calculate the radial distribution function inside the frameworks, considering rc = 20 Å.

To model the interaction with the frameworks, we used the DREIDING force field,[43] where, for the missing elements (Cd, Cu, Mg, Sr, Ti, and Zr), we added the Lennard-Jones parameters from the UFF force field.[44] We applied Lorentz–Berthelot mixing rules to compute the interaction parameters. For all simulations shown in the main text we blocked inaccessible pores.[45,46] Given that some structures from our selection have nonpermeable pores, as determined using the Zeo++ code,[47,48] we also performed simulations without pore blocking and discuss them in the Supporting Information.

Results and Discussion

Error Bounds

In the following, we discuss the error in the potential energy that comes from truncating the pairwise potential. This error is directly related to the radial distribution function for the adsorbed particles in the microporous material. For a pairwise potential of the formthe average potential energy in the canonical ensemble, ⟨U⟩, can be written as a double integral over the product of the pair–correlation function, g(2), given as the ensemble average over pair–correlations of positionsand the pairwise energy term u(|r – r|):[49]Here, ρ = N/V is the number density of particles, and we changed the variables to center-of-mass R = (r + r)/2 and relative coordinates r = r – r. Note that the more popular radial distribution function, g(r), is obtained by integrating the pair-correlation function, g̃(2)(r, R), over the center-of-mass positions R and the spherical angles ϕ and θ

The error in the potential energy due to the truncation, εtrunc., that one introduces by neglecting the interactions beyond a cutoff, i.e., for |r| > rc, isWe can rearrange the equation using the mean-value theorem. To do so, we assume that both the potential u(|r|) and the pair–correlation function g̃(2)(r, R) are continuous over the range of integration, and we introduce a new upper boundary of integration, p, to approximate the indefinite integral with a definite integral. This new upper boundary of integration, p, is large enough such that u(|p|) = 0 (i.e., |p| → ∞). Finally, we define a value ξ ∈ [rc, p] for the relative distances, and we rewrite the error as

The tail corrections are based on the assumption that g̃(2)(r, R) = 1 for r ≥ rc; hence, the error εt.c. on the potential energy can be written as

Comparing eq with eq , we see that the error using tail correction is in general smaller than for the simple truncation if ∫dRg̃(2)(ξ, R) > 0.5, ∀ ξ. In terms of the radial distribution function, this corresponds to g(r) > 0.5 ∀ r > rc; if this is the case, one can be confident that the use of homogeneous tail corrections, instead of simple truncation, is the more accurate choice. We will show in the following sections that this is the case for all the microporous materials that we selected for this study and that have permeable pores, as well as for the simulations of Macedonia and Maginn,[25] who were the first to question the applicability of tail corrections in porous systems.

Proof of Concept on a Simplified Model

As a numeric example of the behavior of tail corrections in crystalline materials, we first start recalling the simple case of a Lennard-Jones fluid with the well-known pair potentialfor which we can use a parametrization of the radial distribution function from Matteolia et al.[34]Figure illustrates that the tail correction gives us an estimate of the energy that is less sensitive to the cutoff radius, i.e., we have εtrunc. > εt.c. for the physical sensible region of rc > σ.

Figure 1

Error for truncated potential (trunc. error) and truncation with tail correction (t. c. error) for a Lennard-Jones fluid.

We can now perturb the radial distribution function from Matteolia et al. to imitate the behavior in a crystal, i.e. by using a smaller exponential decay and adding additional oscillations (see the Supporting Information) as shown in Figure . Note that similar approximations have already been used in the literature to describe porous materials, e.g., aluminum oxide.[50,51]

Figure 2

Error for truncated potential (trunc. error) and truncation with tail correction (t. c. error) for a perturbed Lennard-Jones fluid.

In ref (25), Macedonia and Maginn make the case that for the tail correction to be exact g(r) must converge to one for sufficiently large r. Here, we take a slightly different point of view in the sense that we show that the errors related to different values of the cutoff radius are smaller if we use a tail correction rather than simple truncation for which different groups use different values of the cutoff radius. For periodic systems, where g(r1, r2) = g(r1 + R, r2 + R) ∀ R, and where the R form a Bravais lattice,[52] one might intuitively expect some error compensation due to the oscillatory nature of g(r), and still, it is not clear a priori to what extent this is the case in porous materials. However, we could show in this perturbed model that, even when the periodicity is broken, one can still find that the simulations with tail corrections are less sensitive to the cutoff radius of the truncation (see Figure ).

Framework-CH4 Radial Distribution Functions

In Figure we show the framework-CH4 radial distribution functions in the limit of infinite dilution (i.e., one particle per unit cell) for all the microporous materials in our set that have pores that are accessible to our methane bead.

Figure 3

CH4-framework radial distribution functions. Dashed lines indicate g(r) = 1; dotted lines indicate g(r) = 0.5.

We observe that, within 20 Å, most of the radial distribution functions do not converge to g(r) ≈ 1. Additionally, we also observe for all the structures with accessible pores that g(r) > 0.5 (we discuss the cases with nonaccessible cases in the Supporting Information), which is the threshold we derived by comparing eq with eq and indicates that the use of tail corrections should be a more appropriate choice, if compared to simple truncation.

Convergence of Henry Coefficient at Different Pressures

For the Henry coefficient, KH, we can writeThese expressions show that the Henry coefficients depend exponentially on the tail correction,[23] and therefore, they are a useful probe to study the influence of the tail correction on simulation results.

For all structures, we find that KH – KH,converged is smaller with tail corrections than without (cf. Figure ), especially at cutoff radii around 12 Å, which are commonly used in current molecular simulations. We also find that the converged results of simulations with and without tail corrections are identical within the error of the simulations (cf. details in the Supporting Information). In the Supporting Information we also demonstrate that the more desirable convergence behavior with respect to the cutoff can also be found for the more complex cases of adsorption of water and longer alkanes.

Figure 4

Convergence of the Henry coefficient KH as a function of the cutoff distance with and without tail corrections. We show the relative error due to the large spread of Henry coefficients and give the absolute values in the Supporting Information.

Convergence of Methane Loading at Different Pressures

To further illustrate the applicability of tail corrections, we compared real-case adsorption simulations of methane to assess the effect of tail corrections at different cutoff radii. The actual application property for gas storage, the deliverable capacity DC, is defined as the difference between the amount of methane stored at storage pressure and the methane that remains at depletion pressure, DC65 bar,5.8 bar = θ65 bar – θ5.8 bar, where we have chosen the depletion pressure set by the Advanced Research Projects Agency-Energy (ARPA-E) of the US Department of Energy that is part of the design target to make adsorbed natural gas storage competitive with compressed natural gas storage.[53,54]

Also for this property we find that simulations with tail corrections generally show a more desirable convergence behavior with respect to the cutoff than simulations without tail corrections (cf. Figure ).

Figure 5

Convergence of deliverable capacities DC65 bar,5.8 bar as a function of the cutoff distance with and without tail corrections. ΔDC = DC(rc) – DC(rc = 24 Å).

Notably, in the case of deliverable capacity, the error coming from the truncation (without tail corrections) can lead to both an overestimation or underestimation, while, when considering the total uptake at a certain pressure, the truncation leads to weaker interactions and always lower values (see the Supporting Information).

Because it is common for high-throughput screenings for methane deliverable capacity to employ cutoffs around 12 Å, the findings from Figure indicate the use of tail correction might be a preferable choice for these studies. Indeed, the deliverable capacity is, in all the cases, less sensitive to the choice of the cutoff. In contrast, simulations with truncated potentials lead to nonsystematic deviations of around 10% from the cutoff-converged results.

Conclusions

In this work we investigated, for microporous frameworks, the influence of tail corrections on the convergence behavior of gas adsorption properties with respect to the truncation radius of the pairwise dispersion potential. We found that for most of the systems, the strict condition g(r) > 0.5 ∀ r > rc holds, and homogeneous tail corrections give results that are far less sensitive to the details of the truncation of the potential (e.g., the cutoff radius) and therefore are the preferable choice. These findings are proven to be valid for a test set of diverse structures, including zeolites, metal–organic frameworks, and covalent organic frameworks.

Given that there is no commonly agreed cutoff for the potential, we recommend the use of tail corrections when modeling gas adsorption in microporous materials in order to allow for a more consistent direct comparison of the results from different simulation studies.