Lennard-Jones Force Fields and Potential Energy

Lennard-Jones force field parameters are stored in comma-separated-value format in rc[:paths][:forcefields].

Interaction of an adsorbate with the crystal is modeled as pair-wise additive and with Lennard-Jones potentials of the form:

V(r) = 4 * ϵ * [ x ^ 12 - x ^ 6 ], where x = σ / r

The Lennard-Jones force field input files, e.g. UFF.csv contain a list of pure (i.e. X-X, where X is an atom) sigmas (σ) and epsilons (ϵ) with units Angstrom (Å) and Kelvin (K), respectively. Note that, e.g., in the UFF paper, the Lennard-Jones potential is written in a different form; thus, parameters need to be converted to correspond to the functional form used in PorousMaterials.jl.

Building Blocks of PorousMaterials: Lennard-Jones Force Fields

Loading Force Field Files and Accessing Attributes

Reading in Lennard-Jones force field parameters is made easy with the LJForceField function. Let's load in the parameters from the Universal Force Field file (UFF.csv):

# read in Lennard-Jones force field parameters from the Universal Force Field
ljforcefield = LJForceField("UFF", r_cutoff=14.0, mixing_rules="Lorentz-Berthelot")
# output
Force field: UFF
Number of atoms included: 108
Cut-off radius (Å) = 14.0
   Ra-   Ra ϵ =  203.30076 K, σ =    3.27583 Å
   Cl-   Cl ϵ =  114.23087 K, σ =    3.51638 Å
   Al-   Al ϵ =  254.12595 K, σ =    4.00815 Å
   Be-   Be ϵ =   42.77367 K, σ =    2.44552 Å
   Re-   Re ϵ =   33.21250 K, σ =    2.63171 Å
   Cr-   Cr ϵ =    7.54830 K, σ =    2.69319 Å
   Na-   Na ϵ =   15.09659 K, σ =    2.65755 Å
   Sb-   Sb ϵ =  225.94565 K, σ =    3.93777 Å
   Cf-   Cf ϵ =    6.54186 K, σ =    2.95155 Å
   Kr-   Kr ϵ =  110.70833 K, σ =    3.68921 Å
   Ni-   Ni ϵ =    7.54830 K, σ =    2.52481 Å
    S-    S ϵ =  137.88220 K, σ =    3.59478 Å
  CH4-  CH4 ϵ =  148.00000 K, σ =    3.73000 Å
   Fm-   Fm ϵ =    6.03864 K, σ =    2.92749 Å
   Ru-   Ru ϵ =   28.18030 K, σ =    2.63973 Å
   Tl-   Tl ϵ =  342.18940 K, σ =    3.87274 Å
   re-   re ϵ =    0.00010 K, σ =    5.00000 Å
   Tm-   Tm ϵ =    3.01932 K, σ =    3.00589 Å
