Examples

Below we provide a few simple examples of short Monte Carlo simulations with MoSDeF Cassandra.

NVT simulation of methane

import mbuild
import foyer
import mosdef_cassandra as mc
import unyt as u

# Use mBuild to create a methane molecule
methane = mbuild.load("C", smiles=True)

# Create an empty mbuild.Box
box = mbuild.Box(lengths=[3.0, 3.0, 3.0])

# Load force field
oplsaa = foyer.forcefields.load_OPLSAA()

# Use foyer to apply force field to methane
methane_ff = oplsaa.apply(methane)

# Create box and species list
box_list = [box]
species_list = [methane_ff]

# Use Cassandra to insert some initial number of methane molecules
mols_to_add = [[50]]

# Define the System
system = mc.System(box_list, species_list, mols_to_add=mols_to_add)
# Define the MoveSet
moveset = mc.MoveSet("nvt", species_list)

# Run a simulation at 300 K for 10000 MC moves
mc.run(
    system=system,
    moveset=moveset,
    run_type="equilibration",
    run_length=10000,
    temperature=300.0 * u.K,
)

NPT simulation of methane

import mbuild
import foyer
import mosdef_cassandra as mc
import unyt as u

# Use mbuild to create molecules
methane = mbuild.load("C", smiles=True)

# Create an empty mbuild.Box
box = mbuild.Box(lengths=[3.0, 3.0, 3.0])

# Load force field
oplsaa = foyer.forcefields.load_OPLSAA()

# Use foyer to apply force field
methane_ff = oplsaa.apply(methane)

# Create box and species list
box_list = [box]
species_list = [methane_ff]

# Use Cassandra to insert some initial number of species
mols_to_add = [[5]]

# Define the System
system = mc.System(box_list, species_list, mols_to_add=mols_to_add)
# Define the MoveSet
moveset = mc.MoveSet("npt", species_list)

# Here we must specify the pressure since we are performing a
# NpT simulation. It can be provided in the custom_args dictionary
# or as a keyword argument to the "run" function.
custom_args = {
    "pressure": 1.0 * u.bar,
}

# Run a simulation with at 300 K with 10000 MC moves
mc.run(
    system=system,
    moveset=moveset,
    run_type="equilibration",
    run_length=10000,
    temperature=300.0 * u.K,
    **custom_args,
)

NVT simulation of methane and propane mixture

import mbuild
import foyer
import mosdef_cassandra as mc
import unyt as u

# Use mbuild to create methane and propane molecules
methane = mbuild.load("C", smiles=True)
propane = mbuild.load("CCC", smiles=True)

# Create an empty mbuild.Box
box = mbuild.Box(lengths=[3.0, 3.0, 3.0])

# Load force field
oplsaa = foyer.forcefields.load_OPLSAA()

# Use foyer to apply the force field
typed_methane = oplsaa.apply(methane)
typed_propane = oplsaa.apply(propane)

# Create box and species list
box_list = [box]
species_list = [typed_methane, typed_propane]

# Use Cassandra to insert some initial number of species
mols_to_add = [[100, 50]]

system = mc.System(box_list, species_list, mols_to_add=mols_to_add)
moveset = mc.MoveSet("nvt", species_list)

mc.run(
    system=system,
    moveset=moveset,
    run_type="equilibration",
    run_length=10000,
    temperature=200.0 * u.K,
)

GEMC simulation of methane (united atom)

import mbuild
import foyer
import mosdef_cassandra as mc
import unyt as u

# Use mbuild to create a coarse-grained CH4 bead
methane = mbuild.Compound(name="_CH4")

# Create two empty mbuild.Box
# (vapor = larger, liquid = smaller)
liquid_box = mbuild.Box(lengths=[3.0, 3.0, 3.0])
vapor_box = mbuild.Box(lengths=[4.0, 4.0, 4.0])

# Load force field
trappe = foyer.forcefields.load_TRAPPE_UA()

# Use foyer to apply force field
typed_methane = trappe.apply(methane)

# Create box and species list
box_list = [liquid_box, vapor_box]
species_list = [typed_methane]

mols_to_add = [[350], [100]]

system = mc.System(box_list, species_list, mols_to_add=mols_to_add)
moveset = mc.MoveSet("gemc", species_list)

moveset.prob_volume = 0.010
moveset.prob_swap = 0.11

thermo_props = [
    "energy_total",
    "energy_intervdw",
    "pressure",
    "volume",
    "nmols",
    "mass_density",
]

