Analysing adsorption isotherms with MuSiC¶

Tina Düren <t.duren@bath.ac.uk> and Gaël Donval <g.donval@bath.ac.uk>

The purpose of this tutorial is to see what information can be extracted from gcmc calculations using the example of the study of $\mathrm{CO}_2$ adsorption in two very similar metal-organic frameworks (mofs), Irmof-1 and Irmof-3. You won’t perform any simulation yourself this time around and we won’t describe in detail how such a simulation is practically carried out: this will be the object of a more advanced tutorial.

Objectives

Understanding what mofs are through two examples: Irmof-1 and Irmof-3.
Learning how to analyse gcmc results carried out with MuSiC (the multipurpose simulation program we used to run the gcmc calculation).
Getting full adsorption isotherms for $\mathrm{CO}_2$.
Learning how to compare different mofs with respect to their adsorption properties.

Background¶

Metal-organic frameworks have recently attracted much attention as materials for carbon capture materials due to their high surface area, low material density and the sheer endless possibilities in combining metal centres with organic linkers in a kind of chemical “Lego” set (see Fig. 11). The linkers can be functionalised to tailor the host guest interactions further.

Fig. 11 a) Schematic representation of the self-assembly process of mofs from corner and bridging units. b) Building blocks for Irmof-1 and Irmof-3. c) Structure of Irmof-1, one complete, larger cavity of the framework is shown. The transparent blue sphere was added to demonstrate the size of the cavity and to emphasize the open three-dimensional framework. Most mofs are porous which explains why they are used for gas adsorption and separation.¶

We are interested in Irmof-1 and Irmof-3 for $\mathrm{CO}_2$ adsorption: as you can see in Fig. 12, they have very similar structures: Irmof-3 is basically just Irmof-1 with an extra amine group on its ligands. This similarity can also be seen in Tab. 2 where we compared different physical properties.

Fig. 12 Illustration of Irmof-1 and Irmof-3.¶

Tab. 2 Physical properties of Irmof-1 and Irmof-3.¶
Property	Irmof-1	Irmof-3
density (g/cm$^3$)	0.59	0.67
pore volume (Å$^3$)	13916	13561
molecular weight of unit-cell (g/mol)	6156	6852
volume of unit-cell (Å$^3$)	17237	17067
accessible surface area (m$^2$/g)	3717	3245
accessible surface area (m$^2$/cm$^3$)	2156	2174
pore size (Å)	10.8\|13.9	9.1\|13.9

One the one side, Irmof-1 seems to be more porous but on the other side the amine groups in Irmof-3 could favour interactions with $\mathrm{CO}_2$ and increase the uptake. To determine which framework has better $\mathrm{CO}_2$ adsorption properties, we ran a few gcmc calculations. We obtained adsorption isotherms for full range of pressure points and at different temperatures: $298\,\mathrm{K}$, $313\,\mathrm{K}$ and $328\,\mathrm{K}$.

You are going to copy the result of those calculations to your directory and we are going to analyse them.

Preparing the analysis¶

Connect to Balena the usual way and get on a computing node:

$ sint
$ module load music/tutorial

The module commands are used to make MuSiC available to you (see How do you access specific programs? for more details about the module system).

As we mentioned at the beginning, you won’t be running simulations this time. You will need to copy files from the following location on Balena to your scratch space:

/beegfs/scratch/group/ce-molsim/tutorials/music_analysis.tar.xz

This is obviously an archive file that needs to be uncompressed first.

Running a gcmc calculation in MuSiC¶

Now if you have a look at those files you copied over, there are two folders each containing 3 folders (one for each temperature point). Let us start by going into IRMOF1/298K.

gcmc calculations with MuSiC typically starts from a so-called control file, in this case, gcmc.ctr. If you have a look in it, you should see different sections (starting with ------) with an explanation for each line. Note that # serves as a comment delimiter: everything following a # is intended to be read by human beings and is completely ignored by MuSiC. Most elements should be self-explanatory: you are not expected to understand everything yet and those of you who are going to use MuSiC will follow dedicated tutorials.

In General Information though, you should see a reference to restart files contained in a restart/ folder: these files are written by MuSiC and as the name suggests, they are used to restart/resume a calculation. The same thing goes for config (configuration) files in config/ that are the data files containing the results of the calculations. These are the files we are going to postprocess and analyse.

The section Atomic Types defines the name of each atom used in the simulation followed by a file name (.atm). Molecules (and framework) are defined in a very similar way in the Molecule Types section (but with .mol files instead). These files exist somewhere on the disk but describing them would be going too far for our current purpose.

The temperature of the simulation is set under GCMC Information and the pressure points are defined in terms of fugacity in the same section: they are read from fugacity_298K.dat. If you open that file, you should see a list of 15 fugacity points.

Note

At low pressure, fugacity and pressure can be considered equal so it is fine if you think in terms of pressure.

