The effect of atomic point charges on adsorption isotherms of CO₂ and water in metal organic frameworks

Abstract

The interactions between metal–organic frameworks (MOFs) and adsorbates have been increasingly predicted and studied by computer simulations, particularly by Grand-Canonical Monte Carlo (GCMC), as this method enables comparing the results with experimental data and also provides a degree of molecular level detail that is difficult to obtain in experiments. The assignment of atomic point charges to each atom of the framework is essential for modelling Coulombic interactions between the MOF and the adsorbate. Such interactions are important in adsorption of polar gases like water or carbon dioxide, both of which are central in carbon capture processes. The aim of this work is to systematically investigate the effect of varying atomic point charges on adsorption isotherm predictions, identify the underlying trends, and based on this knowledge to improve existing models in order to increase the accuracy of gas adsorption prediction in MOFs. Adsorption isotherms for CO₂ and water in several MOFs were generated with GCMC, using the same computational parameters for each material except framework point charge sets that were obtained through a wide range of computational approaches. We carried out this work for 6 widely studied MOFs; IRMOF-1, MIL-47, UiO-66, CuBTC, Co-MOF-74 and SIFSIX-2-Cu-I. We included both MOFs with and without open metal sites (OMS), specifically to investigate whether this property affects the predicted adsorption behaviour. Our results show that point charges obtained from quantum mechanical calculations on fully periodic structures are generally more consistent and reliable than those obtained from either cluster-based QM calculations or semi-empirical approaches. Furthermore, adsorption in MOFs that contain OMS is much more sensitive to the point charge values, with particularly large variability being observed for water adsorption in such MOFs. This suggests that particular care must be taken when simulating adsorption of polar molecules in MOFs with open metal sites to ensure that accurate results are obtained.

1 Introduction

Metal–organic frameworks (MOFs) are porous crystalline materials consisting of coordination bonds between transition-metal cations and organic ligands. A wide array of metals and organic linkers have been used to produce MOFs, giving rise to structures with different properties, which can be tailored to specific requirements or applications (Moghadam et al. 2017).The intrinsic properties of these materials include very high porosity (Furukawa et al. 2010) and surface areas (Farha et al. 2012), which can be employed in various applications, such as gas storage and separation (Ma and Zhou 2010), sensor devices (Li et al. 2012) or drug delivery (Farrusseng 2011). In particular, MOFs have been widely investigated as promising adsorbents for carbon dioxide capture, normally by designing new materials with high affinity for CO₂ (Kenarsari et al. 2013; Torrisi et al. 2010). In this context, however, a crucial and often overlooked aspect is the effect of the presence of trace amounts of water vapour in the feed stream on the performance of the process. Water, being a highly polar molecule, can compete with CO₂ for adsorption in strongly adsorbing sites of the MOF, and therefore affects the efficiency of carbon capture (Yazaydin et al. 2009a). Indeed, some MOFs contain open metal sites (OMS), also known as coordinatively unsaturated sites (CUS), which are unsaturated metal centers developed upon solvent removal during the activation step. These sites can form strong coordination bonds with adsorbates such as CO₂, water or unsaturated hydrocarbons by electron donation from π orbitals of the adsorbate to the metal (Zhang et al. 2015). This phenomenon affords greater control over selectivity in gas adsorption, for example increased affinity for gases such as CO₂ or ethene, making these materials excellent candidates for gas separation applications (Campbell et al. 2018). However, water also adsorbs very strongly at the CUS, and the exact mechanism of competitive adsorption at these sites is still not fully understood.

Despite their great potential as adsorbents, the vast number and variety of MOF structures make it extremely time consuming and error prone to experimentally evaluate their properties. This obstacle is being increasingly overcome by using molecular simulation. This computational method allows high throughput screening of existing or hypothetical MOFs by modelling the MOF structure and its interactions with the adsorbate(s). It is therefore possible to efficiently identify materials with properties that are desirable for a specific application, for example their selectivity for a specific gas mixture, by simulating their adsorption behaviour. Once the best performing MOFs have been narrowed down, they can be investigated further by experimental methods. Watanabe and Sholl (2012) used this approach to screen over 30,000 MOFs for CO₂/N₂ separation, and consequently identified 1163 promising materials to be further examined by computer modelling. Importantly, molecular simulation also allows testing of theoretical MOFs that have not been synthesized yet, thus making it possible to identify new hypothetical materials with desirable properties. For example, Wilmer et al. (2012b) used computer modelling to screen 130,000 hypothetical MOFs for their CO₂ adsorption properties, while Haldoupis et al. (2012) screened ~ 500 MOFs for CO₂/N₂ selectivity, shortlisting 11 promising materials to be investigated further. Grand Canonical Monte Carlo (GCMC) simulations are normally used to predict adsorption isotherms in MOFs. The accuracy of the predictions relies on determining the potential energy of the system, and thus it is crucial that the interactions between adsorbates and the MOF are described correctly. This is typically done by using forcefields, a set of parameters describing the interactions between all atoms in the system. In the standard approach, repulsion and dispersion interactions are described by the Lennard-Jones (LJ) potential, while permanent electrostatic interactions are described by a set of atomic point charges. In this work, we assess the effect of framework point charge selection on the resulting adsorption isotherm predictions.

The electron density around a molecule consists of electrons that are constantly in motion, but they are likely to be found in some locations more than others. Electrostatic interactions thus arise from differences in electronegativity between atoms, and they are governed by Coulomb’s law. For simulation purposes, the electrostatic potential is normally represented by assigning point charges to each atom to simplify the computation process. These are therefore important parameters in adsorption simulations, particularly with polar species such as water. CO₂ is a quadrupolar molecule—due to the difference in electronegativities, there is a partial positive charge on the carbon and partial negative charges on the oxygens. Due to this effect, the molecule will orient itself preferentially based on the interactions of these charges; therefore electrostatic interactions also have to be accounted for in simulations of CO₂ adsorption. To clarify, although strictly speaking most classical intermolecular interactions (including repulsion and dispersion) are electrostatic in origin, we use the term “electrostatic” in the context of this paper to describe interactions arising from permanent molecular multipole moments, described through the use of effective atomic point charges.

There is currently no universally accepted system of point charge assignment, because point charges are not experimentally observable properties. As a consequence, multiple methods have been developed to represent these interactions as accurately as possible. In most cases, a quantum mechanical (QM) calculation using methods such as Density Functional Theory (DFT) (Parr and Yang 1994) or Hartree–Fock (HF) (Szabo and Ostlund 2006) is first carried out, followed by applying mathematical analysis methods that assign point charges to each atom in order to best represent the full electron distribution of the system. The approach of assigning point charges is convenient for simulation purposes; however, it has several caveats. For example, it is not possible to exactly reproduce the continuous nature of electron density by a set of fixed point charges, and this can sometimes lead to pronounced artefacts. The charge values may also depend strongly on the QM method (e.g. basis set size), on the configuration of the molecule and on the charge determination approach. Indeed, there are various approaches to divide the electron density of the material, and they can generate widely different charges (Watanabe et al. 2011). They can be broadly classified as atomic orbital assignment using basis sets [e.g., Mulliken (1955)], direct partitioning of electron density [e.g., DDEC (Manz and Sholl 2010), or Hirshfeld (1977)], or fitting the electrostatic potential around the molecule [e.g., CHELPG (Breneman and Wiberg 1990), or REPEAT (Campana et al. 2009)]. Furthermore, semiempirical methods that do not require computationally expensive QM calculations have also been developed, such as charges based on bond connectivity [CBAC (Xu and Zhong 2010)] or electrostatic equalisation methods [e.g., EQeq (Wilmer et al. 2012a)].

