DL_MONTE: a general purpose program for parallel Monte Carlo simulation

Abstract

Monte Carlo (MC) represents a powerful simulation tool that can be usefully applied to calculating thermodynamic data. However, such codes are normally bespoke for a particular problem and not widely applicable. In this paper, we report a new flexible and versatile MC code called DL_MONTE, which builds on the highly successful DL_POLY molecular dynamics code to allow the treatment of polymers, minerals, semiconductors and metals in a range of applications on both workstations and highly parallel supercomputers. In addition, to describe its features, we used a recent work to model the phase diagrams of mixed metal oxide nanoparticles using MgO/MnO as an illustration, adsorption of water at the MgO surface and, finally, the adsorption isotherms of CO₂ in different microporous zeolites. The results demonstrate the flexibility of the methodology and how semi-grand and grand canonical MC can be readily applied.

Keywords::

• Monte Carlo
• parallel architecture
• nanoparticles
• zeolites

Introduction

The academic community covers a diverse range of disciplines, and as new scientific challenges emerge, novel processes, devices and materials need to be developed thus stimulating new software functionality. For example, climate change has influenced research into clean energy materials and carbon dioxide (CO₂) sequestration. Hence, a feature of any computer code is that they should be developed to be flexible and sufficiently generalised that they can be used to tackle problems from different areas. In addition, the growth in computer power has largely been due to vendors increasing the number of cores within a processor and/or the number of processors within a computer. Thus, many of today's supercomputers contain approximately 100,000 cores, while research groups often have access to computers containing several thousand cores.

Monte Carlo (MC) and molecular dynamics (MD) are techniques for sampling configurational space. Both techniques have strengths and weaknesses. The MC method is far more suitable for calculating the thermodynamic properties (e.g. orderings, adsorption, partition coefficients, solvation, self-assembly, agglomeration, melting and crystallisation, and free energy differences). Meanwhile, MD allows the calculation of time-dependent quantities such as atomic/molecular vibrations and transport coefficients. MD algorithms have been developed such that multimillion particle simulations can be executed on parallel computers containing many thousands of cores and a number of open-source packages exist (e.g. LAMMPS [1] and DL_POLY [2]). These MD packages support a wide range of potential models or force fields and include efficient algorithms for calculating long-range electrostatic interactions and, in order to carry out calculations on very large simulation cells, data (memory) are distributed over a large number of computer nodes. The distributed data strategy has been applied to a number of MC simulations. However, the stochastic nature of MC moves inhibits the inclusion of long-range electrostatic interactions and, currently, there is a lack of flexible MC software that can be used to study large complex systems.

The philosophy behind DL_MONTE is to develop a versatile MC package to address these limitations and accelerate simulations on a range of platforms from high-end PCs through to large computer clusters. To date, it has been assumed that both stochastic moves of individual particles and long-range interactions provided by the Ewald sum [3] are essential and this restricts parallelisation schemes.

Program structure

DL_MONTE is designed to be complementary to DL_POLY [4] and the force field set (the functions needed to define the interactions in a molecular system) are almost identical. Consequently, the program can be used to study a wide range of compounds (e.g. molecular systems, ionic materials and minerals, metals and semiconductors). As a huge variety of potential forms are available and DL_MONTE is designed to be adaptable, the program is not supplied with its own force field parameters. However, it contains many of the functions contained in other packages GROMOS,[5] DREIDING [6] and AMBER.[7]

Force fields available include the following:

• Ewald summation

• Two-body (non-bonded and bonded)

• Three-body (non-bonded and bonded)

• Bonded four-body interactions (dihedral and inversion interactions)

The total configuration energy of a molecular system may be written as

where , and are analytical interaction functions representing chemical bonds, valence angles, dihedral angles, inversion angles, pair-body forces and three-body forces, respectively. The first four are regarded by DL_MONTE as intra-molecular interactions and the next three as inter-molecular interactions. The position vectors and refer to the positions of the atoms specifically involved in a given interaction. The numbers , and refer to the total numbers of these respective interactions present in the simulated system, and the indices and uniquely specify an individual interaction of each type. The indices i, j (and k) appearing in the pair-body (and three-body) terms indicate the atoms involved in the interaction. Within the non-bonded interactions, we also include force fields to handle metals via density-dependent functions (Sutton–Chen [8], Gupta [9] and embedded atom method [10]) and semiconductors using Tersoff potentials.[11] The bonded and non-bonded interactions are evaluated in parallel using a replicated data strategy similar to that described in Ref. [12].

The last term in Equation (1) accounts for the long-range electrostatic interactions between the ions. The Ewald sum [3] is an accurate way of summing long-ranged potentials, such as the Coulomb potential, in periodic systems. It is generally regarded as the most appropriate method for simulations containing a few thousand ions (i.e. the size typically used in MC simulations). Nevertheless, the Ewald sum is the most computationally intensive part of the calculation and often takes approximately 75% of the simulation CPU time.

