RASPA: molecular simulation software for adsorption and diffusion in flexible nanoporous materials

Abstract

A new software package, RASPA, for simulating adsorption and diffusion of molecules in flexible nanoporous materials is presented. The code implements the latest state-of-the-art algorithms for molecular dynamics and Monte Carlo (MC) in various ensembles including symplectic/measure-preserving integrators, Ewald summation, configurational-bias MC, continuous fractional component MC, reactive MC and Baker's minimisation. We show example applications of RASPA in computing coexistence properties, adsorption isotherms for single and multiple components, self- and collective diffusivities, reaction systems and visualisation. The software is released under the GNU General Public License.

Keywords::

• molecular simulation
• Monte Carlo
• molecular dynamics
• adsorption
• diffusion
• software

1. Introduction

Molecular sieves are selective, high-capacity adsorbents because of their high intracrystalline surface areas and strong interactions with adsorbates. Molecules of different size generally have different diffusion properties in a given molecular sieve, and molecules can be separated on the basis of their size and structure relative to the size and geometry of the apertures of the sieve. Much progress has been made in understanding the subtle interaction between molecules and the confinement, and much of this understanding comes from computer simulations that are able to analyse the chemistry and physics at the atomistic level.

The two main computational approaches to tackle these systems are (i) quantum mechanical calculations and (ii) force field-based simulations. The first approach is required for studying properties such as bond breakage and formation and is available in many excellent commercial and non-commercial packages. The second approach is useful for studying larger systems and for calculating a wide variety of thermodynamic and dynamic properties. Force field-based approaches include Monte Carlo (MC) simulations, molecular dynamics (MD) simulations and energy minimisations. We introduce here a code, RASPA, that focuses on MC, MD and minimisation of systems described by classical force fields.

The RASPA code was written as a collaboration among Northwestern University (USA), the University of Amsterdam (The Netherlands) and the University Pablo de Olavide (Spain), with recent contributions also from the University of Delft (The Netherlands). The code evolved initially from the post-doc project (2006–2009) of David Dubbeldam at Northwestern University, where the Snurr group had another MUltipurpose SImulation Code (MUSIC), which was written in object-oriented Fortran 90.[1,2] MUSIC provides functionality for performing MD and MC simulations in a number of different ensembles, minimisations and free energy calculations for bulk and adsorbed phases using a variety of force fields, but not for treating flexible adsorbent frameworks and hence RASPA was developed. Version 1.0 of this code has been used internally by the authors and a growing list of collaborators. In this paper, we present version 2.0 available for public use. Its main areas of utility are thermodynamic properties of liquids and gases, and adsorption/diffusion behaviour of adsorbates in crystalline nanoporous materials.

Examples of nanoporous materials are clays, carbon nanotubes, zeolites and metal-organic frameworks (MOFs). MOFs are a relatively new class of materials composed of metal nodes connected by organic linkers. MOFs possess almost unlimited structural variety because of the many combinations of building blocks that can be imagined. The building blocks self-assemble during synthesis into crystalline materials that, after evacuation of the structure, may find applications in adsorption separations, gas storage and catalysis.345 MOFs have crystal structures that exhibit unusual flexibility. An extreme example is the ‘breathing MOF’ MIL-53 that expands or shrinks to admit guest molecules such as CO₂ and water.[6] For simulation of zeolites, it is common practice to keep the positions of the framework atoms fixed, but this assumption is not valid for many large-pore MOFs. New algorithms and a new code were developed to handle these systems.

RASPA is a serial code. A single point of an isotherm can be obtained within hours for a simple system and in days for more complicated systems. MC codes are ideally suited for task-farm parallelism. Here, simulations are independent and are run as batches of serial simulations that differ in temperature, pressure, etc. For example, (assuming no hysteresis), each point of an isotherm can be run independently. Memory requirements of MC codes are modest.

Programs can be written in various ways, but often it is true that the fastest codes are probably the hardest to read, while programs strictly based on readability lack efficiency. RASPA (being a ‘research’ code) chooses the middle-ground and is based on the following ideas:

• Correctness and accuracy. For all techniques and algorithms available in RASPA, we have implemented the ‘best’ ones (in our opinion) available in the literature. For example, RASPA uses configurational-bias Monte-Carlo (CBMC) and continuous fractional component Monte Carlo (CFCMC); it uses the Ewald summation for electrostatics; MD is based on symplectic and measure-preserving integrators.

• Functional design. Examining the source code, one can notice that there are not a large number of files. The program is split up according to its function: ‘grid.c’ contains the code to make and use an energy/force grid of a framework, ‘ewald.c’ handles all the electrostatic calculations, ‘mc_moves.c’ contains all the moves to be used in MC, ‘potentials.c’ contains all the van der Waals potentials, etc.

• Input made easy. The requirements for the input files are kept as minimal as possible. Only for more advanced options are extra commands in the input file needed. Also the format of the input is straightforward. Fugacity coefficients and excess adsorption are automatically computed.

• Integrated simulation environment. The code is built up of many functions and routines which can be easily combined. MD can be used in MC and vice versa. Extension and modification of the code is relatively straightforward.

This article provides an overview of the application areas of the code. For a detailed description of the algorithms themselves and the inner working of MC and MD codes, we refer to Refs.[7,8]

2. Units, input and conventions

A small set of internal units needs to be chosen. A convenient set, which is chosen in DLPOLY [9], RASPA, and many other codes, is the following:

• The unit of length l₀ is chosen as the Ångstrom, i.e.  m

• The unit of time is chosen as the picosecond, i.e.  s

• The unit of mass is chosen as the atomic mass unit, i.e.  kg

• The unit of charge is chosen as the unit of proton charge, i.e.  C.

All other units follow from this choice. For example, one Pascal is in internal units. A pressure input of 10 Pa in the input file is converted to ‘internal units’ by dividing by . Similarly, at output the pressure in internal units is converted to Pa by multiplying by .

RASPA used three generic ‘types’ or ‘groups’ for the particles: (1) ‘framework atoms’, (2) ‘adsorbates’ and (3) ‘cations’. (The classification was done with respect to porous materials, for pure fluids the meaning of ‘adsorbates’ reduces to ‘molecules’.) The advantage is that the different components of the total energy are available and the interactions can be examined (also the energies in the Ewald Fourier part are split [10]). Cations are considered as part of the framework (they are included in the total mass of the framework). Another example is when using thermostats, e.g. in LTA5A, a different thermostat can operate on the framework atoms, the adsorbates and the cations (these all vibrate/move at different length- and time scales). There are no restrictions on the number of molecules or the number of components. This allows for example an adsorption simulation of a mixture of and in LTA5A with Na and Ca ions. In this example, there are four components: two adsorbate components, CO₂ and N₂, and two cation components, Na and Ca. For each component the MC move types and relative attempt frequency can be specified. In this case, CO₂ and N₂ can undergo particle insertion and deletion moves, while Na and Ca only use translation (the cations can be simulated as mobile). A typical input for such a simulation at 298 K, 1 bar (1:3 mixture of CO₂/N₂) looks like:

SimulationType            MonteCarlo

NumberOfCycles            250000

NumberOfInitialisationCycles    100000

PrintEvery 1000

Forcefield               GenericZeolites

ModifyOxgensConnectedToAluminium  yes

Framework 0

FrameworkName LTA5A

RemoveAtomNumberCodeFromLabel yes

UnitCells 1 1 1

ExternalTemperature 298.0

ExternalPressure 10000.0