The importance of inclusion and correct estimation of electrostatic interactions between the MOF and the adsorbate has been studied for various systems. Yang and Zhong (2006b) showed the crucial effect of electrostatic interactions in CO₂/CH₄ adsorption simulations, which showed reversed selectivity depending on whether electrostatics were included or not. The system showed selectivity for methane when the electrostatic interactions were not considered, however, higher selectivity for CO₂ was observed once these interactions were accounted for. This was particularly evident at higher pressures, when the long-range electrostatic interactions exerted their effect over the longer distances of multilayer adsorption. McDaniel et al. (2015) screened 424 MOFs for CO₂ and CH₄ uptake, followed by comparing REPEAT, Qeq (Kadantsev et al. 2013) and no charges to in-house Q_sbu charges for selected materials, finding slight underestimation by the QEq method and significantly lower uptake when no charges were used. Zheng et al. (2009) studied the influence of point charges on CO₂ adsorption on 20 MOFs with different structural properties by simulating the adsorption with and without electrostatics and calculating the contribution of framework point charges to CO₂ uptake. The results showed that the effect of framework charges was most significant at low pressures when the strongest CO₂ adsorption was occurring around the metal centres, reaching up to 40% of the total uptake for CuBTC at low pressure. These observations clearly demonstrate the importance of including framework point charges in the adsorption process and also the potential of selectivity enhancement in mixtures with different electrostatic interactions (Yang and Zhong 2006a).

Walton et al. (2008) carried out a study on IRMOF-1, simulating CO₂ adsorption isotherms and comparing them to experimental data. They found that good agreement with experiment was obtained when adsorbate–adsorbate electrostatics were included in the model but, crucially, when adsorbate-MOF electrostatic interactions were omitted. Inclusion of point charges on the MOF was in fact deemed unnecessary because it led to a deterioration in agreement with experiment. This highlights the potential pitfalls, such as error cancellation, when comparing simulations against experimental adsorption isotherms. Watanabe et al. (2011) compared CO₂ adsorption isotherms using framework charges calculated using DDEC, REPEAT, CBAC and Hirshfeld methods, for a variety of Zn-based MOFs including IRMOF-1 and ZIF-90. They concluded that for some MOFs the adsorption behaviour is consistent even when different point charges are used, while for other materials a slight change in the point charges can cause significant deviations in adsorption predictions. Indeed, their results showed that while IRMOF-1 produced fairly similar adsorption isotherms, they were very different in ZIF-90, thus showing that the effect of framework point charges on the adsorption behaviour is dependent on the material and its structure. While Haldoupis et al. (2015) compared DDEC and LoProp charges in M-MOF-74 (M = Mn, Co, Ni, Cu) resulting in significantly different adsorption isotherms, Borycz et al. (2016) compared DDEC and CM5 charges in M-IRMOF-10, finding that although DDEC charges were significantly larger in magnitude than CM5 charges, they both led to similar adsorption isotherms. Babarao et al. (2011) compared CO₂ adsorption on ZIF-68 and ZIF-69 using CHELPG and Mulliken point charges. The results showed that while Mulliken charges led to slight overestimation of adsorption, the difference was not significant. However, as only two structurally similar materials and two sets of charges were employed, this result does not offer a wide overview of the effect of different framework point charges on adsorption behaviour.

Prolonged exposure to water vapour can result in degradation of the MOF structure. Understanding water adsorption in CuBTC has been a topic of many studies, both experimental and theoretical. It has been shown that the open metal sites play a crucial role in water adsorption. Previous work by Yang et al. (2007) has demonstrated that in CO₂/N₂ mixture adsorption simulations, the Cu metal centres will preferentially adsorb CO₂ due to its larger quadrupole moment. However, their affinity for water is expected to be even higher, as water has a large dipole moment. Yazaydin et al. (2009a) have demonstrated this in their GCMC simulations, along with high sensitivity of the point charge on the Cu atom for modelling water adsorption. All these findings lead to the conclusion that the open metal sites play a crucial role in adsorption in MOFs.

A set of DDEC point charges for over 2000 MOFs was developed by Nazarian et al. (2016), who also compared these charges to those obtained by the CBAC and EQeq methods. Their findings showed that higher charges on the metal centre were produced consistently by the CBAC method. Comparison to the semi-empirical EQeq method showed large differences to the charges calculated by the DDEC method, leading to possible inaccuracies in adsorption modelling. Hamad et al. found that simulations with different charges led to difference in thermal expansion results, thus showing that the choice of framework point charges can also affect modelling of structural properties of MOFs (Hamad et al. 2015).

Although the importance of including and accurately defining electrostatic interactions has been demonstrated, the precise effect of different framework point charges on adsorption isotherm predictions has not yet been systematically studied over a wide enough range of MOF topologies. Instead, the comparison of point charges has mostly been carried out for the same group of materials, such as Zn-based MOFs. In this work, we performed a detailed and systematic investigation of the effect of framework point charges on GCMC adsorption predictions in MOFs. We simulated adsorption of both CO₂ and water, as it has been demonstrated that for these adsorbates the electrostatic interactions are of crucial importance in the adsorption process, and due to their relevance for carbon capture applications. To gain a wider perspective and a deeper understanding, we have included MOF materials with different composition and structural properties, representing the most important MOF families. We intentionally included MOFs with and without OMS, to see whether this structural property would have an impact on the results. Methods for point charge assignment investigated in this work covered the most important classes and were based on QM calculations at different levels of theory and with different structural model approaches (e.g., periodic vs. cluster). In order to isolate the effect of framework point charges and decouple it from other sources of error and variation in simulated adsorption isotherms, we have: (i) selected well-established and widely used adsorbate–adsorbate potentials for CO₂ and water, and kept them identical in all simulations; (ii) selected a widely used model (DREIDING) for the framework LJ contribution, together with Lorentz-Berthelot combining rules, and kept the adsorbate-adsorbent LJ parameters constant in all simulations; (iii) run our simulations with always the same set of technical parameters (run lengths, type and probability of MC moves, etc.) so that they are directly comparable to each other. Because our goal is not to judge the quality of the overall force field, or to fit/adjust force field parameters to match any set of reference data, we have intentionally avoided a direct comparison with experimental data. This way, we also eliminate any uncertainty arising from variability in experimental measurements of adsorption isotherms in MOFs, which has been shown to be quite large (Park et al. 2017) and may have otherwise obscured the analysis we wish to carry out. Our results show that the impact of framework point charges varies significantly according to the type of MOF structure. We also observe that some charge determination methods produce isotherms that are consistent with each other, while other methods lead to isotherms that deviate significantly from the rest.

2 Methods

2.1 Adsorption isotherm calculations

