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_simFunction
results, 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 adsorption
  • molecule_templates::Array{Molecule, 1}: an array of the templates of unique adsorbate molecules of which we seek to simulate
  • temperature::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 adsorption
  • ljff::LJForceField: the molecular model used to describe the
  • molecules::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 between n_burn_cycles and 'n_sample_cycles, where n_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 sampling
  • sample_frequency::Int: during the sampling cycles, sample e.g. the number of adsorbed gas molecules every this number of Markov proposals
  • verbose::Bool: whether or not to print off information during the simulation
  • ewald_precision::Float64: desired precision for the long range Ewald summation
  • eos::Symbol: equation of state to use for calculation of fugacity from pressure
  • write_adsorbate_snapshots::Bool: whether the simulation will create and save a snapshot file
  • snapshot_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 adsorbates
  • density_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 by required_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)
source
PorousMaterials.adsorption_isothermFunction
results = 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.

source
PorousMaterials.stepwise_adsorption_isothermFunction
results = 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.

source
PorousMaterials.μVT_output_filenameFunction
filename = μ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 simulation
  • molecule_templates::Array{Molecule, 1}: template of the adsorbate molecules used in adsorption simulation
  • temperature::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: bar
  • ljff::LJForceField: the molecular model used in adsorption simulation
  • n_burn_cycles::Int: number of cycles to allow the system to reach equilibrium before sampling.
  • n_sample_cycles::Int: number of cycles used for sampling
  • comment::String="": remarks to be included in the filename
  • extension::String=".jld2": the file extension

Returns

  • filename::String: the name of the specific .jld2 simulation file
source
PorousMaterials.isotherm_sim_results_to_dataframeFunction
dataframe = 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 file
  • xtal::Crystal: The porous crystal
  • molecule::Array{Molecule{Cart}, 1}: The adsorbate molecule
  • temperature::Float64: The temperature in the simulation, units: Kelvin (K)
  • pressures::Array{Array{Float64, 1}, 1}: The pressures in the adsorption isotherm simulation(s), units: bar
  • ljff::LJForceField: The molecular model being used in the simulation to describe intermolecular Van der Waals interactions
  • n_burn_cycles::Int64: The number of burn cycles used in this simulation
  • n_sample_cycles::Int64: The number of sample cycles used in the simulations
  • where_are_jld_files::Union{String, Nothing}=nothing: The location to the simulation files. defaults to PorousMaterials.rc[:paths][:simulations].
  • comment::String="": comment appended to outputfilename

Returns

  • dataframe::DataFrame: A DataFrame 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
source