Rapid adsorption enthalpy surface sampling (RAESS) to characterize nanoporous materials

Molecular adsorption in nanoporous materials has many large-scale industrial applications ranging from separation to storage. To design the best materials, computational simulations are key to guiding the experimentation and engineering processes. Because nanoporous materials exist in a plethora of forms, we need to speed up the existing simulation tools to be able to screen databases of hundreds of thousands of structures. Here, we describe a new algorithm that quickly calculates adsorption enthalpies by sampling the surface of the material instead of the whole porous space. This surface sampling has been tested on the CoRE MOF 2019 database and has been proven to be more than 2 orders of magnitude faster than the gold standard method (Widom insertion), with an acceptable level of error on an enthalpy value of 0.34 kJ mol−1, and is therefore proposed as a valuable addition to the high-throughput screening toolbox.


S1 Adsorption enthalpy
The adsorption enthalpy values of structures from CoRE MOF 2019 have been calculated using Widom insertion and the RAESS algorithm. In this section, we first give a visual comparison through some scatter-plots. We can visually see the improvement from the initial algorithm to the final algorithm. The raw data of the convergence plots are also shown to highlight the simulation time needed to reach convergence. Finally, the data of the performance plot that compares every enthalpy calculation methods (Widom insertion, Voronoi node sampling, RAESS). Figure S1: Scatter-plot of the enthalpies calculated by our initial algorithm (λ = 2 1/6 and µ = 0) compared to the enthalpies calculated by a 100k step Widom insertion simulation of xenon in structures of CoRE MOF 2019 at 298 K. Figure S2: Scatter-plot of the enthalpies calculated by our initial algorithm (λ = 2 1/6 and µ = 0) compared to the enthalpies calculated by a 100k step Widom insertion simulation of xenon in structures of CoRE MOF 2019 with LCD ≥ 3.7Å at 298 K. Figure S3: Scatter-plot of the enthalpies calculated by our final algorithm (λ = 1.6 and µ = 0.85) compared to the enthalpies calculated by a 100k step Widom insertion simulation of xenon in structures of CoRE MOF 2019 with LCD ≥ 3.7Å at 298 K. The RAESS method relies on the higher weight of the strong sites close to the surface of the pores. If we increase the temperature, the less attractive sites would play an increasing role and the accuracy of the method would drop. To grasp this limitation of the RAESS algorithm at higher temperature, we compared the results of a screening over the CoREMOF 2019 database. Figure S4: Scatter-plot of the enthalpies calculated by our final algorithm (λ = 1.6 and µ = 0.85) compared to the enthalpies calculated by a 12k step Widom insertion simulation of xenon in structures of CoRE MOF 2019 with LCD ≥ 3.7Å at 600 K.
The method is as expected less accurate but it still gives a reasonable correlation on the performance, with an RMSE of 0.70 kJ mol −1 and an MAE of 0.41 kJ mol −1 . The errors have almost doubled when going from 298 K to 600 K.
However, these limitations of the method are not crippling since adsorption processes are usually not performed at very high temperature. High temperatures are commonly used in temperature swing adsorption (TSA) to desorb the adsorbates rather than to adsorb them.

S3 Henry constant S3.1 Derivation of the formula
Here we give a general derivation of the relation between the Henry constant K H and the interaction potential E int of an adsorbate in a porous material. First, for an ideal gas the amount of adsorbed molecules n ads can be expressed using the bulk density of the gas ρ ads,bulk and the volume of the pores V pore : The pore volume can be seen as the continuous sum of each voxel times the Boltzmann probability of presence, which is represented by the following integral of the Boltzmann factors. This integral can then be changed to the average of the Boltzmann factors: If we apply the equation S2 and the ideal gas equation of state P = ρ ads,bulk RT on the bulk gas in equilibrium, we can change the equation S1 to: If we now consider the gravimetric loading L (in mol kg −1 ), we need to divide the equation by mass density of the framework ρ f : Since the Henry's law is described by L = K H × P , we have the final equation:

S3.2 Data comparison
Here, we present the data on the Henry constants calculated by Widom insertion and the final RAESS algorithm. Their accuracy and time efficiency were compared thoroughly.

S4 Xe/Kr Selectivity
Let us consider the following exchange equilibrium between the free gas and the adsorbed gas to model the competition between the xenon and krypton for adsorption sites.

S4.1 Discussion on the rejection condition
To calculate a selectivity, the rejection condition on xenon can be high since we are interested in the most favorable materials for xenon adsorption. But for krypton, we want to accurately describe very low Henry constants, because a selective material would also be a material unfavorable to krypton. For these reasons, the µ parameter needs to be chosen wisely and needs to be low enough to have accurate selectivities.
Here, we give some values of the error on the selectivity depending on the rejection condition we chose. According to the quick study here, the optimal value is µ = 0.5 since it gives the best accuracy in a minimal amount of time.
S4.2 Correlation analysis for µ = 0.5 The selectivity can be compared directly using the a log-scale plot and log-scale metrics. If  To be able to give a thermodynamic interpretation, we can define an exchange Gibbs free energy associated to this selectivity: Using this exchange Gibbs free energy, we can assess much more easily the performance of the approach. The RMSE is about 0.36 kJ mol −1 . We cannot compare it to the adsorption enthalpy errors, since the ranges and interpretation are very different. Here, the selective materials have a negative ∆ exch G Xe/Kr 0 and it goes to a maximum value of about −12.7 kJ mol −1 . The relative error is of course higher on the Gibbs free energy. This is due to a higher uncertainty on the Henry constant and to the denominator term brought by the krypton.

S5 Surface area
Comparison of the surface areas calculated by RAESS and by RASPA2. In our algorithm, we sample at a further distance from the atoms than a surface area calculation algorithm would do. For this reason, the surface areas calculated by our algorithm are correlated to the ones given by RASPA2 by are very different. Figure S10: Scatter plot comparison of surface area calculated by our final algorithm (sampling at 1.6 × σ) and the ones calculated by Raspa with the same values of sigma (m 2 cm −3 ) However, when modifying our code to match the standard surface sampling algorithm, we have a better agreement between both methods.  We can note that the points with weaker correlations correspond to the ones with an LCD greater than 10Å.

S6.2 Amorphous materials
We tested our algorithm on the amorphous database and retrieved the top 20 materials for xenon (the whole csv file can be found on our group's GitHub). The Raspa software is not able to run on these amorphous structures with our computers, the memory is filled up and the job is killed (didn't investigate where this comes from). Therefore, there is no comparison point with a Widom simulation from Raspa.