All adsorption isotherms reported in this work were calculated by Grand Canonical Monte Carlo (GCMC) simulations using Music 4.0 software (Gupta et al. 2003). This method works in the Grand Canonical ensemble, allowing for adsorption simulations in rigid porous materials by keeping the volume (V), temperature (T) and chemical potential (µ) constant but varying the number of adsorbate molecules. The conditions of GCMC simulation closely replicate experimental thermodynamic conditions of measuring adsorption. The simulations yield the absolute adsorbed amount as a function of µ, whereas in experiments, adsorption is typically measured as excess adsorption at varying pressure of the external fluid (P) (Duren et al. 2009). To be able to compare simulated data to experiment, µ has to be related to the gas phase pressure, which is done here by applying the Peng-Robinson equation of state (Peng and Robinson 1976). Furthermore, absolute adsorption measured by GCMC has to be converted to excess adsorption measured by experiment (Coudert and Fuchs 2016). In this work, however, we are interested purely in the effect of different point charges on the predicted adsorption isotherms without comparison to experimental data; hence, only absolute adsorbed amounts are reported. The GCMC simulation process consists of adsorbate molecules being randomly inserted, deleted, translated or rotated within the system, with all types of trial moves being equally weighted. This procedure generates an array of configurations with different potential energy, which are accepted or rejected depending on their Boltzmann probability, following the Metropolis method (Metropolis et al. 1953). To enhance sampling of the relevant regions of phase space, cavity bias insertion and deletion have been used, based on pre-calculated potential energy maps of each type of adsorbate-adsorbent interaction (Snurr et al. 1993). These maps, with a grid spacing of 0.15 Å, were also used to speed up the calculations by applying interpolation between grid points instead of calculating the potential energy on-the-fly at each step. Each simulation, corresponding to a single pressure point, was set for a total of 5.0 × 10⁷ iterations, composed of a 2.0 × 10⁷ cycle equilibration period, followed by a 3.0 × 10⁷ cycle production run. These lengths are sufficient to ensure convergence (see Supporting Information Figs. S5 and S6 for details) and were kept constant in all runs to enable direct comparison. The sampling period was split into 20 blocks for statistical analysis, and error bars were calculated as a 95% confidence interval—in most cases, they are approximately of the size of the symbols used, so they are not visible in the isotherm plots. For selected cases, we also carried out four repetitions of each isotherm using different random seeds, and found all repetitions to agree within the error bars determined above (see analysis in Supporting Information, Figs. S7–S10). Along an isotherm, simulation of each pressure point started from an equilibrated configuration at the pressure immediately preceding it. However, for selected cases we confirmed that the same uptakes, within statistical error, were obtained when simulations were started from an empty framework. All input data files underpinning this publication are openly available from the University of Strathclyde Knowledge Base at https://doi.org/10.15129/a6ad15b6-0762-4181-ad3a-24c70ca7a208.

2.2 Force fields

Fluid–fluid and fluid–solid interactions were calculated using models composed of Lennard-Jones (LJ) sites to describe repulsion and dispersion interactions, and partial point charges to describe permanent electrostatic interactions. The Lennard-Jones potential for an atomic pairs is calculated as:

$$ V_{\text{LJ}} = 4\varepsilon \left[ {\left( {\frac{\sigma }{r}} \right)^{12} - \left( {\frac{\sigma }{r}} \right)^{6} } \right] $$

where σ is the atomic radius (Å), ε is the depth of the potential well (kJ/mol), and r is the inter-atomic distance (Å). Lorentz–Berthelot mixing rules were used to compute the interactions between unlike atoms (Allen and Tildesley 1986):

$$ \varepsilon_{ij} = \sqrt {\varepsilon_{ii} \times \varepsilon_{jj} }\, \text{and} \,\sigma_{ij} = \frac{{\left( {\sigma_{ii} + \sigma_{jj} } \right)}}{2} $$

The electrostatic interactions between atoms ij can be accounted for using Coulomb’s law defined as

$$ {\text{V}}_{\text{c}} = k_{\text{e}} \frac{q_i \times q_j}{r} $$

where q is the partial point charge on each site, and k_e is the electrostatic constant (Jm/C²).

$$ k_{\text{e}} = \frac{1}{{4\pi \varepsilon_{0} }} $$

A LJ cut-off distance was defined as 13 Å. The electrostatic charges exert their effects over longer distances than LJ interactions, and this was accounted for using Ewald summations (Allen and Tildesley 1986) for MOF-adsorbate interactions and Wolf summations (Wolf et al. 1999) for adsorbate–adsorbate interactions. The Ewald summation separately considers short-range interactions in real space and long-range interactions in Fourier space in order to calculate the overall electrostatic potential. It is a highly accurate method which is used commonly to calculate electrostatic interactions in condensed phase systems. However, the Ewald method is fairly computationally expensive which makes it unsuitable for calculating the electrostatic potential of a large number of adsorbate molecules on-the-fly. The Wolf summation method is derived from the Ewald method, but involves a novel cut-off scheme ensuring charge neutrality of the system, thus increasing the calculation speed.

The Universal Force Field (UFF) (Rappe et al. 1992) and DREIDING (Mayo et al. 1990) are frequently used to model gas adsorption in MOFs and they were used in this work to describe the interactions with MOF framework atoms. More precisely, LJ parameters for the MOF metal atoms were taken from UFF and for non-metal atoms from DREIDING. For the adsorbate models, the Transferable Potentials for Phase Equilibria (TraPPE) (Potoff and Siepmann 2001) model was used for CO₂ and the SPC/E (Berendsen et al. 1987) model for water. Framework point charge sets generated by different methods were either obtained from literature or calculated in-house (see next section for details). To investigate solely the effect of varying framework point charges on the adsorption isotherms, the adsorbate–adsorbate potential (both LJ and charges), as well as the MOF-adsorbate LJ parameters, were kept constant in all calculations, and only the MOF-adsorbate electrostatic interactions were varied. This means that for each isotherm, a separate electrostatic interaction energy grid was calculated, followed by a full GCMC simulation. We studied IRMOF-1 (also known as MOF-5) (Li et al. 1999), MIL-47 (Barthelet et al. 2002), UiO-66 (in the fully hydroxylated form) (Cavka et al. 2008), Cu(dpa)₂SiF₆-i (in the interpenetrated form, abbreviated as SIFSIX) (Nugent et al. 2013), CuBTC (also known as HKUST-1) (Chui et al. 1999) and Co-MOF-74 (also known as CPO-27-Co or Co₂(dobdc)) (Dietzel et al. 2005), with all framework structures obtained from the Cambridge Structural Database.These MOFs were chosen with the aim of covering the most well-known and comprehensively studied “families” of MOF structures, as well as ensuring a large degree of topological diversity.

2.3 Point charge calculation methods

Determining point charges on interaction sites remain one of the most challenging aspects of force field development, mainly because they have no direct relation to experimental observables. As a consequence, a multitude of approaches have been suggested over the years, which can be broadly classified in three categories: (1) empirically fitting the charges to match some set of target experimental properties of the system of interest (e.g. thermodynamic properties in the case of pure fluids, or adsorption isotherms in the case of porous materials)—hereafter termed “Empirical”; (2) extracting the charges from ab initio quantum mechanical (QM) calculations—termed “QM-based”; (3) assigning charges based on chemical properties of atoms or small molecular fragments using a set of simple rules (e.g. fragment-based charges; electronegativity equalization methods)—termed “Semi-empirical”. There is also pronounced variability within each class of approaches: class 1 charges will depend strongly on the target properties; class 2 on the details of the underlying QM calculations and on the mathematical procedure to extract the charges; class 3 charges will depend on the properties and set of rules used to obtain them. It is not the purpose of the present paper to provide a comprehensive description of each charge determination method; instead, we will discuss only the key aspects pertaining to the methods employed here. The reader is referred to previous literature sources for additional technical details (e.g. Hamad et al. 2015; Sigfridsson and Ryde 1998; Verstraelen et al. 2016).