custom_args = {
    "run_name": "equil",
    "charge_style": "none",
    "rcut_min": 2.0 * u.angstrom,
    "vdw_cutoff": 14.0 * u.angstrom,
    "units": "sweeps",
    "steps_per_sweep": 450,
    "coord_freq": 50,
    "prop_freq": 10,
    "properties": thermo_props,
}

mc.run(
    system=system,
    moveset=moveset,
    run_type="equilibration",
    run_length=250,
    temperature=151.0 * u.K,
    **custom_args,
)

# Update run_name and restart_name
custom_args["run_name"] = "prod"
custom_args["restart_name"] = "equil"

mc.restart(
    system=system,
    moveset=moveset,
    run_type="production",
    run_length=750,
    temperature=151.0 * u.K,
    **custom_args,
)

GCMC simulation of methane

import mbuild
import foyer
import mosdef_cassandra as mc
import unyt as u

# Use mbuild to create a methane
methane = mbuild.load("C", smiles=True)

# Create an empty mbuild.Box
box = mbuild.Box(lengths=[10.0, 10.0, 10.0])

# Load force field
oplsaa = foyer.forcefields.load_OPLSAA()

# Use foyer to apply the force field
methane_ff = oplsaa.apply(methane)

# Create box and species list
box_list = [box]
species_list = [methane_ff]

mols_to_add = [[100]]

system = mc.System(box_list, species_list, mols_to_add=mols_to_add)
moveset = mc.MoveSet("gcmc", species_list)

custom_args = {
    "chemical_potentials": [-35.0 * (u.kJ / u.mol)],
    "prop_freq": 100,
}

mc.run(
    system=system,
    moveset=moveset,
    run_type="equilibration",
    run_length=1000,
    temperature=300.0 * u.K,
    **custom_args,
)

GCMC simulation of methane adsorption in a solid framework

import mbuild
import foyer
import mosdef_cassandra as mc
import unyt as u

from mosdef_cassandra.examples.structures import carbon_lattice


# Load a structure created with mbuild
lattice = carbon_lattice()
# Use mbuild to create a methane
methane = mbuild.load("C", smiles=True)

# Load force field
trappe = foyer.forcefields.load_TRAPPE_UA()
oplsaa = foyer.forcefields.load_OPLSAA()

# Use foyer to apply the force fields
typed_lattice = trappe.apply(lattice)
methane_ff = oplsaa.apply(methane)

# Create box and species list
box_list = [lattice]
species_list = [typed_lattice, methane_ff]

# Since we have an occupied box we need to specify
# the number of each species present in the initial config
mols_in_boxes = [[1, 0]]

system = mc.System(box_list, species_list, mols_in_boxes=mols_in_boxes)
moveset = mc.MoveSet("gcmc", species_list)

custom_args = {
    "chemical_potentials": ["none", -30.0 * (u.kJ / u.mol)],
    "rcut_min": 0.5 * u.angstrom,
    "vdw_cutoff": 14.0 * u.angstrom,
    "charge_cutoff": 14.0 * u.angstrom,
    "coord_freq": 100,
    "prop_freq": 10,
}

mc.run(
    system=system,
    moveset=moveset,
    run_type="equilibration",
    run_length=10000,
    temperature=300.0 * u.K,
    **custom_args,
)

NVT simulation of SPC/E water

import mbuild
import foyer
import mosdef_cassandra as mc
import unyt as u
from mosdef_cassandra.utils.get_files import get_example_ff_path, get_example_mol2_path

# Load water with SPC/E geometry from mol2 file
molecule = mbuild.load(get_example_mol2_path("spce"))

# Create an empty mbuild.Box
box = mbuild.Box(lengths=[3.0, 3.0, 3.0])

# Load force field
spce = foyer.Forcefield(get_example_ff_path("spce"))

# Use foyer to apply force field
molecule_ff = spce.apply(molecule)

# Create box and species list
box_list = [box]
species_list = [molecule_ff]

# Use Cassandra to insert some initial number of species
mols_to_add = [[50]]

# Define the System
system = mc.System(box_list, species_list, mols_to_add=mols_to_add)
# Define the MoveSet
moveset = mc.MoveSet("nvt", species_list)

# Note here we need to use the angle_style="fixed" keyword argument
# SPC/E geometry is rigid; default angle style is "harmonic"
custom_args = {"angle_style": ["fixed"]}

# Run a simulation with at 300 K with 10000 MC moveset
mc.run(
    system=system,
    moveset=moveset,
    run_type="equilibration",
    run_length=10000,
    temperature=300.0 * u.K,
    **custom_args,
)