View Single Page Version

A Metal Organic Framework

Setting up a metal organic framework system

60 minutes

Summary

In this example we’ll do some baseline comparison between calculated Bragg scattering from a metal organic framework (MOF) and experimental neutron data measured on the NIMROD instrument in 2013.

We’ll go through the steps of importing a crystal structure from a Crystallographic Information File (CIF) file, and setting up a basic simulation so we can compare the calculated structure against some experimental neutron data. The target framework is Cu-BTC (copper(II)-benzene-1,3,5-tricarboxylate, known also as MOF-199, HKUST-1, or Basolite C300) and we’ll use a structure from the Crystallography Open Database - the corresponding CIF file is also included in the example file data for Dissolve.

Preparation

Download the example data as a zip or tar.gz and unpack it if you haven’t done so already. Alternatively, download the data files listed below and save them to a location of your choice.

CIF file: 7108574.cif
Neutron Data (Empty MOF): Empty-CuBTC.mint01

Save your input file under a sensible name in same directory as the example data

Step 1 - Import the CIF

Our first task is to import the CIF file - Dissolve has a wizard to help you do this, and which also allows some tidying of the imported coordinates, replication of the unit cell to create a supercell etc.

On the top left-hand you will see the Source CIF section, click the icon to open the file selector

Select the relevant CIF file (7108574.cif)

Dissolve will parse the CIF and display some basic information from it such as the authors (Peikert, Hoffmann, and Froba), publication reference, and space group (Fm-3m).

You can change the space group from the drop-down menu.

The viewer will now be showing you the unit cell and its contents as determined from the CIF (you might need to zoom out or rotate the view a bit to see something meaningful). If you open the Structure Cleanup tab, you will have some broad choices on the tolerance for removing overlapping atoms (which arise from defined symmetry generators creating symmetry-identical copies). In the Bonding tab, you have some choices for how connectivity between atoms is calculated (either by Dissolve or from defined bond distances in the CIF) and which structural assemblies / groups (typically arising from disordered moieties in the determined crystal structure) are used.

Basic Cu-BTC unit cell

The default options get us most of the way there, but it’s instructive to take a look at the Assemblies that were detected in the CIF. An assembly is a set of symmetry-unique atoms that have one or more configurational states (i.e. they have different positions) and can represent disorder, different conformers of molecules, or even completely different molecules within the same unit cell. If you look at Peikert, Hoffman, and Froba’s original paper you’ll see that the focus is on an amino derivative of the BTC linkers, and both the original Cu-BTC and the amino-substituted structure are represented in the CIF. The default assembly/group combination of Global/Default, B/Default and A/1 gives the normal Cu-BTC structure, while disabling A/1 and enabling B/2 and C/2 gives you the amino-substituted version.

Before we proceed, make sure that the default Cu-BTC groups (Global/Default, B/Default, and A/1) are the only ones enabled.

Now, if you look closely at the crystal structure you’ll see some oxygen atoms floating around in space…

Dangling oxygens on Cu sites

These are in fact from coordinated water molecules present when the crystal structure was determined - the water hydrogen atoms are not present in the CIF data. We will remove these oxygen atoms since we want a perfectly “dry” unit cell.

Open the Structure Cleanup section

The Structure Cleanup page has several options for cleaning up various aspects of the structure we currently have, for example removing free atoms/ions. There is a specific option for removing water molecules (Remove water and coordinated oxygens) which will also remove terminal oxygen atoms based on the assumption that hydrogen positions were not available in the CIF. Ours are just free atoms, however, so:

Enable the Remove unbound atoms option

You should see those terminal oxygen atoms disappear from the structure, leaving us with a “pure” Cu-BTC framework. There are a couple of pages left in the wizard which allow us to create a supercell from the unit cell,

Click ok to complete the dialog and generate the Cu-BTC species

Our new species not only contains the atoms of the framework, but also the unit cell (or supercell) definition. You will also see that a Configuration has been added which represents the unit cell by referencing the periodic species.

Step 2 - Apply a Forcefield

Time to tell our Cu-BTC how to move! There are lots of forcefields out there for metal organic frameworks, with varying complexity and functional form. Dissolve implements a few of them for specific cases, and we’ll use the one by Zhao et al..