By far the most common approach to obtain point charges for MOFs falls into class 2 above—first a QM calculation is carried out on the whole framework or on suitably selected fragments (or “clusters”), followed by a mathematical analysis to extract point charges. This leads us to the first main distinction related to the nature of the QM calculation—periodic or cluster model representations of the molecular systems. Strictly speaking, one should use QM calculations on the entire framework (periodic model approach) to fully capture all electronic effects. Unfortunately, these calculations can become quite computationally demanding as the system sizes increase, and an alternative approach is to run the QM calculation on a smaller molecular fragment that is representative of the MOF functionalities. This, however, introduces a certain degree of arbitrariness in the choice of the cluster(s) and in the choice of atoms to cap the clusters truncated from a periodic framework. As we will see below, this may have a significant effect on the electrostatic interactions in adsorption simulations. It is also important to notice that several charge calculation procedures were not designed to work with periodic QM calculations, and are therefore restricted to cluster calculations. Finally, the level of theory of the QM calculation (e.g., exchange–correlation functional, basis set size) may also have a pronounced effect on the charge values.

Once the QM calculation has been carried out, the question is then how to extract the set of point charges that best represents the electronic environment of the molecule. The multitude of methods to do this can again be broadly classified into 3 types: (2.1) methods based on analysis of the QM basis set; (2.2) methods based on partitioning the electron density around each atom; (2.3) methods based on fitting to the electrostatic potential. The first type is based on the idea of assigning molecular orbitals to individual atoms using basis sets, and then summing up the population of electrons belonging to each orbital. These “orbital-based” methods then differ according to the choice of how to partition shared orbitals. For example, in the Mulliken method (1955), if a particular density matrix component corresponds to one basis set on atom A and another basis set on atom B, then half the electron population associated with this density matrix component is assigned to each atom. This makes it relatively simple and computationally efficient to obtain charges using this method, which is why it has been widely used. However, Mulliken charges are not only affected by the problems related with the different atomic electronegativities but are also very sensitive to changes in basis sets and molecular geometry. The Löwdin method (1950) suffers from similar problems, while the more recent Natural Population Analysis (Reed et al. 1985) is able to overcome at least some of those limitations, generally leading to more consistent point charges (Verstraelen et al. 2016). The more recent LoProp method of Gagliardi et al. (2004) extends this approach by taking into account multipole moments and polarizabilities. An important drawback of these methods is that they cannot be applied to periodic QM calculations using delocalized basis sets, unless plane-wave to atomic orbitals projection algorithms are employed (Dunnington and Schmidt 2012).

The second type, of which the Bader (1975) and Hirshfeld (1977) methods are the earliest examples, relies on a direct partitioning of the electron density—hence why they are often called “Atoms in Molecules” methods. The key difference with respect to the first type is that electron-density partition methods compute atomic charges as functionals of the electron and spin density distributions. Hence, they approach a mathematical limit as the basis set is improved towards completeness, and do not exhibit the problems of basis set dependence that are inherent in orbital-based methods like Mulliken or NPA. The Hirshfeld method works by assigning the charge density at each location in proportion to that of reference isolated atoms. This normally leads to rather low charges that tend to underestimate the electrostatic potential around each atom, which has led to the development of numerous improved methods over the years (Verstraelen et al. 2016). One example is the Iterative Stockholder Approach (ISA), in which the partitioning process is optimized iteratively leading to a more realistic representation of the electrostatic potential. Another method that has been widely used in adsorption simulations is DDEC (Density-Derived Electrostatic and Chemical charges) (Manz and Sholl 2010). DDEC has several advantages over the original Hirshfeld method, such as providing a much more accurate representation of the electrostatic potential and being generally applicable to both cluster and periodic QM calculations. A rather different approach of this type are Bader charges, which assign the charges based on a numerical analysis of the gradient and Laplacian of the electron density. Although they are quite useful to obtain chemical insight, they tend to strongly overestimate the electrostatic potential, and so have rarely been used in adsorption calculations.

Finally, the third type of QM-based charges are obtained through a direct fitting of the electrostatic potential around the molecule (or fragment) of interest. Several methods of this type are available for molecular clusters, with the most prominent examples being ESP (Electrostatic Potential derived charges) (Momany 1978), RESP (Restricted ESP) (Bayly et al. 1993), CHELPG (Charges from Electrostatic Potential using a Grid-based method) (Breneman and Wiberg 1990), and MSK (Merz-Singh-Kollman) (Singh and Kollman 1984). They differ mostly in the construction of the grid of points in which to fit the electrostatic potential (for example, MSK uses concentric surfaces around each atom, while CHELPG uses a uniform 3-dimensional grid) and in the setting of any constraints to the fitting (for example, RESP imposes constraints on buried atoms to try to ensure chemically realistic charges). This approach has the advantage that the charges are designed to explicitly reproduce the electrostatic potential energy around the molecule/fragment, which is precisely what then goes into the force field calculation. However, the nature of the numerical fitting process means that charges are often quite sensitive to details of the calculation (level of theory, basis set, conformation, etc.). Recently, a generalized version of RESP that is applicable to both cluster and periodic model approaches, called REPEAT, has been developed (Campana et al. 2009).

Determining point charges through approaches 1 (empirical) or 2 (QM-based) can be quite time consuming and computationally intensive. The advent of high-throughput screening of MOFs has brought the need to develop methods that can yield chemically reasonable charges with low computational requirements. One approach that has been widely used for this purpose is based on charge equilibration (Qeq) (Rappe and Goddard 1991). In a nutshell, this works by assigning point charges that minimize an energy function, which is constructed on the basis of measurable properties like electronegativities or ionization potentials. The choice of energy function separates the different varieties of this method, which include the Extended Charge Equilibration (EQeq) method of Wilmer et al. (2012a), developed specifically for MOFs. A different approach, also developed for MOFs, is the Connectivity-Based Atom Contribution (CBAC) (Xu and Zhong 2010). This method assigns charges to representative atoms of MOF building blocks with the same bonding environment, and is thus designed to be very fast and highly transferrable. The atomic charge database is itself constructed from QM calculations on a variety of small clusters representative of the most common MOF functionalities.