The calculation itself was run by executing the script run_gcmc. A script is a list of commands gathered into a text file: it is equivalent to typing each command in turn at the command line without the hassle of actually remembering and typing all of them.

Warning

Don’t run run_gcmc: it would require hours to finish!

You can have a look at it though using less: the export commands define global variables that MuSiC is going to use (in this case to know where to look for atomic and molecular definitions). It is followed by an actual call to MuSiC: music_gcmc gcmc.ctr.

Note

Global variables are like information shared at the terminal level. All the programs run in it can then access and use them. MuSiC needs ATOMSDIR, MOLSDIR, PMAPDIR and EMAPDIR to be set to an appropriate value to work

Getting raw isotherms¶

We need to extract the isotherms from the config files we were mentioning earlier. For that we need to use a new MuSiC subprogram: music_post and a suitable control file dedicated to post-processing. In principle the following should work:

$ music_post post.ctr
Please set your ATOMSDIR environment variable

But as shown above, it doesn’t because the ATOMSDIR global variable is not set. You could reuse the commands from run_gcmc:

$ export ATOMSDIR="../music_files/atoms"
$ export MOLSDIR="../music_files/molecules"
$ export PMAPDIR="../music_files/maps"
$ export EMAPDIR="../music_files/maps"
$ music_post post.ctr

These would work perfectly well! However this is very cumbersome to remember and type every time! Instead we are going to make a script: all these commands can be typed in a simple text file and executed all at once just by calling that file.

To do that, just copy and paste the commands above in a new file (like run_post). Be sure you don’t paste the $ as they only are visual indications. Then add #!/bin/bash on its own at the top of the file. Make it executable:

$ chmod u+x run_post

Finally call it:

$ ./run_post

Neat, isn’t it?

Hint

You could instead just create a copy of run_gcmc and replace the last line, music_gcmc gcmc.ctr, by music_post post.ctr.

MuSiC will create three new files:

isotherm.IRMOF1 that does not contain any gcmc data and can be ignored (Irmof-1 is the framework: it is never changed!).
isotherm.CO2 that contains the raw loading/uptake (in molecules per unit-cell of Irmof-1) vs. fugacity. The file also ends by an attempt to fit the Langmuir isotherm equation to the data (not a good idea in this case, do you see why?). We are not going to use this file directly though because we need an extra processing step to get all the relevant data.
IRMOF1.CO2.298K.post contains information about energy and loading for each pressure point as well.

Converting fugacities to pressures¶

These files provide raw results out of the calculation. However experimentalists measure pressures not fugacities. We need to convert the quantities we calculated to quantities they measured. The program music_analysis does just that:

$ music_analysis < analysis.ctr

The file analysis.ctr is quite self-explanatory: it just tells music_analysis where to find the data it needs (i.e. the output files of music_post) and where to output its results. This will create two files:

energy.IRMOF1.298K.dat which contains various energy-related information; we won’t use it and you can ignore it.

adiso.IRMOF1.298K.dat which contains the isotherm data; this is the file you want to use for your plots.

Plot N_abs vs. pressure and vs. fugacity. Can you see why we said pressure and fugacity are roughly the same for low values? When would you consider them different?

Hint

You can create a run_analysis script as well to automate things or even better: why don’t you add the command above at the end of run_post? Then you’d just need to call ./run_post to get isotherm data and post-process them to get pressures!

Understanding absolute and excess uptake¶

The output of a gcmc simulation is the absolute amount adsorbed, i.e. the total number of adsorbate molecules present in the pore, whereas experimentally the excess amount adsorbed is measured. The difference between the absolute amount adsorbed and the excess amount adsorbed is shown schematically in Fig. 13.

Fig. 13 Schematic illustration the concept of excess mount of adsorbed particles. a) Free particles naturally present in the framework if there were no interactions (i.e. free gas in “empty box” at given pressure and temperature). b) Absolute amount of adsorbed particles: raw number of particles in the box. c) Excess amount of particles: particles that are there because of the interactions with the framework.¶

The density of fluid molecules close to the solid is higher than the bulk density due to the intermolecular interactions between the fluid molecules and the pore walls. Yet, even without adsorption taking place, a certain amount of gas molecules would be present in the pore volume depending on pressure and temperature. The excess number of molecules, $n_\mathrm{ex}$, is therefore related to the absolute number of molecules, $n_\mathrm{abs}$ the following expression:

\[n_\mathrm{ex} = n_\mathrm{abs} - V_g \rho_g\]

Here, $V_g$ is the pore volume and $\rho_g$ is the molar density of the bulk gas phase.

In practice, $V_g$ was calculated with the second virial coefficient as described here. In our case, MuSiC used the Peng-Robinson equation of state to calculate the bulk gas density (case depicted by Fig. 13a). The pore volume of Irmof-1 and Irmof-3 are given in Tab. 2.

