GCMC Class

The GCMC class serves as the main interface for handling Grand Canonical Monte Carlo (GCMC) simulation logic.

class flames.gcmc.GCMC(model, framework_atoms, adsorbates, temperature, pressure, device, vdw_radii, vdw_factor=0.6, framework_energy=None, adsorbate_energy=None, max_translation=1.5, max_rotation=np.float64(1.5707963267948966), max_deltaE=1.555, save_frequency=100, save_rejected=False, output_to_file=True, output_folder=None, debug=False, fugacity_coeff=1.0, random_seed=None, cutoff_radius=6.0, automatic_supercell=True, void_fraction=0.0)

Bases: BaseSimulator

Base class for Grand Canonical Monte Carlo (GCMC) simulations using ASE.

This class employs Monte Carlo simulations under the grand canonical ensemble (\(μVT\)) to study the adsorption of molecules in a framework material. It allows for movements such as insertion, deletion, translation, and rotation of adsorbate molecules within the framework.

Currently, it supports any ASE-compatible calculator for energy calculations.

Parameters:

• model (ase.calculators.calculator.Calculator) – The calculator to use for energy calculations. Can be any ASE-compatible calculator. The output of the calculator should be in eV.

• framework_atoms (ase.Atoms) – The framework structure as an ASE Atoms object.

• adsorbates (Adsorbate) – The adsorbate structure as an Adsorbate object.

• temperature (float) – Temperature of the ideal reservoir in Kelvin.

• pressure (float) – Pressure of the ideal reservoir in Pascal.

• device (str) – Device to run the simulation on, e.g., 'cpu' or 'cuda'.

• vdw_radii (np.ndarray) – Van der Waals radii for the atoms in the framework and adsorbate. Should be an array of the same length as the number of atomic numbers in ASE.

• max_deltaE (float, optional) – Maximum energy difference (in eV) to consider for acceptance criteria. This is used to avoid overflow due to problematic calculations. Default is 1.555 eV (approx. 150 kJ/mol).

• vdw_factor (float, optional) – Factor to scale the Van der Waals radii. Default is 0.6.

• framework_energy (float or None, optional) – Pre-calculated potential energy of the empty framework in eV. If not provided, it will be calculated during initialization.

• adsorbate_energy (float or None, optional) – Pre-calculated potential energy of the adsorbate molecule in eV. If not provided, it will be calculated during initialization.

• max_translation (float, optional) – Maximum translation distance. Default is 1.5.

• max_rotation (float, optional) – Maximum rotation angle (in radians). Default is 90 degrees (converted to radians).

• save_frequency (int, optional) – Frequency at which to save the simulation state and results. Default is 100.

• save_rejected (bool, optional) – If True, saves the rejected moves in a trajectory file. Default is False.

• output_to_file (bool, optional) – If True, writes the output to a file named output_{temperature}_{pressure}.out in the results directory. Default is True.

• output_folder (str or None, optional) – Folder to save the output files. If None, a folder named results_<T>_<P> will be created.

• debug (bool, optional) – If True, prints detailed debug information during the simulation. Default is False.

• fugacity_coeff (float, optional) – Fugacity coefficient to correct the pressure. Default is 1.0. Only used if criticalTemperature, criticalPressure, and acentricFactor are not provided.

• random_seed (int or None, optional) – Random seed for reproducibility. Default is None and will generate a random seed automatically if not provided.

• cutoff_radius (float, optional) – Interaction potential cut-off radius used to estimate the minimum unit cell. Default is 6.0.

• automatic_supercell (bool, optional) – If True, automatically creates a supercell based on the cutoff radius. Default is True.

• max_length (float or None, optional) – Maximum length (in Angstroms) for any side of the supercell. If None, no maximum length is enforced. This can be used to limit the size of the supercell for computational efficiency. Default is None.

• criticalTemperature (float, optional) – Critical temperature of the adsorbate in Kelvin.

• criticalPressure (float, optional) – Critical pressure of the adsorbate in Pascal.

• acentricFactor (float, optional) – Acentric factor of the adsorbate.

• void_fraction (float, optional) – Void fraction of the adsorbate.

equilibrate(equilibration_steps=0, LLM=False, batch_size=False, run_ADF=False, uncertainty='uSD')

