Data structure for a regular [equal spacing between points in each coordinate] grid of points superimposed on a unit cell box (Box). Each grid point has data, data, associated with it, of type T, stored in a 3D array.

Attributes

• box::Box: describes Bravais lattice over which a grid of points is super-imposed. grid points on all faces are included.
• n_pts::Tuple{Int, Int, Int}: number of grid points in x, y, z directions. 0 and 1 fractional coordinates are included.
• data::Array{T, 3}: three dimensional array conaining data associated with each grid point.
• units::Symbol: the units associated with each data point.
• origin::Array{Float64, 1}: the origin of the grid.

grid = energy_grid(crystal, molecule, ljforcefield; resolution=1.0, temperature=298.0, n_rotations=750)

Superimposes a regular grid of points (regularly spaced in fractional coordinates of the crystal.box) over the unit cell of a crystal, with n_gridpts dictating the number of grid points in the a, b, c directions (including 0 and 1 fractional coords). The fractional coordinates 0 and 1 are included in the grid, although they are redundant. Then, at each grid point, calculate the ensemble average potential energy of the molecule when its mass is centered at that point. The average is taken over Boltzmann-weighted rotations.

The ensemble average is a Boltzmann average over rotations: - R T log ⟨e⁻ᵇᵁ⟩

Arguments

• crystal::Crystal: crystal in which we seek to compute an energy grid for a molecule. grid.box will be framework.box.
• molecule::Molecule: molecule for which we seek an energy grid
• ljforcefield::LJForceField: molecular model for computing molecule-crystal interactions
• resolution::Union{Float64, Tuple{Int, Int, Int}}=1.0: maximum distance between grid points, in Å, or a tuple specifying the number of grid points in each dimension.
• n_rotations::Int: number of random rotations to conduct in a Monte Carlo simulation for finding the free energy of a molecule centered at a given grid point.

This is only relevant for molecules that are comprised of more than one Lennard Jones sphere.

• temperature::Float64: the temperature at which to compute the free energy for molecules where rotations are required. Lower temperatures overemphasize the minimum potential energy rotational conformation at that point.
• units::Symbol: either :K or :kJ_mol, the units in which the energy should be stored in the returned Grid.
• center::Bool: shift coords of grid so that the origin is the center of the unit cell crystal.box.
• verbose::Bool=true: print some information.

Returns

• grid::Grid: A grid data structure containing the potential energy of the system

write_cube(grid, filename, verbose=true)

Write grid to a .cube file format. This format is described here: http://paulbourke.net/dataformats/cube/ The atoms of the unit cell are not printed in the .cube. Instead, use .xyz files to also visualize atoms.

Arguments

• grid::Grid: grid with associated data at each grid point.
• filename::AbstractString: name of .cube file to which we write the grid; this is relative to rc[:paths][:grids].
• verbose::Bool: print name of file after writing.
• length_units::String: units for length. Bohr or Angstrom.

grid = read_cube(filename)

Read a .cube file and return a populated Grid data structure. It will detect and skip over atomic information if it is present in the file.

Arguments

• filename::AbstractString: name of .cube file to which we write the grid; this is relative to rc[:paths][:grids]
• has_units::Bool=true: flag for function to read units from file header

Returns

• grid::Grid: A grid data structure

n_pts = required_n_pts(box, dx)

Calculate the required number of grid pts in a, b, c unit cell directions required to keep distances between grid points less than dx apart, where dx is in units of Angstrom.

voxel_id = xf_to_id(n_pts, xf)

Returns the indices of the voxel in which it falls when a unit cube is partitioned into a regular grid of n_pts[1] by n_pts[2] by n_pts[3] voxels. Periodic boundary conditions are applied.

Arguments

• n_pts::Tuple{Int, Int, Int}: The number of points for each axis in the Grid
• xf::Array{Float64, 1}: The fractional coordinates to be converted to an id

Returns

• id::Array{Int, 1}: The array indices for storing this point in space