We have collected from literature reports a large number of framework point charge sets for all the MOFs under study, spanning all the charge calculation approaches described above, as well as different QM levels of theory (in the case of approach 2) including cluster and periodic calculations (computational details are collected in Tables S1, S2 of the Supporting Information). These were then used, together with the remaining force field parameters described in section 2.2, to generate individual adsorption isotherms for each point charge set. In some cases (notably for SIFSIX MOF), we have calculated our own framework charges from QM calculations, as follows. DFT calculations using periodic models of selected MOF structures were carried out with both VASP (Kresse and Hafner 1993; Kresse and Furthmüller 1996a, b) and CP2K (Laino et al. 2006) software. VASP calculations used the PBE functional (Perdew et al. 1996), without spin polarization, Projected Augmented Wave (PAW) potentials (Kresse and Joubert 1999) for core electrons, a cutoff of 415 eV for plane-wave basis sets and a grid of 1 × 1 × 2 k-points. CP2K calculations also used the PBE functional with a double zeta plus polarization (DZVP) basis set (Godbout et al. 1992) and optimised Goedecker-Teter-Hutter pseudopotentials (Goedecker et al. 1996). The energy cut-off selected was 500 Ry, the calculations used Γ-point sampling and spin polarization was accounted for. For the SIFSIX MOF, we also carried out cluster calculations using Gaussian 09 (Frisch et al. 2009), with the M06-L functional (Zhao and Truhlar 2006) and a 6-31G** basis set (Hariharan and Pople 1973). DDEC charges were computed using the Chargemol code and the DDEC6 variant (Limas and Manz 2018). Further details of the in-house QM calculations and charge determinations are provided in the Supporting Information (Figs. S1–S4 and Tables S3 and S4).

3 Results and discussion

3.1 IRMOF-1

IRMOF-1 (Fig. 1) belongs to a class of MOFs called the isoreticular MOFs (IRMOFs), which are characterised by their cubic topology. It was first reported in 1999. It forms a cubic network that consists of Zn₄O units joined by linear 1,4-benzenedicarboxylate links (Tranchemontagne et al. 2008). Even though it does not contain OMS, it has been shown to be very water unstable, similarly to other structures with Zn as a central atom, which limits the material’s potential applications (Castillo et al. 2008).

Fig. 1

a IRMOF-1 building block showing the different uniquely charged atoms. The corresponding charges are listed in Table 1. b IRMOF-1 framework structure

IRMOF-1 is one of the most widely studied MOFs in the scientific literature. As such, we were able to find a very large number of charge sets for this framework, reported in Table 1. We also carried out our own calculations using DDEC based on a periodic QM calculation (see section 2.3). It is clear from Table 1 that the point charges on each individual atom show a very wide variability (see also the averages and standard deviations at the bottom of Table 1). This is particularly the case for buried atoms, like Zn and O2, which tend to have a small effect on the electrostatic potential. We proceed to compare the adsorption isotherms for CO₂ and water at 298 K in IRMOF-1, obtained from GCMC simulations using each set of charges (Figs. 2, 3 and 4). A single plot comparing all charge sets is shown in Fig. S12.

Table 1 Charge sets for IRMOF-1 calculated by different methods, obtained from literature sources

Fig. 2

Adsorption isotherms of: a CO₂; b water in IRMOF-1 at 298 K using point charge sets obtained by periodic methods. Isotherms calculated without any framework charges are shown as a black line. Error bars are the size of the symbols used

Fig. 3

Adsorption isotherms of: a CO₂; b water in IRMOF-1 at 298 K comparing DDEC point charges (thick red line) to charges obtained from QM cluster calculations. The charge calculation method for each set is reported in the legend. Isotherms calculated without any framework charges are shown as a black line. Error bars are the size of the symbols used (Color figure online)

Fig. 4

A snapshot of CO₂ (green) adsorption in IRMOF-1 at 298 K and 250 kPa with a DDEC and b Fischer charges

Figure 2a compares adsorption isotherms of CO₂ in IRMOF-1 using framework point charges which were obtained from periodic QM calculations. First of all, our results corroborate previous conclusions that framework charges have a significant effect on adsorption simulations of CO₂ in MOFs (Kadantsev et al. 2013; McDaniel et al. 2015; Yang and Zhong 2006b; Zheng et al. 2009), since it is clear that the isotherm calculated without framework point charges lies significantly below all the other isotherms. All the Manz charge sets were generated from the same underlying DFT calculation, but differ in the method used to extract the point charges, so they make for a particularly interesting comparison. The Hirshfeld method produced much lower magnitude charges than the other methods, but this is consistent over all atoms and the relative difference between atoms is similar (see Table 1). Adsorption predictions using the Hirshfeld method are slightly lower than the rest, but the difference is not statistically significant. The ISA charges show very good agreement with DDEC except for C2 and C3 atoms, which are half and double in magnitude compared to DDEC. However, when translated into adsorption isotherm calculations, these differences in point charge magnitudes turn out to be negligible. In fact, it is quite remarkable that all isotherms obtained using framework charges from periodic calculations (periodic point charges) are statistically consistent with each other, with the exception of the Bader charge set. The latter is markedly different from the rest, both in shape and capacity, indicating much stronger framework/CO₂ interactions at low pressure. The isotherm implies the strength of this interaction is strong enough to overcome the intermolecular CO₂ interactions at low pressure, hence it does not follow the slight inflection of the other isotherms. The fact that experimental isotherms of CO₂ on IRMOF-1 exhibit this pronounced inflexion (Walton et al. 2008), points to the inadequacy of the Bader method for obtaining point charges for adsorption simulations. The tendency of Bader charges to overestimate adsorption in MOFs was previously reported by Liu et al. (2009). Overall, with the exception of Bader, all the periodic point charge sets should provide consistent adsorption predictions of CO₂in IRMOF-1. In the following, we will use the DDEC charges from the group of Manz and Sholl as reference points for comparison in all MOFs.

The water isotherms (Fig. 2b) show a type V shape, where the interaction between the water molecules is stronger than their interaction with IRMOF-1. Very slow uptake is initially observed, which increases once the surface coverage is sufficient to induce water clustering, causing a step-wise increase in the adsorbed amount. This is consistent with the previously reported hydrophobic nature of this material (Walton et al. 2008). We would expect that water, being a much more polar molecule than CO₂ due to a strong permanent dipole moment, would show a stronger effect of framework point charges on adsorption isotherms. It is therefore quite interesting that also in the case of water, all the adsorption isotherms obtained from periodic point charges, with the exception of Bader and Hirshfeld, show very good agreement with each other. The Hirshfeld method showed slightly lower adsorption in the case of CO₂ and this is magnified in the water adsorption isotherm, such that differences are now statistically significant. It is important to note that simulations of water adsorption in hydrophobic materials may suffer from convergence problems (Zhang and Snurr 2017). In this paper, we opted to run all isotherms with exactly the same number of MC trials since our aim was to compare simulated isotherms against each other, rather than against experimental data. In any case, we have tested the convergence of individual simulation points, shown in Figs. S5, S6 in Supporting information. These results, together with the fact that several repetitions of the calculations show close agreement with each other (Fig. S8), give us confidence that the isotherms represent converged equilibrium uptakes.

Figure 3 compares framework point charges obtained by cluster methods, using periodic DDEC charges as reference. The Belof charge sets use the same charge determination method and are both based on the same underlying Hartree–Fock calculation, but applied different basis sets. The fact that they produce very consistent results suggests that the basis set may not have a strong effect, at least with the MSK charge calculation method. Several data sets (e.g. most CHELPG sets) coalesce around the reference DDEC charges, yielding statistically indistinguishable isotherms for both water and CO₂, despite the fact that individual atomic point charges, such as on the Zn atom, vary greatly among them (see Table 1). This implies that the values of individual point charges are somewhat underdetermined, and that the important figure-of-merit for assessment of point charge sets should be how well they reproduce the underlying electrostatic potential. It also emphasises the usefulness of comparing point charge sets on the basis of adsorption isotherms, which is the ultimate goal of adsorption simulations, rather than focusing on the individual charges.