Component 0 MoleculeName            sodium

  MoleculeDefinition         Cations

  TranslationProbability      1.0

  RandomTranslationProbability  1.0

  ExtraFrameworkMolecule      yes

  CreateNumberOfMolecules    32

Component 1 MoleculeName calcium

  MoleculeDefinition Cations

  TranslationProbability 1.0

  RandomTranslationProbability  1.0

  ExtraFrameworkMolecule yes

  CreateNumberOfMolecules 32

Component 2 MoleculeName CO2

  MoleculeDefinition TraPPE

  MolFraction 0.25

  BlockPockets yes

  BlockPocketsFilename LTA

  TranslationProbability 1.0

  RotationProbability 1.0

  ReinsertionProbability 1.0

  SwapProbability 2.0

  ExtraFrameworkMolecule no

  CreateNumberOfMolecules 0

Component 3 MoleculeName N2

  MoleculeDefinition TraPPE

  MolFraction 0.75

  BlockPockets yes

  BlockPocketsFilename LTA

  TranslationProbability 1.0

  RotationProbability 1.0

  ReinsertionProbability 1.0

  SwapProbability 2.0

  ExtraFrameworkMolecule no

  CreateNumberOfMolecules 0

Numbering of frameworks, components, etc., is based on the C-convention, i.e. starting from zero. The ‘SimulationType’ line sets the simulation type (e.g. MC, MD, transition state theory (TST), minimisation). RASPA ‘cycles’ for MD is the number of integration steps. For MC a cycle is max(20, N)-move-attempts with N being the number of molecules. In one cycle each molecule experiences on average one MC move attempt (either accepted or rejected). More molecules require more MC moves and the use of a ‘cycle’ allows for a specification of the simulation length that is relatively insensitive to the number of molecules. Once a line like ‘Framework 0’ has been read, the lines below it refer to that framework until another line like ‘Framework 1’ is encountered (it is possible to define multiple frameworks per system). Similarly, the MC-move probabilities are set for a specific component (set with the ‘Component [number] MoleculeName [string]’ line). For each component you can specify mole fraction, fugacity coefficient, MC moves, whether the component is an adsorbate or a cation, the initial number of molecules, etc. The MC-move probabilities are appropriately normalised by the code and therefore only have to be given relative to each other. In this example, twice the number of ‘swap’-moves (i.e. insertion/deletion) will be used compared to translation, rotation and reinsertion. If no fugacity coefficients are given in the input, then the Peng–Robinson equation of state will be used to convert pressure to fugacity. Therefore, if a fugacity coefficient of unity is specified, then adsorption will be computed as a function of fugacity instead of pressure. Any specified initial number of molecules will be created at start up using the CBMC algorithm. This avoids the need to create the initial positions by hand.

In addition to this input file, force-field files need to be created. If one uses a generic force field, then a simple ‘Forcefield GenericZeolites’ is sufficient. If one creates force field files in the current directory, then these files are read instead of the generic files. This is convenient for force field fitting where one needs to change the parameters frequently. The first two force field files are ‘force_field_mixing_rules.def’ and ‘force_field.def’. The first is read to construct an initial force field based on the parameters for each atom-type and using mixing rules. The second file allows you to overwrite an interaction pair directly. Both ways of specifying a force field occur in the zeolite and MOF literature. A third file making up the force field is ‘pseudo_atoms.def’ which defines atom properties such as the name, atomic weight, charge and so on.

Next, for each component a definition must be provided. Many molecules have been already been created for general use and for these one can simply specify, e.g. ‘MoleculeDefinition TraPPE’. The atom names from this file are also those used in the files for defining the molecules. These list all the bond, bend and torsion interactions. For defining a flexible framework, the bond, bend and torsions are specified by type and an algorithm searches automatically for all occurrences of these. Details and examples are given in the manual accompanying the source code.

3. Vapour–liquid coexistence

RASPA was initially developed to simulate porous materials. However, it can also be used to model vapour–liquid equilibrium (VLE) as illustrated in this section.

3.1 Coexistence properties

The enthalpy of vapourisation (or heat of vapourisation) is the enthalpy change required to transform a given quantity of a substance from a liquid into a gas at a given pressure. The enthalpy of vapourisation is given without approximation by

(1)

where is the internal energy per molecule, is the pressure and is the volume per molecule. can be conveniently computed in the canonical Gibbs ensemble. In the canonical Gibbs ensemble, the two fluid phases (i.e. vapour and liquid) are explicitly simulated in two separate simulation boxes. Martin and Biddy noted that, for properties such as enthalpy of vapourisation that involve the pressure, it is preferable to calculate the pressure from the vapour phase.[15] The observed error bars on liquid box pressures are quite large in a molecular simulation and an equilibrated Gibbs ensemble simulation has the same pressure in both boxes. The Gibbs ensemble method, which is implemented in RASPA, is ideally suited to compute properties such as the vapour pressure, gas and liquid densities, compressibility, heat of vapourisation, second virial coefficients, boiling and critical points.

3.2 NVT Gibbs ensemble for vapour–liquid equilibrium

The Gibbs ensemble MC simulation technique allows direct simulation of phase equilibria in fluids.[16,17] NVT (also called ‘canonical’) Gibbs ensemble simulations are performed in two separate microscopic regions, each with periodic boundary conditions. The temperature T is held constant, as a well as the total volume V and the total number of particles N. The equilibrium conditions are (i) equal temperature, (ii) equal pressure and (iii) equal chemical potential for each species in the two boxes. The equal temperature in both boxes is imposed via the MC scheme for configurational equilibration (MC moves like translation and rotation). Condition (ii) is enforced by a volume move, and condition (iii) by a particle transfer move. The volume move makes one box larger and the other box smaller, which leads to pressure equilibration. The transfer of particles between the boxes leads to equal chemical potential.

Van der Waals parameters are very difficult to obtain from experiment or from ab initio calculations. However, the VLE curves are very sensitive to, e.g. the strength parameter and size parameter of the Lennard-Jones potentials. Martin and Siepmann developed the transferable potentials for phase equilibria (TraPPE) force field for a large variety of molecules. It includes (and historically started with) united-atom linear and branched alkanes.[12,13] Figure 1 shows simulation data calibrated to experimentally available VLE data. The simulation results of RASPA agree very well with the simulation data of Martin and Siepmann. By fitting sequentially to methane (CH₄), ethane (CH₃), propane (CH₂), isobutane (CH) and neopentane (C), the five atom types in the model can be uniquely fitted. Other data for alkanes can then be used to validate the force field. These types of simulations are also useful for re-optimising force fields for a different cut-off or change from tail-correction to a shifted potential type.[18,19] The chief advantage of VLE-fitted force fields is that over a large range of pressures and temperatures, the density of the fluid is accurately reproduced.

Figure 1 (Colour online) Vapour–liquid coexistence curves for methane, ethane, propane, isobutane and neopentane computed in the Gibbs ensemble. Line, experimental data taken from the NIST database [11]; closed symbols, previous simulation data of Martin and Siepmann [12,13] using the Towhee code [14], open symbols, this work using RASPA. The order of the data from top-to-bottom is the same as the order in the legend.

For pure component systems, the Gibbs phase rule requires that only one intensive variable (usually the temperature) can be independently specified when two phases coexist. The density and pressure are obtained from the simulation. By contrast, for multi-component systems, pressure can be specified in advance, with the total system being considered at constant NpT. The only change necessary is that the volume changes in the two regions are now independent.

4. Adsorption

4.1 Adsorption in the NpT Gibbs ensemble