Click on the Cu48O192C288H96 species tab

From the available forcefields choose Zhao2010 and click Next

You can filter forcefields by keywords in name and description by using the filter box at the top-right of the forcefield selection controls.

We will use the default Determine atom types for all atoms option to add atom types for every atom in the species, so click Next

Here you’ll see that there are some conflicts with atom types assigned during the CIF import - we will simply overwrite them with those from the forcefield, so click Next

For the intramolecular terms we want to assign them and reduce to master terms (the default settings) so click Next to proceed

There will be no conflicts with existing master terms, so click Finish to exit the wizard

Our forcefield for the MOF is now ready, but because we’re only treating the unit cell (rather than a supercell) we need to reduce the pair potential range since the default (15.0 Å) is longer than the half-cell length of the Cu-BTC unit cell ($l$ = 26.3336 Å):

Go to the Forcefield tab, Pair Potentials section

In the controls, change the Range from 15.0 to 12.0

Step 3 - Set Up the Simulation

Create a Configuration

When importing the original CIF file a configuration was automatically created, so in this example we’re ready to go. It might be instructive to look at the nodes in the generator, specifcally the Add node for the framework species. If you look at its control parameters you will see that the BoxAction parameter is Set, meaning that the configuration’s box is set from that defined in the species.

Set up Layers

Since our framework is one big periodic molecule there is no point in trying out Monte Carlo translation and rotation moves on it, so we’ll evolve it with just molecular dynamics.

The MD module in this layer runs every step and for a larger number of steps compared to the one in the standard MC/MD layer.

Finally, we need to calculate our correlation functions:

Go to the RDF / Neutron S(Q) tab

Click on the SQ module to display its options

Set the QDelta to 0.01

Set the QBroadening to OmegaDependentGaussian with a FWHM of 0.02

Click on the NeutronSQ module to display its options

For the Reference keyword in the Reference Data settings group, select the file Empty-CuBTC.mint01 and set the format of the data to mint

We’re changing to a finer binning in $Q$ so our Bragg features have better definition.

Add Bragg Scattering

One thing missing from the standard layers is the calculation of Bragg scattering, so we must add a Bragg module in to this layer.

Show the module palette for the current layer by clicking the Show Available Modules button at the bottom of the module list on the left

Drag a Bragg module from the Correlations category up to the module list, placing it in-between the existing GR and SQ modules

Select the new Bragg module to show its options

Change the QMax to 5.0

The Bragg module is responsible only for calculating the reflection intensities from a target configuration - we must tell the SQ module to take those intensities and incorporate them into the structure factors:

Go to the RDF / Neutron S(Q) tab

Click on the SQ module to display its options

Find the Bragg Scattering options group and set the IncludeBragg option to reference the only Bragg module we have (called “Bragg”!)

Step 4 - Assessing the F(Q)

Start the simulation running:

Ctrl-R

After just a few tens of steps you should see that the MD is doing what it’s supposed to and giving our framework a little bit of disorder as it flexes and vibrates:

Cu-BTC at 300 K after several iterations of molecular dynamics

If we take a look at the output of the NeutronSQ we can see how our simulation compares to the experimental data:

Go to the G(r) / Neutron S(Q) layer tab

Select the NeutronSQ module and go to its Output tab

Comparison of simulated (black) and reference (red) total structure factors

We’ve zoomed in to quite a narrow $Q$ range (0 < $Q$ < 2.0 Å^-1) in the above graphic in order to show in detail the Bragg peaks from the sample, and you’ll see that there is pretty good agreement between the simulated and experimental data in terms of positions, but the simulated data are way too sharp and need some broadening applied. The exact degree of broadening to be applied varies from sample to sample and depends on the crystallite size and uniformity, and the resolution of the instrument on which the data were measured. Without a calibrant sample - i.e. a crystalline sample of specific quality, crystallite size, and uniformity - against which to assess instrument resolution, a general broadening must be applied which accounts for all these possible factors.

We’ll make an initial guess as to the right broadening here - first we’ll try to get the peak widths at the lowest $Q$ values a bit better by applying a Q-independent factor:

Go to the RDF / Neutron S(Q) tab

Click on the SQ module to display its options