Use pyMSER to get the equilibrated statistics of the simulation.

Parameters:

• equilibration_steps (int) – Number of steps to use for equilibration. Default is 0.

• LLM (bool) – If True, use the Leftmost Local Minima (LLM) on the determination of the equilibration point by pyMSER. this can underestimate the equilibration point in some situations, but generate good averages for well-behaved scenarios. Default is False.

• batch_size (int) – Batch size to use for speedup the equilibration process. Default is 100.

• run_ADF (bool) – If True, run the Augmented Dickey-Fuller (ADF) test to confirm for stationarity. Default is False.

• uncertainty (str) – The type of uncertainty to use for the equilibration process. Default is “uSD”. Options are: - “uSD”: uncorrelated Standard Deviation - “uSE”: uncorrelated Standard Error - “SD”: Standard Deviation - “SE”: Standard Error

Return type:

None

get_adsorbates_index(tag=None)

Get a list of index of the adsorbate molecules in the current system. If a tag is provided, only returns the indices for adsorbates matching that tag.

Returns:

A list of lists containing the indices of the adsorbate molecules.

Return type:

list[list]

get_average_ads_energy()

Compute the average adsorption energy per adsorbed molecule.

The adsorption energy is calculated as:

E_ads_avg = [E_total - sum_i(N_ads_i * E_adsorbate_i) - E_framework] / N_ads

where:

• E_total: current total energy of the system (simulation)

• N_ads: number of adsorbed molecules

• E_adsorbate: energy of an isolated adsorbate molecule

• E_framework: energy of the empty framework

The result is converted from simulation units to kJ/mol per adsorbate.

Returns:

average_adsorption_energy – The average adsorption energy per molecule in kJ/mol. Returns 0.0 if no molecules are adsorbed (N_ads == 0).

Return type:

float

get_framework_mass()

Calculate the mass of the framework in kg.

Returns:

The mass of the framework in kg.

Return type:

float

get_number_of_adsorbates(system=None)

Get the number of adsorbates in the current system. It considers the possibility of having multiple adsorbate species in the simulation.

Parameters:

system (ase.Atoms | None) – The ASE Atoms object representing the current system.

Returns:

n_adsorbate – A dictionary mapping adsorbate names to their respective counts.

Return type:

dict[str, int]

get_total_ads_energy()

Compute the total adsorption energy of the system.

The adsorption energy is calculated as:

E_ads_total = E_total - sum_i(N_ads_i * E_adsorbate_i) - E_framework

Returns:

total_adsorption_energy – The total adsorption energy of the system in eV. Returns 0.0 if no molecules are adsorbed (N_ads == 0).

Return type:

float

insert_adsorbates(n_adsorbates, ads_name, max_attempts=1000)

Insert a given number of adsorbate molecules into the framework at random positions without overlap.

Parameters:

• n_molecules (int) – Number of adsorbate molecules to insert.

• ads_name (str) – Name of the adsorbate type to insert.

• max_attempts (int) – Maximum number of attempts to insert the molecules without overlap. Default is 1000.

Return type:

None

load_state(state_file)

Load the state of the simulation from a file.

Parameters:

state_file (str) – Path to the file containing the saved state of the simulation.

Return type:

None

md(nsteps, atoms=None, time_step=0.5, ensemble='NVT', thermostat='NoseHoover', output_interval=100, movie_interval=100, calculator=None, update_state=True, trajectory_file=None, set_momenta=True, **kwargs)

Run a molecular dynamics simulation using the specified ensemble and thermostat.

Parameters:

• nsteps (int) – Number of steps to run the MD simulation.

• time_step (float, optional) – Time step for the MD simulation in femtoseconds (default is 0.5 fs).

• ensemble (str, optional) – The ensemble to use for the MD simulation (default is “NVT”). Can be one of “NVT” or “NPT”.

• thermostat (str, optional) – The thermostat to use for the MD simulation (default is “NoseHoover”). For NVT, can be one of “Berendsen”, “NoseHoover”, or “Langevin”. For NPT, can be one of “Berendsen”, “NoseHoover”, or “MTKNPT”.

• output_interval (int, optional) – The interval for logging output (default is 100 steps).

• movie_interval (int, optional) – The interval for saving trajectory frames (default is 100 steps).