The basic model for a neutral periodic system is a system of charged point ions mutually interacting via the Coulomb potential. Ewald's method replaces a potentially infinite sum in real space by two finite sums: one in real space and one in reciprocal space (in addition, a self-energy correction is also required).[13] Thus, the Ewald method can be evaluated via two steps. First, each ion is effectively neutralised (at long range) by the superposition of a spherical Gaussian cloud of opposite charge centred on the ion. The combined assembly of point ions and Gaussian charges becomes the real space part of the Ewald sum, which is now short ranged and can be evaluated with the non-bonded interactions. The second step is to superimpose a second set of Gaussian charges, this time with the same charges as the original point ions and again centred on the point ions (so nullifying the effect of the first set of Gaussians). The potential due to these Gaussians is obtained from Poisson's equation and is solved as a Fourier series in reciprocal space. The complete Ewald sum requires an additional correction, known as the self-energy correction, which arises from a Gaussian acting on its own site.

For molecular systems, as opposed to systems comprised simply of point ions, additional modifications are necessary to correct for the excluded Coulombic interactions. These are intra-molecular bonded interactions and it is possible to specify these interactions for any given molecular type on input. In the real-space sum, these are simply omitted. In reciprocal space however, the effects of individual Gaussian charges cannot easily be extracted, and the correction is made in real space. It amounts to removing terms corresponding to the potential energy of an ion due to the Gaussian charge on a neighbouring charge m (or vice versa). This correction appears as the final term in the full Ewald formula below.

where N is the number of ions in the system and the same number discounting any excluded (intramolecular) interactions. represents the number of excluded atoms in a given molecule and includes the atomic self-correction. is the simulation cell volume and is a reciprocal lattice vector defined by

where are integers and are the reciprocal space basis vectors. With these definitions, the Ewald expression given above is applicable to general periodic systems, and in practice, its convergence is controlled by three variables: the real-space cut-off, the convergence parameter α and the largest reciprocal space vector used in the reciprocal space sum. As the calculation of the reciprocal space part of the Ewald sum is the most computationally intensive part of the calculation, much of the optimisation has been achieved in this section of the code (as noted above, the real-space summation can be included within the calculation of the non-bonded two-body interactions). The calculation is organised so that the outer loop of the summation is over the reciprocal lattice vectors and the inner summation over the number of atoms. This has two advantages. First, many MC ‘moves’ involve the translation/rotation/insertion of a single atom or a group of atoms in a molecule (such as CO₂). Therefore, the number of reciprocal vectors is much greater than the number of atoms being moved, this allows for parallelisation of the largest vector in the summation. As with many other MC codes, we store the components of the reciprocal space summation to prevent excessive recalculation (the new energy is evaluated by subtracting the old value and adding the new value).[14] However, for larger simulation cells, this can lead to large memory requirements. Thus, a second advantage of placing the reciprocal lattice vectors in the outer loop is that it is possible to distribute these over different nodes of a cluster thus significantly reducing the memory overhead of the calculation.

The DL_MONTE program was used throughout this study, and the calculations were undertaken on either the SCARF cluster or the UK national supercomputer (HECToR). The computational demand for these calculations and especially including the effects of long-range interactions are demanding and, therefore, required a parallel computational strategy to be developed. This may include running several similar calculations on different nodes (i.e. task farming). However, this was found not to be the most beneficial strategy due to the long start-up (equilibration) times. The calculations presented in this paper were performed using between 4 and 128 cores (depending on the number of atoms/molecules in the simulation cell) and used the message passing interface (MPI) interface. In Figure 1, the time taken for two example simulations (these are typical calculations rather than artificial simulations chosen to demonstrate beneficial parallel scaling) as a function of the number of cores has been plotted. From these plots, one can observe that, for even simple calculations such as MgO (512 ions in total and within the canonical ensemble), it is beneficial to use approximately 128 cores. For a larger, more complex example, such as the grand canonical Monte Carlo (GCMC) simulation of CO₂ in zeolite, good scaling is achieved up to approximately 256 cores. The difference in scaling between these two examples is largely due to the size of the simulation (both the number of atoms and the number of reciprocal lattice vectors). Thus, our strategy for parallelisation for MC simulations on Beouwolf style clusters can be applied to a diverse range of simulations (i.e. NVT, NPT, grand canonical, semi-grand canonical and Gibbs ensembles).

Figure 1 The time taken for two representative simulations. The first is for GCMC simulation of CO₂ in zeolite-A (approximately 2400 zeolite ions and 100 CO₂ molecules) and the second for the NPT simulation of bulk MgO (512 ions). Both calculations used the Ewald summation technique. The zeolite simulations were carried out on the UK supercomputing facility, HECToR, a Cray XE6 and the MgO simulation on the IBM Idataplex, Blue Wonder, at Daresbury Laboratory. All times are referenced to those taken on a single core.

Scientific applications

Nanoparticles