C_CO2-C_CO2 ϵ =   27.00000 K, σ =    2.80000 Å
    W-    W ϵ =   33.71572 K, σ =    2.73417 Å
    O-    O ϵ =   30.19318 K, σ =    3.11815 Å
   Nd-   Nd ϵ =    5.03220 K, σ =    3.18496 Å
   Tb-   Tb ϵ =    3.52254 K, σ =    3.07449 Å
   Th-   Th ϵ =   13.08371 K, σ =    3.02549 Å
   Zr-   Zr ϵ =   34.72216 K, σ =    2.78317 Å
    F-    F ϵ =   25.16099 K, σ =    2.99698 Å
   Co-   Co ϵ =    7.04508 K, σ =    2.55866 Å
   Fr-   Fr ϵ =   25.16099 K, σ =    4.36540 Å
   Gd-   Gd ϵ =    4.52898 K, σ =    3.00055 Å
   Rh-   Rh ϵ =   26.67064 K, σ =    2.60944 Å
   Pu-   Pu ϵ =    8.05152 K, σ =    3.05044 Å
   Lw-   Lw ϵ =    5.53542 K, σ =    2.88295 Å
   Ar-   Ar ϵ =   93.09564 K, σ =    3.44600 Å
   Ca-   Ca ϵ =  119.76629 K, σ =    3.02816 Å
   Cm-   Cm ϵ =    6.54186 K, σ =    2.96313 Å
    N-    N ϵ =   34.72216 K, σ =    3.26069 Å
   As-   As ϵ =  155.49489 K, σ =    3.76850 Å
   Yb-   Yb ϵ =  114.73409 K, σ =    2.98897 Å
   Se-   Se ϵ =  146.43693 K, σ =    3.74623 Å
    Y-    Y ϵ =   36.23182 K, σ =    2.98006 Å
   Am-   Am ϵ =    7.04508 K, σ =    3.01213 Å
   Pt-   Pt ϵ =   40.25758 K, σ =    2.45354 Å
    I-    I ϵ =  170.59148 K, σ =    4.00904 Å
   Fe-   Fe ϵ =    6.54186 K, σ =    2.59430 Å
   Ba-   Ba ϵ =  183.17197 K, σ =    3.29900 Å
   Hf-   Hf ϵ =   36.23182 K, σ =    2.79831 Å
   Es-   Es ϵ =    6.03864 K, σ =    2.93907 Å
   Po-   Po ϵ =  163.54640 K, σ =    4.19524 Å
   Eu-   Eu ϵ =    4.02576 K, σ =    3.11191 Å
    C-    C ϵ =   52.83807 K, σ =    3.43085 Å
   Zn-   Zn ϵ =   62.39924 K, σ =    2.46155 Å
   Cs-   Cs ϵ =   22.64489 K, σ =    4.02419 Å
   Mn-   Mn ϵ =    6.54186 K, σ =    2.63795 Å
   Rn-   Rn ϵ =  124.79849 K, σ =    4.24513 Å
   Bk-   Bk ϵ =    6.54186 K, σ =    2.97471 Å
   Ir-   Ir ϵ =   36.73504 K, σ =    2.53015 Å
   Rb-   Rb ϵ =   20.12879 K, σ =    3.66516 Å
   In-   In ϵ =  301.42860 K, σ =    3.97608 Å
   Hg-   Hg ϵ =  193.73958 K, σ =    2.40988 Å
   Te-   Te ϵ =  200.28144 K, σ =    3.98232 Å
   At-   At ϵ =  142.91439 K, σ =    4.23177 Å
   Bi-   Bi ϵ =  260.66780 K, σ =    3.89323 Å
   Cu-   Cu ϵ =    2.51610 K, σ =    3.11369 Å
   Tc-   Tc ϵ =   24.15455 K, σ =    2.67091 Å
   Sn-   Sn ϵ =  285.32557 K, σ =    3.91283 Å
   Pa-   Pa ϵ =   11.07083 K, σ =    3.05044 Å
   Lu-   Lu ϵ =   20.63201 K, σ =    3.24287 Å
   Mo-   Mo ϵ =   28.18030 K, σ =    2.71902 Å
   Ac-   Ac ϵ =   16.60625 K, σ =    3.09855 Å
    U-    U ϵ =   11.07083 K, σ =    3.02460 Å
   Li-   Li ϵ =   12.58049 K, σ =    2.18359 Å
   Er-   Er ϵ =    3.52254 K, σ =    3.02104 Å
   Ta-   Ta ϵ =   40.76080 K, σ =    2.82415 Å
S_H2S-S_H2S ϵ =  122.00000 K, σ =    3.60000 Å
   Cd-   Cd ϵ =  114.73409 K, σ =    2.53728 Å
   Os-   Os ϵ =   18.61913 K, σ =    2.77960 Å
   Ti-   Ti ϵ =    8.55473 K, σ =    2.82860 Å
    B-    B ϵ =   90.57955 K, σ =    3.63754 Å
    V-    V ϵ =    8.05152 K, σ =    2.80099 Å
H_H2S-H_H2S ϵ =   50.00000 K, σ =    2.50000 Å
   Si-   Si ϵ =  202.29432 K, σ =    3.82641 Å
   Ga-   Ga ϵ =  208.83618 K, σ =    3.90481 Å
   Au-   Au ϵ =   19.62557 K, σ =    2.93373 Å
   Mg-   Mg ϵ =   55.85739 K, σ =    2.69141 Å
    K-    K ϵ =   17.61269 K, σ =    3.39611 Å
   Ag-   Ag ϵ =   18.11591 K, σ =    2.80455 Å
   Sc-   Sc ϵ =    9.56117 K, σ =    2.93551 Å
   Ge-   Ge ϵ =  190.72027 K, σ =    3.81305 Å
   Nb-   Nb ϵ =   29.68996 K, σ =    2.81969 Å
   Ce-   Ce ϵ =    6.54186 K, σ =    3.16804 Å
   Pm-   Pm ϵ =    4.52898 K, σ =    3.16002 Å
   Pd-   Pd ϵ =   24.15455 K, σ =    2.58272 Å
   Dy-   Dy ϵ =    3.52254 K, σ =    3.05400 Å
   Sr-   Sr ϵ =  118.25663 K, σ =    3.24376 Å
   Ho-   Ho ϵ =    3.52254 K, σ =    3.03707 Å
   No-   No ϵ =    5.53542 K, σ =    2.89364 Å