• calculator (ase.calculators.calculator.Calculator or None, optional) – The calculator to use for energy calculations. If None, the default model will be used.

• update_state (bool, optional) – If True, updates the current state of the simulation with the final state after MD. If False, returns the final state without updating. This difference is important for using the same method in the GCMC simulations when a MD step is used as a move. Default is True.

• trajectory_file (str or None, optional) – If provided, saves the trajectory to the specified file instead of the default trajectory.

• set_momenta (bool, optional) – Whether to set the atomic momenta to a Maxwell-Boltzmann distribution of the simulation temperature.

• **kwargs (optional) –

  Additional parameters passed directly to the specific MD thermostat (e.g., taut, tdamp, friction).

  NVT Berendsen:

  • tautfloat, optional

    Time constant for the Berendsen thermostat in fs (default is 1.0 fs).

  NVT Nose-Hoover:

  • tdampfloat, optional

    Time constant for the Nose-Hoover thermostat in fs (default is 50.0 fs).

  • tchainint, optional

    Number of thermostats in the Nose-Hoover chain (default is 3).

  • tloopint, optional

    Number of loops for the Nose-Hoover chain (default is 1).

  NVT Langevin:

  • frictionfloat, optional

    Friction coefficient for the Langevin dynamics (default is 0.01).

  NPT Berendsen:

  • isotropicbool, optional

    Whether to use isotropic pressure coupling (default is True).

  • compressibilityfloat, optional

    Compressibility of the material in bar^-1 (default is 1e-4 bar^-1).

  • tautfloat, optional

    Time constant for the Berendsen thermostat in fs (default is 10.0 fs).

  • taupfloat, optional

    Time constant for the Berendsen barostat in fs (default is 500.0 fs).

  NPT Nose-Hoover:

  • ttimefloat, optional

    Time constant for the Nose-Hoover thermostat in fs (default is 25.0 fs).

  • ptimefloat, optional

    Time constant for the Parrinello-Rahman barostat in fs (default is 75.0 fs).

  • bulk_modulusfloat, optional

    Bulk modulus of the material in GPa (default is 30.0 GPa).

  NPT MTK:

  • tdampfloat, optional

    Time constant for the Nose-Hoover thermostat in fs (default is 50.0 fs).

  • pdampfloat, optional

    Time constant for the Nose-Hoover barostat in fs (default is 500.0 fs).

  • tchainint, optional

    Number of thermostats in the Nose-Hoover chain (default is 3).

  • pchainint, optional

    Number of barostats in the Nose-Hoover chain (default is 3).

  • tloopint, optional

    Number of loops for the Nose-Hoover thermostat chain (default is 1).

  • ploopint, optional

    Number of loops for the Nose-Hoover barostat chain (default is 1).

  • vol_constraintbool, optional

    If True, the (N, V, sigma_a = 0, T)-ensemble is sampled, which allows for full cell fluctuations while keeping the cell volume fixed (default is False).

Return type:

None | Atoms

npt(nsteps, time_step=0.5, mode='iso_shape', driver='MTKNPT', set_momenta=True, output_interval=100, movie_interval=100, calculator=None, **kwargs)

Run a NPT simulation using the Berendsen thermostat and barostat.

Parameters:

• nsteps (int) – Number of steps to run the NPT simulation.

• time_step (float, optional) – Time step for the NPT simulation (default is 0.5 fs).

• mode (str, optional) – The mode of the NPT simulation (default is “iso_shape”). Can be one of “iso_shape”, “aniso_shape”, or “aniso_flex”.

• driver (str, optional) – The driver to use for the NPT simulation. Can be Berendsen, NoseHoover or MTKNPT (default is “MTKNPT”).

• set_momenta (bool, optional) – Whether to set the atomic momenta to a Maxwell-Boltzmann distribution of the simulation temperature.

• output_interval (int, optional) – The interval for logging output (default is 100 steps).

• movie_interval (int, optional) – The interval for saving trajectory frames (default is 100 step).

• calculator (ase.calculators.calculator.Calculator or None, optional) – The calculator to use for energy calculations. If None, the default model will be used.

• kwargs (optional) – Arguments passed to the ase molecular dynamics class.

Return type:

None | Atoms