We have investigated the thermodynamic properties of ternary nanoparticle mixtures using DL_MONTE. The objective here is to identify the distribution of cations within the nanoparticle and how this effects the thermodynamic properties. In our pilot study, we examined the MgO–MnO system as we have used this previously in studies of bulk MnO–MgO solid solutions and the resulting phase diagram, calculated using exchange MC, is in good agreement with experiment with a consolute temperature of approximately 1150 K.[15] In addition, we have successfully used this model to calculate segregation energies in thin films.[16] The surface of a crystal always provides a different elastic and electrostatic environment from that of the bulk, and so there is a free energy difference between the energy associated with any defect, including an impurity, in the bulk and at the surface of a nanoparticle. Impurities are thus driven to or from the surface. The basis for all our calculations is the well-known MC method modified for nanoparticles of variable composition as described below. All nanoparticles were cubic and constructed from 8 atom units identical to the cubic unit cell (u.c.) of the bulk material (i.e. face centered cubic (FCC)). Calculations were undertaken on MnO–MgO mixtures comprising 1728 ions (6 × 6 × 6), 2744 ions (7 × 7 × 7), 4096 ions (8 × 8 × 8), 5832 ions (9 × 9 × 9), 8000 ions (10 × 10 × 10), 13,824 ions (12 × 12,12) and 21,952 ions (14 × 14 × 14). The Coulombic interactions were calculated directly and the simulation carried out within the NVT ensemble. We have considered only the {1 0 0} termination of the nanoparticle as the {1 0 0} surfaces of both MgO and MnO are considerably more stable than other terminations. The role of both point defects (e.g. vacancies and interstitials) and linear defects (e.g. steps and dislocations) on the thermodynamics properties of the nanoparticles are neglected. We do not report results for smaller nanoparticles as the estimated consolute temperature is less than that at which the exchange of ions in the simulations is very inefficient (approximately 500 K).

In simulations of non-stoichiometric ionic materials and solid solutions using ‘standard’ MC, kinetic barriers prevent sampling the whole of the configurational space as almost always only one cation arrangement – the initial configuration – is sampled. We have described elsewhere MC exchange simulations in which both the atomic configuration and the atomic coordinates of all the atoms are changed.[17] A single cycle allows for one of three possible moves to be chosen at random with equal probability:

• Random displacement of an ion. Allowing random moves of randomly selected atoms takes account of relaxation throughout the nanoparticle and vibrational effects. To determine whether the change is accepted or rejected, the usual Metropolis algorithm is applied. The maximum change in the atomic displacement for each ionic species has a variable , and its magnitude is adjusted automatically during the simulation to maintain an acceptance/rejection ratio of approximately 0.37.

• An exchange of ionic configuration. Changes in ionic configuration are attempted by exchanging the position of a Mn²⁺ ion and a Mg²⁺ ion, both chosen at random. Again, the Metropolis algorithm is used to accept or reject any attempted exchange.

• A test change in identity of a cation. We evaluate the potential energy change which would result if one species, B, were to be converted into another, A. This change in energy is related to the corresponding change in chemical potential [15] by

Thus for MgO–MnO, we evaluate the energy associated with the conversion of a randomly chosen Mg ion to a Mn ion, and the reverse transformation. We emphasise that the change of Mg into Mn, and vice versa, is only a temporary substitution to determine the free energy difference for a given and ; after Mg/Mn is calculated, the MC simulation of the nanoparticle continues with the initial composition. The test-identity change method is used in preference to the grand canonical ensemble as good statistics on can easily be achieved. To calculate the phase diagram, the calculated values of are fitted to the equation:

By integrating Equation (4) with respect to composition, we obtain the variation in free energy with respect to x (x is the mole fraction of a component, in this instance the mole fraction of Mg²⁺) at each temperature.

Each calculation consisted of an initial equilibration of 5 × 10⁸ cycles followed subsequently by 5 × 10⁸ production cycles. The larger clusters (>4096 ions) pose a significant computational challenge, and for this reason, the calculations were performed using 8–128 cores in the parallel replicated data strategy described above.

One of the principle thermodynamic quantities is the enthalpy required to create a mixed nanoparticle from the stoichiometric binary nanoparticles. We call this the enthalpy of forming, which is measured from the simulation energies, and define the quantity by the following reaction:

The enthalpy of forming for nanoparticles ranging in size between 1728 and 21,952 ions is presented in Figure 2 and can be interpreted from the change in bond lengths of the component oxides on the formation of a mixed oxide. The Mg–O bond length in the bulk oxide is approximately 2.10 Å compared with the longer Mn–O bond length of 2.17 Å. The strain of the ‘bonds’ gives rise to a positive enthalpy. However, in the nanoparticles starting with a MgO nanoparticle, when Mn²⁺ ions are initially added the enthalpy of forming is negative because the larger Mn ions progressively replace Mg ions at corners, edges and faces (Figure 3). This is energetically favourable as the Mn ions have weaker bonds and thus lower interface energies. The opposite is true for Mg ions replacing Mn ions in the MnO nanoparticle. The relative importance of the strain energy increases with the size of the nanoparticle (Figure 2).

Figure 2 The enthalpy of forming, ΔH (kJ mol^–1), for different size nanoparticles. All calculations were carried out at 800 K.