O_CO2-O_CO2 ϵ =   79.00000 K, σ =    3.05000 Å
   Sm-   Sm ϵ =    4.02576 K, σ =    3.13596 Å
   Br-   Br ϵ =  126.30814 K, σ =    3.73197 Å
   Pb-   Pb ϵ =  333.63466 K, σ =    3.82819 Å
   Xe-   Xe ϵ =  167.06894 K, σ =    3.92352 Å
   Np-   Np ϵ =    9.56117 K, σ =    3.05044 Å
    P-    P ϵ =  153.48201 K, σ =    3.69456 Å
   La-   La ϵ =    8.55473 K, σ =    3.13775 Å
   Md-   Md ϵ =    5.53542 K, σ =    2.91680 Å
    H-    H ϵ =   22.14167 K, σ =    2.57113 Å
   Pr-   Pr ϵ =    5.03220 K, σ =    3.21258 Å
   He-   He ϵ =   28.18319 K, σ =    2.10430 Å

This also prints all of the atoms included in the loaded forcefield with their given ϵ and σ. This was excluded because it would use too much space on this page.

We can access attributes LJForceField such as pure_σ, pure_ϵ, and interaction values:

# access the Lennard-Jones epsilon & sigma for Xe 
ljforcefield.pure_ϵ[:Xe] # K
ljforcefield.pure_σ[:Xe] # Å

# access the Lennard-Jones epsilon & sigma for Xe-C interactions
ljforcefield.ϵ[:Xe][:C]  # K
ljforcefield.σ²[:Xe][:C] # Å (store σ² for faster computation)

Checking Force Field Coverage

When running simulations, it is necessary to have the force field terms for all of the atoms. This can be checked using forcefield_coverage:

# check is the atoms in a crystal are covered
xtal = Crystal("SBMOF-1.cif")
forcefield_coverage(xtal, ljforcefield)

# check if the atoms in a molecule are covered
molecule = Molecule("CO2")
forcefield_coverage(molecule, ljforcefield)

Simulation Box and the Cutoff Radius

Find the replication factors needed to make a supercell big enough to fit a sphere with the specified cutoff radius.In PorousMaterials.jl, rather than replicating the atoms in the home unit cell to build the supercell that serves as a simulation box, we replicate the home unit cell to form the supercell (simulation box) in a for loop.The replication_factors function ensures enough replication factors such that the nearest image convention can be applied.

r_cutoff = 14.0 # Å
repfactors = replication_factors(xtal.box, r_cutoff) 
# output
(3, 6, 2)

Potential Energies: Van der Waals

What is the van der Waals potential energy of a Xe adsorbate inside SBMOF-1 at Cartesian coordinates [0.0, 1.0, 3.0] using the UFF as a molecular model?

# load molecule and convert it to fractional
molecule = Molecule("Xe")
molecule = Frac(molecule, xtal.box) 
translate_to!(molecule, Cart([0.0, 1.0, 3.0]), xtal.box) # need box b/c we're in Cartesian
energy = vdw_energy(xtal, molecule, ljforcefield) # K
# output
5.73882798944654e6

detailed docs

Forcefields

PorousMaterials.LJForceField — Type

Data structure for a Lennard Jones forcefield.

Attributes

name::String: name of forcefield; correponds to filename
pure_σ::Dict{Symbol, Float64}: Dictionary that returns Lennard-Jones σ of an X-X interaction, where X is an atom. (units: Angstrom)
pure_ϵ::Dict{Symbol, Float64}: Dictionary that returns Lennard-Jones ϵ of an X-X interaction, where X is an atom. (units: K)
σ²::Dict{Symbol, Dict{Symbol, Float64}}: Lennard Jones σ² (units: Angstrom²) for cross-interactions. Example use is sigmas_squared[:He][:C]
ϵ::Dict{Symbol, Dict{Symbol, Float64}}: Lennard Jones ϵ (units: K) for cross-interactions. Example use is epsilons[:He][:C]
r²_cutoff::Float64: The square of the cut-off radius beyond which we define the potential energy to be zero (units: Angstrom²). We store σ² to speed up computations, which involve σ², not σ.

source

PorousMaterials.replication_factors — Function

repfactors = replication_factors(unitcell, r_cutoff)

Find the replication factors needed to make a supercell big enough to fit a sphere with the specified cutoff radius. In PorousMaterials.jl, rather than replicating the atoms in the home unit cell to build the supercell that serves as a simulation box, we replicate the home unit cell to form the supercell (simulation box) in a for loop. This function ensures enough replication factors such that the nearest image convention can be applied.