Four of the charge sets (Dubbeldam, Mu and the two Belof sets) lead to slightly stronger adsorption of CO₂, but these effects are magnified for water—as expected, the latter shows a more pronounced dependence on the electrostatic interactions. The Fischer charges produced the highest adsorption isotherms, and also showed the largest magnitude of charges on the aromatic ring atoms (see Table 1), possibly caused by over-polarization. CO₂ molecules adsorb quite strongly in the vicinity of the aromatic groups (see Fig. 4), which means that this enhancement of electrostatic interactions leads directly to higher adsorbed amounts. In fact, when we plot the CO₂ adsorbed amount against the electrostatic component of the fluid–solid interaction energy for both periodic and cluster charge sets (see Fig. S11), a linear relationship is observed, while the LJ component of the interaction energy is mostly the same (see Table S5). Overall, although most simulations using cluster-based framework charges lead to isotherms in agreement with periodic charges, the larger degree of variability observed in the cluster methods suggests that some degree of care should be taken when applying this approach to predict adsorption.

In Fig. 5, we compare framework point charges obtained by semi-empirical methods, again using periodic DDEC charges from Manz and Sholl as a reference. The original Qeq method yields charges that lead to a significantly lower adsorption isotherm for CO₂, and hence are likely underestimating the electrostatic potential. Interestingly, however, the same effect is not manifested in the water adsorption isotherm, which agrees closely with the reference one. The precise reason for this is unclear at this point. The improved EQeq method of Wilmer et al. (2012a), designed to reproduce QM charges in MOFs, yields isotherms in excellent agreement with the reference method. It is interesting to see, however, that an earlier version of EQeq actually leads to significantly higher adsorption than the reference isotherm. Finally, the CBAC method produces a CO₂ isotherm that is within statistical error of the DDEC isotherm, which is somewhat expected since this system was used as a test case for the development of the method. However, the CBAC water isotherm is significantly shifted to the left, emphasising that small differences in charges can lead to pronounced changes in adsorption of polar molecules. Due to the larger variability and uncertainty of these semi-empirical methods, care should be taken when accurate predictions of individual adsorption isotherms are required. However, they can provide quite useful alternatives to more computationally intensive methods for high-throughput screening of a large number of materials.

Fig. 5

Adsorption isotherms of: a CO₂; b water in IRMOF-1 at 298 K comparing DDEC point charges to charges obtained by semi-empirical approaches. Isotherms calculated without any framework charges are shown as a black line. Error bars are the size of the symbols used

3.2 MIL-47

MIL-47 (Fig. 6) is classed as a metal dicarboxylate. It forms a three-dimensional network consisting of chains of corner-sharing metal octahedra interlinked by benzene-dicarboxylate groups. Within this network there are one dimensional diamond-shaped pore tunnels. This material has been reported to have hydrophobic character and shows good water stability (Salles et al. 2011).

Fig. 6

a MIL-47 building block showing the different uniquely charged atoms. The corresponding charges are listed in Table 2. b MIL-47 framework structure

Several framework charge sets for MIL-47 obtained from literature reports are provided in Table 2. It is clear that the charges on the metal atom, which is only weakly exposed to adsorbate molecules, show the greatest degree of variation. Figure 7 shows adsorption isotherms for CO₂ and water in MIL-47 obtained using those point charge sets. As for IRMOF-1, most charge sets calculated from QM calculations, either cluster or periodic, yield CO₂ isotherms within statistical error of each other. In particular, the two sets of isotherms from periodic QM calculations show excellent consistency, despite variation in the values of the individual charges. The isotherm obtained with the EQeq method is slightly higher than the rest, which again suggests an inherent loss of accuracy in favour of computational speed. These differences are already statistically significant for CO₂ and are further amplified in the water isotherms (Fig. 6b), with the Mulliken charge set now showing a somewhat larger difference relative to the reference case. All simulated isotherms show a characteristic sigmoidal shape, again reflecting the relatively hydrophobic nature of this material.

Table 2 Charge sets for MIL-47 calculated by different methods, obtained from literature sources

Fig. 7

Adsorption isotherms of: a CO₂; b water in MIL-47 at 298 K using different point charge sets for the framework atoms. Isotherms calculated without any framework charges are shown as a black line. Error bars are the size of the symbols used

3.3 UiO-66

UiO-66, a zirconium-based MOF, consists of hexa-nuclear Zr₆O₄(OH)₄ inorganic nodes which form lattices via 1,4-benzene-dicarboxylate (BDC) linker, forming robust 3D structures (Fig. 8). This MOF is characterised by a high surface area and very high thermal stability. The inorganic centre contains polar -OH groups and this is reported to be the reason for this MOF’s exceptional stability, specifically its ability to undergo a reversible change from its hydroxylated structure to a dehydroxylated form without evoking any changes in the linked carboxylate ligands (Kandiah et al. 2010). DeCoste et al. (2013) investigated the water stability of carboxylate-containing MOFs, including CuBTC and Mg-MOF-74, and found that UiO-66 was the most water stable. This stability is likely to be the result of narrow pores and sterically hindered metal carboxylate sites, making these less accessible to water. The fully saturated metal centres don’t have the ability to coordinate with other molecules, unlike structures with OMS.UiO-66 has been known for its tendency to contain a significant number of defects in the form of missing linkers, which would affect the adsorption properties considerably (Ghosh et al. 2014). In this work, a fully hydroxylated structure with no defects was used for simulation purposes, so as to more reliably assess the effect of framework point charges on adsorption isotherms.

Fig. 8

a UiO-66 building block showing the different uniquely charged atoms. The corresponding charges are listed in Table 3. b UiO-66 framework structure

We were only able to identify four distinct sets of framework point charges for this MOF in literature (see Table 3), but these span all classes of charge calculation methods (semi-empirical, periodic QM-based and cluster QM-based). As for the two previous MOFs under study, the greatest variability in charge values is observed on the buried Zr and O3 atoms, with more exposed atoms exhibiting more consistent charge values. Predicted CO₂ isotherms (Fig. 9a) are mostly self-consistent with the exception of the CBAC charge set, which shows somewhat lower uptake than the remaining isotherms. The other three isotherms are statistically indistinguishable.Greater variation is seen in the water adsorption isotherms (Fig. 9b). Although the three QM-based isotherms show the same sigmoid shape and only slight variations in amount adsorbed, the CBAC isotherm shows a strangely high uptake at low pressure. Figure 10 compares snapshots of low pressure adsorption using DDEC and CBAC charges; the latter (Fig. 10b) shows clustering of water molecules in the secondary pores, which could lead to higher adsorption due to increased water–water interactions. Interestingly, while the CO₂ isotherm computed without framework charges shows only slightly lower adsorption than the reference DDEC isotherm, it very strongly underestimates water adsorption in the entire pressure range. Once again, this is caused by the pronounced sensitivity of water to details of the electrostatic interactions with the MOF framework.

Table 3 Charge sets for UiO-66 calculated by different methods, obtained from literature sources

Fig. 9

Adsorption isotherms of: a CO₂; b water in UiO-66 at 298 K using different point charge sets for the framework atoms. Isotherms calculated without any framework charges are shown as a black line. Error bars are the size of the symbols used

Fig. 10

A snapshot of water (blue) adsorption at 298 K and 250 kPa in UiO-66 with a DDEC and b CBAC charges (Color figure online)

3.4 CuBTC