Figure 3 (Colour online) Structure of 21,952 ion nanoparticle containing 5% Mn²⁺ ions. The Mn²⁺ ions preferentially decorate corners, edges and faces. Mg²⁺, Mg²⁺ and O^2– are coloured green, purple and red, respectively.

The phase diagram of the ternary nanoparticles relative to those of the binary nanoparticles can be determined using free energy differences given by Equations (4) and (5) (a detailed description of the methodology is given in [15]). Calculated phase diagrams are presented in Figure 4. As anticipated from the enthalpy of forming, the Mn²⁺ ions are much more soluble in the binary MgO nanoparticle than vice versa.

Figure 4 The calculated phase diagram of the 1728 nanoparticle. The solid line is the binodal and the dashed line the spinodal.

DL_MONTE was used to calculate the miscibility of the binary components for a range of temperatures and nanoparticle size. From this information, the consolute temperature (the critical temperature above which the components of a mixture are miscible, T_c) can be determined as a function of nanoparticle size. The variation with the total number of ions is nonlinear (Figure 5). Limitations in the efficiency of the MC method and the range in the nanoparticle size have prevented us from determining whether the curve is sigmoidal as observed in previous studies.

Figure 5 The consolute temperature T_c as a function of nanoparticle size (total number of ions).

Water on MgO{100} surface

The interaction of water with surfaces has received sustained interest from both experimental and theoretical scientists for a considerable period of time due to its importance in biology, geosciences and chemistry. Consequently, there have been numerous experimental studies (e.g. [18,19]) and theoretical studies (e.g. 2021222324).

GCMC is a standard simulation technique for computing the adsorption isotherms of molecules on surfaces. As it is rather well known (see e.g. [13,25]), we give here only a brief summary of GCMC. GCMC provides a way of sampling configurations with a weighting corresponding to the grand canonical ensemble, in which the number of molecules is allowed to fluctuate. Naturally, in the present situation, only the number of H₂O molecules varies. Four kinds of MC moves are used: translations, rotations, insertions and deletions. At each step, a random decision is made as to which of these moves will be attempted, the probabilities of each being fixed. The probability of acceptance or rejection of each trial move is given by the usual Metropolis rule, based on the change of total energy caused by the move. A similar procedure is used for rigid rotations. There are specified probabilities of deciding to attempt insertion or deletion, these two probabilities being equal. The average number of water molecules, , is controlled over a number of simulations by adjusting the B parameter of Adams [26]. In the case of insertion, a single attempt is made to insert an H₂O molecule, and the probability of accepting the insertion is given by the usual formula:

Likewise for deletion, a single attempt is made to delete an H₂O molecule chosen at random, the probability of this move being given by

The parameter B is related to the chemical potential, , by the equation:

where T and V are the temperature and volume of the system, respectively, and is the Boltzmann constant. is the thermal de Broglie wavelength for the water molecule

where h is the Planck's constant and m is the mass of a water molecule. For the conditions that interest us in this work, the spatial probability distribution of H₂O molecules is strongly peaked in the surface region, and the probability of finding molecules away from this region is negligible, so that we can regard essentially all the molecules in the system as being adsorbed on the surface. The simulations are performed in periodic boundary conditions, so that the MgO crystal is represented by a slab, with vacuum layers separating the periodically repeated images of the slabs. If we denote A as the area of one side of a slab in each simulation cell, then the mean number of adsorbed molecules per unit surface area is . Sometimes, it will be convenient to refer to the ‘coverage’, i.e. the mean number of adsorbed molecules per adsorption site. For the purpose of defining this, we choose to regard each surface Mg ion as an adsorption site. If the area per surface Mg ions is denoted by , then the coverage is defined as .

In the present work, we treat the MgO lattice as rigid, with ions on their bulk positions, the lattice parameter (twice the Mg–O nearest-neighbour distance) being a₀ = 4.27 Å at all temperatures. The interaction potential between each H₂O molecule and the ions of the lattice is a sum of atom–atom potentials having the Buckingham form:

The parameters A_ij and of the short-range repulsions and the parameter C_ij of the dispersion interactions are taken from the study of McCarthy et al. [21], who obtained them by fitting to the results of Hartree–Fock calculations on the interaction of an H₂O molecule with a MgO slab, with the Colle–Salvetti scheme used to correct for electron correlation.[27] The charges on the Mg and O ions of the lattice also have the values given by McCarthy et al. The water molecules are treated as rigid, and the interactions between them are described by the SPC model, the molecular geometry also being taken from this model. The cut-off used in the simulations was 10 Å. We conclude this section by noting some technical details. The MgO (1 0 0) surface was orientated parallel to the XY plane and consisted of eight MgO layers. The vacuum gap was approximately 16 Å and the number of surface Mg ions was 144 (72 at each surface and corresponds to a 6 × 6 u.c.). The influence of cell size on the simulation results were checked against a larger cell (an 8 × 8 surface cell with 128 Mg ions at each surface). In addition to insertion/deletion events, the water molecules were allowed to move (translate) and rotate as a rigid body using the methods described in Ref. [13]. The translation, rotation and insertion/deletion steps were undertaken with a probability of 0.20, 0.20 and 0.60, respectively. The magnitude of the translation and rotation were controlled by parameters that were adjusted automatically during the simulation to maintain an acceptance-to-attempted move ratio of 0.37.[13]