The fundamental concept in adsorption science is the adsorption isotherm. It is the equilibrium relationship between the quantity of the molecules adsorbed and the pressure or concentration in the bulk fluid phase at constant temperature.[21] The Gibbs ensemble method can be used to compute adsorption isotherms in nanoporous materials.[16,22] One of the boxes contains the framework, while the other box contains the fluid phase (either gas or liquid) that is in equilibrium with the adsorbed phase. For adsorption of a system of n components, the Gibbs phase rule requires that n+1 intensive variables be set, if you consider the adsorbent as an additional component. These n+1 variables are conveniently taken as the temperature, the pressure of the fluid phase, and n–1 mole fractions in the fluid phase. The system is then simulated using the NpT Gibbs ensemble. The fluid-phase box is maintained at constant pressure (and temperature) by applying volume moves. For adsorption in a flexible framework, the adsorbed-phase box is also maintained at constant pressure, but volume moves in a multi component Gibbs ensemble are performed independently for each box. For the simulation of adsorption in a rigid framework, the volume moves on the adsorbed-phase box are switched off; there is no requirement for mechanical equilibrium.[16] The equilibrium constraints are equal temperature in both systems and equal chemical potentials in the bulk and in the interior of the framework (similar to the VLE, the chemical potential equilibrium is enforced by particle swap moves between the boxes).

Figure 2(a) shows adsorption isotherms of xylenes in MIL-47 using the NpT-Gibbs ensemble. Although we use a slightly different model (all-atom OPLS/DREIDING/UFF model,[23]) the results generally agree with the previous simulation data of Castillo et al. [24]. Note that experiments measure ‘excess adsorption’ while simulations calculate ‘absolute adsorption’. Excess adsorption is the number of molecules in the nanopores in excess of the amount that would be present in the pore volume at the equilibrium density of the bulk gas.[25,26] RASPA calculates absolute adsorption and, if the pore-volume fraction is given (can be computed separately), it also calculates the excess adsorption for convenience. With this type of modelling, the experimental results are well reproduced. The snapshots from the simulation (Figure 3) explain the ortho-selectivity: because ortho-xylene is commensurate with the size of the channel, it forms two layers of molecules that stack very efficiently.[24,23] The advantage of the Gibbs adsorption method is that the reservoir is explicitly simulated and hence the conversion from pressure to fugacity is consistently computed with the same force field. This avoids having to find and use an accurate equation of state for the adsorbates. Downsides include having to explicitly simulate the fluid phase (which can be expensive, especially in the liquid phase), and also the computed fugacity coefficient depends on the quality of the chosen force field and representation of the adsorbates.

Figure 2 (Colour online) Pure component isotherms of o-, m- and p-xylene in MIL-47 at 423 K: closed symbols represent experimental data [20], open symbols represent simulation data.

Figure 3 (Colour online) Snapshot of o-xylene in MIL-47 at 433 K, (left) view along the channel, (right) side view with the channel 45° rotated around the channel axis (the line in the centre is an edge of the unit cell). The 1D channels of MIL-47 are about 8.5 Å in diameter, which optimally stacks molecules that are commensurate with this dimension (i.e. o-xylene). Figure courtesy of Ariana Torres Knoop.

4.2 Adsorption in the grand-canonical ensemble

In the limit of low pressure, fugacity and pressure are equal (i.e. the fugacity coefficient is unity). There is, therefore, no need to explicitly simulate the fluid phase. But also if an accurate equation of state is available, or if the fugacity coefficient is known experimentally, or if one is simply interested in adsorption as a function of fugacity, then the reservoir computation is not necessary. In the grand-canonical (GC) ensemble ( ensemble), the chemical potential is imposed at fixed temperature in a fixed volume (determined in this case by the crystallographic definition of the host framework). Insertion and deletion moves are used in the ensemble to equilibrate the system at the fixed value of the chemical potential which is directly related to the fugacity by (where is the reference chemical potential).

Figure 2(b) shows the results for the xylene-MIL-47 system using the GC ensemble. The results of the Gibbs ensemble and grand canonical Monte Carlo (GCMC) simulations agree very well. In the GCMC simulations, the pressure was converted to fugacity using the Peng–Robinson equation of state, which uses the critical temperature, critical pressure and the ‘acentric factor’ that has been tabulated for many compounds. As output it gives whether under these conditions the fluid is a gas or liquid (or metastable), as well as properties such as the fugacity coefficient, the compressibility and the density of the bulk fluid phase. The density is needed to convert absolute adsorption to excess adsorption and vice versa.

In theory and simulation, it is common to plot loading as a function of fugacity because these plots are unaffected by the gas–liquid transition. For validation of the GCMC capabilities of RASPA, we show in Figure 4 adsorption results for benzene in MFI-type zeolite compared to the previous simulation results of Hansen [27] using the BIGMAC code.[28] Both simulations use fugacity and absolute loadings, and the agreement is excellent.

Figure 4 (Colour online) Pure component adsorption isotherms of benzene in MFI at 603  and 703 K. Closed symbols, previous simulation result of Hansen [27] using BIGMAC [28]; open symbols, this work using RASPA.

Simulations of mixture adsorption require the specification of chemical potentials for each component. For gas phase adsorption, an empirical equation of state can be employed. The same treatment is often not readily generalisable to liquid mixtures because of the lack of accurate activity models.[29] A convenient simulation set-up for such systems is to use NpT Gibbs-ensemble simulations using three boxes: (i) the adsorbed phase with the host framework, (ii) the solution phase and (iii) a vapour phase transfer medium. The molecules are not swapped directly between the adsorbed phase and the liquid phase but instead rely on the vapour phase as an intermediate transfer medium.[29]

Measuring mixed-gas adsorption experimentally is difficult. The ideal adsorption solution theory (IAST) of Myers and Prausnitz [30] is often used to estimate the mixture loading from the pure component isotherms. The validity of IAST can be checked using simulations. Figure 5 shows single component isotherms, and results for an equimolar four-component mixture of para-, meta-, ortho-xylene and ethylbenzene in MAF-x8 at 433 K.[23] The IAST prediction is validated with explicit mixture simulations and for this system the IAST is applicable. Another reason to validate IAST is because it is convenient to use IAST as the input for breakthrough simulations.[31]

Figure 5 (Colour online) Xylene separation at 433 K using MAF-x8,[23] (a) single components fitted with dual-site Langmuir–Freundlich isotherms, (b) equimolar mixture simulations and IAST prediction based on single component isotherms.

4.3 Adsorption in the ensemble (flexible frameworks)

