10
Sorption
The Sorption module simulates adsorption of molecules into porous crystal structures such as zeolites. More than one type of sorbate molecule can be used in the simulation, and pressure and temperature conditions can be varied.
Introduction
Information obtained
Information obtained from sorption simulations includes adsorption isotherms, binding sites, binding energies, loading capacity as a function of pressure, and Henry's constant.
Output forms
Output can be viewed in four forms:
Sections in this chapter
- Using Sorption
- Running a sorption simulation
- Settings for the energy calculation
- Output during the simulation
- Setting up and running the simulation
- Analysis of sorption trajectory files
- Theory
- References
Using Sorption
Sorption is used to simulate sorption of small molecules, called sorbates, into porous 3D frameworks. The frameworks are typically microporous inorganic structures such as zeolites and alumino-phosphates. Characterizing the behavior of these materials has important application in fields of catalysis and separations technology.
The Sorption module has been designed to be used by experimentalists, as well as computational chemists. Output-analysis features, such as automatic calculation and display of isotherms, make it easy to directly compare Sorption data with bench results.
Sorption implements a rapid Monte Carlo statistical mechanics method (Metropolis et al. 1953, Allen and Tildesley 1987); three different types of simulation can be performed:
Running a sorption simulation
The general steps involved in a sorption simulation are outlined below. Detailed procedural directions on how to use the module follow in subsequent sections.
1. Prepare models.
a. Create or load from a file one or more sorbate molecules. The
sorbates must be both nonperiodic and of neutral charge. If
the sorbate is flexible, a trajectory file containing representative
conformations of the sorbate may be generated, using
molecular dynamics or conformational analysis.
b. Create or load from a file the framework model. The framework
model must be a 3D periodic structure of both primitive
(P1) symmetry and neutral charge. If a surface
calculation is desired, a vacuum-slab model must be generated.
- Some commonly used sorbates and framework structures can be found in the Cerius2-Models /sorbates, /catalysts, and /zeolites directories.
2. Specify the simulation type and set up various simulation
options (that is, temperature, length of run, move sizes, move
probabilities, and output).
3. Choose a forcefield and other energy calculation options to be
used during the simulation. A number of forcefields developed
specifically for Sorption energy calculations are provided in the
Cerius2-Resources/FORCE-FIELD directory. For more information
on the sor_yashonath1.01, sor_demontis1.01, sor_
pickett1.01, and watanabe-austin1.01 forcefields, see the Forcefield-Based
Simulations book.
4. Run the simulation.
- First, Sorption generates random configurations by translating and rotating, and, depending on the ensemble type, creating and destroying sorbate molecules in the framework. The maximum size of translation and rotation moves is set through the Sorption Step Sizes control panel. (See the online help for more detailed information on control panels.)
- Next, Sorption either accepts or rejects the configurations on the basis of energy and, in the grand canonical ensemble, pressure. Configurations with lower energy are more likely to be accepted.
- As the simulation proceeds, a large set of configurations is generated. Typically, 105-107 configurations are generated in a simulation.
5. Upon completion of the simulation, configurations saved to a
sorption trajectory file can be analyzed: energy distribution,
average interaction energies, positional distributions, and so on
can be extracted from these trajectories.
General methodology
Settings for the energy calculation
Energy calculation in Sorption is confined to intermolecular or nonbond energies; that is, sorbate-framework energies and sorbate-sorbate energies. Because both the framework and the sorbates are held internally rigid during the simulation, no valence energy terms are required in the sorption energy calculation.
Sorption relies on the Open Force Field (OFF) for loading the forcefield, calculating and assigning forcefield atom types, and, if the Universal forcefield is used, assigning bond order. However, unlike modules such as Minimizer and Dynamics, Sorption doesn't use the energy terms selection or any of the calculation preferences from OFF.
You probably want to be somewhat familiar with OFF before running sorption simulations. For information about OFF, see the Cerius2 Simulation Tools book.
Bump checks
Before energy is calculated, Cerius2 looks for bad contacts between sorbates and the framework. Configurations in which sorbate and framework atoms are very close together have high energies, and conformations with high energies are very likely to be rejected. The bad contact check (sometimes called a bump check) simply saves calculation time by rejecting the high-energy configurations without having to do the energy calculation.
The default value for the bad contact rejection factor is 0.5; that is, the configuration is rejected if any sorbate and framework atoms are closer to each other than half of their van der Waals radii. In general, the rejection factor should be between 0.5 and 0.8. The aim is to find the right balance between needlessly rejecting acceptable configurations and wasting time calculating high-energy configurations.
Excluded volumes
Certain areas within framework structures may be physically inaccessible to a sorbate molecule diffusing through the lattice. Due to the random nature of the insertion method, these areas might still be sampled during a Sorption calculation, leading to an over-estimation of loading. To counter this effect, you can calculate an accessibility grid using the Free Volumes panel, accessed by selecting the Simulation Controls/Excluded Volumes item from the SORPTION card, and clicking the Free Volumes... button. You may then use this grid to prohibit sampling of inaccessible regions.
van der Waals radius
reduction
Before running a sorption simulation, you may want to reduce the van der Waals radii of some elements, particularly atoms with large positive charges where the actual radii are significantly smaller than the van der Waals radii of the neutral atoms (see the Cerius2 Modeling Environment book).
Sorption calculates the van der Waals and Coulomb energy of the sorbate/framework system. Hydrogen bonds are not explicitly calculated, but hydrogen-bonding effects can be included through careful choice of van der Waals forcefield parameters and charges.
van der Waals energy
VDW energy is calculated by summing all pair interactions within a specified sphere, the radius of the sphere being determined by the interaction cutoff distance.
Minimum image convention
The minimum-image convention is used to approximate van der Waals energy in the periodic framework. In the minimum image convention, an atom is considered to interact only with its closest neighbor atoms in a periodic box around it. For more about the minimum image convention, see the Forcefield-Based Simulations book and Allen and Tildesley (1987).
When the minimum-image convention is applied to the evaluation of inter-sorbate energy, the sorbate molecules are treated as whole unit, and not split across minimum-image borders. That is, if the center of a molecule lies within the border, the energy interactions of the whole sorbate molecule are considered, not just those of atoms within the border. This differs from the treatment of framework/sorbate interactions, because sorbates are not considered to repeat periodically in the way that the framework does.
Note
Coulomb energy
Including Coulomb interactions in the sorption-energy expression is optional. The Coulomb energy calculation takes time and, thus, you should check if the Coulomb energy contribution is significant compared to the VDW energy before beginning a long simulation. For weakly charged systems, omitting the Coulomb energy can greatly reduce simulation time with minimal loss of accuracy.
Ewald summation method
Sorbate/framework electrostatics are evaluated using the Ewald summation method. This method accelerates the long-range Coulomb calculation by splitting the slowly converging real-space sum into a quickly converging real-space sum and a reciprocal-space sum. If a long calculation, or series of calculations, is to be performed using Ewald summation, it may be expeditious to pre-calculate the reciprocal space component over a grid. Thereafter, lookup and interpolation will be used rather than explicit calculation. The duration of the calculation must be sufficiently long to justify the time required to calculate the grid. For more information on the Ewald sum, see Forcefield-Based Simulations book.
Note
Setting energy calculation options for a sorption simulation
1. Open the Sorption Energy control panel (see the online help for
more information) by selecting the Energy Calculation item
from the SORPTION card.
Forcefield
2. If the default forcefield and atom typing rules are adequate for
the simulation, check the Use Automatic Force-Field Options
box and skip to step 5, below.
- The default forcefield is usually the Universal 1.02 forcefield. To see which forcefield is set as the default, open the Load Preferences control panel in the Open Force Field module. (Select the Energy Expression/Automate Setup item on the OPEN FORCE FIELD card to open the Automate Energy Setup control panel, then select the Preferences... button to open the Load Preferences control panel.)
3. To use a forcefield other than the default, but have atom types
for that forcefield automatically calculated:
a. Check the Use Automatic Force-Field Options box.
b. Load your choice of forcefield. For more information on
loading forcefields, see the Cerius2 Simulation Tools book.
Atom types
4. To override or alter calculated atom types:
a. Uncheck the Use Automatic Force-Field Options check box.
b. Load your choice of forcefield.
c. Assign and alter forcefield atom types (see the Cerius2 Simulation
Tools book).
Bad contacts
5. If you want, alter the bad contact rejection factor from its
default value of 0.5.
- If there are cations in the framework or strongly positive atoms in the sorbate molecules, you may want to reduce the VDW radii of these elements to values approaching the ionic radii.
- VDW radii can be edited through the Edit Elements control panel accessed by choosing Element Defaults... from the Build pulldown.
van der Waals energy
6. If you want, alter the interaction cutoff from its default value.
Because the interaction cutoff must be consistent with the minimum
image convention, its size is limited by the cell dimensions.
7. If charges have been assigned to the framework and sorbate
models, then you may want to check the Include Coulomb
Energy check box.
Coulomb energy
8. If Coulomb energy is included, set the Ewald parameters for
the electrostatic energy summation:
a. Open the Sorption Ewald control panel (see the online help)
by selecting the Ewald Parameters... button on the Sorption
Energy control panel.
b. The real-space cutoff is the same as the interaction cutoff on
the Sorption Energy control panel.
c. Under the Automatic Parameter Estimation heading, enter
the required accuracy for the Coulomb energy calculation in
kcal/mol.
d. Click the Estimate Ewald Parameters action button to
update both the reciprocal-space cutoff and the Ewald sum
constant to suitable values.
Using Ewald Grids
9. If long calculations involving electrostatic interactions are to be
performed, it may be desirable to pre-compute part of the
Ewald summation and use interpolation. The steps outlined
above should first be followed for specifying the Ewald accuracy.
a. Open the Sorption Ewald Grid control panel by selecting the
Ewald Grid... button on the Sorption Energy control panel.
b. Ensure that the current model is set to the host lattice. Specify
a name for the grid file, then click the Calculate Ewald
data grid action button.
c. To ensure that the grid is used in later calculations, check the
Use Ewald data grid checkbox, and select the grid from
which to read the Ewald data. If the grid already exists from
a previous calculation, step b above can be omitted.
Output during the simulation
There are several ways to monitor and record the progress of the simulation. Before beginning the simulation, use the Sorption Output control panel to set these forms of output.
- Trajectory files can get very large. For example, in a simulation cell of 624 atoms with a maximum of 250 water sorbate molecules run for 10 million configurations with output every 10,000 steps, the trajectory file would be about 11 MB.
- Sorption trajectory files are non-ASCII and as such cannot be viewed with a text editor.
Note
Although control panels cannot be opened or closed during a simulation, the text, model, and graph windows can be resized, layered, shrunk to icons, and expanded out from icons without interrupting the simulation. This allows you to inspect output in all three windows, even if they are overlapping at the start of the simulation.
Setting output options for the simulation
1. Open the Sorption Output control panel by selecting Simulation
Controls/Output from the SORPTION card.
Model window output
2. If you want to see the position of sorbates in the framework
updated in the model window during the course of the simulation,
check the Update Model During Simulation box and
specify a model update frequency. This option is not recommended
for long sorption runs.
3. If you want to see the position of sorbates in the model window
at the conclusion of the simulation, check the Update Model at
End of Simulation box.
- If the model is updated in this way, the sorbates become part of the framework model, and, if another sorption simulation is run on the same framework, you must first delete the old sorbate molecules.
Plotted output
4. If you want progress of the simulation to be sent to the graph
window during the simulation, choose the Use Graphical Simulation
Monitor option and specify an update frequency. Plots
sent to the graph window depend upon the type of simulation
chosen.
Trajectory file output
5. If you want to save the simulation to a trajectory file so that it
can be analyzed later, choose the Write Trajectory File During
the Simulation box and specify the trajectory write frequency.
Remember that trajectory files can get quite large. Also, specify
the root name of the trajectory file name in the Filename Seed
entry box.
Snapshot file output
6. If you want to save sample conformations to a file so that it can
be analyzed later, choose the Write .msi snapshots box and
specify the snapshot frequency. Remember that these files can
get quite large. Also, specify the root name of the snapshot file
name in the Filename Seed entry box.
Note
Text output
7. If you want output sent to the text window during the course
of the simulation, choose the Output Text During the Simulation
box and specify the write frequency for the text. This text
will also be echoed to a log file.
Setting up and running the simulation
Basic run parameters for the simulation (that is, run length, temperature, and pressure) can be found on the Run Sorption control panel. The more advanced options, which rarely need altering, can be found on smaller control panels accessed by selecting the Simulations Controls item.
Framework
The sorbing framework for the simulation must be in the current model (in the model window) at the start of the simulation. It must be a 3D periodic model with primitive symmetry, overall neutral in charge, and of appropriate size; that is, not too small (see Note on page 234).
Surface
If a surface constrained calculation is to be performed, the 2D surface model must first have been converted to a vacuum-slab model. This can be achieved easily using the CRYSTAL BUILDER module.
a. Having generated a 2D surface, via the SURFACE
BUILDER module (see the Builder manual on line), open
the Preferences... subpanel on the Crystal Building panel.
b. Specify the Vacuum Thickness. This should be at least twice
as big as the direct cutoff which will be subsequently used.
c. Specify the surface Orientation.
d. Convert the surface to a 3D vacuum-slab model by clicking
the BUILD CRYSTAL button.
Surface Setup
In a surface constrained calculation, the sorption placements are restricted to lie within a specified distance above a defined surface. The surface height plus short-range cutoff distance should not exceed the width of the vacuum region. Prior to calculation, a check is made to determine if the specified surface orientation appears to be correct, and a warning will be issued if the surface orientation seems inconsistent with the model.
Sorbate
Sorbate molecules are typically small (neutral) molecules. Though they can reside in any model space, to participate in the simulation, they must appear on the Run Sorption control panel.
Move probabilities
Move probabilities control the choice of move and, during a simulation, are used to generate a new trial configuration from the previous configuration.
In the fixed-loading simulation, the only choice of move is between a translation or a rotation.
In the fixed-pressure or Henry's constant simulations, the move can be the creation of a new sorbate, the destruction of an existing one, a rotation, or a translation. Creation and destruction probabilities are always set equal, ensuring microscopic reversibility.
By default, each type of move is given equal probability. Nonetheless, you can bias the choice by editing the probabilities for each move type in the Move Probabilities control panel. For example, if the sorbate is a single atom or a spherically shaped molecule, reducing the rotation probability shortens the simulation with no loss of information.
Maximum step sizes
The size of a random translation or rotation move for a given sorbate is limited by the maximum step size. Typical default starting values are 1 Å for the maximum translation step and 50° for the maximum rotation.
This maximum step size can be a significant parameter in the simulation. If it is high, the risk of the new sorbate making bad contacts with the framework is large, thus generation of valid new configurations may be slow. However, if the maximum step size is low, the simulation takes a long time to generate configurations far from the starting point; in effect, too many similar configurations are generated. Both these cases cause the simulation to take a long time to equilibrate; thus, it is desirable to strike a balance between them.
The maximum step-size values can be held constant throughout a simulation, if you are satisfied that they give an acceptable equilibration period. By default, however, they are used as starting values for the simulation and are automatically rescaled during the run to try to minimize the equilibration period.
Rescaling step sizes
How to determine the correct bias for the step size, such that equilibrium is achieved most quickly? The equilibration rate is conventionally held to be highest when the success rate for the sorbate and move type in question is approximately 30%-50%. By default, the target rate is set for acceptance of half the steps sampled in the simulation.
Three kinds of control parameters can be specified on the Sorption Step Scaling control panel:
A single rescale frequency is set per simulation, but a unique rescale factor and target success rate can be set for each move type for each sorbate.
To set up and run the simulation
Once the sorbate and framework models have been loaded (see "Running a sorption simulation" on page 231), energy options set (see "Setting energy calculation options for a sorption simulation" on page 235), and output options set (see "Setting output options for the simulation" on page 238), you are almost ready to begin the simulation.
1. Make the model that contains the framework use the current
model.
2. Open the Run Sorption control panel (see the online help) from
the Run item on the SORPTION card.
3. Choose the type of simulation from the popup.
4. Enter the simulation temperature in the entry box.
5. Enter the length of run. This is the total number of steps, including
both accepted and rejected configurations.
6. Enter sorbate molecules into the list box by entering the model
space number of each sorbate into the Toggle Model in Sorbate
List entry box.
- To remove a sorbate from the list, simply re-enter the model space number in the entry box.
- If you use multiple sorbates for the fixed loading calculation, the program first attempts to add the requested number of species. If it cannot add these due to energetic considerations, the program issues a warning message and the calculation continues.
- If a fixed pressure simulation, enter the partial pressure (or, for a non-ideal gas, the fugacity) for the sorbate into the list box.
7. If flexible sorbate molecules are to be used, a trajectory file containing
representative conformations must be associated with
the sorbate. This trajectory file may have been generated by
using molecular dynamics in OFF, or by conformation analysis.
To associate the trajectory file with the sorbate, open the Simulation
Controls/Sorbate Conformers panel. Select the desired
sorbate from the Sorbate List, which will be highlighted. Select
a trajectory file from the Conformer Files browser. The selected
sorbate and conformer file will be added to the list. To specify
that the calculation should use the structures from the conformer
file, rather than the sorbate structure, as loaded, click
the Use checkbox adjacent to the conformer file of interest.
- The program will then extract a conformer with random probability for each specified conformer file/sorbate combination. If a conformer file is not provided, the single conformation of the sorbate, as loaded, will be used.
8. Open the Move Probabilities control panel (see the online help)
by selecting the Simulation Controls/Move Probabilities item
from the SORPTION card. By default, all moves are set with
equal probability; if you want to alter the probability, enter new
values in the entry boxes.
9. To set the maximum step size per move, open the Sorption Step
Sizes control panel (see the online help) by choosing Simulation
Controls/Step Sizes from the SORPTION card. The
Henry's constant simulation does not use step-size control; for
this simulation type, skip to Step 11.
- Default step sizes are 1 Å for translations and 50° for rotations. If you want to specify other step sizes, edit the entry boxes.
10. To modify the procedure for step-size rescaling, open the Sorption
Step Scaling control panel (see the online help) by selecting
the Rescale Controls... button on the Sorption Step Sizes control
panel:
a. By default, step sizes are automatically rescaled to help the
simulation reach equilibrium most efficiently. If you want
the maximum step size to remain constant throughout the
simulation, deselect the Automatically Rescale Step Sizes
check box.
b. Specify the frequency of the rescale and, for each sorbate
and move type, both a rescale factor and a target success
rate. The maximum step size is multiplied or divided by this
rescale factor if the requested acceptance rate is not met. The
acceptance rate is based on the number of accepted to total
moves during the previous rescale.
11. If a surface constrained calculation is being run, the surface
region must first be defined.
a. Open the Simulation Controls/Surface Setup panel. Check
the Surface Constrained Calculation.
b. Specify the Surface Height and Surface Orientation.
Ensure that the orientation is consistent with that used to
generate the vacuum-slab representation of the model.
Begin the simulation
12. Return to the Run Sorption control panel and select the Run
Simulation button.
13. During the simulation, you can move and rescale windows,
enabling you to look at the output in the model, graph, and text
windows.
Stopping the simulation
14. The simulation can be stopped before completion by selecting
the Interrupt button in the dialog box that appears during the
run. Choose the Stop Current Process ASAP option. This ends
the simulation at the end of the current step, writes the normal
information summary to the text window, and updates the
model window and trajectory file (if appropriate).
Restarting the simulation
15. If the calculation has not converged over the specified number
of steps, or if it has ended prematurely (perhaps due to a problem
with disk space availability), you can restart or resume the
calculation using data stored in the restart data file. To restart a
calculation, first load the following saved coordinates and status
files back into Cerius2, for use in the Sorption calculation:
- These files contain the coordinates of the framework and sorbates molecules, as well as the status of the calculation. Loading the files before restarting ensures that the framework structure and sorbate models exactly match the models stored in the .msi files.
16. Be sure to add the sorbate models to the sorbate list in the exact
same order as originally, that is, in their correct numerical
sequence, as specified by the file name.
17. When you restart a calculation using the Restart panel, the statistics
from the terminated calculation are read from the .rst file,
and the calculation proceeds from the last stored configuration.
The calculation statistics are written to the restart file at the
same frequency as the textport output.
Analysis of sorption trajectory files
Sorption analysis tools help you analyze the simulation by presenting the trajectory file data in several ways. Five sets of plots can be created:
Each is described on the following pages.
Trajectory file plots
Trajectory file data are plotted as energy vs. step number, and, for the fixed-pressure simulation, as cell loading vs. step number.
These plots present the raw data of the trajectory file; that is, the energy and loading are instantaneous values rather than the cumulative average values that are displayed in other plots.
Plotting the trajectory file in this way is often the first step in simulation analysis. These plots give an overall view of the trajectory file and are particularly helpful in estimating where the simulation equilibrated.
Mass distribution plots
The mass distribution of all sorbate molecules in the analysis range is plotted. The distribution is projected down a particular zone axis onto a periodic cell. You specify the zone axis and the two basis vectors for the plot cell.
Each cell axis is divided evenly into bins, and the number of molecules having centers that fall within each bin square is counted. This value is plotted on the graph.
Energy distribution plots
A histogram of energy distribution for each sorbate is plotted. You specify the maximum energy, minimum energy, and resolution. 
Loading-curve plots
Loading-curve plots combine average cell-loading values from a set of fixed sorption simulations. The framework and sorbate models must be identical in each set. In the example below, three loading curves from sets of simulations are presented on the same graph. Six fixed-pressure simulations are made for each temperature.
Because the same analysis range is chosen for each simulation, be sure to use the range based on the simulation that was slowest to equilibrate.
Mass cloud plots
In mass-cloud analysis, the center of mass of each sorbate molecule in each configuration is displayed as a dot in the model space. The mass cloud can be displayed concurrently with the framework model or alone in the model window. The mass cloud results in a 3D picture that shows the preferred positions of the sorbates in the lattice. The dots can be colored according to sorbate type or configuration energy.
Plotting sorption trajectory files
All graphs, with the exception of mass cloud plots, created in the procedures outlined below can be saved, reloaded, and edited using the Graphs module.
Note
Plotting trajectory files
Plotting trajectory files creates a raw plot of instantaneous energy as a function of configuration and, for the fixed pressure simulation, a plot of instantaneous loading versus configuration.
1. Open the Sorption Analysis control panel (see the online help)
by selecting Analysis from the SORPTION card.
2. Load the trajectory file of the simulation using the browser box
in the Sorption Analysis control panel.
3. Specify the range of steps you want to include in the analysis.
4. Click the Plot Trajectory File button.
Plotting mass distribution of sorbates in framework
This procedure produces a 2D plot in which sorbate positions are projected onto a specified crystal plane.
1. Open the Sorption Analysis control panel by selecting Analysis
from the SORPTION card.
2. Load the trajectory file of the simulation using the browser box
in the Sorption Analysis control panel.
3. Specify the range of steps you want to include in the analysis
(probably excluding the initial pre-equilibrium steps of the simulation).
4. Open the Sorption Mass Plot control panel by choosing the first
Preferences... button of the Sorption Analysis control panel.
5. Choose the axis of projection for the mass distribution.
6. Cerius2 suggests suitable first and second axes vectors based on
the projection vector, but you can override these suggestions by
editing the entry boxes.
7. The default of 20 bins on each axis gives a mass distribution
graph of 21 by 21 points. Increase these bin numbers for a graph
of higher resolution.
8. Select the Plot Mass Distribution button on the Sorption Analysis
control panel to create the plot.
- If the simulation contains more than one sorbate type, a separate plot for each sorbate is sent to the graph window. Please see the online help for more information about the control panels.
Plotting energy distribution of sorbates
This procedure plots a graph showing the number of configurations as a function of energy.
1. Open the Sorption Analysis control panel (see the online help)
by selecting Analysis from the SORPTION card.
2. Load the trajectory file of the simulation using the browser box
in the Sorption Analysis control panel.
3. Specify the range of steps you want to include in the analysis.
4. Click the Plot Energy Distribution button.
- This produces a default energy distribution plot with an energy range of -50 to 50 kcal/mol divided into 100 bins (10 kcal/mol per bin).
5. To rescale the energy graph:
a. Open the Sorption Energy Plot control panel (see the online
help) by choosing the Preferences... button next to the Plot
Energy Distribution button on the Sorption Analysis control
panel.
b. Edit the entry boxes to set more appropriate minimum and
maximum energy limits and bin groupings.
c. Display the rescaled graph by clicking Plot Energy distribution
on the Sorption Analysis control panel.
Plotting a loading curve
Plots of this type are available only if you have a series of fixed-pressure trajectory files. The framework and sorbates must be identical in all of the files; file names must all have the same seed and be numbered in consecutive order starting from one, that is:
<run_name>_1.sor, <run_name>_2.sor, <run_name>_3.sor...
If you have the correct files, but incorrect names, simply rename the files to form an appropriately numbered set (use the UNIX mv command).
1. Open the Sorption Analysis control panel (see the online help)
by selecting Analysis from the SORPTION card.
2. Load any sorption trajectory file in the set using the browser
box on the Sorption Analysis control panel.
3. Specify the range of steps you want to include in the analysis
(probably excluding the initial pre-equilibrium steps of the simulation).
This range applies to all simulations in the set.
4. Select the Plot Loading Curve button.
- The loading curve plot of sorbate molecules per cell as a function of pressure is sent to the text window. A separate data set for each sorbate type is generated, but all data sets are presented on the same graph.
Plotting mass clouds (dot density maps) of sorbate positions
Mass clouds are displayed in the Model window. They allow you to view preferred sorbate locations in 3D. (Mass clouds are most often shown in the model space that contains the framework.)
1. Open the Sorption Analysis control panel (see the online help)
by selecting Analysis from the SORPTION card.
2. Load the trajectory file of the simulation using the browser box
in the Sorption Analysis control panel.
3. Specify the range of steps you want to include in the analysis.
4. Open the Mass Cloud control panel (see the online help) by
choosing the second Preferences... button on the Sorption
Analysis control panel.
5. Choose to color the cloud by sorbate type or by energy. Cerius2
automatically sets the property range to that of the trajectory
file; edit this range if you want.
- The default color is red. In the Color Cloud by Energy selection, dots representing the highest and lowest energy sorbates are red, and dots representing intermediate energies are a rainbow of colors in between. On some monitors, red dots can appear very faint; if the red dots are hard to see, try a different color setting (for example, white (0) or a yellow shade (45 through 85)).
6. With the framework model in the current model space or in an
empty model space, select the Plot Mass Cloud button on the
Sorption Analysis control panel.
- Mass clouds can take some time to calculate and, once calculated, can take up a lot of memory.
Theory
Simulation methods
As mentioned above, three simulation methods are available in the Sorption module:
Details of each method are discussed below.
Fixed loading (canonical ensemble) simulation
In the fixed loading (or canonical ensemble) simulation, the number of sorbates in the framework is held constant throughout the simulation.
Note
The initial configuration is generated by placing the sorbate at an arbitrary position in the framework cell. The initial coordinates of the sorbate in the framework are the same as the coordinates of the sorbate in its model window. In this way, you can choose to begin a simulation by placing the sorbate in a pore of the framework.
Each subsequent configuration is generated by either a random translation or rotation of the sorbate molecule. The choice of move is governed by move probabilities and limited by the maximum translation and rotation step size. The maximum step size values can be adjusted automatically during the course of the simulation in order to achieve the target acceptance rate (equilibrium). For more information about move probabilities and acceptance rates, see "Setting up and running the simulation" on page 239 and online help for the Move Probabilities control panel.
Each generated configuration is accepted or rejected using a Metropolis (Metropolis 1953) algorithm based on the configurational energy change:
Eq. 1 
Where:
P = Probability of the move being accepted.
ýE = Energy change between the new configuration and the previous configuration.
k = Boltzmann's constant.
T = Temperature of the simulation.
If the change in energy (ýE) resulting from the move is negative, the new configuration is accepted. If the move results in an energy increase (ýE>0), the Boltzmann factor (exp(-ýE/kT)) is computed and compared to a randomly generated number between zero and one. If the random number is less than the Boltzmann factor, the move is selected; otherwise it is rejected.
The simulation takes several steps to equilibrate from its initial position; this equilibration can be recognized by the convergence of the configurational energy. Simulation statistics are valid only for results taken from a sufficiently large number of configurations after this equilibration period.
Fixed pressure (grand canonical ensemble) simulation
The second simulation method available in Sorption uses a grand canonical ensemble; that is, the number of particles in the system is not fixed, but the chemical potential of each species is fixed. Sorption translates the chemical potential into the partial pressure (or fugacity) of each component. The basis for the simulation is as follows:
Note
The initial configuration is generated by one of four moves, for which the acceptance criteria are different. Each type of move is described below:
Eq. 2 
Where:
P = Probability of the new configuration being accepted.
ýE = Energy change between the new configuration and the previous configuration.
k = Boltzmann's constant.
T = Temperature of the simulation.
Ni = Current number of molecules of component i in framework.
fi = Fugacity of component i in the gas phase.
V = Cell volume. If a surface constrained calculation is performed, the volume used corresponds only to the framework and surface regions -- i.e., the vacuum layer above the specified surface region is not included.
Eq. 3 
For fixed loading, the move is chosen according to given move probabilities. For translations and rotations, a maximum value is placed on the allowed movement. This value can be adjusted automatically during the simulation to achieve a target acceptance rate (see "Setting up and running the simulation" on page 239 and online help for the Sorption Step Scaling control panel).
The result of this simulation is a set of configurations that converge towards the specified fugacity and temperature. The simulation takes several steps to equilibrate from its original random position. For accurate statistical results, the steps made prior to equilibration should be excluded in the analysis, and a large number of configurations should be generated after the equilibration period.
Henry's constant simulation
A third type of sorption simulation calculates Henry's constant, which is defined as the simulation-cell loading divided by the sorbate pressure in the limit of vanishing pressure:
Eq. 4 
Where:
K = Henry's constant with units of molecules/cell/kPa.
P = Pressure in units of kPa.
Loading = Number of sorbate molecules in the cell.
Henry's constant can be calculated by performing a series of fixed pressure simulations, and extrapolating the ratio to zero pressure. More conveniently, it can be found from the expression (Bezus et al. 1978):
Eq. 5 
Where:
K = Henry's constant.
r = Positional degrees of freedom for the sorbate molecule.
= Rotational degrees of freedom for the sorbate molecule.
U(r,) = Energy of a sorbate molecule in the cell at position r
orientated by .
The integral of Eq. 5 can be approximated by a finite sum:
Eq. 6 
Where:
N = Number of steps in the simulation.
KN = Henry's constant approximated after N steps.
Vcell = Volume of the cell
In the simulation, the positions ri and the orientations i are chosen at random. No step size control is used, and the Metropolis algorithm is not applied. A check for bad contacts is made to identify high energy positions that give negligible contribution to the integrand (see "Settings for the energy calculation" on page 232).
If the graphical-simulation monitor is selected, KN is displayed as N increases. This can be used to assess the number of steps, N, needed to obtain a well-converged estimate of K.
References
Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J. Chem. Phys., 21 1087 (1953).
Allen, M.P.; Tildesley, D.J. Computer Simulation of Liquids, Clarendon Press, Oxford, (1987).
Karavias, F.; Myers, A.L. Molecular Simulation, 8 23 (1991).
Stapleton, M. R.; Tildesley, D. J.; Quirke, N. J. Chem. Phys., 92 4456 (1990).
Bezus, A. G.; Kiselev, A. V.; Lopatkin, A. A.; Pham Quang Du J. Chem.Soc., Faraday Trans.2, 74 367 (1978).
Ewald, P.P. Ann. d. Physik, 64 253 (1921).
Last updated December 14, 1998 at 01:02PM PST.
Copyright © 1998, Molecular Simulations, Inc. All rights
reserved.
µgas = mframe
(where