We have performed GCMC simulations on the H₂O/MgO system at a series of temperatures T = 500, 400, 350, 300, 275, 250 and 225 K. These temperatures were chosen because it became clear from exploratory simulations that the critical temperature was somewhere between 200 and 300 K. At each temperature, we performed simulations at a rather closely spaced set of values chosen so that the resulting mean values of coverage spanned the range from 0 to approximately − 80 kJ mol^− 1. As an example of the number of steps required, we note that, to achieve this statistical error at 350 K, runs of 1 × 10⁸ steps were needed. However, the required length of the runs needed increases strongly as T is reduced, and at 225 K, up to 5 × 10⁸ steps were needed.

The first result presented is the entire isotherm calculated for 350 K in Figure 6. At low values of (approximately − 72.0 kJ mol^− 1), very few water molecules adsorb to the surface.

Figure 6 The average number of adsorbed water molecules, as a function of the chemical potential, . The calculation of is described in the text. All the calculations were performed at 350 K.

As the chemical potential is increased, there is a rapid increase in the average number of water molecules until a plateau is reached at  kJ mol^− 1. This point corresponds to identical numbers of adsorption sites, as defined above, and water molecules. The water coverage increases slightly with further increases in the chemical potential until  kJ mol^− 1. When the chemical potential is set above this value, water molecules are observed both on the surface and in the vacuum. We present in Figure 7 the energy distribution of water molecules at two points along the isotherm ( and  kJ mol^–1), i.e. the step and plateau regions.

Figure 7 (Colour online) The distribution of water molecule energies determined at  kJ mol^–1 and  kJ mol^− 1. The energies displayed are taken from representative simulations at 275 K.

The distribution of energies of the water molecules is both a function of temperature and cover. For the step region, the distribution of energies for the water molecules forms a well-defined peak. This is due to the structure of the water molecules on the surface. The average number of water molecules, , has not reached a monolayer and the water oxygen atom is able to sit almost above an Mg²⁺ cation with the H atoms pointing towards the O^2– ions (Figure 8(a)). We also note that the peak in energy is approximately 60 kJ mol^− 1, which is similar to the value of 73.3 kJ mol^− 1 determined by energy minimisation.[21] In the plateau region of the isotherm, the distribution of molecule energies both broadens and has an extra peak at approximately 50 kJ mol^− 1 (Figure 7). This is due to a greater diversity in molecular orientations to the surface. In Figure 8(b), water molecules can be observed with both H atoms coordinated to the oxide ions and with only a single H atom coordinated to the oxide ion.

Figure 8 (Colour online) The structure of water molecules at the surface for (a)  = 69 and (b)  = 86. The Mg²⁺ and O^2– ions in the substrate are coloured green and red, respectively. The O and H atoms in the water molecule are coloured blue and white, respectively.

CO₂ adsorption in zeolites

The adsorption and separation of CO₂ from gas streams is an important area of research, and there are many published works, both experimental and computational, on adsorption and diffusion in bulk zeolite structures 28293031323334353637383940 using a variety of potential models.41424344 Surface effects are very important for adsorption and much less work has been carried out on the adsorption and transport of CO₂ at zeolite surfaces, so this is where we chose to focus our work. However, first we needed to determine an accurate potential model. One of the most sensitive experimental techniques for testing adsorption models is the adsorption isotherm; at low pressures, the uptake will depend almost exclusively on the gas–framework interaction, then as the pressure increases the uptake will depend on the relative strength of the gas–gas and gas–framework interactions, and finally at high pressures the potential model at closer intermolecular distances is probed.

Thus, we used the GCMC feature within DL_MONTE to determine the equilibrium number of CO₂ molecules within a given zeolite structure and, by comparison with experimental isotherms, probe the relative strengths of the interactions. Iterative adjustment of the inter-molecular potentials for different zeolites is used to derive a set of potential parameters for application in a range of materials. This is particularly important for zeolite systems as we require flexibility in being able to model both siliceous zeolites and aluminosilicates with different cations, particularly Na⁺ and K⁺. Another important feature is that, for these systems, there are a large number of experimental adsorption isotherms available along with heats of adsorption to test the effectiveness of the model as a function of composition.

The CO₂ potential that we chose to use in this work was the EPM2 model by Harris and Yung.[45] The EPM2 model has been widely used for the simulation of CO₂; it accurately reproduces the liquid–vapour coexistence curve and the critical properties. For the zeolites, we needed a model that included a variety of different cations and also hydroxyl groups for the modelling of surfaces, and was compatible with CO₂ and water to cover all of the simulations that we intended to carry out. We chose the CLAYFF force field,[46] originally parameterised for the simulation of clays, which can be thought of as chemically similar to zeolites because of their aluminosilicate composition. It includes parameters for a range of metals including Li, Na and Ca and includes water and the hydroxyl group. The CLAYFF force field has been used successfully for the simulation of zeolites.[47] These authors used the CLAYFF parameters for the zeolite, and added O–Si–O and Si–O–Si bend terms to accurately reproduce the zeolite structures. Moreover, Kerisit et al. [48] successfully used the EPM2 model with the CLAYFF model in simulations of scCO₂/H₂O mixtures above a forsterite surface.