The ensemble [35,36] is the natural ensemble to compute adsorption for flexible frameworks. The system is considered as two components, where the chemical potential of component 1 (the guest species) is kept constant (and has variable particle number), while component 2 (the framework) has constant particle number. As in GCMC, only the adsorbed phase is simulated, but now the volume moves are included to hold the pressure constant. For a single component system, it is not possible to vary three intensive variables independently because of the Gibbs–Duhem relation (from which Gibbs' phase rule follows) which relates them. However, for two (or more) species systems, it is possible to derive, rigorously, a statistical ensemble in which , and and are held fixed. For this ensemble, is the chemical potential of the adsorbate and is the fixed number of atoms of the framework (host). This is a hybrid statistical ensemble which has some properties similar to the single species () and () ensembles. In () MC simulation, one carries out (at least) three distinct types of trial procedures [35,36] (i) the conventional configuration change moves, (ii) the change of volume and/or size of the system and (iii) a creation or deletion move.

Figure 6 shows simulated adsorption results of CO₂ in a flexible IRMOF-1 compared to simulations using a rigid structure (the energy-minimised structure with the same force field), and also compared to experimental data. The results for the rigid and flexible model are very similar for this system and in excellent agreement with experimental data. The computation of adsorption in the flexible structure was feasible because the IRMOF-1 structure stays relatively close to its equilibrium structure. The framework motions are efficiently sampled using the MC/MD-hybrid move.[37,38]

Figure 6 (Colour online) CO₂ adsorption in IRMOF-1 showing step-like isotherms [32]. Lines, experimental data [33];, filled coloured points, rigid framework model [34];, open symbols, flexible framework model.[34]

4.4 Efficient algorithms for open ensembles

A system where the number of molecules varies is called an open system. All open-ensemble methods suffer from a major drawback: the insertion and deletion probabilities become vanishingly low at high densities. This problem is particularly severe for long chain molecules. For adsorption simulations, the fraction of successful insertions into the pores becomes too low. To increase the number of successfully inserted molecules, the CBMC technique was developed in the early 1990s.7394041 Instead of generating ideal gas configurations and trying to insert the molecule as a whole, the CBMC method inserts chains part by part, biasing the growth process towards energetically favourable configurations, and therefore significantly reduces overlap with the framework and other particles.

An alternative scheme to remedy the insertion problem is the recently developed CFCMC method of Shi and Maginn 424344. The system is expanded with an additional particle whose interactions with the other atoms in the system are scaled by a parameter , where . Note that only the inter-molecular energy is scaled (not the intra-molecular energy). Many variations on the algorithm are possible. For example can be changed per molecule or per atom. Both methods slowly ‘inflate’ and ‘deflate’ the molecule like a balloon but differently.

CFCMC uses conventional NVT MC for thermalisation (such as translation, rotation and/or MC–MD hybrid moves), but in addition attempts to change of the fractional molecule using . is chosen uniformly between and and scaled to achieve around 50% acceptance. However, many systems show a behaviour where -changes are hard. An additional bias on can be used. This bias will be removed by the acceptance rules. A careful calibration of can make histograms flat and hence can avoid that the system gets stuck in a certain range of . There are three possible outcomes of a change of to : (i) remains between 0 and 1; there is no change in the number of particles, nor in the positions, nor in the intra-molecular energies. Only and the inter-molecular energy have changed. (ii) becomes larger than 1; when exceeds unity, , the current fractional molecule is made fully present (), and a new fractional molecule is randomly inserted with . Shi and Maginn used a methodology where a rigid conformation is chosen from a ‘reservoir’ of ideal gas molecules generated before the simulation. (iii) becomes smaller than 0: when falls below 0, , the current fractional molecule is removed from the system , and a new fractional molecule is chosen with .

RASPA implements both CBMC and CFCMC, but also a combination (named CB/CFCMC) of the two developed by Torres-Knoop et al. [45]. The basic CFCMC algorithm is used with -biasing, but the insertion and deletion moves are performed using configurational biasing. Figure 2(c) shows that smoother curves are obtained by using this method. The method leads to very reliable results. Other implemented methods in RASPA to improve the efficiency of MC simulations are parallel-tempering and mole fraction replica exchange.[46,47]

5. Screening

Continued research and investments in high-performance computing have produced computing platforms that are now fast enough to permit predictive simulations and large-scale screening studies. Simulation (virtual) screening is significantly cheaper than experimental screening and can be used to increase the successful hit rate. A hierarchical or step-wise approach is often used.

• Initial screening (millions of structures). Screening on the basis of properties that can be computed very quickly. Example properties are pore volume, surface area, pore-size distribution, Henry coefficients and heats of adsorption at infinite dilution.

• High throughput screening (hundreds of thousands of structures). In the pressure range of practical interest, the heat of adsorption and loadings are simulated (usually using relatively short runs). This allows a comparison of structures and an elucidation of structure–property relationships.4950515248

• Detailed analysis (tens or hundreds of structures). Detailed analysis of the most promising structures could include simulations of single-component and mixture isotherms, diffusivities and efficiency estimates of the performance in a fixed-bed adsorber using breakthrough simulations.[23,53]

Figure 7 shows two examples of screening. Gómez-Gualdrón et al. [48] investigated physical limits for methane storage and delivery in nanoporous materials, with a focus on whether it is possible to reach a methane deliverable capacity of 315 cm³(STP)/cm³ in line with the adsorption target established by the ARPA-E agency. Using GCMC simulations, methane adsorption and delivery properties were studied in a population of 122,835 hypothetical pcu MOFs and 39 idealised carbon-based porous materials. From the simulation results, an analytical equation was obtained that delimits the necessary material properties to reach specific methane deliverable capacity targets. This high-throughput analysis elucidates how both geometric and chemical properties, such as void fraction, volumetric surface area and heat of adsorption, impact methane deliverable capacity.

Figure 7 (Colour online) Screening results for (a) physical limits for methane storage and delivery in about hundred thousand nanoporous materials,[48] (b) fixed bed performance of para-selective MOFs.[23] Figure (a) courtesy of Diego A. Gómez-Gualdrón.

The second example is a detailed analysis by Torres-Knoop et al. on separation of benzene, toluene, meta-, ortho-, para-xylene and ethylbenzene (BTEX process).[23] Many ortho-xylene selective structures have been found, but finding para-selective structures is much harder. Small pore structures are able to separate para-xylene using ‘sieving’ (the smaller molecules fit into the structure, but the larger molecules are excluded), but these are unable to separate para-xylene from ethylbenzene (same smallest dimension) and are usually diffusion limited. Torres-Knoop et al. studied about 30 structures in full detail and elucidated why some structures are ortho-selective and others are para-selective. Using snapshots, the reason for a selectivity was explained. Snapshots in Figure 3 show that strong ortho-selectively can be obtained by a two-layer molecular stacking in the MIL-47 structure. The ortho-xylene fits in perfectly, while meta-xylene fits less well. Para-xylene and ethylbenzene are too long and are forced to align obliquely. In their work, a para-selective structure was sought. Using the same mechanism, it is then required that the channel dimension are perfectly commensurate with the para-xylene dimensions. The screening found a strongly para-selective structure (MAF-x8). The single component isotherms were computed, the IAST prediction was validated with mixture simulations and the IAST solution was the input for simulating breakthrough curves. These breakthrough simulations have the ‘cycle-time’ as output, i.e. the time needed before needing to start the expensive desorption process. As can be seen from Figure 7, the MAF-x8 would be better than the currently used technology (BaX). The study also revealed that other para-xylene-selective structures, such as MIL-125 and JUC-77, would be diffusion limited (decreasing their performance).

RASPA provides perl scripts to submit jobs for screening purposes. One specifies the list of adsorbates, structures, temperatures, pressure range and number of pressure points, whether fugacity or pressure is used, whether the points are equally spaced in normal or log-scale, etc. The script then generates all the necessary input files and has as output the job scripts needed to submit it to a cluster (with a single command).

6. Reactive Monte Carlo

The RxMC method allows computation of equilibrium properties for chemically reacting or associating fluids.[54,55] The method samples forward and backward reactions using MC moves directly (without going through the transition states). No chemical potentials or chemical potential differences need to be specified for the reaction steps, just the stoichiometry of the reactions. Essentially, the method enhances the GCMC with a ‘forward’ and ‘backward’ reaction step, which ensures that the chemical reaction equilibria between the reactants and the products are maintained.

As an example, the industrially important propene metathesis is described by three equilibrium reactions [57]

• 2C₃H₆ ↔ C₂H₄ + trans-C₄H₈

• 2C₃H₆ ↔ C₂H₄ + cis-C₄H₈

• cis-C₃H₆ ↔ trans-C₄H₈

Only two reactions are independent and need to be included. In addition to the MC moves associated with simulating a chosen ensemble, also ‘reaction’ moves are performed:

• randomly choose a reaction,

• randomly choose whether to do a forward or backward reaction (this determines the ‘reactant’ and ‘product’ molecule types),

• randomly select the reactant molecules and remove them from the system,

• insert the product molecules at random positions,

• accept or reject the reaction step with the appropriate acceptance probability.

Inserting molecules at high densities is difficult and even more so when one needs to insert several molecules at the same time. To overcome these difficulties, Hansen et al. [56] and Jakobtorweihen et al. . Recently, Rosch and Maginn combined the CFCMC method with the RxMC method.[44] For the propene metathesis reactions simulated with reactions 2 and 3, six fractional reaction molecules are added to the system. Each reaction has an associated reaction (between 0 and 1) and MC moves are performed trying to change to . When a forward reaction is performed and, if accepted, is set to . If a backward reaction is performed and, if accepted, is set to . The insertions and deletions are biased in which allows the method to efficiently overcome insertion and deletion difficulties.

Figure 8 shows the results from RASPA using the CFCMC–RxMC method compared to the CBMC–RxMC simulation results of Hansen et al. [56]. The RASPA results are in excellent agreement. Previously, Rosch and Maginn [44] validated their implementation with the results of Hansen et al. and also found excellent agreement.

Figure 8 (Colour online) Selectivity of the propene metathesis reaction system in the temperature range between 300 and 600 K. Closed symbols, previous simulation result of Hansen et al. [56]; open symbols, this work using RASPA.

7. Diffusion

7.1 Molecular dynamics

In many applications of nanoporous materials, the rate of molecular transport inside the pores plays a key role in the overall process. The size of the pores is usually of the same order as the size of the adsorbates. Diffusion properties of guest molecules in nanoporous materials can therefore be quite sensitive to small differences between different host materials. Molecular-level modelling has become a useful tool for gaining a better understanding of diffusion in nanoporous materials. Many different diffusion coefficients can be defined for guest molecules in nanoporous materials, but it is useful to put them into two general classes: transport diffusivities and self-diffusivities. The transport (or Fickian) diffusivity describes the transport of mass and the decay of density fluctuations in the system, while self-diffusion describes the diffusive motion of a single particle. By omitting the thermodynamic contribution in the transport diffusion, the so-called ‘corrected diffusivity’ (also called ‘collective diffusivity’) is obtained. The self-, corrected and transport diffusivities are equal only in the limit of zero loading.

In an equilibrium molecular dynamics simulation, the self-diffusion coefficient of component is computed by taking the slope of the mean-squared displacement (MSD) at long times

(2)

where is the number of molecules of component , is the spatial dimension of the system, is the time and is the centre-of-mass of molecule of component . The order- algorithm for measuring MSDs conveniently and efficiently captures correlations over short, medium and long times.[59] In crystalline materials the MSDs become linear beyond , where is the repeating distance (usually the unit cell distance). Since the MSD accuracy rapidly decreases over increasing times, a good practice is to fit the diffusivities from a few data points after the MSD has become linear.

For a single adsorbed component, the transport diffusion coefficient is given by

(3)

Note that here first the distances are summed and then squared. Hence collective diffusion can be considered a ‘centre-of-mass’ diffusion. The thermodynamic factor is

(4)

where denotes the concentration (adsorbate loading in the framework), and can be obtained from the adsorption isotherm or from the fluctuation formula.[60] The omission of the thermodynamic factor in Equation (3) leads to the ‘corrected diffusivity’ (also called ‘collective diffusivity’) . The concept of collective diffusivity can be extended to multi-component systems using

(5)

where the elements for components and are known as the Onsager elements. Maxwell–Stefan diffusivities are related to the elements of and can be obtained by matrix inversion.[61,62] Equations relating to the Fickian diffusion coefficients can also be derived.[63] From a phenomenological point of view, there are three different approaches to setting up the flux–driving force relationship for diffusion in nanoporous materials under non-equilibrium conditions. The Fickian-, Maxwell–Stefan- and Onsager formulations are strictly equivalent, and all three viewpoints are needed for different purposes.

If the framework is kept fixed, the potential energy surface induced by the framework can be pre-computed.[64,65] Instead of looping over all framework atoms in order to compute the host-adsorbate energy at each time step, one can construct a 3D grid before the simulation and then obtain the energy by interpolation during the simulation. The more points in the grid the higher the accuracy. RASPA implements the triclinic grid interpolation scheme in three dimensions of Lekien and Marsden [45,66]. The algorithm is based on a specific matrix that provides the relationship between the derivatives at the corners of the elements and the coefficients of the tricubic interpolant for this element. The cubic interpolant and its first three derivatives are continuous and consistent. The same grids can, therefore, be used for both MC and MD with no additional energy drift besides the drift due to the integration scheme, i.e. the energy gradients are the exact derivatives of the energy at each point in the element.

RASPA uses symplectic or measure-preserving and reversible integrators. Symplectic integrators tend to preserve qualitative properties of phase space trajectories: trajectories do not cross, and although energy is not exactly conserved, energy fluctuations are bounded. Since symplectic integrators preserve the topological structure of trajectories in phase space, they are more reliable for long-term integration than non-symplectic integrators. The implemented NVE integrator is the symplectic and time reversible integrator for molecules with an arbitrary level of rigidity, developed by Miller et al. [67] based on a novel quaternion scheme. Thermo- and barostats can be combined with the Miller scheme to control the temperature and pressure.[68,70,71] Figure 9 shows that this type of scheme has excellent energy conservation over many nanoseconds. The CO₂ is modelled as rigid and integrated using the quaternion integration scheme of Miller et al. [67]. Separate thermostats are used to thermostat the translation and the rotation of the molecules. The framework–molecule interactions are computed using a grid interpolation scheme.[45,66] The grid spacing was 0.1 Å. A separate grid is used for each van der Waals interaction of the O and the C of the adsorbate molecule, and another grid is used for the real part of the Ewald summation. In the Fourier part of the Ewald summation, the contribution of the rigid framework atoms is pre-computed at the start of the simulation.[10]

Figure 9 (Colour online) Molecular dynamics of CO₂ at 20 molecules per IRMOF-1 unit cell (2 × 2 × 2 system) at room temperature using the quaternion integration scheme of Miller et al. [67]: (a) the individual contributions to the conserved quantity, (b) the instantaneous and average energy drift. The temperature is maintained by using a Nose´–Hoover chain.[68] The molecule–framework interactions are computed using a grid-interpolation scheme [45,66].

Figure 10 shows the self- and collective diffusivities of small gases (H₂, N₂, Ar, CH₄ and CO₂) in IRMOF-1 at 298 K as a function of loading using the force field of Skoulidas and Sholl [69]. The results of RASPA agree quantitatively with previous simulation results of Skoulidas and Sholl. Self-diffusivities can be computed very accurately because it is a single particle property (which can be averaged over all particles). The collective diffusivity is much more difficult to compute because it is a system property. The self-diffusivities are more strongly influenced by correlation effects (kinetic and vacancy correlations) than the collective diffusivities. Correlations between the particles increase with loading.

Figure 10 (Colour online) Simulated diffusivities of small gases in IRMOF-1 at 298 K. Closed symbols, previous simulation result of Skoulidas and Sholl [69]; open symbols, this work using RASPA.

7.2 Dynamically corrected transition state theory

For some systems, the molecules move too slowly and the diffusion coefficients cannot be calculated reliably using MD. An alternative approach is to use transition state theory (TST). In the TST approximation, one computes a rate constant for hopping between states A and B by computing the equilibrium particle flux through the dividing surface. The dividing surface should partition the system into two well-defined states along a reaction coordinate, which describes the progress of the diffusion event from state to state . In many nanoporous materials, the reaction coordinate follows directly from the geometry of the confinement. For example, in Figure 11 the reaction coordinate for methane in LTL (Linde Type L)-type zeolite is shown: the projection of the position on the channel-axis. The location of the dividing surface is denoted by . In ‘dynamically corrected’ TST, one computes the hopping rate over the barrier in two steps [72,73]:

• the relative probability is computed to find a particle at the dividing surface relative to finding it in state ,

  Figure 11 (Colour online) A typical snapshot of a tagged methane particle (green) in LTL-type zeolite restrained to the barrier surface at an average loading of three methane molecules per unit cell (there are two parallel channels per unit cell) at 300 K. Four unit cells each of 7.474 Å in length are shown. The constrictions are caused by the 12-T-membered rings, which form free energy barriers impeding diffusion. The free energy profile in dimensionless units at this average loading is plotted in white, where the reaction coordinate is chosen parallel to the channel direction. If the free energy barriers are high enough, diffusion can be considered a hopping process from minimum to minimum (, , etc).

  

• the average velocity at the top of the barrier is computed as (assuming that the particle velocities follow a Maxwell–Boltzmann distribution), and the probability (dynamical correction) that the system ends up in state is obtained by running short MD trajectories from the dividing surface.

The transmission rate from cage to cage is then given by

(6)

where , is the Boltzmann constant, the temperature, the mass involved in the reaction coordinate and the Helmholtz free energy as a function of . Calculating TST rate constants is therefore equivalent to calculating free energy differences. The exact rate can be recovered by running short MD trajectories from the dividing surface to compute a dynamical correction.[72] The extension to non-zero loading (or to a flexible framework) simply involves sampling these effects ‘in the background’.[73,74] In Figure 11 the free energy profile for methane in LTL is plotted for an average loading of three molecules per unit cell (a unit cell contains two channels). The barrier of this free energy profile is denoted as and corresponds to an entropic constriction of the channel. The dynamic correction is computed from many snapshots with one particle constrained to and particles free (the snapshots are easily sampled using MC). Each snapshot is used to start an MD path with initial velocities sampled from a Maxwell–Boltzmann distribution. The velocity of the barrier particle is pointing towards cage . For all of these snapshots, MD paths are simulated and the flux at the top of the barrier is computed. Figure 12 shows data for methane, ethane and propane in LTL-type zeolite with dcTST compared to MD (for this system both are feasible). It can be seen that the two methods give identical results. The dcTST method, however, is also applicable for slow diffusion ( m) that is (currently) impossible to compute with MD.

Figure 12 (Colour online) Diffusion of methane (C₁), ethane (C₂) and propane (C₃) at 300 K as a function of loading in LTL-type zeolite computed by TST, dcTST and MD.

The dcTST sampling is an example where it is convenient to be able to constrain MC moves. RASPA includes ways to constrain the movement of a component to a line, plane or sub-volume (box, cylinder, etc.). Any attempt to move a particle outside the sub-volume is rejected. Another useful feature is ‘blocking’ volumes that are large enough to contain a molecule, but where that volume is not accessible from the main channel. Yet another common example is a channel system where one would like to always have equal particles for each channel. This can be achieved by creating a different component for each channel, and only allow MC moves to a cylinder that encompasses that channel.

8. Material properties

8.1 Surface area, void fraction and pore-size distribution

Surface area is the most basic property of porous materials. Along with pore volume, surface area has become the main benchmark characterisation method for any porous material. The surface area is usually determined for experimental samples by measuring a nitrogen isotherm at 77 K and then applying the Brunauer–Emmett–Teller (BET) model. Walton and Snurr [75] examined the consistency of the surface areas obtained from the BET model with those calculated geometrically from the crystal structure for several prototypical MOFs with varying pore sizes. Geometric surface areas can be calculated by using a simple MC integration technique in which a nitrogen probe (3.681 Å) molecule is rolled along the surface of the framework.[26,76,77] Walton and Snurr provided compelling evidence for the importance of calculating the BET surface area from the proper region of the adsorption isotherm.[75] Commercial ‘BET’ instruments are typically set to automatically choose a fixed range for BET fitting. The operator must ensure that the range results in consistent BET model parameters.

The pore-size distribution (PSD) can be calculated geometrically in RASPA using the method of Gelb and Gubbins [77,78]. For every point in the void volume, the largest sphere is found that encloses the point but does not overlap with any framework atoms. This yields the cumulative pore volume curve. Let be the volume of the void space ‘coverable’ by spheres of radius r or smaller; a point x is in only if we can construct a sphere of radius that overlaps and does not overlap any substrate atoms. This volume is equivalent to that enclosed by the pore's ‘Connolly surface’. is a monotonically decreasing function of and is easily compared with the ‘cumulative pore volume’ curves often calculated in isotherm-based PSD methods. The derivative is the fraction of volume coverable by spheres of radius but not by spheres of radius and is a direct definition of the pore size distribution. The function can be calculated by a MC volume integration in RASPA.

Knowledge of the density of porous materials is critical for full characterisation of the adsorbents and for fixed-bed adsorber design studies. Talu and Myers [25] proposed a simulation methodology that mimics the experimental procedure. For consistency with experiment, the helium void fraction is determined by probing the framework with a helium molecule using the Widom particle insertion method:

(7)

Widom insertion uses a probe particle that is inserted at random positions to measure the energy required for or obtained by insertion of the particle in the system. Usually a reference temperature of 25°C (298 K) is chosen for the determination of the helium void volume. Computationally, one can also use the limit of the pore size distribution to evaluate the void fraction. The helium void fraction is needed to convert absolute loadings to excess values (or vice versa).

8.2 Thermal and mechanical properties

Structural flexibility is a well-known property of MOFs.[79] For example, MIL-53 exhibits breathing [6] and IRMOFs exhibit negative thermal expansion.[34] RASPA allows a wide variety of flexible models for the framework. Flexible models are needed to obtain properties such as thermal expansion of the framework itself. Thermal and mechanical transport properties are calculated either from equilibrium Green–Kubo relations or by setting up a small non-equilibrium flux across the system. Thermal expansion can be calculated from MD simulations using a flexible framework.

Figure 13 plots the unit cell size of IRMOF-1 as a function of temperature.[34] The structure become smaller with increasing temperature, in contrast to most materials which expand when you heat them. The classical models, quantum mechanical results and experimental data are qualitatively consistent. The quantum results at zero Kelvin depend very much on the used level of theory and basis sets. The negative thermal expansion of IRMOFs is a direct result of the inherent structure of MOFs: linker molecules connected via metal corners. The wiggling of the linkers becomes larger with increasing temperature leading to a reduced ‘projected’ length. Viewed oppositely, when the temperature is lowered the linkers ‘stretch out’. Therefore, negative thermal expansion is very likely a generic property of many MOFs.

Figure 13 (Colour online) Unit cell size of IRMOF-1 as a function of temperature; experimental data, quantum simulations and three classical models. For further details see text and Ref. [34].

9. Zero Kelvin modelling

9.1 Unit cells

To study the shape and size of unit cells, one needs a model for the framework itself. Core-shell models are very suitable for minerals, metal-oxides and zeolites. The core-shell model introduces charged, massless shells around each core. For minimisation of core-shell models, the generalised Hessian matrix contains both the cores and the shells because the shells need to be minimised with respect to the cores, too. In the shell model, the short-range repulsion and van der Waals interactions are taken to act between the shell particles.

In Table 1 we compare our results to those obtained with GULP [80,81] for the minimisation of various zeolites and minerals using Baker's mode-following technique.[83] Using mode-following minimisation, the gradients on the cell and particles can be lowered arbitrarily close to zero and a true minimum energy is obtained (i.e. all positive eigenvalues of the Hessian matrix). The results of the RASPA and GULP codes are identical. The structures are minimised in space group P1 (no symmetry) with a cut-off of 12 Å using the mode-following technique with full unit cell fluctuation, i.e. all cell lengths and angles are allowed to change. The GULP simulations are computed using GULP 3.1 and experimental data are taken from Schröder and Sauer [82]. The RASPA simulations were fully converged to forces smaller than 10 ^− 8 K/Å and 10 ^− 8 K/strain (1 degree Kelvin is  eV).

Table 1 Comparison of RASPA vs GULP for two core-shell models: the model of Catlow and Gale as provided in GULP [80,81] and Schröder and Sauer [82].

Structure	RASPA						GULP
	[Å]	[Å]	[Å]	[]	[]	[]	[Å]	[Å]	[Å]	[]	[]	[]
FAU (Catlow)	24.226	24.226	24.226	90	90	90	24.226	24.226	24.226	90	90	90
FAU (Schröder)	24.631	24.631	24.631	90	90	90	24.632	24.632	24.632	90	90	90
FAU (exp.)	24.26	24.26	24.26	90	90	90
MFI (Catlow)	19.979	19.739	13.320	90.814	90	90	19.980	19.740	13.320	90.813	90	90
MFI (Schröder)	20.425	20.205	13.634	90	90	90	20.425	20.205	13.634	90	90	90
MFI (exp.)	20.11	19.88	13.37	90.7	90	90
CHA (Catlow)	9.198	9.198	9.198	94.752	94.752	94.752	9.198	9.198	9.198	94.751	94.751	94.751
CHA (Schröder)	9.335	9.335	9.335	94.524	94.524	94.524	9.335	9.335	9.335	94.524	94.524	94.524
CHA (exp.)	9.326	9.326	9.326	94.7	94.7	94.7
TON (Catlow)	13.814	17.387	5.002	90	90	90	13.815	17.388	5.002	90	90	90
TON (Schröder)	14.127	17.779	5.155	90	90	90	14.127	17.779	5.155	90	90	90
TON (exp.)	13.86	17.42	5.04	90	90	90
SOD (Catlow)	8.767	8.767	8.767	90	90	90	8.767	8.767	8.767	90	90	90
SOD (Schröder)	8.951	8.951	8.951	90	90	90	8.951	8.951	8.951	90	90	90
SOD (exp.)	8.83	8.83	8.83	90	90	90
Coesite (Catlow)	7.025	12.290	7.115	90	122.485	90	7.026	12.290	7.115	90	122.485	90
Coesite (Schröder)	7.187	12.538	7.256	90	122.803	90	7.187	12.538	7.256	90	122.803	90
Coesite (exp.)	7.14	12.37	7.17	90	120.3	90
-quartz (Catlow)	4.835	4.835	5.346	90	90	120	4.836	4.836	5.346	90	90	120
-quartz (Schröder)	4.988	4.988	5.506	90	90	120	4.988	4.988	5.506	90	90	120
-quartz (exp.)	4.92	4.92	5.41	90	90	120

9.2 Elastic constants

Elastic constants express the degree to which a material possesses elasticity and mechanical stability (Born criteria). The elasticity tensor is the second derivative of the energy with respect to strain , and can be described in terms of fluctuations in the stress tensor [84]

(8)

where the first term on the right is the so-called Born term

(9)

and the second and third terms are the stress-fluctuations term and ideal gas term, respectively. The is the Kronecker's delta, the function is 1 if the variables are equal, and 0 otherwise.

At zero Kelvin, the elastic constants reduce to the Born term minus a ‘relaxation term’ [85]

(10)

Note that the derivative needs to be evaluated at constant zero gradient , which is an algebraic relation between the coordinates at zero temperature. That is, the state before and after a strain is applied must be in a state of zero net force. When more than one particle is present in the system, this requires a ‘relaxation’ of the atoms relative to one another when the system is strained.[85] The zero temperature limit of the stress fluctuation term in Equation (8) is the relaxation term (and the ideal gas term vanishes in this limit). All expressions in Equation (10) are contained in the generalised Hessian matrix, which is the central quantity used in Baker's minimisation scheme. The elastic constants at 0 K can therefore be computed with very high accuracy. In Table 2 we show the elastic constants at 0 K computed from RASPA and GULP for several zeolites and other silicates. The results are identical. Note that RASPA uses the upper triangular matrix for the simulation cell with the direction of the lattice always aligned with the -axis. This means that during the minimisation, the cell does not change orientation. This is convenient when computing elastic constants (which are directional) because the elastic constants are computed along the Cartesian axes (so the crystal should be aligned with these axes).

Table 2 Comparison of elastic constants at zero Kelvin in units of GPa as computed by GULP and RASPA.

	MFI (RASPA)	MFI (GULP)	CHA (RASPA)	CHA (GULP)	Coesite (RASPA)	Coesite (GULP)	-quartz (RASPA)	-quartz (GULP)
	97.71	97.75	124.05	124.07	124.69	124.75	94.55	94.59
	88.99	89.03	122.24	122.26	210.93	210.95	94.55	94.59
	79.36	79.40	119.85	119.88	160.29	160.35	116.04	116.06
	28.65	28.67	17.51	17.52	28.28	28.32	49.97	50.00
	26.27	26.29	17.91	17.92	71.00	71.00	49.97	50.00
	23.09	23.10	18.32	18.32	43.55	43.58	38.06	38.08
	12.09	12.10	57.70	57.71	84.15	84.16	18.43	18.41
	26.49	26.49	56.82	56.84	98.09	98.09	19.67	19.67
	− 2.15	− 2.16	− 5.53	− 5.52	–	–	− 14.48	− 14.49
	–	–	− 5.53	− 5.52	− 34.12	− 34.13	–	–
	–	–	− 5.11	− 5.10	–	–	–	–
	25.44	− 25.45	56.05	56.06	45.24	45.25	19.67	19.67
	− 7.93	− 7.93	− 6.32	− 6.31	–	–	14.48	14.49
	–	–	− 5.48	− 5.47	9.73	9.73	–	–
	–	–	− 7.24	− 7.23	–	–	–	–
	− 5.78	− 5.78	− 8.71	− 8.70	–	–	–	–
	–	–	− 8.10	− 8.09	− 62.03	− 62.02	–	–
	–	–	5.31	− 5.30	–	–	–	–
	–	–	− 2.73	− 2.73	–	–	–	–
	–	–	− 2.69	− 2.69	0.63	0.65	–	–
	− 0.78	− 0.78	− 2.53	− 2.53	–	–	− 14.48	− 14.49

9.3 Approach angles

In catalysis, insight into chirality transfer from catalyst to reactant can be gained by performing constrained minimisations of the reactant–catalyst complex. Because the mechanism of asymmetric induction for epoxidation of olefins by (salen)Mn catalysts is thought to involve steric interactions between the olefin and the catalyst, the direction of olefin approach has been studied by Oxford et al. [86] using hybrid MC simulations combined with classical optimisations. Four main directions of approach to the Mn-oxo moiety have been proposed in the literature (Figure 14(a)); see Ref. [86] and references therein. To examine the approach of 2,2-dimethyl-2H-chromene to the active site of (salen)Mn = O, the potential energy surface for rotation around the dihedral angle of approach was mapped using constrained classical optimisations. In these optimisations using RASPA, the distance between the oxo ligand and C₁ atom was constrained to 2.0 Å, and the angle defined by the manganese atom, oxoligand and C₁ atom was constrained to 122 (see Figure 14(b)) because this geometry is similar to that expected in the transition state.[87] Hard constraints were employed, using the r ²-SHAKE, cos²-SHAKE and the -SHAKE algorithms [88] for the bond, bend and torsion angle constraints, respectively. The minimisation method used guaranteed that the minimum found was a true minimum (all eigenvalues of the Hessian matrix are positive).

Figure 14 (Colour online) Asymmetric induction for epoxidation of olefins by (salen)Mn catalysts. (a) Proposed directions of approach to the active Mn-oxo moiety of (salen)Mn. is the approach angle defined by the midpoint between the oxygen atoms in the salen ligand, the manganese atom, the oxo ligand and the carbon of the reactant forming a bond with the oxo ligand. (b) The bond-, bend- and dihedral constraints of the saddle point. The inset shows the two enantiofaces (Si and Re) of 2,2-dimethyl-2H-chromene. Figure courtesy of G.A.E. Oxford.

Figure 15 shows adsorption energies of 2,2-dimethyl-2H-chromene on the homogeneous catalyst as a function of the approach angle. The approach angle was the dihedral angle defined by the midpoint of the oxygen atoms in the salen catalyst, the manganese atom, the oxo ligand and the C atom of the reactant that would be forming a bond with the oxo ligand in the first transition state to the radical intermediate. The C atom of chromene was chosen assuming that the phenyl group would stabilise the radical on the C atom. The Si and Re enantiofaces of chromene are defined by the chirality of the C atom in the reactant–catalyst complex. The reactant–catalyst complex was optimised at 10 intervals in . The simulations predict the Re enantioface to be the preferred enantioface, in agreement with experiments. The Re enantioface favours the approach from (approximately equivalent to the approach from direction C in Figure 14), while the Si enantioface prefers to approach from direction D, with a minimum in energy at .

Figure 15 (Colour online) Adsorption energies of 2,2-dimethyl-2H-chromene on the homogeneous salen catalyst as a function of the approach angle.

10. Visualisation

It is difficult to determine the connectivity, shape and size of a channel/cage system just by examining the atomic positions of the framework. Early simulation work therefore used visualisations of energy contour plots and 3D density distributions, e.g. for benzene in silicalite,[90] to obtain siting information. The visualization toolkit (VTK) is an open-source, freely available software system for 3D computer graphics, image processing and visualisation.[91] RASPA has a stand-alone utility, written in C++, which visualises output files written by RASPA using VTK. RASPA produces 3D VTK files for visualising channel structures and 3D VTK files of the histograms of molecule positions during adsorption. For mixtures, a 3D histogram is produced for each component. This allows one to check and study ‘segregation’ of molecules.[92] Figures 3 and 11 have been made using RASPA and VTK.

Figure 16(a) shows the MFI-type zeolite. The orthorhombic unit cell has edge lengths  Å,  Å and  Å. The visualisation shows two straight channels, two ‘zig–zag’ channels and four intersections per unit cell. The depicted surface is how a methane molecule would feel the adsorption surface. To visualise molecules inside the structure, the pore walls can be rendered transparent. RASPA generates 3D energy landscapes using the free energy obtained from the Widom insertion method. The simulation cell is divided up into, e.g. 150 × 150 × 150 voxels. The adsorbate is randomly inserted millions of times and the voxels corresponding to the atom positions of the adsorbate are updated with the current Boltzmann weight. The resulting data-set has regions with value , which correspond to overlap with the structure. The ratio of the non-zero values to the total number of voxels is the void fraction. Multiplying by the volume of the unit cell, we can compute the pore volume.

Figure 16 (Colour online) Energy landscape of MFI. The MFI unit cell has edge lengths a = 20.022 Å, b = 19.899 Å and c = 13.383 Å, with cell angles . The MFI pore system (a) consists of straight channels running in the c-direction, which are connected via ‘zig–zag’ channels. About 29% of the structure is void. Colour code: oxygen (red), silicon (yellow). The snapshot (b) and density plot (c) are at 433 K and 100 kPa. Pictures adapted from Ref. [89].

Figure 16(b) shows a snapshot of 2,3-dimethylbutane in MFI. To see the molecules themselves, the framework has to be either cut open or rendered transparent. The combination of the snapshot and the transparent framework allows for an analysis of adsorption sites, molecular positions and orientations, and molecule–molecule correlations. Snapshots are useful to detect differences in adsorption sites of the various species. For example, in this system the linear alkanes predominantly adsorb in the channels while the dibranched molecules adsorb first in the intersections.

Snapshots are very useful, but sometimes one needs to examine a large number of them to start to see a pattern. By keeping track of the atomic positions using 3D histograms the ‘density’ can be visualised. Figure 16(c) is the average of many snapshots. Therefore, the density is very convenient to obtain the siting information at the unit cell level. The picture is made using atomic positions (you could also use the centre of mass position) and, therefore, gives information on the average configuration (position and orientation). Using this type of approach we previously showed that average positions and occupations of the adsorption sites of argon and nitrogen in IRMOF-1 match well with experiments.[93]

11. Conclusions

We have provided an overview of the algorithms that RASPA implements and showed examples of its application in computing coexistence properties, adsorption isotherms for single and multiple components, self- and collective diffusivities and reaction systems. RASPA is provided as source code under the GPL. The login information for the ‘git’-server can be obtained by emailing one of the authors of this manuscript. RASPA is provided without any kind of support or warranty. It should be viewed as an educational ‘research-code’ that could be useful for researchers working in the field.

Acknowledgements

The partition function values for the RxMC propene metathesis were computed by Sayee Prasaad Balaji; Diego A. Gómez-Gualdrón provided Figure 7a; Ariana Torres Knoop provided Figure 3; R. Krishna provided the IAST computation in Figure 5. We would like to thank the following people for their help and input to improve the program and for the very helpful discussions about the algorithms: Sayee Prasaad Balaji, Youn-Sang Bae, Xiaoying Bao, Rocío Bueno Pérez, Nicholas C. Burtch, Tom Caremans, Ana Martín Calvo, Yamil Colon, Juan Manuel Castillo Sanchez, Allison Dickey, Tina Düren, Titus van Erp, Denise Ford, Houston Frost, Rachel Getman, Pritha Ghosh, Elena García Pérez, Gloria Oxford, Sudeep Punnathanam, Almudena Garcia Sanchez, Juan Jose Gutierrez Sevillano, John J. Low, Patrick Merkling, Patrick Ryan, Lev Sarkisov, Ben Sikora, Ariana Torres Knoop, Krista S. Walton, Chris Wilmer, Ozgur Yazaydin and Decai Yu. Very special thanks to Thijs Vlugt.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This material is supported by the Netherlands Research Council for Chemical Sciences (NWO/CW) through a VIDI grant (David Dubbeldam), by the European Research Council through an ERC Starting Grant [grant number ERC-StG-279520] (Sofia Calero), and by the U.S. National Science Foundation Grant [grant number DMR-1308799] (Randall Snurr).

Notes

^1. Email: [email protected]

^2. Email: [email protected]

^3. Email: [email protected]