nvt(nsteps, time_step=0.5, set_momenta=True, driver='NoseHoover', output_interval=100, movie_interval=100, calculator=None, **kwargs)

Run a NVT simulation using the Berendsen thermostat.

Parameters:

• nsteps (int) – Number of steps to run the NVT simulation.

• time_step (float, optional) – Time step for the NVT simulation in femtoseconds (default is 0.5 fs).

• set_momenta (bool, optional) – Whether to set the atomic momenta to a Maxwell-Boltzmann distribution of the simulation temperature.

• driver (str, optional) – The driver to use for the NVT simulation. Can be “NoseHoover” (default), Langevin, or “Berendsen” (Not recommended, for testing purposes only).

• output_interval (int, optional) – The interval for logging output (default is 100 steps).

• movie_interval (int, optional) – The interval for saving trajectory frames (default is 100 step).

• calculator (ase.calculators.calculator.Calculator or None, optional) – The calculator to use for energy calculations. If None, the default model will be used.

• kwargs (optional) – Arguments passed to the ase molecular dynamics class.

• NoseHoover (For)

• tdamp (-) – Time constant for the thermostat in fs (default is 50.0 fs).

• tchain (-) – Length of the Nose-Hoover chain (default is 3).

• tloop (-) – Number of loops for the Nose-Hoover chain (default is 1).

• Langevin (For)

• friction (-) – Friction coefficient for the Langevin thermostat in fs^-1 (default is 0.1 fs^-1).

• Berendsen (For)

• tau (-) – Time constant for the Berendsen thermostat in fs (default is 1.0 fs).

optimize_adsorbate(max_steps=1000, max_force=0.05)

Optimize the adsorbate structure using the provided calculator.

Parameters:

• max_steps (int, optional) – Maximum number of optimization steps (default is 1000).

• symm_tol (float, optional) – Tolerance for symmetry (default is 1e-3).

• max_force (float, optional) – Maximum force tolerance for convergence (default is 0.05 eV/Å).

Returns:

The optimized adsorbate structure.

Return type:

ase.Atoms

optimize_framework(max_steps=1000, opt_cell=True, fix_symmetry=True, hydrostatic_strain=True, symm_tol=0.001, max_force=0.05)

Optimize the framework structure using the provided calculator.

Parameters:

• max_steps (int, optional) – Maximum number of optimization steps (default is 1000).

• tol (float, optional) – Tolerance for convergence (default is 1e-5).

Returns:

The optimized framework structure.

Return type:

ase.Atoms

restart()

Restart the simulation from the last state.

This method loads the last saved state from the trajectory file and restores the simulation to that state. It also loads the uptake, total energy, and total adsorbates lists from the saved files if they exist.

Return type:

None

run(N)

Run the Grand Canonical Monte Carlo simulation for N iterations.

Return type:

None

save_results(file_name=None, batch_size=False, run_ADF=False, uncertainty='uSD', LLM=False)

Save a json file with the main results of the simulation.

Parameters:

• file_name (str) – Name of the output file. Default is ‘GCMC_Results.json’.

• LLM (bool) – If True, use the Leftmost-Local Minima (LLM) method to determine the equilibration time. This is only recommended for high-throughput simulations, and sometimes can underestimate the true equilibration point. Default is True.

• batch_size (int) – Batch size to use for speedup the equilibration process. Default is False, which means 2% of the total number of steps.

• run_ADF (bool) – If True, run the Augmented Dickey-Fuller (ADF) test to confirm for stationarity. Default is False.

• uncertainty (str) – The type of uncertainty to use for the equilibration process. Default is “uSD”. Options are: - “uSD”: uncorrelated Standard Deviation - “uSE”: uncorrelated Standard Error - “SD”: Standard Deviation - “SE”: Standard Error

• LLM – If True, use the Leftmost-Local Minima (LLM) method to determine the equilibration time. This is only recommended for high-throughput simulations, and sometimes can underestimate the true equilibration point. Default is False.

Return type:

None

set_adsorbates(adsorbates, adsorbate_energy=None)

Set the adsorbate structure for the simulation.

Parameters:

• adsorbates (Adsorbate) – The list of adsorbate structure(s) as an Adsorbate object.

• adsorbate_energy (float or None, optional) – The energy of the adsorbate in eV. If None, the energy will be calculated using the provided model.