The zeolite structures were optimised using energy minimisation with the METADISE code.[49] The zeolite structures (Figure 9) minimised with the CLAYFF potentials have comparable lattice parameters to experiment, as listed in Table 1.

Figure 9 (Colour online) Structure of (a) FAU viewed along (1 1 0) and (b) MFI viewed along (0 1 0) and (1 0 0). Grey colour represents silicon and red colour represents oxygen.

Table 1 Simulated zeolite lattice parameters compared to experiment.

	Calculated lattice parameters (Å)			Experimental lattice parameters (Å)
	a	b	c	a	b	c
Siliceous MFI^a	20.047	19.935	13.378	20.051	19.876	13.368
Siliceous FAU^b	24.056	24.056	24.056	24.258	24.258	24.258
Na-FAU^c	24.416	24.412	24.421	24.638	24.638	24.638
K-FAU^d	24.543	24.455	24.424	24.998	24.998	24.998

^a Experimental data from Artioli et al. [50].

^b Experimental data from Hriljac et al. [51].

^c Experimental data from Wise [52].

^d Experimental data from Lim et al. [53].

The MFI structure compares very favourably with experiment, while some of the faujasite structures exhibit slightly contracted lattice parameters. However, there are a couple of possible explanations for this discrepancy. First, the simulated and experimental systems are not exactly comparable, as our Na-FAU simulations have a Si:Al ratio of 1, and only Na⁺ as a cation, whereas the experimental composition is Ca_0.95Mg_0.22Na_1.04K_0.02[Al_3.40Si_8.60O₂₄] · 6H₂O – note that as well as having a different composition, there is also water present in this structure which will affect the lattice parameters; secondly, energy minimisation calculations do not take into account any crystal vibrations and hence neglect temperature effects. The structures of siliceous faujasite (FAU) and silicalite (MPI) are shown in Figure 9.

The adsorption isotherms were generated using the GCMC features in the DL_MONTE code. A different selection procedure was used for the insertion/deletion of molecules from that described above. This allows for the partial pressure to be used as a simulation parameter (the two methods are equivalent). The insertion of a molecule is governed by the selection procedure:

While the selection procedure for deletion was

where P_g is the partial pressure of the gas. The simulations were run using the SCARF cluster on 32 cores. The zeolite framework was held rigid but we allowed the cations to move. All simulations included a period of equilibration followed by at least 4 × 10⁶ steps, until equilibrium was reached. The potential cut-off was 12 Å. CO₂ molecules could be inserted, deleted, rotated and translated. In siliceous zeolites the translation, rotation and insertion/deletion of CO₂ were undertaken with a probability of 0.2, 0.2 and 0.6, respectively. For aluminosilicates, the cations were allowed to move as well, so the probabilities were 0.05 for cation translation, 0.25 for CO₂ translation, 0.2 for CO₂ rotation and 0.5 for CO₂ insertion/deletion. For our model system of faujasite, there are regions of the zeolite structure (the sodalite cages) that CO₂ cannot enter due to the narrow cage windows. We found that during the GCMC simulations, CO₂ was being inserted into these regions, so we introduced dummy atoms into the sodalite cages to prevent any CO₂ from being adsorbed there.

Each isotherm is the result of a series of simulations at different pressures. The number of CO₂ molecules in the simulation cell (NCO₂) at each pressure is plotted to give the isotherm. The heat of adsorption is calculated from five isotherms generated over a narrow temperature range, usually 263, 268, 273, 278 and 283 K (Figure 10(a)) and the data are used to plot adsorption isosteres, which are graphs of log p versus 1/T at constant coverage (Figure 10(b)). The heat of adsorption (isosteric heat) (q), as a function of coverage (Figure 10(c)), can be calculated from the gradient using the Clausius–Clapeyron relationship (Equation (14)).

Figure 10 (Colour online) Adsorption of CO₂ in Na-FAU (a) adsorption isotherms, (b) adsorption isosteres and (c) heats of adsorption.

Development of model

Initially we focused on the potentials for siliceous zeolites and produced isotherms of CO₂ in siliceous MFI and FAU using potentials generated by mixing the Lennard-Jones 12-6 CLAYFF zeolite parameters with the EPM2 CO₂ potentials using the Lorentz–Berthelot mixing rules but we found this to overestimate at low pressure and underestimate at high pressure. To correct the inaccuracy, we altered the parameters for the CO₂–zeolite interaction, while keeping CO₂–CO₂ interactions purely as EPM2. The new CO₂–zeolite interaction mixed the CLAYFF zeolite with CLAYFF parameters for the oxygen of CO₂ and EPM2 with a slightly modified r (Lennard-Jones size parameter) value for C. Figure 11 shows isotherms for CO₂ in silicalite (MFI) using this modified potential, which can be seen to reproduce the experimental isotherms of Yamazaki et al. [54] very well. Using these isotherms, we calculate the heat of adsorption to be 17 kJ/mol which is slightly lower than the value published by Yamazaki et al. of 23.4 kJ/mol; however, their calculation was over the temperature range 273–473 K, while our calculation was between 263 and 283 K, which may account for the different values.