A non-replicated supercell has 1 as the replication factor in each dimension (repfactors = (1, 1, 1)).

Arguments

unitcell::Box: The unit cell of the crystal
r_cutoff::Float64: Cutoff radius beyond which we define the potential energy to be zero (units: Angstrom)

Returns

repfactors::Tuple{Int, Int, Int}: The replication factors in the a, b, c directions

source

PorousMaterials.forcefield_coverage — Function

forcefield_coverage(atoms, ljforcefield)
forcefield_coverage(molecule, ljforcefield)
forcefield_coverage(crystal, ljforcefield)

Check that the force field contains parameters for every species in atoms::Atoms. Will print out which atoms are missing.

Arguments

atoms::Atoms: a set of atoms
ljforcefield::LJForceField: A Lennard Jones forcefield object containing information on atom interactions

Returns

all_covered::Bool: returns true if all species in the atoms are covered by the force field.

source

Potential Energy

Missing docstring.

Missing docstring for PotentialEnergy. Check Documenter's build log for details.

Missing docstring.

Missing docstring for SystemPotentialEnergy. Check Documenter's build log for details.

Nearest Image Conventions

Missing docstring.

Missing docstring for nearest_image!. Check Documenter's build log for details.

Electrostatics Energy

PorousMaterials.Eikr — Type

eikr = Eikr(eikar, eikbr, eikcr)

mutable struct for holding the eikr vectors

Attributes

eikar::OffsetArray{Complex{Float64}}: array for storing e^{i * ka ⋅ r}; has indices 0:kreps[1] and corresponds to recip. vectors in a-direction
eikbr::OffsetArray{Complex{Float64}}: array for storing e^{i * kb ⋅ r}; has indices -kreps[2]:kreps[2] and corresponds to recip. vectors in b-direction
eikcr::OffsetArray{Complex{Float64}}: array for storing e^{i * kc ⋅ r}; has indices -kreps[2]:kreps[1] and corresponds to recip. vectors in c-direction

source

Missing docstring.

Missing docstring for total. Check Documenter's build log for details.

Missing docstring.

Missing docstring for electrostatic_potential. Check Documenter's build log for details.

PorousMaterials.electrostatic_potential_energy — Function

ϕ = electrostatic_potential_energy(crystal, molecule, eparams, eikr)

Compute the electrostatic potential energy of a molecule inside a crystal.

The electrostatic potential is created by the point charges assigned to the crystal atoms in crystal.charges. Periodic boundary conditions are applied through the Ewald summation. The spurious self-interaction term is neglected here because we are looking at differences in energy in a Monte Carlo simulation.

Warning: it is assumed that the crystal is replicated enough such that the nearest image convention can be applied for the short-range cutoff radius supplied in eparams.sr_cutoff_r.

Arguments

crystal::Crystal: Crystal structure (see crystal.charges for charges)
molecule::Molecule: The molecule being compared to the atoms in the crystal.
eparams::EwaldParams: data structure containing Ewald summation settings
eikr::Eikr: Stores the eikar, eikbr, and eikcr OffsetArrays used in this calculation.

Returns

pot::EwaldSum: Electrostatic potential between crystal and molecule (units: K)

source

ϕ = electrostatic_potential_energy(molecules, eparams, box, eikr)

Compute the electrostatic potential energy of a system comprised of an array of Molecules.

The EWald summation is used here in a double for loop; do not use this function for Monte Carlo simulations because it is computationally expensive.

Returns an EwaldSum type containing short-range and long-range contributions to the Ewald sum as well as the spurious self-interaction and intramolecular interactions. Access via (ϕ.sr, ϕ.lr, ϕ.self, ϕ.intra).

Units of energy: Kelvin

Arguments

molecules::Array{Molecules, 1}: array of molecules comprising the system.
eparams::EwaldParams: data structure containing Ewald summation settings
box::Box: the box the energy is being computed in
eikr::Eikr: Stores the eikar, eikbr, and eikcr OffsetArrays used in this calculation.

Returns

ϕ::GGEwaldSum: The total electrostatic potential energy

source

PorousMaterials.precompute_kvec_wts — Function

kvectors = precompute_kvec_wts(kreps, box, α, max_mag_k_sqrd=Inf)

For speed, pre-compute the weights for each reciprocal lattice vector for the Ewald sum in Fourier space. This function takes advantage of the symmetry: cos(-k⋅(x-xᵢ)) + cos(k⋅(x-xᵢ)) = 2 cos(k⋅(x-xᵢ))

If max_mag_k_sqrd is passed, k-vectors with a magnitude greater than max_mag_k_sqrd are not included.

Arguments