Find the Bragg Scattering options group and set the BraggQBroadening to GaussianC2 with a FWHM of 0.05 and a FWHM(x) of 0.02

Now run the simulation for twenty steps or so:

Ctrl-F

Enter 20 and click OK

Let’s go back to the NeutronSQ module and see what the comparison looks like now:

First broadening of simulated structure factors

As you can see this looks much better! We can still observe that as we go to Bragg peaks at higher $Q$ values are too intense, so we’ll increase the $Q$-dependent broadening as well:

Set the BraggQBroadening of the SQ module to GaussianC2 with a FWHM of 0.055 and a FWHM(x) of 0.03

Second broadening of simulated structure factors

Not looking too bad considering we haven’t taken account of the NIMROD instrument’s “interesting” peak shape, and are representing our Bragg peaks as simple Gaussians.

Step 5 - Assessing the G(r)

Let;s now turn our attention to the total radial distribution function:

Go to the G(r) / Neutron S(Q) layer tab

Select the NeutronSQ module and go to its Output tab

Click on Total G(r)

You should see something similar to this:

Cu-BTC G(r) at 300 K from simulation (black line) and experiment (red line)

Again, not bad - the forcefield is doing a decent job of reproducing the experimental structure. While many of the peaks are present and occur at the correct distances there are a couple that are slightly off the mark. We can assign those peaks to families of structural correlations (or specific ones if you’re lucky) by cross-referencing the peak positions with either the equilibrium values for forcefield bond terms (for the lowest $r$ values) or the distances between atoms at the termini of angle interactions. For the latter these values must be worked out by hand since the angle terms typically specify equilibrium angles instead than distances, so get your trigonometry sharpened up (or look at the individual g(r) for those atom types and read off the peak distance).

The following table does this for the first few peaks.

Experimental $r$, Å	Simulated $r$, Å	(Possible) Peak Assignment
0.9 - 1.06	0.94	C–H (benzene ring)
1.34	1.33	C–C (benzene ring) and C–O
1.92	1.96	Cu–O and C(2)–H (benzene ring)
2.4	2.39	C(2)–O (angle) and C(3)–C(3) (angle)
2.77	2.83	C(3)–O (proximity)

One note of caution - we are assuming here that the experimental structure is as “clean” as that of the model. Adsorbed molecules or impurities / inconsistencies in the MOF structure may also be responsible for the observed differences between model and reality.

It is prudent to adjust bond distances in the forcefield at this point in order to get better agreement. While this is possible to an extent for small molecular systems, for a periodic framework it is considerably more difficult, but we’ll try anyway! Let’s start with the Cu–O interaction occurring in the experimental at 1.92 Å, and which our forcefield is over-estimating.

Go to the Forcefield tab, Master Terms section

Find the Cu-O bond interaction and change the value of the eq parameter from 1.969 to 1.94

We need to run the simulation long enough to allow the GR module to re-average its data, which is accumulated over five datasets and collected every five iterations.

Ctrl-F

Enter 25 and click OK

Once these steps are completed, the simulation should now show much better agreement with the experimental data at 1.92 Å.

Exercises

Now that we have illustrated the basic set up of a crystalline framework simulation, here are a few additional exercises:

1. Even better G(r) - Further improvements to the total G(r)

Exercise 1 - Better Total G(r)

We saw how to improve agreement between the experiment and simulation by modifying the equilibrium distance of the Cu–O master bond in the forcefield. How much better can we make the agreement?

Here we’ll focus on the peak around 2.8 Å, whose assignment is as follows:

Experimental $r$, Å	Simulated $r$, Å	(Possible) Peak Assignment
2.77	2.83	C(3)–O (proximity)

The peak appears to be associated to the proximity of the C(3) and O atoms in the structure, and which essentially are connected by a torsion interaction C3–C2–C1–O. The central interaction C2–C1 is the bond connecting the carboxylate group to the benzene ring. So how can we improve (reduce) the distance between the C3 and O atoms at the termini of this interaction?

One possibility is to change the equilibrium bond distance of the middle C2–C1 in order to push the C3 and O atoms apart. But what effect will this have on other parts of the structure as represented by the G(r) and, indeed, the F(Q)?