Figure 11 (Colour online) Adsorption isotherms of CO₂ in siliceous MFI, coloured according to temperature. Data points are our simulations, line is experimental data from Yamazaki et al. [54] at 273 K.

For aluminosilicate zeolites, cations are introduced into the structure, which creates a net charge in the system. We therefore altered the partial charges to maintain a neutral cell. For a Si:Al ratio of 1 in Na-FAU, the O_z (oxygen of zeolite) partial charge was changed by 11% from 1.05 to 1.16875 to keep the cell charge neutral. To negate the effect of this change, we ran a series of tests comparing to experimental values at pressures of 5 and 25 bar by Maurin et al. [41] of N = 93 and 107 molecules CO₂/u.c., respectively. Table 2 lists some of the parameter changes that we trialled. The first test was to increase ϵ O_z (O_z ϵ value) by a series of percentages: 5%, 11%, 20% and 25% (Table 2, tests 1–4). We found that a higher ϵ O_z increased the loading at low pressure more than at high pressure, to give a loading that was too low at higher pressure. We then tried these same percentage increases in ϵ O_z, while also increasing r_Al (Al r value) by 0.15 Å (tests 5–8). We chose to increase r_Al because, in the CLAYFF model, Si and Al have the same potential parameters, the only difference being the partial charges. Al³⁺ is a larger ion than Si⁴⁺, so we thought that increasing r_Al would account for this and reduce some of the increased loading at low pressure. The best result with r_Al increased by 0.15 Å was when ϵ O_z increased by 11% (the same percentage that the charge had changed by), although this was still giving too low adsorption at high pressure. We therefore kept these parameters and tried reducing r O_z (O_zr value) by 0.05, 0.1 and 0.15 Å (1.4%, 2.8% and 4.2%) but found that this increased the loading at low pressure, but not as high pressure as we wanted (Table 2, tests 9–12). Finally, we tried reducing r_Na (Na r value) by 10%, 15% and 20%, finding that this increased the loading as we had hoped (tests 13–15). As the subtleties of these parameter changes was hard to discern from just two data points, however, we generated complete isotherms with all of the r_Na parameter changes and calculated the heats of adsorption for each of them.

Table 2 Sensitivity on potential parameters of the number of CO₂ molecules adsorbed at two test pressures, 5 and 25 bar, to compare with experimental data from Maurin et al. [41].

Potential modifications					N, molec. CO2/u.c.
Test	rAl (Å)	ϵ Oz	(Å)	rNa (Å)	P = 5 bar (target = 93)	P = 25 bar (target = 107)
Base	3.706	0.006739	3.5532	2.6378
1		+5%				91.7
2		+11%			88.6	94.9
3		+20%				93.8
4		+25%			92.9	93.7
5	+4%	+5%				93.5
6	+4%	+11%			90.6	93.3
7	+4%	+20%				92.4
8	+4%	+25%			90.5	92.0
9	+4%	+11%	− 1.4%		90.9	103.3
10	+4%	+11%	− 2.8%		91.7	100.7
11	+4%	+11%	− 4.2%		93.8	101.9
12	+4%	+25%	− 4.2%		94.0	102.4
13	+4%	+11%		− 10%	95.8	103.2
14	+4%	+11%		− 15%	96.9	103.6
15	+4%	+11%		− 20%	99.4	105.2

The DL_MONTE method was able to distinguish well between these small changes in potential sets, changing both the isotherms and the slope of the heats of adsorption graph (Figure 12). As the isotherms were all relatively similar and it was hard to choose the best, the final decision on the parameters that we chose to use was based on being able to reproduce the heat of adsorption, providing the isotherm was satisfactory. We therefore chose: ϵ O_z increased by 11%, r_Al increased by 0.15 Å and r_Na decreased by 15% (Table 2 and Figure 11, test 14). The full potential set for siliceous and aluminosilicate zeolites is listed in Table 3.

Figure 12 (Colour online) Heat of adsorption of CO₂ in Na-FAU, Si:Al ratio = 1. Green: test 13, blue: test 14, red: test 15, orange: experimental data from Maurin et al. [41].

Table 3 Lennard-Jones potential parameters and partial charges.

	ϵ (eV)	r(Å)	Charge
Zeolites
Si	7.98817 × 10^− 8	3.7061	2.1
Al	7.98817 × 10^− 8	3.8564	1.575
Ozs	0.006739	3.5532	− 1.05
Oza	0.007501	3.5532	− 1.16875
Na	0.005641	2.2421	1.0
K	0.004336	3.1810	1.0
CO2–CO2
Oc	0.006930	3.4044	− 0.3256
C	0.002421	3.0946	0.6512
CO2-zeolite
Oc	0.006739	3.5532	− 0.3256
C	0.002421	3.3191	0.6512
Three-body terms	(eVÅ^2)	θijk (°)
O–Si–O	12.1	109.47
O–Al–O	12.1	109.47