kreps::Tuple{Int, Int, Int}: number of k-vector replications required in a, b, c
box::Box: the simulation box containing the reciprocal lattice.
α::Float64: Ewald sum convergence parameter (units: inverse Å)
max_mag_k_sqrd::Float64: cutoff for |k|² in Fourier sum; if passed, do not include

k-vectors with magnitude squared greater than this.

Returns

kvectors::Array{Kvector, 1}: array of k-vectors to include in the Fourier sum and their

corresponding weights indicating the contribution to the Fourier sum.

source

PorousMaterials.setup_Ewald_sum — Function

eparams = setup_Ewald_sum(box, sr_cutoff_r; ϵ=1e-6, verbose=false)

Given the short-range cutoff radius and simulation box, automatically compute Ewald convergence parameter and number of k-vector replications in Fourier space required for a given precision. Constructs and returns Ewald parameters data type with this information.

Also, pre-compute weights on k-vector contributions to Ewald sum in Fourier space.

Also, allocate OffsetArrays for storing e^{i * k ⋅ r} where r = x - xⱼ and k is a reciprocal lattice vector.

Arguments

box::Box: the simulation box containing the reciprocal lattice.
sr_cutoff_r::Float64: cutoff-radius (units: Å) for short-range contributions to Ewald
ϵ::Float64: desired level of precision. Typical value is 1e-6, but this does not
verbose::Bool: If true will print results

Returns

eparams::EwaldParams: data structure containing Ewald summation settings

corresponding weights indicating the contribution to the Fourier sum.

source

PorousMaterials.total_electrostatic_potential_energy — Function

total_ϕ = total_electrostatic_potential_energy(molecules, eparams, box, eikr)

Calculates the total electrostatic potential energy of an array of Molecules using a Grand Canonical Monte Carlo (GCMC) algorithm. #TODO add to this

Arguments

molecules::Array{Molecule, 1}: The molecules comprising the system.
eparams::EwaldParams: data structure containing Ewald summation settings
box::Box: The box the energy is being computed in.
eikr::Eikr: Stores the eikar, eikbr, and eikcr OffsetArrays used in this calculation.

Returns

ϕ::GGEwaldSum: The total electrostatic potential energy

source

total_ϕ = total_electrostatic_potential_energy(crystal, molecules, eparams, eikr)

Explanation of totalelectrostaticpotential_energy that uses crystal

Arguments

crystal::Crystal: Crystal structure (see crystal.charges for charges)
molecules::Array{Molecule, 1}: The molecules comprising the system.
eparams::EwaldParams: data structure containing Ewald summation settings
eikr::Eikr: Stores the eikar, eikbr, and eikcr OffsetArrays used in this calculation.

source

Van der Waals Energy

PorousMaterials.lennard_jones — Function

energy = lennard_jones(r², σ², ϵ)  (units: Kelvin)

Calculate the lennard jones potential energy given the square of the radius r between two lennard-jones spheres. σ and ϵ are specific to interaction between two elements. Return the potential energy in units Kelvin (well, whatever the units of ϵ are).

Arguments

r²::Float64: distance between two (pseudo)atoms in question squared (Angstrom²)
σ²::Float64: sigma parameter in Lennard Jones potential squared (units: Angstrom²)
ϵ::Float64: epsilon parameter in Lennard Jones potential (units: Kelvin)

Returns

energy::Float64: Lennard Jones potential energy

source

PorousMaterials.vdw_energy — Function

energy = vdw_energy(crystal, molecule, ljforcefield)

Calculates the van der Waals interaction energy between a molecule and a crystal. Applies the nearest image convention to find the closest replicate of a specific atom.

WARNING: it is assumed that the framework is replicated sufficiently such that the nearest image convention can be applied. See replicate and replication_factors.

source

ggenergy = vdwenergy(moleculeid, molecules, ljforcefield, simulationbox)

Calculates van der Waals interaction energy of a single adsorbate molecules[molecule_id] with all of the other molecules in the system. Periodic boundary conditions are applied, using the nearest image convention.

Arguments

molecule_id::Int: Molecule ID used to determine which molecule in molecules we wish to calculate the guest-guest interactions
molecules::Array{Molecule, 1}: An array of Molecule data structures
ljforcefield::LJForceField: A Lennard Jones forcefield data structure describing the interactions between different atoms
simulation_box::Box: The simulation box for the computation.

Returns

gg_energy::Float64: The guest-guest interaction energy of molecules[molecule_id] with the other molecules in molecules

source

PorousMaterials.vdw_energy_no_PBC — Function

potenergy = vdwenergynoPBC(atomsi, atomsj , ljff)

compute vdw potential energy without periodic boundary conditions

source