| Property Prediction |
In this tutorial you will:
First you must load a model of the zeolite.
| Select the File/Load Model... command from the menu bar in the Visualizer window. Go to the Cerius2-Models/zeolites directory, select MFI.msi, and click the LOAD button. |
You will next make this into a superlattice.
| Select the Crystal Building item, and click the Crystalline Superlattice button on the resulting panel. |
This generates a 1x1x2 supercell.
| Close the Crystal Building panel. |
2. Sketching a pentane molecule
Next you sketch a pentane molecule in a new model window.
| Click the + button in the Visualizer window to create a new, empty model space. |
| Open the Sketcher panel by selecting the Build/3D-Sketcher command from the Visualizer menu. |
| Select the Draw with Hydrogens option, then select the Sketch with option. Click the left mouse button in the model window five times in a zig-zag fashion to construct the pentane molecule. |
| Click the arrow button at the top left of the Sketcher panel. |
| Click and hold down the CLEAN button until the atoms are no longer moving. |
| Rename the model by going to the Visualizer panel, highlighting the name of the model with the mouse, and typing pentane. Close the Sketcher panel. |
Next you will load in the appropriate forcefield for pentane sorption.
With the forcefield loaded, you can input the parameters for your sorption study.
| In the Run Sorption panel, change the Length of run from 100000 to 25000. |
| Make the MFI model active by clicking the diamond next to it in the Visualizer panel. Now go back to the Run Sorption panel and click the RUN SIMULATION button. |
Note that the number of accepted creation attempts is very low, and the calculation is far from convergence. This problem progressively worsens at higher pressures, for larger molecules, and for mixed sorbate systems.
| Select File/Load Model... from the menu in the Visualizer panel. Go to the Cerius2-Models/zeolites directory and select MFI.msi. Click LOAD. Close the Load Model panel. |
| Go to the BUILDERS 1 deck and select the CRYSTAL BUILDER card. Select the Crystal Building item and click the Crystalline Superlattice button in the resulting panel. Close the Crystal Building panel. |
Next you assign atom types to the crystal lattice, using the default cbmc.frc forcefield file.
| Go to the SORPTION TOOLS deck and select the FLEXISORB card. Select the Forcefield/Typing item to open the panel and click the Calculate Atom Types button. |
Now you must check that atom types have been assigned.
| Use the Atom Labels popup in the Visualizer panel to change the atom labeling from NO LABEL to FFTYPE. |
The Si atoms should be assigned sx and the O assigned oz. Now you need to set the atom types for the pentane sorbate.
Now you will add the pentane sorbate to the sorbate list.
| Select the Run item from the FLEXISORB card and use the pulldown selection tool in the Sorbate List to select the pentane sorbate molecule. |
7. Creating the potential energy maps
| Make sure that the Task is set to Create Energy Maps, and Set the File Prefix to MFI_C5. |
| Make the MFI_1 model the current model. Change the FFTYPE atom labeling back to NO LABEL using the popup in the Visualizer panel. Now go to the Run Flexisorb panel and click the RUN button. |
Once the calculation has completed, you can view the output.
| Select the Analysis/View Output item from the FLEXISORB card and select the MFI_C5.log file. |
Now you will load in one of the potential energy grids using the Analysis isosurface commands.
| Select the Analysis/Surfaces item from the FLEXISORB card and set the Isosurface Value to 10. Now click the Create New Surface button. |
| Delete the isosurface by clicking the Delete Surface radio button in the Isosurfaces panel. |
9. Map properties onto surfaces
| Select the Geometry/Free Volumes... item from the Visualizer menu bar, and calculate a surface for the zeolite by clicking the CALCULATE button in the Free Volume panel. |
| Select the surface named CAVITY from the Edit Surface list in the Isosurfaces panel. Click the Property Map... button to bring up the Property Map panel. |
| Go back to the Isosurfaces panel and click the Delete Surface button. |
10. Perform the ideal gas calculation
| In the Run Flexisorb panel, toggle the Task from Create Energy Maps to Ideal Gas. |
Now you will set up the simulation parameters.
In more realistic simulations much larger numbers of steps should be run to ensure convergence.
| Click the RUN button in the Run Flexisorb panel. |
While the calculation is running, you can monitor its progress using the Job Control item on the FLEXISORB card. Click the UPDATE button periodically to check whether a calculation has completed.
Statistics will not begin to be gathered until this number of steps has been reached.
| Select the Analysis/Cloud Points item from the FLEXISORB card, and load the MFI_C5_1.grd grid file. Click the Create Cloud Points button to display the data cloud showing the sorbate distributions. |
| When you are finished, click the Delete Could Points button. |
| Select the File/Load Model... command from the Visualizer menu bar. Change the File Format in the resulting panel from MSI to CAR, select the MFI_C5_4.car file, and click LOAD. |
| Go to the BUILDERS 1 card deck and select the CRYSTAL BUILDER card. Select the Visualization item. On the resulting panel, change the Crystal Cell Display Range to 2 2 3 and click the ENTER button. |
This completes the tutorial lesson.
General Methodology
Preparing a model
The first step in preparing a model is to sketch, or load from a file, one or more sorbate molecules. If you desire a united-atom model for the sorbate, then you may need to remove the hydrogen atoms from the sorbate (see below). Once you have done this, you should load the framework structure. Some commonly used sorbates and framework structures can be found in the Cerius2-Models directories: /sorbates, /catalysts, and /zeolites. Before you can run a calculation, you must reduce the host lattice in symmetry to triclinic, P1, symmetry. You can accomplish this using the Crystal Superlattice command in the CRYSTAL BUILDER card (found in the BUILDERS 1 card deck).
| The sorbate molecules must not contain rings, and the host framework must be orthorhombic. |
Selecting a forcefield
You must assign forcefield parameters to the host lattice and the sorbates. First you must select a forcefield, using the Forcefield/Select item on the FLEXISORB card (found in the SORPTION TOOLS card deck). The default forcefield selection is cbmc.frc. This forcefield is based on the familiar cvff.frc forcefield, but contains additional parameters for simulating alkanes in zeolites. If you wish, you may select other CVFF format forcefields using the browser (click the SELECT button), though their ability to simulate sorption data may be limited.
| This calculation does not use OFF. Selecting forcefields in the OFF SETUP card deck has no bearing on Flexisorb forcefield calculations. |
Atom typing
Once you select an appropriate forcefield, you must assign atom types to the molecules. Selecting the Forcefield/Typing item from the FLEXISORB card opens the Atom Typing panel, which allows you to calculate atom types, or to have them be automatically determined when you launch calculation. If you select the Perform Automatic Typing option when you have manually set any atom types, these manual settings will be overridden (if they do not coincide).
| No charges are used in the simulations. |
United-atom versus all-atom
The cbmc.frc forcefield contains some united atom parameters (ch1, ch2, ch3, ch4) for the interaction of hydrocarbons with zeolites. Using these potentials versus all atom potentials will yield significantly shorter calculation and convergence times. These potentials are automatically assigned if you have constructed a united-atom sorbate molecule (that is, a sorbate molecule with the hydrogens stripped).
| It is best not to mix united-atom and all-atom models within the same simulation. |
Calculating energy maps
The Flexisorb program uses energy maps in conjunction with lookup and interpolation for computational expediency, rather than calculating the energies explicitly. Prior to running any calculation, you must first generate these energy maps.
Job control
You can monitor the progress of the calculations by selecting the FLEXISORB/Job Control item and using the options in the resulting Flexisorb Job Control panel. The default behavior is to run all calculations in the background on the local machine. Since many files required for the calculation are binary, it is not currently possible to launch calculations on remote hosts.
Calculating gas phase chemical potential
The Grand Canonical Monte Carlo method relies upon determination of the equilibrium between the chemical potential for a sorbate molecule in the gas phase, and its chemical potential within a host lattice. You must determine the chemical potential in the gas phase before running any further calculations.
Using existing map files
This calculation uses the energy information stored in the previously computed map files. The program will try to determine the map files to use based upon the names of the host lattice and the atom types. However, if the host name is not the same as that you used in the generation of the map files, you can explicitly specify the maps using the Energy maps panel (accessed by selecting the Flexisorb/Run item and clicking Input... on the panel, then clicking Map Files... on the second panel). You should be careful when explicitly specifying the maps, that they were generated for the same system. Information stored within the maps, such as the atom types and forcefield, can be queried as an additional check, using the INFO buttons on the Selected Maps list within the Energy maps panel.
Calculating the isosteric heat
The setup for the calculation of the isosteric heat of adsorption is similar to that for the Ideal Gas calculation. You must first ensure that the Task is set to Isosteric Heat, and that the sorbates are listed in the exact same order used for the Ideal Gas computation. In the panel accessed by clicking the Input... button (Isoteric Heat), you should verify that the appropriate Ideal Gas File is selected, or select it using the browser. Set the Temperature (K) to the same value as that used for the Ideal Gas simulation.
Predicting the uptake isotherm
To begin this calculation, you must once again select the Flexisorb/Run item from the FLEXISORB card, then set the Task (on the Run Flexisorb panel) to Uptake Isotherm. The calculation of the uptake isotherm requires that you know the temperature and pressure of the system. You can set these parameters, together with other control parameters, using the commands in the panel accessed by clicking the Input... button. As with the calculation of the isosteric heat (see Calculating the isosteric heat), the Ideal Gas File must be specified on the Uptake Isotherm panel, and be consistent with the current simulation. If you plan to select multiple sorbate molecules, you must also specify the mole fraction of each, from which the partial pressures of each sorbate can be calculated. The total mole fraction must sum to unity.
The GCMC calculation attempts four different moves: translation, insertion, deletion, and re-growth. You can modify the relative ratios of these moves to try to improve convergence of the system. However, the total probabilities must always sum to unity. In certain cases, it may be a good idea to bias the insertion attempts to low energy regions of the system, since there may be a greater probability of acceptance. However, at higher pressures this may be a disadvantage, since the low energy sites may already be occupied. It may also be undesirable in certain frameworks, where a low energy sites may be physically inaccessible to a larger molecule. Therefore, you should exercise care in using this feature.
Non-Ideal gases
Many gases exhibit non-ideal behavior, particularly at higher pressures. At low pressures, the observed behavior is often close to ideal. This deviation from ideality is described using the equation of state of the gas, via the BWR formalism. You can load a data file for each sorbate in question by clicking the BWR Files... button on the Uptake Isotherm panel. Several files are provided, and you may select others using the file browser. If this data is unavailable, ideal gas behavior is assumed.
Output files
In additional to an output file, <file_prefix>.log, the uptake isotherm calculation may produce several other output files, for use in analysis. Files containing conformational snapshots of the system may be written. Number density grids giving information on the preferential siting of the sorbates may also be generated. You can prevent the creation of these files using the commands in the panel opened by clicking the Output... button on the Run Flexisorb panel. Table files containing details of statistics and energies are also be produced by default. You may plot these using the panel accessed by going to the GRAPHS card in the TABLES & GRAPHS card deck, and selecting the File/Import Insight TBL item.
Analysis
Depending on the type of calculation you perform, several output files may be produced. The Analysis pullright on the FLEXISORB card contains items that allow you to view the output log files and visualize the generated grid data.
You can delete an isosurface by selecting it from the Edit Surface list and clicking the Delete Surface button.
| Once an isosurface is colored by property, you cannot subsequently recolor it unless you recreate the isosurface. |
Creating a slice plane
A slice plane is a two-dimensional property map. The data is colored by value, using a spectral range, at points on a user-defined plane. You can create multiple slice planes, using the panel accessed by clicking the Analysis/Slices item. You must first load an appropriate grid file by clicking the Files... button, selecting the desired file from the Grid Files panel, and clicking LOAD. You should designate a Position and Direction for the slice (in Cartesian coordinates). The Position specifies a point in space through which the plane will pass, while the Direction specifies a plane normal. To calculate the average position of the selected atoms, select a group of atoms, then click the radio button next to the Position parameter. Similarly, you can calculate the best fit plane for a selected group of atoms by clicking the radio button next to the Direction parameter. You can then create the slice plane by clicking the Create New Slice button on the Slices panel.
Displaying a cloud point diagram
The cloud point diagram displays grid data points colored by their data value. If a grid is very dense, the data points are automatically sparsed. To view a cloud point diagram, you should first select the Analysis/Cloud Points item on the FLEXISORB card. You must then load a grid file, using the panel accessed when you click the Grid Files... button. You the should specify the Property range to plot and the Color mapping, then click the Create Cloud Points button.
| Data points lying outside the specified range are not plotted. Also, you may only display a single cloud point per model. |
Note that each of the grid display options has its own grid browsing and loading capability, so different display options may be simultaneously used on different grids. You must be careful to select the appropriate grid for the display style of interest. For example, you could load a grid for Isosurfaces, but neglect to load a grid for Cloud Points.
Monte Carlo methods are widely used to compute phase equilibria. Common to all these methods is the insertion and deletion of molecules to relate the chemical potential, a "statistical" property, to "mechanical" properties such as density and pressure. In standard canonical ensemble Monte Carlo, this is accomplished through the use of Widom insertions (Widom, 1963). In grand canonical Monte Carlo (GCMC) simulations, molecules are created and destroyed to maintain a fixed chemical potential, µ. In Gibbs ensemble Monte Carlo (GEMC) (Panangiotopoulos, 1987), molecules are transferred between two systems representing homogenous regions of the different phases.
Theory
General CB-GCMC Acceptance Rules
Here we derive the form of the Monte Carlo acceptance rules that remove the bias introduced by preferentially generating low-energy conformations. Let 
mn be the one-step transition probability of going from state m to state n in a Monte Carlo simulation. Let 
mn be the probability of attempting such a move, and let 
m and n be the grand canonical ensemble probability distribution functions of states m and n having Nm and Nn molecules, respectively. Microscopic reversibility requires that:
Eq. 14
For the grand canonical ensemble,
Eq. 15
where µ is the chemical potential of the system and 
is the potential energy of the system which is a function of r, the multi-dimensional vector containing the Cartesian coordinates of all the atoms in the Nm molecules. The acceptance rule for moves which guarantees that Eq. 14 is met is:
Eq. 16
If conformations are generated randomly, then 
mn = nm and standard Monte Carlo acceptance rules are recovered. In the Flexisorb code, the form of the attempt probabilities mn and nm are altered such that the total compute time required to take proper thermodynamic averages is minimized. To see how this method works, consider the case of moving between states m and n, where Nn = Nm + 1. Through substitution of Eq. 15, into Eq. 16, the acceptance rule for insertions is:
Eq. 17
and for deletions:
Eq. 18
where 
is the difference in energy between states m and n. Acceptance rules for thermal equilibrations can also be obtained by considering the case when states m and n contain the same number of molecules. The general acceptance rule for these moves is:
Eq. 19
Pressure-explicit forms of the insertion and deletion acceptance rules can be obtained as follows. The chemical potential of a species in the gas phase is the sum of an ideal gas and residual term:
Eq. 20
Note that for mixtures at low to moderate pressures, it is reasonable to assume that the gas phase forms an ideal solution, and thus the chemical potential of each component can be computed from its partial pressure by Eq. 20. Calculation of µRES can be accomplished through integration of an equation of state at the simulation pressure and temperature via the following relation:
Eq. 21
where 
is the fugacity coefficient, Zc is the compressibility factor, and is the molar density of the gas phase. In the Flexisorb module, the empirical (but accurate) Benedict-Webb-Rubin (BWR) equation of state is used to evaluate Zc. If BWR coefficients do not exist for the species under consideration, then the pressure set at the start of a simulation will actually be the fugacity. For sufficiently low pressures, 
1 and µRES can be safely neglected. The ideal gas chemical potential is given by:
Eq. 22
where the first term is a kinetic energy term comprised of known quantities. The second term in Eq. 22 is zero for simple spheres, but is nonzero for molecules having internal degrees of freedom. ZIG is the ideal gas configurational integral of the molecule:
Eq. 23
and 
is the volume associated with the generalized coordinates, q, given by:
Eq. 24
The cbmc.frc forcefield employs a classical flexible model in the limit of infinitely stiff bond lengths, which results in 2M + 1 generalized coordinates per molecule of M atoms. It is shown in Maginn et al. (1995) that for a known sampling distribution, 
(q):
Eq. 25
While 
(q) can be any probability distribution, here it represents the configurational-biased distributions described in later sections. Evaluation of ZIG requires the simulation of a single molecule in the ideal gas phase. Notice that ZIG is only a function of temperature and sorbate identity. Since it is possible to obtain ZIG for a range of temperatures during a single integration, this calculation needs to be performed only once for each sorbate examined. Since the calculation is fast, however, in the present implementation of Flexisorb, ZIG is evaluated for each sorbate and temperature with an independent simulation.
Eq. 26
while for deletions, substitution of Eq. 20 into Eq. 18 gives:
Eq. 27
where f is the fugacity of the gas phase. Eq. 26 and Eq. 27, along with Eq. 19, are the acceptance rules used in the Flexisorb module. One is at liberty to use any form of the attempt probabilities 
mn and nm that optimize the efficiency of the algorithm. Care must be taken, however, to ensure that the form of the attempt probabilities that appear in these equations actually corresponds to the attempt probabilities used in the simulations, so that the biasing is properly removed. A combination of energy and conformational biasing moves are used in the Flexisorb module. These moves exploit the fact that the zeolite lattice is rigid, so certain locations within the pores are more favorable than others. The code also uses the fact that the internal degrees of freedom are much "stiffer" than external degrees of freedom, therefor biasing with these degrees of freedom is more aggressive. Complete details of the form of the attempt probabilities can be found in Macedonia and Maginn (1998).
Following Kiselev and coworkers (1985), the zeolite is modeled as a rigid lattice of oxygen atoms, with lattice atoms taken from the MSI databank. Adsorbate molecules are modeled using either an all-atom (AA) forcefield, or a united-atom (UA) forcefield with Lennard-Jones parameters optimized for zeolite adsorption (Vlugt et al., 1998). The UA forcefield parameters are summarized in Table 2.
| Non-Bond |
 /kB[K]