• n_adsorbates (int) – Number of adsorbate molecules in the framework.

Return type:

None

set_framework(framework_atoms, framework_energy=None)

Set the framework structure for the simulation.

Parameters:

• framework_atoms (ase.Atoms) – The new framework structure as an ASE Atoms object.

• framework_energy (float or None, optional) – The energy of the framework in eV. If None, the energy will be calculated using the provided model.

Return type:

None

set_state(state)

Set the current state of the simulation.

Parameters:

state (ase.Atoms) – The current state of the simulation as an ASE Atoms object.

Return type:

None

step(iteration)

Perform a single Grand Canonical Monte Carlo step. It will randomly select a move based on the move weights and attempt to perform it. The uptake, total energy, and total adsorbates lists are updated accordingly.

Parameters:

iteration (int) – The current iteration number.

Return type:

None

try_deletion(adsorbate_tag)

Try to delete an adsorbate molecule from the framework. This method randomly selects an adsorbate molecule and try to apply the deletion.

If there are no adsorbates, it returns False.

Parameters:

adsorbate_tag (int) – The tag of the adsorbate molecule to be deleted.

Returns:

True if the deletion was accepted, False otherwise.

Return type:

bool

try_identity_swap(adsorbate_tag)

Try to swap the identity of two adsorbate molecules in the system. This method randomly selects two adsorbate molecules and attempts to swap their identities.

Parameters:

adsorbate_tag (int) – The tag of the first adsorbate molecule to be swapped.

Returns:

True if the swap was accepted, False otherwise.

Return type:

bool

try_insertion(adsorbate_tag)

Try to insert a new adsorbate molecule into the framework. This method randomly places the adsorbate in the framework and checks for van der Waals overlap. If there is no overlap, it calculates the new potential energy and decides whether to accept the insertion based on the acceptance criteria.

Parameters:

adsorbate_tag (int) – The tag of the adsorbate molecule to be inserted.

Returns:

True if the insertion was accepted, False otherwise.

Return type:

bool

try_npt_md(n_teps=500, **kwargs)

Try to perform a NPT-MD move of the system using the MTK thermostat.

Parameters:

n_teps (int) – Number of time steps to run the NPT-MD simulation.

Returns:

True if the NPT-MD move was accepted, False otherwise.

Return type:

bool

try_nve_md(n_teps=50, **kwargs)

Try to perform a NVE-MD move of the system using the Nose-Hoover thermostat.

Parameters:

n_teps (int) – Number of time steps to run the NVE-MD simulation.

Returns:

True if the NVE-MD move was accepted, False otherwise.

Return type:

bool

try_nvt_md(n_teps=500, **kwargs)

Try to perform a NVT-MD move of the system using the Nose-Hoover thermostat.

Parameters:

n_teps (int) – Number of time steps to run the NVT-MD simulation.

Returns:

True if the NVT-MD move was accepted, False otherwise.

Return type:

bool

try_reinsertion(adsorbate_tag)

Try to delete and reinsert an adsorbate molecule. This method randomly selects an adsorbate molecule, deletes it, and tries to reinsert it at a new random position within the framework.

If there are no adsorbates, it returns False.

Parameters:

adsorbate_tag (int) – The tag of the adsorbate molecule to be reinserted.

Returns:

True if the reinsertion was accepted, False otherwise.

Return type:

bool

try_rotation(adsorbate_tag)

Try to rotate an adsorbate molecule within the framework. This method randomly selects an adsorbate molecule and applies a random rotation. It checks for van der Waals overlap and calculates the new potential energy.

Parameters:

adsorbate_tag (int) – The tag of the adsorbate molecule to be rotated.

Returns:

True if the rotation was accepted, False otherwise.

Return type:

bool

try_translation(adsorbate_tag)

Try to translate an adsorbate molecule within the framework. This method randomly selects an adsorbate molecule and applies a random translation. It checks for van der Waals overlap and calculates the new potential energy.

Parameters:

adsorbate_tag (int) – The tag of the adsorbate molecule to be translated.

Returns:

True if the translation was accepted, False otherwise.

Return type:

bool

property base_iteration: int

Get the base iteration for the GCMC simulation.

Returns:

The base iteration count.

Return type:

int