CuBTC (also known as MOF-199 and HKUST-1) was first discovered by researchers at the Hong Kong University of Science and Technology in 1990 (Fig. 11). It is sold commercially as Basolite C300 and it is one of the most frequently studied MOFs (Moghadam et al. 2017). CuBTC consists of central copper ions linked with 1,3,5-benzenetricarboxylate (BTC) acid ligands. These linkages form a porous crystalline structure with two central copper ions which are bound to 4 BTC ligands via 2 oxygens on each ligand and to solvent (usually water) via 2 oxygens. Activation of CuBTC results in removal of solvent molecules and leaves two coordinatively unsaturated copper ions with available binding sites. The copper ions form open metal sites where the metal atom is exposed, and these sites have a high affinity and selectivity for electron-donating adsorbates. Consequently, CuBTC has a high affinity for water and will readily adsorb moisture (Castillo et al. 2008).

Fig. 11

a CuBTC building block showing the different uniquely charged atoms. The corresponding charges are listed in Table 4. b activated CuBTC framework structure

Due to the large number of studies carried out on CuBTC, several distinct framework charge sets could be obtained from literature, as shown in Table 4. It is interesting to see that in this case, the charge variability is practically uniform among all atoms of the framework. This is because there are no buried atoms in CuBTC, and even the metal sites are exposed to the surface, hence they all contribute significantly to the electrostatic potential. Nevertheless, as we will see below, the observed variation in the charge magnitude at the OMS is likely to have a pronounced effect on adsorption in this MOF.

Table 4 Charge sets for CuBTC calculated by different methods, obtained from literature sources

As for IRMOF-1, due to the large number of charge sets considered, we have opted to split the comparison in two groups (a single plot showing all isotherms is available in Fig. S13). In Fig. 12 we collect isotherms obtained with periodic QM charges and with semi-empirical charges, while in Fig. 13 we compare isotherms obtained with cluster QM charges. For CO₂ (Figs. 12a and 13a) all the QM-based charge sets (both periodic and cluster) yield adsorption isotherms in good statistical agreement with each other, with a variability that is certainly well within the observed uncertainty in experimental isotherms reported by Park et al. (2017). The EQeq charges also produce an isotherm in agreement with the QM-based ones. However, the CBAC method and the empirically adjusted Castillo charges lead to a statistically significant increase in adsorbed amounts, albeit not by a large extent. In water adsorption (Figs. 12b and 13b), the variation in adsorption isotherms is much more significant. As for the previous MOFs, the isotherms based on periodic QM charges are still in agreement with each other, reinforcing the consistency of this charge determination approach, even for MOFs with OMS. The QM cluster-based isotherms show a much more significant degree of variability, with one of the charge sets leading to a two-fold decrease in the pressure at the isotherm inflection point (from ~ 1250 to ~ 600 kPa). The semi-empirical sets also show rather extreme differences from the reference periodic calculations, even for the EQeq charges, which had yielded CO₂ isotherms in good agreement with DDEC.

Fig. 12

Adsorption isotherms of: a CO₂; b water in CuBTC at 298 K using point charge sets obtained by periodic QM and semi-empirical methods. Isotherms calculated without any framework charges are shown as a black line. Error bars are the size of the symbols used

Fig. 13

Adsorption isotherms of: a CO₂; b water in CuBTC at 298 K comparing DDEC point charges to charges obtained by cluster methods. Isotherms calculated without any framework charges are shown as a black line. Error bars are the size of the symbols used

The large variety of water adsorption isotherms indicates that the adsorption process of water in CuBTC is different from that of CO₂, which could be due to the large dipole moment of water molecule, as opposed to CO₂ which has only a permanent quadrupole, as well as due to specific interactions with the open Cu site. In fact, one can observe a broad correlation between the steepness of the slope of the water isotherms at low pressure and the point charge on the Cu atom (see Table 4). For example, the highest Cu charge of + 1.25 is in the Castillo set, which shows one of the strongest water adsorption isotherms. The Liu1, Yazaydin2 and Yang sets show almost identical slopes and all have Cu charges around + 1.1; the Babarao, Fischer and Liu2 sets have slightly lower slopes and Cu charges around + 1.05; finally, the periodic sets of Nazarian, Wilmer and Zang have the smallest slopes and Cu charges around + 0.9. The correlation is not perfect, however, with the CBAC, Yazaydin1 and EQeq sets showing steep slopes and relatively low Cu charges. Nevertheless, the sensitivity of water adsorption at low pressure on the Cu point charge indicates that the presence of OMS plays a key role in the water adsorption behaviour. This suggests that great care must be taken when selecting framework point charges for simulating adsorption of strongly polar molecules such as water in MOFs with open metal sites. In this context, it is worth noting that the Castillo charges were specifically designed to match available water adsorption data on CuBTC, yet yield isotherms for both water and CO₂ that are not consistent with QM-based approaches. The same could be said, to some extent, of the CBAC charge set, which attempts to empirically account for the OMS interaction and incorporate it in the charges, resulting in over-polarization and increased adsorption prediction. Based on previous work on MOFs with OMS, a more physically-grounded approach would involve using consistent QM-based charges while separately treating the specific interactions between water and the metal site through a bespoke interaction model (Campbell et al. 2017, 2018; Fischer et al. 2012, 2014).

3.5 Co-MOF-74

Co-MOF-74 (Fig. 14) belongs to a family of MOFs designated as M-MOF-74 (M = Zn, Ni, Co, Fe, Mg, Mn, Ca, or Sr). It consists of 5-coordinated metal ions linked to 2,5-dioxoterephthalate, forming wide 1-dimensional hexagonal pores around 1.1 nm in diameter. This MOF contains a high concentration of OMS that are formed upon removal of the solvent molecules attached to the metal. Co-MOF-74 and its analogues have been reported to have very low water stability, with their structure degrading after exposure to even small amounts of moisture (DeCoste et al. 2013).

Fig. 14

a Co-MOF-74 building block showing the different uniquely charged atoms. The corresponding charges are listed in Table 5. b Co-MOF-74 framework structure

Six distinct framework point charge sets for this MOF were gathered from literature reports (Table 5). Most charges show relatively low variability, with the exception of Co. In fact, comparing with CuBTC (Table 4), for which many more data sets were available, the fluctuations in the charge of the unsaturated metal are quite significant for Co-MOF-74. It should be noted, however, that this is mainly due to the unphysically low charge produced by the EQeq method, as well as the rather large charge from LoProp (for which no data was found on CuBTC).

Table 5 Charge sets for Co-MOF-74 calculated by different methods, obtained from literature sources

Figure 15 shows simulated adsorption isotherms on Co-MOF-74 using different framework point charge sets. The point charges obtained from the LoProp method carry the highest charge on the metal centre and hydrogen and lowest charge on most of the other atoms. In both CO₂ (Fig. 15a) and water (Fig. 15b), this isotherm deviates the most from the rest and predicts the highest uptake. Conversely, the EQeq method has by far the lowest charge on the metal, and consequently predicts the lowest adsorbed amount for water. However, the same is not observed in the CO₂ isotherm, for which the EQeq isotherm is the second largest. It is clear that the effect of electrostatic potential around the metal site has a much more pronounced effect on water adsorption, as observed above for CuBTC. Co-MOF-74 has a large number of OMS facing the hexagonal pores, which explains the high sensitivity of the adsorption prediction on the metal centre charge. Unlike the rest of the studied MOFs, the type V isotherm due to weak water-adsorbent interactions is not observed here. Again, this is likely to result from the high concentration of OMS pointed directly into the pore channel, facilitating the interaction between these sites and the water molecules. Finally, it is important to note that once again most of the QM-based charge sets, with the exception of LoProp as discussed above, lead to consistent isotherms for both gases.