|
[Å]
| ||
|---|---|---|---|---|
| CH3 | 98.1 | 3.77 | ||
| CH2 | 47.0 | 3.93 | ||
| CH | 12.0 | 4.10 | ||
| CH3-O | 80.0 | 3.60 | ||
| CH2-O | 58.0 | 3.60 | ||
| CH-O | 58.0 | 3.60 | ||
| Bond Angle |
k [kcal rad-2]
|
o [degrees]
| ||
CH -CH2-CH
| 62.1 | 114 | ||
|
CH-CH-CH
| 62.1 | 112 | ||
| Torsion | ao [kcal] | a1 [kcal] | a2 [kcal] | a3 [kcal] |
|
CH-CH2-CH2-CH
| 0 | 0.70544461 | -0.13549350 | 1.57235284 |
The non-bonded interactions are calculated for interaction sites in different molecules or those separated by more than three bonds in the same molecule by the 12-6 Lennard-Jones potential function (Allen and Tildesley, 1987), truncated at a radius of 13.8 Å. Due to the periodicity of the system, tail corrections are not applied. Bond lengths are fixed at a value of 1.54 Å. Bond angle bending interactions are calculated by a harmonic bending potential given by:
Eq. 28
Dihedral angle potentials are treated with the form (Jorgenson et al., 1984) given by:
Eq. 29
| Non-Bond |
/kB[K]
|
 [Å]
