Grand-canonical Monte Carlo Simulations
Simulate the adsorption of CO$_2$ in FIQCEN_clean_min_charges (CuBTC) at 298 K at 1 bar using the Universal Force Field:
xtal = Crystal("FIQCEN_clean.cif") # load a crystal structure
strip_numbers_from_atom_labels!(xtal) # clean up the atom labels
ljforcefield = LJForceField("UFF", r_cutoff=12.8) # load the UFF forcefield
molecule = Molecule("CO2") # load the CO2 molecule
temperature = 298.0 # K
pressure = 1.0 # bar
# conduct Grand-Canonical Monte Carlo simulation (VERY short, should use thousands of cycles!)
results, molecules = μVT_sim(xtal, molecule, temperature, pressure, ljforcefield,
n_burn_cycles=50, n_sample_cycles=50)
Or, compute the entire adsorption isotherm at once, parallelized across many cores (this works by cleverly queuing a μVT_sim
for each pressue across the specified number of cores for optimal efficiency):
pressures = [0.2, 0.6, 0.8, 1.0] # bar
# loop over all pressures and compute entire adsorption isotherm in parallel
results = adsorption_isotherm(xtal, molecule, temperature, pressures, ljforcefield,
n_burn_cycles=50, n_sample_cycles=50)
Or, compute the adsorption isotherm in a step-wise manner, loading the molecules from the previous simulation to save on burn cycles:
results = stepwise_adsorption_isotherm(xtal, molecule, temperature, pressures, forcefield,
n_burn_cycles=50, n_sample_cycles=50)
detailed docs
Grand-Canonical Monte Carlo Simulations
PorousMaterials.μVT_sim
— Functionresults, molecules = μVT_sim(xtal, molecule_templates, temperature, pressure,
ljff; molecules=Array{Molecule, 1}[], settings=settings)
Runs a grand-canonical (μVT) Monte Carlo simulation of the adsorption of a molecule in a xtal at a particular temperature and pressure using a Lennard Jones force field.
Markov chain Monte Carlo moves include:
- deletion/insertion
- translation
- reinsertion
- identity change (if multiple components) see here
Translation stepsize is dynamically updated during burn cycles so that acceptance rate of translations is ~0.4.
A cycle is defined as max(20, number of adsorbates currently in the system) Markov chain proposals.
Arguments
xtal::Crystal
: the porous xtal in which we seek to simulate adsorptionmolecule_templates::Array{Molecule, 1}
: an array of the templates of unique adsorbate molecules of which we seek to simulatetemperature::Float64
: temperature of bulk gas phase in equilibrium with adsorbed phase in the porous material. units: Kelvin (K)pressures::Array{Float64, 1}
: pressure of bulk gas phase in equilibrium with adsorbed phase in the porous material for each adsorbate. units: bar the adsorptionljff::LJForceField
: the molecular model used to describe themolecules::Array{Array{Molecule{Cart}, 1}, 1}
: a starting configuration of molecules in the xtal with an array per species.n_cycles_per_volume::Int
: the number of MC cycles per volume, split evenly betweenn_burn_cycles
and'n_sample_cycles
, wheren_burn_cycles
is the number of cycles to allow the system to reach equilibrium before sampling; and,n_sample_cycles
is the number of cycles used for samplingsample_frequency::Int
: during the sampling cycles, sample e.g. the number of adsorbed gas molecules every this number of Markov proposalsverbose::Bool
: whether or not to print off information during the simulationewald_precision::Float64
: desired precision for the long range Ewald summationeos::Symbol
: equation of state to use for calculation of fugacity from pressurewrite_adsorbate_snapshots::Bool
: whether the simulation will create and save a snapshot filesnapshot_frequency::Int
: the number of cycles taken between each snapshot (after burn cycle completion)calculate_density_grid::Bool
: whether the simulation will keep track of a density grid for adsorbatesdensity_grid_dx::Float64
: The (approximate) space between voxels (in Angstroms) in the density grid. The number of voxels in the simulation box is computed automatically byrequired_n_pts
.density_grid_molecular_species::Symbol
: the adsorbate for which we will make a density grid of its position (center).density_grid_sim_box::Bool
:true
if we wish for the density grid to be over the
entire simulation box as opposed to the box of the crystal passed in. false
if we wish the density grid to be over the original xtal.box
, before replication, passed in.
autosave::Bool
:true
if we wish to automatically save the simulation results to the standard path/filename.results_filename_comment::AbstractString
: An optional comment that will be appended to the name of the saved file (if autosaved)
PorousMaterials.adsorption_isotherm
— Functionresults = adsorption_isotherm(xtal, molecule_templates, temperature,
partial_pressures,
ljff; n_sample_cycles=5000,
n_burn_cycles=5000, sample_frequency=1,
verbose=true, ewald_precision=1e-6, eos=:ideal,
write_adsorbate_snapshots=false,
snapshot_frequency=1, calculate_density_grid=false,
density_grid_dx=1.0, density_grid_species=nothing,
density_grid_sim_box=true,
results_filename_comment="", show_progress_bar=false)
Run a set of grand-canonical (μVT) Monte Carlo simulations in parallel. Arguments are the same as μVT_sim
, as this is the function run in parallel behind the scenes. The only exception is that we pass an array of pressures, and we only consider a single species. To give Julia access to multiple cores, run your script as julia -p 4 mysim.jl
to allocate e.g. four cores. See Parallel Computing.
PorousMaterials.stepwise_adsorption_isotherm
— Functionresults = stepwise_adsorption_isotherm(xtal, molecule_templates, temperature, pressures,
ljff; n_sample_cycles=5000,
n_burn_cycles=5000, sample_frequency=1,
verbose=true, ewald_precision=1e-6, eos=:ideal,
write_adsorbate_snapshots=false,
snapshot_frequency=1, calculate_density_grid=false,
density_grid_dx=1.0, density_grid_species=nothing,
density_grid_sim_box::Bool=true,
results_filename_comment="", show_progress_bar=false)
Run a set of grand-canonical (μVT) Monte Carlo simulations in series. Arguments are the same as μVT_sim
, as this is the function run behind the scenes. An exception is that we pass an array of pressures. The adsorption isotherm is computed step- wise, where the ending configuration from the previous simulation (array of molecules) is passed into the next simulation as a starting point. The ordering of pressures
is honored. By giving each simulation a good starting point, (if the next pressure does not differ significantly from the previous pressure), we can reduce the number of burn cycles required to reach equilibrium in the Monte Carlo simulation. Also see adsorption_isotherm
which runs the μVT simulation at each pressure in parallel.
PorousMaterials.μVT_output_filename
— Functionfilename = μVT_output_filename(xtal, molecule_templates, temperature,
pressures, ljff, n_burn_cycles,
n_sample_cycles; comment="", extension=".jld2")
This is the function that establishes the file naming convention used by μVT_sim.
Arguments
xtal::Crystal
: porous xtal used in adsorption simulationmolecule_templates::Array{Molecule, 1}
: template of the adsorbate molecules used in adsorption simulationtemperature::Float64
:temperature of bulk gas phase in equilibrium with adsorbed phase in the porous material. units: Kelvin (K)pressures::Array{Float64, 1}
: partial pressures of bulk gas phase in equilibrium with adsorbed phase in the porous material. units: barljff::LJForceField
: the molecular model used in adsorption simulationn_burn_cycles::Int
: number of cycles to allow the system to reach equilibrium before sampling.n_sample_cycles::Int
: number of cycles used for samplingcomment::String=""
: remarks to be included in the filenameextension::String=".jld2"
: the file extension
Returns
filename::String
: the name of the specific.jld2
simulation file
PorousMaterials.isotherm_sim_results_to_dataframe
— Functiondataframe = isotherm_sim_results_to_dataframe(desired_props, xtal,
molecule, temperature,
pressures, ljff,
n_burn_cycles, n_sample_cycles;
where_are_jld_files=nothing)`
convert the .jld2
results output files auto-saved from μVT_sim
into a DataFrame
. each row of the DataFrame
corresponds to a pressure in the adsorption isotherm. desired_props
is an array of desired properties from the results. to locate the requested output files, this function calls μVT_output_filename
to retrieve the file names.
Arguments
desired_props::Array{String, 1}
: An array containing names of properties to be retrieved from the.jld2
filextal::Crystal
: The porous crystalmolecule::Array{Molecule{Cart}, 1}
: The adsorbate moleculetemperature::Float64
: The temperature in the simulation, units: Kelvin (K)pressures::Array{Array{Float64, 1}, 1}
: The pressures in the adsorption isotherm simulation(s), units: barljff::LJForceField
: The molecular model being used in the simulation to describe intermolecular Van der Waals interactionsn_burn_cycles::Int64
: The number of burn cycles used in this simulationn_sample_cycles::Int64
: The number of sample cycles used in the simulationswhere_are_jld_files::Union{String, Nothing}=nothing
: The location to the simulation files. defaults toPorousMaterials.rc[:paths][:simulations]
.comment::String=""
: comment appended to outputfilename
Returns
dataframe::DataFrame
: ADataFrame
containing the simulated data along the adsorption isotherm, whose columns are for the specified properties
Note
A range of pressures can be used to select a batch of simulation files to be included in the DataFrame
.
Example
xtal = Crystal("SBMOF-1.cif")
molecule = Molecule("Xe")
ljff = LJForceField("UFF", mixing_rules="Lorentz-Berthelot")
temperature = 298.0 # K
pressures = 10 .^ range(-2, stop=2, length=15)
dataframe = isotherm_sim_results_to_dataframe(["pressure (bar)", "⟨N⟩ (mmol/g)"],
xtal, molecule, temperature,
pressures, ljff, 10000, 10000)
dataframe[Symbol("pressure (bar)")] # pressures
dataframe[Symbol("⟨N⟩ (mmol/g)")] # adsorption at corresponding pressures