Examples

Below are some examples of how to use msibi. These are examples of the API only, and all files shown are place holders. See the Tutorials for instructions on accessing a minimal provided dataset needed to work through an msibi tutorial.

Quick Example: Optimizing bond-stretching

Here is a simple example using MSIBI to learn a bond-stretching force from a single state point:

import hoomd
from msibi import MSIBI, State, Bond

# This is the context/management class for MSIBI
# Set simulation parameters, call `add_state` and `add_force` methods to store other MSIBI objects.
optimizer = MSIBI(
  nlist=hoomd.md.nlist.Cell,
  integrator_method=hoomd.md.methods.ConstantVolume,
  thermostat=hoomd.md.methods.thermostats.MTTK,
  thermostat_kwargs={"tau": 0.1},
  method_kwargs={},
  dt=0.0001,
  gsd_period=int(1e4)
)

# Create a State instance, pass in a path to the target trajectory
stateA = State(name="A", kT=5.0, traj_file="cg_trajectory.gsd", alpha0=0.7, n_frames=100)

# For each force you want to optimize, create an instance, set optimize=True
AA_bond = Bond(type1="A", type2="A", optimize=True, nbins=80)
AA_bond.set_polynomial(x_min=0.0, x_max=0.5, x0=0.22, k2=5000, k3=0, k4=0)
AB_bond = Bond(type1="A", type2="B", optimize=True, nbins=80)
AB_bond.set_polynomial(x_min=0.0, x_max=0.5, x0=0.22, k2=5000, k3=0, k4=0)
# Add all states and forces to the optimization class (MSIBI)
optimizer.add_state(stateA)
optimizer.add_force(AA_bond)
optimizer.add_force(AB_bond)
optimizer.run_optimization(n_iterations=10, n_steps=2e5)

# See distribution comparison
AA_bond.plot_distribution_comparison(state=stateA)
AB_bond.plot_distribution_comparison(state=stateA)

# Save potentials
AA_bond.save_potential("AA_bond.csv")
AB_bond.save_potential("AB_bond.csv")

Example 2: Multiple states, multiple forces

In this example, we set fixed bond and angle potentials that are included during iteration simulations.
The bond potential will set a fixed harmonic force, while the angle potential will be set from a table potential file.
This illustrates a use case of stringing together multiple MSIBI optimizations.
For example, one MSIBI optimization can be used to learn and obtain a coarse-grained angle potential table file which can then be set and held fixed while learning pair potentials in a subsequent MSIBI optimization.

import hoomd
from msibi import MSIBI, State, Pair, Bond, Angle

optimizer = MSIBI(
  nlist=hoomd.md.nlist.Cell,
  integrator_method=hoomd.md.methods.ConstantVolume,
  thermostat=hoomd.md.methods.thermostats.MTTK,
  thermostat_kwargs={"tau": 0.1},
  method_kwargs={},
  dt=0.0001,
  gsd_period=int(1e4)
)

# Create 3 State instances, pass in a path to the target trajectory
stateA = State(name="A", kT=2.0, traj_file="stateA.gsd", alpha=0.2, n_frames=50)
stateB = State(name="B", kT=4.0, traj_file="stateB.gsd", alpha=0.5, n_frames=50)
stateC = State(name="C", kT=6.0, traj_file="stateC.gsd", alpha=0.3, n_frames=50)
optimizer.add_state(stateA)
optimizer.add_state(stateB)
optimizer.add_state(stateC)

# Add bond and set a harmonic force (e.g. fit to Boltzmann inverse of the distribtion)
bondAA = Bond(type1="A", type2="A", optimize=False)
bondAA.set_harmonic(r0=1.4, k=800)
optimize.add_force(bondAA)

# Add angle and load previously obtained table potential
angleAA = Angle(type1="A", type2="A", type3="A", optimize=False)
angleAA.set_from_file("angleAA.csv")
optimize.add_force(angleAA)

# Create a Pair instance to be optimized.
pairAA = Pair(type1="A", type2="A", optimize=True, r_cut=3.0, nbins=100)
# Call the set_lj() method to set an initial guess potential
pairAA.set_lj(r_min=0.001, r_cut=3.0, epsilon=1.0, sigma=1.0)
optimizer.add_force(pairAA)

# Run 20 MSIBI iterations
optimizer.run_optimization(n_steps=2e6, n_iterations=20)
pairAA.save_potential("pairAA.csv")