Applying Interactions

Click to focus. When focused, press m to open the options menu. Refresh the page to restart the simulation.

Overview

The previous tutorials set the Hamiltonian $H=0$. In this one, let’s place particles in an external gravitational field and add a pairwise interaction that penalizes particle overlaps:

$$H=\sum _i\alpha {\vec{r}}_i\cdot \hat{y}+\sum _i\sum _{j>i}{U}_{\mathrm{step}}\left(\left|{\vec{r}}_j-{\vec{r}}_i\right|\right)$$ where $$U_{\mathrm{step}}(r)=\{\begin{array}{ll}\epsilon & r<\sigma \\ 0& r\ge \sigma \end{array}$$

• Objective: Demonstrate the use of external and pairwise potentials in MC simulations.
• File: hoomd-rs/examples/mc-tutorial/applying-interactions.rs
• Run (interactively):

  cargo run --release --features "bevy" --example applying-interactions
  

• Run (in batch mode):

  cargo run --release --example applying-interactions
  

Bodies and Sites

The previous tutorials introduced the microstate as a collection of point bodies. That was an oversimplification. In hoomd-rs, a body is an ordered collection of sites that are defined in the local reference frame of the body. The microstate gains all of its degrees of freedom from the bodies it contains.

You get to choose the types of the body properties and the site properties in your simulation model. At a minimum, both must have position. The body properties (denoted as the B generic) and the site properties (denoted by the S generic) may be the same or different. However, the body properties type must be able to transform a site properties from the body reference frame to the system reference frame.

All interactions on bodies are a function only of its sites and are computed in the system reference frame. Understanding this will help as you review the documentation for the types used later in this tutorial: External, Linear, Boxcar, Isotropic, and PairwiseCutoff. For a complete reference on bodies, sites, and all their related traits, read the hoomd-microstate documentation.

In this tutorial, the bodies will still be points. Specifically, that means each body has Point<Cartesian<2>> for its body properties type (B), and a single site at the origin (in the body reference frame) which has Point<Cartesian<2>> for its site properties (S) type. The constructor Body::point() used in these tutorials is a convenient way to create point bodies.

The next tutorial will demonstrate a case where bodies and sites have different properties and show how to create bodies with more than one site.

Use Declarations

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;

The Simulation Model

Here is the type that holds the simulation model:

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;

In order, Microstate’s generic types are the body properties, the site properties, and the boundary condition.

Construct the Simulation Model

The new() method constructs a new simulation model:

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;

Parameters

Assign all the model parameters in one code block:

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;

box_length is the side length of the square simulation box, maximum_distance is the largest distance a translation trial move can take, alpha is the strength of the gravitational potential, epsilon is the strength of the pairwise potential, sigma is the range of the pairwise potential, and macrostate holds the temperature set point (in units of energy).

External Potential

This code implements the external potential term in the Hamiltonian: $$\sum _i\alpha {\vec{r}}_i\cdot \hat{y}$$

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;

Linear computes $\alpha \vec{r}\cdot \hat{y}$ in its energy() method. Linear also implements the SiteEnergy trait whose method site_energy takes a single site properties argument: $U(s)$. Building on that, External wraps any type that implements SiteEnergy and sums over the energy contributed by each site in the microstate: ${\sum }_{i}U(s_i)$. External implements the DeltaEnergyOne trait which Sweep will use to evaluate the change in energy $\Delta E$ of a trial that moves one body.

Pairwise Potential

This code implements the pairwise potential term in the Hamiltonian: $$\sum _i\sum _{j>i}{U}_{\mathrm{step}}\left(\left|{\vec{r}}_j-{\vec{r}}_i\right|\right)$$

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;

The Boxcar function implements $U_{\mathrm{step}}(r)$ via the UnivariateEnergy trait. $$U_{\mathrm{step}}(r)=\{\begin{array}{ll}\epsilon & r<\sigma \\ 0& r\ge \sigma \end{array}$$ Isotropic is a wrapper that computes $U(\left|{\vec{r}}_j-{\vec{r}}_i\right|)\left[\left|{\vec{r}}_j-{\vec{r}}_i\right|<r_{\mathrm{cut}}\right]$ in its implementation of SitePairEnergy. The site_pair_energy() method is a more general form that depends on the full set of properties of the two interacting sites: $U(s_i,s_j)$. The PairwiseCutoff type sums over all pairs of sites that are within a distance of $r_{\mathrm{cut}}$ and do not belong to the same body: $$\sum _i\sum _{j>i}U\left(s_i,s_j\right)\left[b_i\ne b_j\right]$$ Finally, PairwiseCutoff implements the DeltaEnergyOne trait which Sweep will use to evaluate the change in energy $\Delta E$ of a trial move.

Important

In hoomd-rs, it is YOUR responsibility to determine the appropriate maximum $r_{\mathrm{cut}}$ for your choice of SitePairEnergy. You might be used to other simulation codes, HOOMD-blue for example, that automatically determine this maximum for you based on the parameters of the inner types. That is not possible in hoomd-rs as your site_pair_energy could be any arbitrary code.

The Hamiltonian

To sum the external and pair energies, place them in a struct: $$H=U_{\mathrm{external}}+U_{\mathrm{pair}}$$

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;

See The Hamiltonian Type for the definition of this struct and instructions on how to derive implementations of traits like DeltaEnergyOne as a sum over the struct’s fields. In this example, translate_sweep.apply() calls hamiltonian.delta_energy_one() to evaluate $\Delta E$ when needed.

You can use hamiltonian to compute properties of the system:

• hamiltonian.total_energy(&microstate) - The total energy of the system.
• hamiltonian.linear.total_energy(&microstate) - The total external energy term.
• hamiltonian.linear.0.site_energy(&site.properties) - The contribution of a single site to the external energy.
• hamiltonian.pairwise_cutoff.total_energy(&microstate) - The total pair energy term.
• hamiltonian.pairwise_cutoff.site_pair_energy(&site_i, &site_j) - The contribution of a pair of sites to the pair energy.

The types External and Isotropic are single element tuples. To access the parameters of the inner types, you need access the elements of these tuples:

• hamiltonian.linear.0.alpha - Strength of the linear external potential.
• hamiltonian.pairwise_cutoff.0.r_cut - Maximum cutoff radius of of the pair potential.
• hamiltonian.pairwise_cutoff.0.interaction.epsilon - Strength of the pairwise step potential.

Due to Rust’s ownership model, you cannot use names like boxcar.epsilon to refer to parameters after constructing hamiltonian. You can read more about ownership in The Rust Programming Language.

Trial Moves

Apply translation trial moves to the bodies:

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;

Microstate

Confine the bodies and sites inside of a closed square. While the previous tutorial showed how you could implement custom boundary conditions, this one uses the built in Rectangle type:

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;

Construct the VecCell spatial data structure and pass it to the Microstate when using PairwiseCutoff. PairwiseCutoff will use the spatial data structure to efficiently compute the pairwise interactions.

Important

If you forget this step, your simulation will run very slowly!

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;

Set the nominal search radius to the same value as the r_cut used with PairwiseCutoff. To aid you in choosing the right search radius, the maximum_interaction_range() computes the maximum interaction range needed for all terms in the Hamiltonian.

Construct the Microstate with the square boundary and vec_cell spatial data:

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;

Initialize the Struct

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;

The Hamiltonian Type

The Hamiltonian of the system must be represented by a single type. To sum terms from multiple types (External and PairwiseCutoff here), define a struct with one field for each term in the Hamiltonian:

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;

derive([TotalEnergy, DeltaEnergyOne]) implements those traits for the Hamiltonian type where the result is the sum of the corresponding trait over all the struct’s fields. derive([MaximumInteractionRange)] similarly determines the maximum of all the maximum interaction ranges over all the struct’s fields.

Tip

When there is only one term in your Hamiltonian (e.g. pairwise_cutoff), you can use it directly as hamiltonian without the need for a new struct.

Implement Simulation

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;

Advance the Simulation

The advance() method moves the simulation forward one step:

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;

Add New Bodies

Every 100 steps, add a new body near the top of the simulation box:

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;

Apply Trial Moves

Attempt one translation trial move for each body in the microstate:

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;

The previously unused temperature $kT$ now has meaning in this simulation model. A pair of overlapping disks in this model results in $U=1000kT$. The probability of accepting a trial move that adds an overlap is $e^{\frac{\Delta E}{kT}}=e^{-1000}$ which is identically $0$ in f64 arithmetic. Therefore, translate_sweep.apply() will never add new overlaps.

However, the unconditional add_body() above can place overlapping bodies. When a pair of overlapping disks is placed, translate_sweep.apply() will accept trial moves that keep the same number of overlaps because $\Delta E=1000-1000=0$.

Reset the Simulation

Eventually, the boundary will completely fill with particles. The overlapping disks are visually disconcerting in the interactive example, so let’s reset it when needed.

The pair potential in this example adds 1000 for every pair of disks that overlap. Remove all bodies from the microstate when total pairwise energy exceeds a threshold:

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;

Get the Simulation Step

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;

Execute the Simulation in Batch Mode

Visual, interactive simulations are great for teaching and quickly testing new models. They are less practical when you scale your workflow to run thousands of simulations on HPC resources. For that, you need to run in batch mode and write the simulation results to one or more files.

Construct the Simulation Model

To run the simulation, construct the Fill simulation model:

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;

Create a GSD Trajectory

A GSD file stores the positions of the sites in frames. Each frame is a snapshot of the simulation state at a specific step in the simulation. Start by creating a new GSD file:

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;

Advance the Simulation

Call advance() many times in the main loop:

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;

Write Frames to the GSD File

Call append_microstate to write to the GSD file periodically:

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;

If you would like to monitor the simulation progress, you can use println! (or log and env_logger) to write any output you like.

Note

This main() function runs in batch mode. There is a different main() (not shown here) used in the interactive example. The interactive example does not write the GSD file.

Run the Simulation

In a terminal, execute the following command to run the simulation in batch mode:

cargo run --release --example applying-interactions

Visualize the Simulation Results

Open the generated applying-interactions.gsd in Ovito or another visualization tool to see the simulation results. Ovito will render the disks as spheres with the expected diameter of 1 by default.

Render with Tachyon, and you should see something like:

Conclusion

This tutorial showed you how to define interactions in your simulation model via the Hamiltonian.

Navigate to the top of the page and refresh to see the simulation in action again. Notice how the disks fall to the bottom of the boundary and do not overlap, except when newly added. Wait long enough and you will see the simulation clear the bodies.

Complete Code

use hoomd_geometry::shape::Rectangle;
use hoomd_interaction::{
    DeltaEnergyOne, External, MaximumInteractionRange, PairwiseCutoff,
    TotalEnergy, external::Linear, pairwise::Isotropic, univariate::Boxcar,
};
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;

// Remove the cfg_attr(...) line when using this code outside the hoomd-rs/examples directory.
#[cfg_attr(feature = "bevy", derive(Resource))]
struct Fill {
    /// Positions of all the bodies in the simulation.
    microstate: Microstate<
        Point<Cartesian<2>>,
        Point<Cartesian<2>>,
        VecCell<SiteKey, 2>,
        Closed<Rectangle>,
    >,
    /// How sites interact with other sites and fields.
    hamiltonian: Hamiltonian,
    /// Trial moves to apply.
    translate_sweep: Sweep<Translate<Cartesian<2>>>,
    /// Temperature set point.
    macrostate: Isothermal,
}

impl Fill {
    /// Construct a new fill simulation.
    fn new() -> anyhow::Result<Fill> {
        let box_length = 30.0;
        let maximum_distance = 0.15;
        let alpha = 10.0;
        let epsilon = 1000.0;
        let sigma = 1.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let linear = External(Linear {
            alpha,
            plane_origin: Cartesian::default(),
            plane_normal: [0.0, 1.0].try_into()?,
        });

        let boxcar = Boxcar {
            epsilon,
            left: 0.0,
            right: sigma,
        };
        let pairwise_cutoff = PairwiseCutoff(Isotropic {
            interaction: boxcar,
            r_cut: sigma,
        });

        let hamiltonian = Hamiltonian {
            linear,
            pairwise_cutoff,
        };

        let translate =
            Translate::with_maximum_distance(maximum_distance.try_into()?);
        let translate_sweep = Sweep(translate);

        let square = Rectangle::with_equal_edges(box_length.try_into()?);
        let vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        Ok(Fill {
            microstate,
            hamiltonian,
            translate_sweep,
            macrostate,
        })
    }
}

#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange)]
struct Hamiltonian {
    linear: External<Linear<Cartesian<2>>>,
    pairwise_cutoff: PairwiseCutoff<Isotropic<Boxcar>>,
}

impl Simulation for Fill {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        let boundary = self.microstate.boundary();
        let y = boundary.0.edge_lengths[1].get() / 2.0 - 0.5;
        if self.microstate.step().is_multiple_of(100) {
            self.microstate.add_body(Body::point([0.0, y].into()))?;
        }

        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();

        if self.hamiltonian.linear.total_energy(&self.microstate) > 20_000.0 {
            self.microstate.clear();
        }

        Ok(())
    }

    /// Get the current simulation step.
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory.
#[cfg(not(feature = "bevy"))]
fn main() -> anyhow::Result<()> {
    use hoomd_gsd::hoomd::HoomdGsdFile;
    use hoomd_microstate::AppendMicrostate;

    let mut simulation = Fill::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("applying-interactions.gsd")?;

    for _ in 0..100_000 {
        simulation.advance()?;

        if simulation.step().is_multiple_of(1_000) {
            hoomd_gsd_file
                .append_microstate(&simulation.microstate)?
                .end()?;
        }
    }

    Ok(())
}

#[cfg(feature = "bevy")]
mod applying_interactions_interactive;
#[cfg(feature = "bevy")]
use applying_interactions_interactive::main;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;