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.LJForceFieldType

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_factorsFunction
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_coverageFunction
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.EikrType
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_energyFunction
ϕ = 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_wtsFunction
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_sumFunction
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_energyFunction
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_jonesFunction
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_energyFunction

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