Note: Z, zeolite; C, CO₂; s,siliceous zeolite; a, aluminosilicate zeolite.

To change the Si:Al ratio, the partial charges must also be changed slightly, so potentials need to be adjusted to match. The technique that we have found to work well is to adjust the O_z charges to make the cell neutral and then adjust the ϵ O_z by the same percentage from the original CLAYFF parameters; all other potential parameters are unchanged from those in Table 3. An isotherm for Na-FAU with a Si:Al ratio of 1.18 is shown in Figure 13, with excellent agreement with an experimental isotherm.[55] The potential parameters have an O_z charge of − 1.15885 and ϵ O_z of 0.007437 eV.

Figure 13 (Colour online) Adsorption isotherm of CO₂ in Na-FAU with Si:Al ratio = 1 (green) and Si:Al ratio = 1.18 (red and blue) at 303.2 K. Data points are our simulations, blue line is represents experimental data from Pillai et al. [55].

Now that we have a reliable model, we can use it as a predictive tool to probe the effect of parameters such as composition on adsorption. Figure 13 shows the effect on CO₂ adsorption of changing the Si:Al ratio from 1 to 1.18 in Na-FAU, thereby having 96 and 88 Na⁺/u.c., respectively. This small change in composition has a significant effect on the CO₂ uptake, and our model is sensitive enough to be able to discern these differences. The effect that this change of composition has on the CO₂ uptake also highlights how difficult it can be to compare simulated to experimental data as experimentally the samples may not be ideal. For example, if ion exchange is used to change the cation species in a zeolite, it may not be possible to get complete exchange, thus leaving residual amounts of the other cation species which would change the adsorption characteristics.

In summary, in this study, we have used the GCMC functionality in DL_MONTE to derive a set of potentials that can reproduce adsorption isotherms and heats of adsorption in a variety of zeolites. The technique is sensitive enough to distinguish between subtle changes in the potential parameters, to allow optimisation of a model for CO₂ adsorption in zeolites which can now be used as a predictive tool.

Conclusion

We have introduced a new flexible MC code called DL_MONTE which is now available to the simulation community. It can be readily applied to a wide range of materials, from covalently bonded molecules through insulators, semiconductors to metals because of the extensive range of force fields available. The inclusion of the Ewald summation scheme allows for the explicit and reliable calculation of the electrostatic contribution to the total interaction energy. Further important developments have included, first, that the input is largely compatible with the widely used DL_POLY MD code to both provide comparisons and reduce the time required for starting simulations, and secondly, unusually for MC simulation codes, there is a measure of parallelisation via MPI so that the user can exploit large parallel computer nodes.

An important capability of DL_MONTE is allowing the user to go beyond the usual NVE and NVT ensembles, and in this paper, we have illustrated this by describing three applications where we apply semi-grand and grand canonical ensembles. These include the use of the semi-grand canonical ensemble to model the phase stability of mixed oxide nanoparticles and the GCMC option to investigate an oxide water interface. In the final application, we again use GCMC, in this case to model CO₂ sorption in microporous alumina silicates, but use the approach to improve the quality of the force fields as a pre-cursor to a full investigation.

Our results on mixed MgO/MnO nanoparticles show that the semi-grand ensemble is a highly efficient approach for condensed phase systems, particularly when the atom sizes are similar. We report the phase diagrams and our prediction that there is a marked variation of consolute temperature with temperature and it is only when the nanoparticle has in excess of 25,000 atoms that it converges to bulk behaviour.

The GCMC results of water on MgO surfaces have enabled us to evaluate a more reliable estimate of the critical point and address how this is linked to the adsorption behaviour. In contrast, our work on the sorption of CO₂ in the microporous zeolites enables us to distinguish and rationalise the sorption behaviour as a function of framework structure and composition. Furthermore, the presence of detailed experimental isotherms for certain structures and compositions has allowed us to investigate the sensitivity of the results to the details of the intermolecular interactions and indeed to derive new and more reliable parameters.

In summary, the newly available DL_MONTE MC code now allows the academic community to access and use this highly flexible simulation engine. To obtain a copy of the code, please contact J.A.P. and is available through the web link (http://www.ccp5.ac.uk).

Acknowledgements

The authors acknowledge useful discussions with Professor Mike Gillan (UCL) and Professor Neil Allan (Bristol). Part of this work, funded by NERC, made use of the facilities of HECToR, the UK's national high-performance computing service, which is provided by UoE HPCx Ltd at the University of Edinburgh, Cray Inc. and NAG Ltd, and funded by the Office of Science and Technology through EPSRC's High End Computing Programme. The SCARF computing resources were provided by STFC's e-Science facility. DL_MONTE was developed under the auspices CCP5.

Notes

^1. Email: [email protected]

^2. Email: [email protected]