You may have noticed that adiso.IRMOF1.298K.dat contains a quantity conspicuously called N_ex: this is the excess amount of adsorbed $\mathrm{CO}_2$ we were talking about. Compare N_abs and N_ex vs. pressure. Does it make a large difference? Would N_abs be enough to compare adsorption properties of Irmof-1 and Irmof-3?

Producing snapshots¶

As for the molecular dynamics simulations, snapshots are useful for visualisation. However unlike in molecular dynamics simulations, full movies aren’t very informative. Do you see why?

Here instead, we are going to look at single configurations, one for each pressure point. We have conveniently provided another control file for you to use: snapshot.ctr. The only line you need to focus on is this one (near the end of the file):

RESTARTFILE restart/CO2.IRMOF1.298K.1

It is telling MuSiC which restart file to use. To create the snapshot, copy run_gcmc to run_snap and modify the last line to look like this:

music_gcmc snapshot.ctr

Then save the file and execute it:

$ ./run_snap

MuSiC will then create a new file, finalconfig.xyz, which you can visualise in Vmd. This snapshot was taken for the first config file which corresponds to the first pressure value in pressure.dat: $10\,\mathrm{kPa}$.

Give that file a relevant name (finalconfig.xyz will be overwritten each time you take a snapshot):

$ mv finalconfig.xyz snap.10kPa.298K.xyz

Note

When you visualise a snapshot with Vmd, you will notice that you have $\mathrm{AgNe}_2$ molecules in your framework. Fear not, this is just a visualisation trick to differentiate the carbon and oxygen of $\mathrm{CO}_2$ and the ones present in the framework! Said differently: the calculations were performed with $\mathrm{CO}_2$ but we changed them to $\mathrm{AgNe}_2$ afterwards to have nice distinct colors.

Hint

It is easier to look at the structures if they are shown in the orthographic mode (Display ‣ Orthographic)

To have a look at snapshots from other pressure points, go into pressure.dat, choose the pressure point you want to visualise, determine which restart file number it corresponds to, select the correct restart file in snapshot.ctr and execute your script.

Tasks¶

You should now have all the tools (and files) to perform the following tasks and answer the corresponding questions.

Hint

Though we don’t provide any extra control file or script in the other folders, keep in mind that the ones you used until now are perfectly transferable to those folders. So you don’t have to really restart from scratch every time.

Explore and discuss the adsorption mechanism of $\mathrm{CO}_2$ in Irmof-1.

For this task, use the simulation in Irmof-1 at $298\,\mathrm{K}$ and investigate the pore filling mechanism with the snapshots. Where does adsorption begin? As pressure increases, how do the fluid–fluid and fluid–framework interactions change? Discuss the shape of the isotherm.
Explore and discuss the influence of temperature on $\mathrm{CO}_2$ adsorption in Irmof-1. How is the isotherm affected by the increase in temperature?
Calculate and discuss the isosteric heat of adsorption as function of loading using the Clausius-Clapeyron equation. What is its sign? Is it expected?
Compare the adsorption properties of Irmof-1 and Irmof-3. Discuss if these materials show any potential as materials for carbon capture.

Going further¶

MuSiC requires only three basic elements to work:

atom files (.atm) containing information about atom species;
molecule files (.mol) containing the structures of the framework you want to study as well as the guest molecules;
force-field information contained in aa.interaction (for atom–atom interactions), mm_sim.interactions (for molecule–molecule interactions) and intra.interactions (for intramolecular interactions).

Atom and structural details can be found in tables and in the literature. There is little room for subjectivity in their choice usually.

The choice of proper force-fields is much more complicated though: there are many of them derived using different methods and focusing on reproducing a more or less limited range of physical properties.

In our case, we modelled the fluid–fluid and fluid–adsorbent interactions using Lennard-Jones potentials… If you consider particle $i$ and particle $j$ at distance $r_{ij}$ from each other then the interaction between $i$ and $j$ is written:

\[U(r_{ij}) = 4 \varepsilon_{ij} \left[ \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{12} - \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{6} \right]\]

That equation is defined by two parameters: $\sigma$ and $\varepsilon$. We chose our parameters from the uff database while the parameters for $\mathrm{CO}_2$ were taken from this reference.

If you have a look at them, these references only provide parameters for atoms of the same kind: for instance, you would get $\sigma_{ii}$ and $\varepsilon_{ii}$ for a C–C pair or an O–O pair but not for a C–O pair. To get those cross-interactions back, we used the usual Lorentz–Berthelot’s combination rule which tends to work fairly well:

\[\sigma_{ij} = \frac{\sigma_{ii}+\sigma_{jj}}{2} \quad \mathrm{and} \quad \varepsilon_{ij} = \sqrt{\varepsilon_{ii} \varepsilon_{jj}}\]

Electrostatic interactions were calculated using Ewald summation method.

The simulations were run for $1500000$ steps where each gcmc step consists of either randomly translating, rotating, inserting or deleting a molecule.