xf = id_to_xf(voxel_id, n_pts)

Given a voxel_id in a Grid, return the fractional coordinates to which this voxel corresponds.

Arguments

• n_pts::Tuple{Int, Int, Int}: The number of voxels along each axis in the Grid
• voxel_id::Array{Int, 1}: the voxel coordinates in grid.data

Returns

• xf::Array{Float64, 1}: The fractional coordinates corresponding to the grid voxel

update_density!(grid, molecule, species)

updates the density grid based on an array of molecules. If a molecule doesn't match the specified species it won't be added to the density grid. This function doesn't calculate the actual densities, it will need a ./ = num_snapshots at the end of the GCMC simulation.

Arguments

• grid::Grid: the grid to be updated
• molecules::Array{Molecule, 1}: An array of molecules whose positions will be added to the grid
• species::Symbol: The species of atom that can be added to this density grid

accessibility_grid, nb_segments_blocked, porosity = compute_accessibility_grid(crystal,
probe_molecule, ljforcefield; resolution::Union{Float64, Tuple{Int, Int, Int}}=1.0, energy_tol=10.0, energy_unit=:kJ_mol,
verbose=true, write_b4_after_grids=false, block_inaccessible_pockets=true)

Overlay a grid of points about the unit cell. Compute the potential energy of a probe molecule at each point. If the potential energy is less than energy_tol, the grid point is declared as accessible to an adsorbate; otherwise inaccessible.

If block_pockets is true: Then perform a flood fill algorithm to label disparate (unconnected) segments in the grid.

Then build a graph whose vertices are the unconnected segments in the flood-filled grid and whose edges are the connections between the segments across the periodic boundary.

Then find any simple cycles in the grid. Any vertex that is involved in a simple cycle is considered accessible since a molecule can travel from that segment in the home unit cell to the same segment but in a different unit cell. If any vertex is not involved in a cycle, the segment is declared as inaccessible and all grid points in this segment are re-labeled as inaccessible.

Returns accessibility_grid::Grid{Bool} and nb_segments_blocked, the latter the number of segments that were blocked because they were determined to be inaccessible.

Arguments

• crystal::Crystal: the crystal for which we seek to compute an accessibility grid.
• probe_molecule::Molecule a molecule serving as a probe to determine whether a given

point can be occupied and accessed.

• LJForceField::LJForceField: the force field used to compute the potential energy of

the probe molecule

• • resolution::Union{Float64, Tuple{Int, Int, Int}}=1.0: maximum distance between grid points, in Å, or a tuple specifying the number of grid points in each dimension.
• energy_tol::Float64: if the computed potential energy is less than this, we declare the

grid point to be occupiable. Also this is the energy barrier beyond which we assume the probe adsorbate cannot pass. Units given by energy_units argument

• energy_units::Symbol: units of energy (:kJ_mol or :K) to be used in determining

threshold for occupiability and whether molecule can percolate over barrier in channel. (see energy_tol)

• write_b4_after_grids::Bool: write a .cube file of occupiability for visualization both

before and after flood fill/blocking inaccessible pockets

accessible(accessibility_grid, xf)
accessible(accessibility_grid, xf, repfactors)

Using accessibility_grid, determine if fractional coordinate xf (relative to accessibility_grid.box is accessible or not. Here, we search for the nearest grid point. We then look at the accessibility of this nearest grid point and all surroudning 9 other grid points. The point xf is declared inaccessible if and only if all 10 of these grid points are inaccessible. We take this approach because, if the grid is coarse, we can allow energy computations to automatically determine accessibility at the boundary of accessibility e.g. during a molecular simulation where inaccessible pockets are blocked.

If a tuple of replication factors are also passed, it is assumed that the passed xf is relative to a replicated accessibility_grid.box so that xf is scaled by these rep. factors. So xf = [0.5, 0.25, 0.1] with repfactors=(1, 2, 4) actually is, relative to accessibility_grid.box, fractional coordinate [0.5, 0.5, 0.4].