|
|---|---|---|
|
CH
| 19.6276 | 3.88 |
| H | 19.1243 | 2.45 |
| O | 114.7422 | 2.86 |
| Bond-Angle |
k [kcal rad-2]
|
o [degrees]
|
| H-C3-H | 39.5 | 107.67 |
| H-C2-H | 39.5 | 106.4 |
|
C-C3-H
| 44.4 | 111.26 |
|
C-C2-H
| 44.4 | 109.91 |
|
C-C1-H
| 44.4 | 109.22 |
|
C-C2-C
| 46.6 | 110.55 |
|
C-C2-C
| 46.6 | 109.8 |
| Torsion |
k [kcal]
|
n-o [degrees]
|
|
[*] - C - C - [*]
| 1.4225 | 3-0 |
As with the UA forcefield, non-bonded interactions are treated with the 12-6 Lennard Jones function, truncated at 13.8 Å with no tail corrections. Bond lengths are fixed at 1.105 Å for C-H and 1.526 Å for C-C. Bond angle bending interactions are described by Eq. 28, while the dihedral angle potential function is given by:
Eq. 30
Unlike Eq. 29, this function does not favor the trans state over the gauche state. Yet, through combination of Eq. 30 and intramolecular non-bond interactions, the appropriate barrier of rotation is ultimately realized.
Prior to the start of a simulation, the simulation box size is automatically computed to ensure that the minimum image condition is not violated. A potential energy map is tabulated on a three-dimensional grid for fast potential interpolation (June et al., 1990). An ensemble of molecular fragments are also generated for use in the construction of molecules during insertion and thermal equilibration steps. In addition, ideal gas configuration integrals (and thus chemical potentials) for each sorbate are evaluated via MC integration. If these pre-calculation items have already been computed, the code does not execute these calculations. Each MC cycle consists of an attempted molecule insertion, deletion, translation, or cut and regrow thermal equilibration move. For good statistics, it is recommended that the system be allowed to equilibrate for anywhere from two million to ten million moves, prior to taking averages over the last two million moves. Based on past experience, however, the biasing is so effective that within 500,000 moves, all but the bulkiest molecules at the highest loadings will be essentially equilibrated. The magnitude of molecule translations is adjusted automatically throughout the simulation such that 50% of translations are accepted.
Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids, Clarendon Press: Oxford (1987).
References