Fig. 15

Adsorption isotherms of: a CO₂; b water in Co-MOF-74 at 298 K using different point charge sets for the framework atoms. A plot of the water isotherms using a larger pressure range is provided in Fig. S14. Isotherms calculated without any framework charges are shown as a black line. Error bars are the size of the symbols used

3.6 Sifsix-2-Cu-I

This MOF belongs to the SIFSIX family and consists of a copper centre and two different ligands, one organic linked to the metal via a nitrogen atom, the other based on silicon surrounded by six fluorine atoms (Fig. 16). MOFs of this family have orthorhombic unit cells and are prone to interpenetration. In fact, we focused specifically on an interpenetrated member of this family, in order to assess the impact of this feature on the electrostatic interactions.

Fig. 16

a Sifsix-2-Cu-I building block showing the different uniquely charged atoms. The corresponding charges are listed in Table 6. b Sifsix-2-Cu-I framework structure

For this MOF, to the best of our knowledge only one published charge set was available, from Pham et al. (2013). As such, we calculated several sets of framework point charges from both cluster and periodic QM calculations, as described in section 2.3 and in the Supporting Information. The full sets of charges are provided in Table 6. As observed for other MOFs, the largest variation in charge values is for the buried Cu and Si atoms, which are both fully coordinated and barely accessible to adsorbate molecules.

Table 6 Different charge sets for Sifsix-2-Cu-I calculated using different methods

Adsorption isotherms of CO₂ on SIFSIX MOF are reported in Fig. 17. The first observation to make is that, as was observed in all the other MOF structures studied in this paper, all isotherms calculated using periodic QM charges are consistent with each other, despite the fact that the charges themselves show significant fluctuations. This observation is independent of the charge determination method (e.g. compare REPEAT VASP and DDEC VASP, which were both obtained from the same underlying QM calculation) and of the type of basis set employed (e.g. DDEC VASP used plane waves while DDEC CP2K used Gaussian plus plane waves). The framework charges extracted from cluster calculations by Pham et al. show substantially lower adsorption isotherms than any of the periodic charge sets. Our own cluster-based charges, obtained from rather large molecular fragments using two different charge determination methods, yield isotherms that are on either side of the periodic ones, albeit much closer than the isotherm obtained from Pham charges. These differences are larger than the statistical error of the simulations. This reinforces the earlier observations that cluster-based charges lead to much more significant variations in adsorbed amounts, and care should be taken when applying them directly without any prior consistency check. In this context, the very good agreement between simulations and experiments observed in the work of Pham et al. (2013) appears rather fortuitous. A more detailed analysis is needed to clarify this issue, but the fact that this particular MOF has very narrow pore spaces (due to interpenetration) could mean that strong adsorption sites become much more sensitive to details of the electrostatic potential energy surface.

Fig. 17

Adsorption isotherms of CO₂ in Sifsix-2-Cu-I at 273 K using different point charge sets for the framework atoms. Error bars are the size of the symbols used

4 Conclusions

In this paper, we reported a detailed and systematic analysis of the effect of the choice of framework point charges on adsorption isotherms predicted by molecular simulations. Point charges obtained by different approaches, covering all types of methodology, were used in GCMC simulations of CO₂ and water adsorption in six different MOF structures. The latter covered some of the most widely studied MOF families, and included frameworks with significantly different topologies and surface functionalities, as well as MOFs with and without open metal sites. Our results allow us to draw the following general conclusions:

• The variation in the values of the framework point charges, for any given MOF, is much larger than the variation in the adsorption isotherms themselves. This is particularly the case for QM-derived charges, which suggests that the main property controlling adsorption predictions is the overall electrostatic potential induced by the framework on the adsorbate molecules. As a consequence, assessing the suitability of point charge determination methods solely on the basis of the charges themselves may lead to erroneous conclusions. We recommend diagnosing the suitability of charge sets by comparing predicted isotherms against reference data, whenever possible.

• Charges derived from periodic QM calculations yielded isotherms that were consistent with each other for all charge determination methods, with the exception of Bader and Hirshfeld. This was the case regardless of the details of the underlying QM calculation. In particular, it is noteworthy that consistent charges were obtained from both DDEC (an electron density partitioning method) and REPEAT (an electrostatic potential fitting method) for all MOFs analysed here. It is hard to draw definitive conclusions for other charge determination methods (e.g. LoProp, ISA) due to the small number of instances analysed. As such, we would recommend this approach as the most reliable for obtaining framework charges for adsorption simulations, at least when the size of the unit cell does not make such calculations prohibitive.

• Framework charges derived from QM cluster calculations were comparatively less reliable—while several sets yielded isotherms in agreement with each other and with periodic charges, there was a significant degree of variation in some cases. Moreover, with one or two exceptions, it was not possible to predict which sets would lead to discrepant isotherms simply by examining the charge values. This means that care should be taken when using charges calculated by this approach, and consistency checks should be sought whenever possible.

• Charges calculated using methods that fit the electrostatic potential (such as CHELPG or REPEAT) or that partition the electron density (such as DDEC or ISA) generally yield isotherms that are statistically consistent with each other. A clear exception are Bader charges, which strongly overestimate the electrostatic potential, and hence the adsorbed amounts. Charges determined from population analysis (e.g. Mulliken) appear to be less reliable, although there are not enough examples in our study to confirm this observation. The decreased reliability of Mulliken charges is reinforced by their significant dependence with the basis set size as reported in the literature.

• Semi-empirical approaches, like EQeq or CBAC, can provide reasonable alternatives to QM-based charges when computational expense is an important limitation (e.g. large-scale screening). However, in some cases, predictions from this class of methods can lead to rather unexpected results (e.g. water in UiO-66). Given the wide variability in MOF framework structures and surface chemistries (including functionalization), we recommend that any charges obtained from semi-empirical methods be validated against reference isotherms obtained from periodic QM charges for prototypical MOFs.

• Water isotherms are much more sensitive to details of the electrostatic interactions than CO₂ isotherms. This was rather expected, due to the significant difference in polarity between these two adsorbates. This effect is emphasised in MOFs that contain open metal sites, as these provide rather strong adsorption sites for water molecules. In such MOFs, attempts to adjust framework point charges (or, indeed, LJ parameters) to implicitly account for coordination bonds at the unsaturated metal site are unlikely to lead to physically consistent adsorption behaviour.

We believe our work will help to improve the reliability and reproducibility of adsorption simulations by providing useful guidelines to calculate or choose point charges for MOF frameworks. By providing a consistent set of adsorption isotherms on reference materials, it should now be possible to systematically assess the effect of other force field parameters (e.g. Lennard-Jones) on the performance of the simulations, and we intend to report on this issue in due course. Together with recent systematic analysis of uncertainty in experimental adsorption measurements (Park et al. 2017), this work enables a less arbitrary assessment of the suitability of molecular models to provide adsorption predictions that can be used in an industrial context.