Introduction

hoomd-rs is a collection of Rust crates that implement particle simulations and related methods. It performs Monte Carlo simulations of hard shapes and interacting particles (both isotropic and anisotropic) as well as molecular dynamics simulations with a variety of particle interaction. hoomd-rs provides public APIs for vector math, geometric primitives, spatial data structures, energy calculations, and all other components of the simulation that users can employ in their own analysis and simulation methods. You can use hoomd-rs to create real-time interactive visualizations of simulations, execute long-running simulations in batch mode on high performance computing resources, and analyze the results of those simulations.

hoomd-rs is the spiritual successor to the Python package HOOMD-blue. While the two share many common features, hoomd-rs provides many capabilities that HOOMD-blue cannot, such as:

• Custom per-particle attributes.
• Custom particle interactions that can depend on custom per-particle attributes.
• Custom vector representations, including curved spaces.
• Custom MC trial moves and acceptance criteria.
• Custom simulation box geometries (including non-periodic simulation boxes).
• Custom visual representations of simulation elements.
• Build native command line applications on Linux, macOS, and Windows.
• Run real-time interactive simulations on desktop platforms or embedded in a web page.

This documentation is written in mdBook. The Reading Books chapter explains how to search, change display settings, and navigate this book. Press ? to see a list of keyboard shortcuts.

Rust

The first thing you need to understand about hoomd-rs is that it is not an application or package that you install. Rather it is a collection of Rust crates that you can use to implement your simulation models.

If you are not familiar with Rust, it is a relatively new programming language that combines the best features of generic languages like Python with the best features of strongly-typed compiled languages like C++. hoomd-rs takes full advantage of these capabilities to provide a simulation framework that is fully customizable while still compiling down to machine code. The Rust Programming Language explains everything you need to know about the language itself.

Before continuing, you will need to install Rust. On Linux, Mac, and WSL you can install Rust with a single command:

$ curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh

For more details, including instructions for a native Windows installation, see the Rust installation documentation (make sure you install the 64-bit build).

Tip

hoomd-rs works very well with native builds on all platforms.

Hello, hoomd-rs!

You don’t need to install anything other than Rust to use hoomd-rs. Rust comes packaged with Cargo. Use Cargo to create a hello-hoomd-rs crate:

$ cargo new hello-hoomd-rs
$ cd hello-hoomd-rs

This example uses the hoomd-interaction crate. Add it to your dependencies with:

$ cargo add hoomd-interaction

Cargo creates a “Hello, world” application by default. Replace src/main.rs with:

use hoomd_interaction::univariate::{LennardJones, UnivariateEnergy};

fn main() {
    let lennard_jones: LennardJones = LennardJones::default();
    println!("lennard_jones(1.5): {}", lennard_jones.energy(1.5));
}

Compile and run the code:

$ cargo run

You should see a number of compiling… messages and then:

lennard_jones(1.5): -0.32033659427857464

Congratulations for making it this far! You have successfully compiled hoomd-rs and used it to compute the Lennard-Jones potential at a given r value. Read the following tutorials to learn how to perform Monte Carlo and molecular dynamics simulations using hoomd-rs.

Tip

cargo run builds in debug mode by default. Debug builds make it easier to troubleshoot problems with your code, but they run very slowly. When your code is working, build in release mode and it will run much faster:

$ cargo run --release

Note

On HPC platforms, you should run (preferably in an interactive compute job):

$ cargo build --release

You can then use the executable placed in the target/release directory in your submitted job scripts. For more details, see the workflow tutorial.

Unlike scripting languages (e.g. Python), saving changes to main.rs will NOT take effect until you run cargo build --release again.

Overview

hoomd-rs is not an application that you install and use. It is a coupled set of Rust crates. Each individual crate serves a narrow purpose and, to the extent possible, is usable on its own without the others. The collection as a whole enables Monte Carlo and molecular dynamics simulations of soft matter systems.

To use hoomd-rs, you need to write a Rust application that links to one or more hoomd-rs crates. What you do in your application is completely up to you. You can learn about what each crate does in its documentation.

Note

Skip to the Monte Carlo tutorials for an introduction to hoomd-rs based on concrete example simulations.

Interaction models

hoomd-interaction implements many interaction models commonly used in soft matter research, such as the Lennard-Jones pair, linear external, and hard shape overlap potentials. When you need something different, implement your own interaction model via one or more of the traits defined in hoomd-interaction.

See the hoomd-interaction documentation for a more detailed introduction, a complete list of all interaction models, and example code.

Representing Shapes

hoomd-geometry defines many common 2D and 3D shapes and implements a number of operations on them. You can use these shapes to define your simulation boundary and site interactions. When you need a different class of shape, implement your own via one or more of the traits defined in hoomd-geometry.

See the hoomd-geometry documentation for a more detailed introduction, a complete list of all shapes, and example code.

Tracking the System’s Microstate

hoomd-microstate defines the Microstate type which tracks all the degrees of freedom for the bodies in your simulation and provides efficient methods to access their sites. Many of the hoomd-interaction traits operate on a system microstate, such as TotalEnergy which computes the total energy of a microstate for a given interaction model.

See the hoomd-microstate documentation for a more detailed introduction, a complete list of all methods, and example code.

Monte Carlo Simulations

[hoomd-md] implements typical Monte Carlo trial moves used in soft matter research, such as local translation and rotation moves. Each trial move is applied to a Microstate (hoomd-microstate) and accepted or rejected based on the given interaction model (hoomd-interaction). When you need a different type of trial move for your simulation model, implement your own via one or more of the traits defined in hoomd-mc.

Writing GSD Files

Use hoomd-gsd to read and write GSD files that store the position and orientation of the sites in your microstate as a function of time. Use GSD files for offline analysis and visualization after your simulation job completes.

All crates

There are more hoomd-rs crates not highlighted in this user guide. Here are links to the documentation for every hoomd-rs crate:

• hoomd-bevy: Facilitate visual, interactive simulations.
• hoomd-gsd: Read and write GSD files.
• hoomd-geometry: Define geometric shapes and operate on them.
• hoomd-interaction: Define interaction models and apply them to the system’s microstate.
• hoomd-linear-algebra: Linear algebra methods used by other hoomd-rs crates.
• hoomd-manifold: Define points in curved space and operate on them.
• hoomd-mc: Perform Monte Carlo simulations.
• hoomd-microstate: Represent the degrees of freedom in your simulation model.
• hoomd-rand: Pseudorandom number generation methods used by other hoomd-rs crates.
• hoomd-simulation: Traits that define the macrostate and the simulation model as a whole.
• hoomd-spatial: Spatial data structures for efficient site-site interaction computations.
• hoomd-utility: Miscellaneous data types used throughout hoomd-rs.
• hoomd-vector: Vector and rotation data types and operations.

Running the Examples

Most of the following tutorials show the example executing in the browser. Click on the display to focus it. Click on the rest of the page to cancel that focus.

While focused, you can interact with the simulation. All the examples have a few common controls. For example, you can press <space> to pause the simulation and then the n to advance it one step at a time. Press m to show/hide the options menu. Use the GUI to control the simulation speed, camera sensitivity, and other options. The GUI also shows the keyboard shortcuts on many controls.

Some examples provide a relatively static view while others give you more opportunity to interact with the simulation. See the on-screen help messages of each for details.

Building Desktop Applications

If you want to change an example’s code and see what happens, you will need compile your changes.

First, clone the hoomd-rs repository if you have not done so already:

$ git clone git@github.com:glotzerlab/hoomd-rs

Then change to the repository directory:

$ cd hoomd-rs

You will find the example code in examples.

To compile and run an example in batch mode (no display, GSD file output), execute:

$ cargo run --release --example {example}

where {example} is the name of the example without the path or extension (e.g. --example applying-interactions, will build and run examples/mc-tutorial/applying-interactions.rs).

To compile and run an example as an interactive desktop application, execute:

$ cargo run --release --features=bevy --example {example}

The interactive examples use the Bevy engine. If you get compile errors when building bevy crates, you may need to install additional software. MacOS is the simplest platform to configure as you only need XCode. On Windows, you will need the Microsoft C++ build tools. On Linux, you will need to install a number of system packages depending on your distribution.

Note

These additional software dependencies are only needed to build examples with --features=bevy and your own code that uses the hoomd-bevy crate. You do not need to install additional software to build and run command line applications with hoomd-rs.

Building Browser Applications

If you are having problems building or running the examples on the desktop, try building the examples for the web.

First install a few additional tools to add web support to your Rust installation:

$ rustup target install wasm32-unknown-unknown
$ cargo install wasm-bindgen-cli wasm-server-runner

To build and run an example, execute:

$ cargo run --features bevy --target wasm32-unknown-unknown --example {example}

then open the printed http://localhost URL in your browser. This works in Linux, MacOS, Windows, and WSL.

Monte Carlo Tutorial

The hoomd-mc crate implements Monte Carlo simulations.

In a Monte Carlo (MC) simulation, you chose the Hamiltonian ($H$) and what types of trial moves to propose. hoomd-rs will:

1. Propose random trial moves that change the microstate: $S\to S^{\prime }$.
2. Compute the change in energy made by each trial move: $\Delta E=H(S^{\prime })-H(S)$.
3. Then accept or reject the trial move based on criteria specific to that trial move.

For each of these steps, you can use built-in types provided by hoomd-mc or write your own. Read the following tutorials to learn how to run and customize MC simulations using hoomd-rs.

Random Walk

Overview

A random walk describes the motion of a point over a series of steps. At each step, the point translates by a random vector. This tutorial shows you how to implement a random walk using hoomd-rs. There are certainly easier ways to write a random walk code, but the purpose of this tutorial is to explain how you can express the components of a MC simulation using hoomd-rs.

• Objective: Demonstrate a minimal MC simulation.
• File: hoomd-rs/examples/mc-tutorial/random-walk.rs
• To build and run:

  cargo run --release --example random-walk
  

Bring Used Names Into Scope

The first lines of any Rust code typically bring all the used names into scope:

use hoomd_interaction::Zero;
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{Body, Microstate};
use hoomd_simulation::macrostate::Isothermal;
use hoomd_vector::Cartesian;

fn main() -> anyhow::Result<()> {
    let mut microstate = Microstate::builder()
        .bodies([Body::point(Cartesian::from([0.0, 0.0]))])
        .try_build()?;

    let translate = Translate::with_maximum_distance(0.15.try_into()?);

    let mut translate_sweep = Sweep(translate);

    let hamiltonian = Zero;

    let macrostate = Isothermal { temperature: 1.0 };

    for _ in 0..100 {
        translate_sweep.apply(&mut microstate, &hamiltonian, &macrostate);
        microstate.increment_step();

        println!("{}", microstate.sites()[0].properties.position);
    }

    Ok(())
}

Rust’s use is similar to Python’s import. See The Rust Programming Language for more information.

hoomd-interaction, hoomd-mc, hoomd-microstate, and hoomd-vector are crates that each implement a part of the simulation. See the User Guide for an overview of the hoomd-rs crates with links to their documentation.

The main() Function

When executed, every Rust application starts in the main() function:

use hoomd_interaction::Zero;
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{Body, Microstate};
use hoomd_simulation::macrostate::Isothermal;
use hoomd_vector::Cartesian;

fn main() -> anyhow::Result<()> {
    let mut microstate = Microstate::builder()
        .bodies([Body::point(Cartesian::from([0.0, 0.0]))])
        .try_build()?;

    let translate = Translate::with_maximum_distance(0.15.try_into()?);

    let mut translate_sweep = Sweep(translate);

    let hamiltonian = Zero;

    let macrostate = Isothermal { temperature: 1.0 };

    for _ in 0..100 {
        translate_sweep.apply(&mut microstate, &hamiltonian, &macrostate);
        microstate.increment_step();

        println!("{}", microstate.sites()[0].properties.position);
    }

    Ok(())
}

An error might occur while running the simulation, so main() returns a Result type. See The Rust Programming Language for more information. When you use anyhow::Result<()> as the return value for main(), the anyhow crate will nicely print any errors that occur in a human-readable format.

The Simulation Model

In the random walk simulation, the microstate contains the position of the point, the sweep applies a trial move to each point, the Hamiltonian is always 0 and the temperature is not relevant.

Microstate

The microstate describes all of the degrees of freedom in the simulation. In this example, it consists of one body represented by a single point. This code builds a microstate and adds one point at the origin:

use hoomd_interaction::Zero;
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{Body, Microstate};
use hoomd_simulation::macrostate::Isothermal;
use hoomd_vector::Cartesian;

fn main() -> anyhow::Result<()> {
    let mut microstate = Microstate::builder()
        .bodies([Body::point(Cartesian::from([0.0, 0.0]))])
        .try_build()?;

    let translate = Translate::with_maximum_distance(0.15.try_into()?);

    let mut translate_sweep = Sweep(translate);

    let hamiltonian = Zero;

    let macrostate = Isothermal { temperature: 1.0 };

    for _ in 0..100 {
        translate_sweep.apply(&mut microstate, &hamiltonian, &macrostate);
        microstate.increment_step();

        println!("{}", microstate.sites()[0].properties.position);
    }

    Ok(())
}

try_build() may fail with an error. The ? will return an error or unpack a valid Result. Read more about ? in The Rust Programming Language.

Trial Moves

A random walk is defined entirely by the trial moves applied to each body. The Translate trial move describes a displacement by a random vector drawn uniformly inside the sphere with radius maximum_distance:

use hoomd_interaction::Zero;
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{Body, Microstate};
use hoomd_simulation::macrostate::Isothermal;
use hoomd_vector::Cartesian;

fn main() -> anyhow::Result<()> {
    let mut microstate = Microstate::builder()
        .bodies([Body::point(Cartesian::from([0.0, 0.0]))])
        .try_build()?;

    let translate = Translate::with_maximum_distance(0.15.try_into()?);

    let mut translate_sweep = Sweep(translate);

    let hamiltonian = Zero;

    let macrostate = Isothermal { temperature: 1.0 };

    for _ in 0..100 {
        translate_sweep.apply(&mut microstate, &hamiltonian, &macrostate);
        microstate.increment_step();

        println!("{}", microstate.sites()[0].properties.position);
    }

    Ok(())
}

The try_into()? ensures that the given f64 value is a positive real value.

translate describes how trial moves should be applied to individual bodies. Now you need to describe how to apply these trial moves to the whole microstate using Sweep:

use hoomd_interaction::Zero;
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{Body, Microstate};
use hoomd_simulation::macrostate::Isothermal;
use hoomd_vector::Cartesian;

fn main() -> anyhow::Result<()> {
    let mut microstate = Microstate::builder()
        .bodies([Body::point(Cartesian::from([0.0, 0.0]))])
        .try_build()?;

    let translate = Translate::with_maximum_distance(0.15.try_into()?);

    let mut translate_sweep = Sweep(translate);

    let hamiltonian = Zero;

    let macrostate = Isothermal { temperature: 1.0 };

    for _ in 0..100 {
        translate_sweep.apply(&mut microstate, &hamiltonian, &macrostate);
        microstate.increment_step();

        println!("{}", microstate.sites()[0].properties.position);
    }

    Ok(())
}

Sweep applies the given trial move to each of the microstate’s bodies in sequence and accepts or rejects each move based on the Metropolis criterion.

Hamiltonian

Set the Hamiltonian $H=0$, so that Sweep will accept every random walk move trial move. The hoomd-interaction crate provides the Zero type that expresses $H=0$:

use hoomd_interaction::Zero;
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{Body, Microstate};
use hoomd_simulation::macrostate::Isothermal;
use hoomd_vector::Cartesian;

fn main() -> anyhow::Result<()> {
    let mut microstate = Microstate::builder()
        .bodies([Body::point(Cartesian::from([0.0, 0.0]))])
        .try_build()?;

    let translate = Translate::with_maximum_distance(0.15.try_into()?);

    let mut translate_sweep = Sweep(translate);

    let hamiltonian = Zero;

    let macrostate = Isothermal { temperature: 1.0 };

    for _ in 0..100 {
        translate_sweep.apply(&mut microstate, &hamiltonian, &macrostate);
        microstate.increment_step();

        println!("{}", microstate.sites()[0].properties.position);
    }

    Ok(())
}

Macrostate

Local trial moves are accepted or rejected based on the temperature of the simulation. You need to supply a macrostate even though this simulation has no interactions. Set the temperature to 1.0 as a placeholder:

use hoomd_interaction::Zero;
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{Body, Microstate};
use hoomd_simulation::macrostate::Isothermal;
use hoomd_vector::Cartesian;

fn main() -> anyhow::Result<()> {
    let mut microstate = Microstate::builder()
        .bodies([Body::point(Cartesian::from([0.0, 0.0]))])
        .try_build()?;

    let translate = Translate::with_maximum_distance(0.15.try_into()?);

    let mut translate_sweep = Sweep(translate);

    let hamiltonian = Zero;

    let macrostate = Isothermal { temperature: 1.0 };

    for _ in 0..100 {
        translate_sweep.apply(&mut microstate, &hamiltonian, &macrostate);
        microstate.increment_step();

        println!("{}", microstate.sites()[0].properties.position);
    }

    Ok(())
}

Advancing the Simulation

Use a for loop to execute the simulation a given number of steps:

use hoomd_interaction::Zero;
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{Body, Microstate};
use hoomd_simulation::macrostate::Isothermal;
use hoomd_vector::Cartesian;

fn main() -> anyhow::Result<()> {
    let mut microstate = Microstate::builder()
        .bodies([Body::point(Cartesian::from([0.0, 0.0]))])
        .try_build()?;

    let translate = Translate::with_maximum_distance(0.15.try_into()?);

    let mut translate_sweep = Sweep(translate);

    let hamiltonian = Zero;

    let macrostate = Isothermal { temperature: 1.0 };

    for _ in 0..100 {
        translate_sweep.apply(&mut microstate, &hamiltonian, &macrostate);
        microstate.increment_step();

        println!("{}", microstate.sites()[0].properties.position);
    }

    Ok(())
}

Apply Translation Moves

The translate_sweep.apply() method step applies a translate move to each body in the microstate:

use hoomd_interaction::Zero;
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{Body, Microstate};
use hoomd_simulation::macrostate::Isothermal;
use hoomd_vector::Cartesian;

fn main() -> anyhow::Result<()> {
    let mut microstate = Microstate::builder()
        .bodies([Body::point(Cartesian::from([0.0, 0.0]))])
        .try_build()?;

    let translate = Translate::with_maximum_distance(0.15.try_into()?);

    let mut translate_sweep = Sweep(translate);

    let hamiltonian = Zero;

    let macrostate = Isothermal { temperature: 1.0 };

    for _ in 0..100 {
        translate_sweep.apply(&mut microstate, &hamiltonian, &macrostate);
        microstate.increment_step();

        println!("{}", microstate.sites()[0].properties.position);
    }

    Ok(())
}

Increment the Step

hoomd-rs makes no assumptions about your simulation model. One step in your model may involve many types of MC trial moves, or a mixture of MD and MC calculations. Therefore, you must explicitly call microstate.increment_step() to indicate that this step is complete:

use hoomd_interaction::Zero;
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{Body, Microstate};
use hoomd_simulation::macrostate::Isothermal;
use hoomd_vector::Cartesian;

fn main() -> anyhow::Result<()> {
    let mut microstate = Microstate::builder()
        .bodies([Body::point(Cartesian::from([0.0, 0.0]))])
        .try_build()?;

    let translate = Translate::with_maximum_distance(0.15.try_into()?);

    let mut translate_sweep = Sweep(translate);

    let hamiltonian = Zero;

    let macrostate = Isothermal { temperature: 1.0 };

    for _ in 0..100 {
        translate_sweep.apply(&mut microstate, &hamiltonian, &macrostate);
        microstate.increment_step();

        println!("{}", microstate.sites()[0].properties.position);
    }

    Ok(())
}

Produce Output

Print the coordinates of the point at each step:

use hoomd_interaction::Zero;
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{Body, Microstate};
use hoomd_simulation::macrostate::Isothermal;
use hoomd_vector::Cartesian;

fn main() -> anyhow::Result<()> {
    let mut microstate = Microstate::builder()
        .bodies([Body::point(Cartesian::from([0.0, 0.0]))])
        .try_build()?;

    let translate = Translate::with_maximum_distance(0.15.try_into()?);

    let mut translate_sweep = Sweep(translate);

    let hamiltonian = Zero;

    let macrostate = Isothermal { temperature: 1.0 };

    for _ in 0..100 {
        translate_sweep.apply(&mut microstate, &hamiltonian, &macrostate);
        microstate.increment_step();

        println!("{}", microstate.sites()[0].properties.position);
    }

    Ok(())
}

Future tutorials will show interactive graphical representations of the simulations and/or write simulation trajectory files.

Return from the main() Function

Return Ok(()) to exit the application with no errors:

use hoomd_interaction::Zero;
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{Body, Microstate};
use hoomd_simulation::macrostate::Isothermal;
use hoomd_vector::Cartesian;

fn main() -> anyhow::Result<()> {
    let mut microstate = Microstate::builder()
        .bodies([Body::point(Cartesian::from([0.0, 0.0]))])
        .try_build()?;

    let translate = Translate::with_maximum_distance(0.15.try_into()?);

    let mut translate_sweep = Sweep(translate);

    let hamiltonian = Zero;

    let macrostate = Isothermal { temperature: 1.0 };

    for _ in 0..100 {
        translate_sweep.apply(&mut microstate, &hamiltonian, &macrostate);
        microstate.increment_step();

        println!("{}", microstate.sites()[0].properties.position);
    }

    Ok(())
}

Conclusion

Now you have learned how to create a microstate and apply random translation trial moves to the points in it. Change to the hoomd-rs directory and execute cargo run --release --example random-walk to run this tutorial. You should see output similar to:

[0.10107592706479421, 0.03942734821995694]
[0.0655227510362725, -0.02209445404554057]
[0.04227546664585716, 0.10644967595431494]
[0.1150968339945419, 0.07985960466968363]
[0.09117191688791704, 0.20862494597769088]
[0.19179261557730518, 0.15143176807959824]
[0.08739701966674704, 0.17504549059226845]
[0.16662973228106845, 0.16482677512825672]
...

Notice the particle moves a little bit on every step. Run the simulation for many steps and notice that the particles can move without bounds. By default, a Microstate has open boundary conditions.

Complete Code

use hoomd_interaction::Zero;
use hoomd_mc::{Sweep, Translate, Trial};
use hoomd_microstate::{Body, Microstate};
use hoomd_simulation::macrostate::Isothermal;
use hoomd_vector::Cartesian;

fn main() -> anyhow::Result<()> {
    let mut microstate = Microstate::builder()
        .bodies([Body::point(Cartesian::from([0.0, 0.0]))])
        .try_build()?;

    let translate = Translate::with_maximum_distance(0.15.try_into()?);

    let mut translate_sweep = Sweep(translate);

    let hamiltonian = Zero;

    let macrostate = Isothermal { temperature: 1.0 };

    for _ in 0..100 {
        translate_sweep.apply(&mut microstate, &hamiltonian, &macrostate);
        microstate.increment_step();

        println!("{}", microstate.sites()[0].properties.position);
    }

    Ok(())
}

Custom Random Walk

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

Overview

hoomd-rs allows you to customize your simulation model at all levels. This tutorial you how to implement a custom boundary condition and a custom trial move. It sets a closed boundary condition that confines bodies inside a circle, and applies trial moves that take discrete steps left, right, down, or up.

• Objective: Demonstrate the customization of a MC simulation.
• File: hoomd-rs/examples/mc-tutorial/custom-random-walk.rs
• Run (interactively):

  cargo run --release --features "bevy" --example custom-random-walk
  

• Run (in batch mode):

  cargo run --release --example custom-random-walk
  

Use Declarations

use rand::{Rng, seq::IndexedRandom};
use std::iter;

use hoomd_geometry::IsPointInside;
use hoomd_interaction::Zero;
use hoomd_mc::{LocalTrial, Sweep, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::AllPairs;
use hoomd_vector::{Cartesian, Metric};

/// Closed circular boundary condition.
struct Circle {
    radius: f64,
}

impl IsPointInside<Cartesian<2>> for Circle {
    fn is_point_inside(&self, point: &Cartesian<2>) -> bool {
        point.distance(&[0.0, 0.0].into()) < self.radius
    }
}

/// Take fixed steps left, right, down, or up.
struct Discrete;

impl LocalTrial<Point<Cartesian<2>>> for Discrete {
    fn propose<R: Rng>(
        &self,
        rng: &mut R,
        body_properties: Point<Cartesian<2>>,
    ) -> Point<Cartesian<2>> {
        let steps = [
            [0.0, -1.0].into(),
            [0.0, 1.0].into(),
            [-1.0, 0.0].into(),
            [1.0, 0.0].into(),
        ];

        let mut trial = body_properties;
        trial.position += *steps
            .choose(rng)
            .expect("steps should have at least 1 element");
        trial
    }
}

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

impl CustomRandomWalk {
    /// Construct a new random walk simulation.
    fn new() -> anyhow::Result<CustomRandomWalk> {
        let n = 1000;
        let radius = 50.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let circle = Circle { radius };

        let microstate = Microstate::builder()
            .boundary(Closed(circle))
            .bodies(iter::repeat_n(Body::point(Cartesian::default()), n))
            .try_build()?;

        let translate_sweep = Sweep(Discrete);

        let hamiltonian = Zero;

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

impl Simulation for CustomRandomWalk {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();
        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<()> {
    let mut simulation = CustomRandomWalk::new()?;
    for _ in 0..1_000_000 {
        simulation.advance()?;
    }

    Ok(())
}

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

hoomd-rs uses rand crate to generate pseudorandom numbers.

Custom Boundary Condition

The Random Walk tutorial used the default open boundary condition. This tutorial shows how you can create a custom closed boundary condition.

Define the Circle Type

First, define a type that describes the boundary. In this case, the boundary is a Circle that has a radius:

use rand::{Rng, seq::IndexedRandom};
use std::iter;

use hoomd_geometry::IsPointInside;
use hoomd_interaction::Zero;
use hoomd_mc::{LocalTrial, Sweep, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::AllPairs;
use hoomd_vector::{Cartesian, Metric};

/// Closed circular boundary condition.
struct Circle {
    radius: f64,
}

impl IsPointInside<Cartesian<2>> for Circle {
    fn is_point_inside(&self, point: &Cartesian<2>) -> bool {
        point.distance(&[0.0, 0.0].into()) < self.radius
    }
}

/// Take fixed steps left, right, down, or up.
struct Discrete;

impl LocalTrial<Point<Cartesian<2>>> for Discrete {
    fn propose<R: Rng>(
        &self,
        rng: &mut R,
        body_properties: Point<Cartesian<2>>,
    ) -> Point<Cartesian<2>> {
        let steps = [
            [0.0, -1.0].into(),
            [0.0, 1.0].into(),
            [-1.0, 0.0].into(),
            [1.0, 0.0].into(),
        ];

        let mut trial = body_properties;
        trial.position += *steps
            .choose(rng)
            .expect("steps should have at least 1 element");
        trial
    }
}

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

impl CustomRandomWalk {
    /// Construct a new random walk simulation.
    fn new() -> anyhow::Result<CustomRandomWalk> {
        let n = 1000;
        let radius = 50.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let circle = Circle { radius };

        let microstate = Microstate::builder()
            .boundary(Closed(circle))
            .bodies(iter::repeat_n(Body::point(Cartesian::default()), n))
            .try_build()?;

        let translate_sweep = Sweep(Discrete);

        let hamiltonian = Zero;

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

impl Simulation for CustomRandomWalk {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();
        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<()> {
    let mut simulation = CustomRandomWalk::new()?;
    for _ in 0..1_000_000 {
        simulation.advance()?;
    }

    Ok(())
}

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

Implement IsPointInside

Then, implement the IsPointInside trait for Circle. The only required method is is_point_inside() which should return true for points inside the boundary and false for points outside the boundary:

use rand::{Rng, seq::IndexedRandom};
use std::iter;

use hoomd_geometry::IsPointInside;
use hoomd_interaction::Zero;
use hoomd_mc::{LocalTrial, Sweep, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::AllPairs;
use hoomd_vector::{Cartesian, Metric};

/// Closed circular boundary condition.
struct Circle {
    radius: f64,
}

impl IsPointInside<Cartesian<2>> for Circle {
    fn is_point_inside(&self, point: &Cartesian<2>) -> bool {
        point.distance(&[0.0, 0.0].into()) < self.radius
    }
}

/// Take fixed steps left, right, down, or up.
struct Discrete;

impl LocalTrial<Point<Cartesian<2>>> for Discrete {
    fn propose<R: Rng>(
        &self,
        rng: &mut R,
        body_properties: Point<Cartesian<2>>,
    ) -> Point<Cartesian<2>> {
        let steps = [
            [0.0, -1.0].into(),
            [0.0, 1.0].into(),
            [-1.0, 0.0].into(),
            [1.0, 0.0].into(),
        ];

        let mut trial = body_properties;
        trial.position += *steps
            .choose(rng)
            .expect("steps should have at least 1 element");
        trial
    }
}

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

impl CustomRandomWalk {
    /// Construct a new random walk simulation.
    fn new() -> anyhow::Result<CustomRandomWalk> {
        let n = 1000;
        let radius = 50.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let circle = Circle { radius };

        let microstate = Microstate::builder()
            .boundary(Closed(circle))
            .bodies(iter::repeat_n(Body::point(Cartesian::default()), n))
            .try_build()?;

        let translate_sweep = Sweep(Discrete);

        let hamiltonian = Zero;

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

impl Simulation for CustomRandomWalk {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();
        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<()> {
    let mut simulation = CustomRandomWalk::new()?;
    for _ in 0..1_000_000 {
        simulation.advance()?;
    }

    Ok(())
}

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

This implementation of IsPointInside is only for the Cartesian<2> vector type (V in the IsPointInside definition). You can read The Rust Programming Language to learn more about generic types and traits.

Tip

Many shapes in hoomd_geometry implement IsPointInside and can be used for closed boundary conditions. This tutorial creates a new Circle type to teach you about customization, but hoomd_geometry::shape::Circle would work just as well.

Custom Trial Move

The Random Walk tutorial moved points by a random distance (up to a maximum) and in a random direction. Look up “random walk” in a text book and you find will a model that makes random hops on a square lattice. You can implement that in hoomd-rs with a custom trial move.

Define the Discrete Type

Similar to the custom boundary, you need to implement the trait LocalTrial for a type. Here is the type:

use rand::{Rng, seq::IndexedRandom};
use std::iter;

use hoomd_geometry::IsPointInside;
use hoomd_interaction::Zero;
use hoomd_mc::{LocalTrial, Sweep, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::AllPairs;
use hoomd_vector::{Cartesian, Metric};

/// Closed circular boundary condition.
struct Circle {
    radius: f64,
}

impl IsPointInside<Cartesian<2>> for Circle {
    fn is_point_inside(&self, point: &Cartesian<2>) -> bool {
        point.distance(&[0.0, 0.0].into()) < self.radius
    }
}

/// Take fixed steps left, right, down, or up.
struct Discrete;

impl LocalTrial<Point<Cartesian<2>>> for Discrete {
    fn propose<R: Rng>(
        &self,
        rng: &mut R,
        body_properties: Point<Cartesian<2>>,
    ) -> Point<Cartesian<2>> {
        let steps = [
            [0.0, -1.0].into(),
            [0.0, 1.0].into(),
            [-1.0, 0.0].into(),
            [1.0, 0.0].into(),
        ];

        let mut trial = body_properties;
        trial.position += *steps
            .choose(rng)
            .expect("steps should have at least 1 element");
        trial
    }
}

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

impl CustomRandomWalk {
    /// Construct a new random walk simulation.
    fn new() -> anyhow::Result<CustomRandomWalk> {
        let n = 1000;
        let radius = 50.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let circle = Circle { radius };

        let microstate = Microstate::builder()
            .boundary(Closed(circle))
            .bodies(iter::repeat_n(Body::point(Cartesian::default()), n))
            .try_build()?;

        let translate_sweep = Sweep(Discrete);

        let hamiltonian = Zero;

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

impl Simulation for CustomRandomWalk {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();
        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<()> {
    let mut simulation = CustomRandomWalk::new()?;
    for _ in 0..1_000_000 {
        simulation.advance()?;
    }

    Ok(())
}

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

Discrete trials always take one step in one direction. Therefore, the Discrete struct needs no fields.

Implement LocalTrial

LocalTrial has one method, propose():

use rand::{Rng, seq::IndexedRandom};
use std::iter;

use hoomd_geometry::IsPointInside;
use hoomd_interaction::Zero;
use hoomd_mc::{LocalTrial, Sweep, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::AllPairs;
use hoomd_vector::{Cartesian, Metric};

/// Closed circular boundary condition.
struct Circle {
    radius: f64,
}

impl IsPointInside<Cartesian<2>> for Circle {
    fn is_point_inside(&self, point: &Cartesian<2>) -> bool {
        point.distance(&[0.0, 0.0].into()) < self.radius
    }
}

/// Take fixed steps left, right, down, or up.
struct Discrete;

impl LocalTrial<Point<Cartesian<2>>> for Discrete {
    fn propose<R: Rng>(
        &self,
        rng: &mut R,
        body_properties: Point<Cartesian<2>>,
    ) -> Point<Cartesian<2>> {
        let steps = [
            [0.0, -1.0].into(),
            [0.0, 1.0].into(),
            [-1.0, 0.0].into(),
            [1.0, 0.0].into(),
        ];

        let mut trial = body_properties;
        trial.position += *steps
            .choose(rng)
            .expect("steps should have at least 1 element");
        trial
    }
}

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

impl CustomRandomWalk {
    /// Construct a new random walk simulation.
    fn new() -> anyhow::Result<CustomRandomWalk> {
        let n = 1000;
        let radius = 50.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let circle = Circle { radius };

        let microstate = Microstate::builder()
            .boundary(Closed(circle))
            .bodies(iter::repeat_n(Body::point(Cartesian::default()), n))
            .try_build()?;

        let translate_sweep = Sweep(Discrete);

        let hamiltonian = Zero;

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

impl Simulation for CustomRandomWalk {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();
        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<()> {
    let mut simulation = CustomRandomWalk::new()?;
    for _ in 0..1_000_000 {
        simulation.advance()?;
    }

    Ok(())
}

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

propose() takes in the properties (body_properties) of a body in the current microstate of the system. It returns the trial body properties randomly generated using the given rng. In this tutorial, the only body property is position.

Enumerate Possible Steps

First, place the possible moves (down, up, left, right) in an array:

use rand::{Rng, seq::IndexedRandom};
use std::iter;

use hoomd_geometry::IsPointInside;
use hoomd_interaction::Zero;
use hoomd_mc::{LocalTrial, Sweep, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::AllPairs;
use hoomd_vector::{Cartesian, Metric};

/// Closed circular boundary condition.
struct Circle {
    radius: f64,
}

impl IsPointInside<Cartesian<2>> for Circle {
    fn is_point_inside(&self, point: &Cartesian<2>) -> bool {
        point.distance(&[0.0, 0.0].into()) < self.radius
    }
}

/// Take fixed steps left, right, down, or up.
struct Discrete;

impl LocalTrial<Point<Cartesian<2>>> for Discrete {
    fn propose<R: Rng>(
        &self,
        rng: &mut R,
        body_properties: Point<Cartesian<2>>,
    ) -> Point<Cartesian<2>> {
        let steps = [
            [0.0, -1.0].into(),
            [0.0, 1.0].into(),
            [-1.0, 0.0].into(),
            [1.0, 0.0].into(),
        ];

        let mut trial = body_properties;
        trial.position += *steps
            .choose(rng)
            .expect("steps should have at least 1 element");
        trial
    }
}

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

impl CustomRandomWalk {
    /// Construct a new random walk simulation.
    fn new() -> anyhow::Result<CustomRandomWalk> {
        let n = 1000;
        let radius = 50.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let circle = Circle { radius };

        let microstate = Microstate::builder()
            .boundary(Closed(circle))
            .bodies(iter::repeat_n(Body::point(Cartesian::default()), n))
            .try_build()?;

        let translate_sweep = Sweep(Discrete);

        let hamiltonian = Zero;

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

impl Simulation for CustomRandomWalk {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();
        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<()> {
    let mut simulation = CustomRandomWalk::new()?;
    for _ in 0..1_000_000 {
        simulation.advance()?;
    }

    Ok(())
}

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

Randomly Select a Step

Then choose one of the steps randomly and add it to the body’s position:

use rand::{Rng, seq::IndexedRandom};
use std::iter;

use hoomd_geometry::IsPointInside;
use hoomd_interaction::Zero;
use hoomd_mc::{LocalTrial, Sweep, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::AllPairs;
use hoomd_vector::{Cartesian, Metric};

/// Closed circular boundary condition.
struct Circle {
    radius: f64,
}

impl IsPointInside<Cartesian<2>> for Circle {
    fn is_point_inside(&self, point: &Cartesian<2>) -> bool {
        point.distance(&[0.0, 0.0].into()) < self.radius
    }
}

/// Take fixed steps left, right, down, or up.
struct Discrete;

impl LocalTrial<Point<Cartesian<2>>> for Discrete {
    fn propose<R: Rng>(
        &self,
        rng: &mut R,
        body_properties: Point<Cartesian<2>>,
    ) -> Point<Cartesian<2>> {
        let steps = [
            [0.0, -1.0].into(),
            [0.0, 1.0].into(),
            [-1.0, 0.0].into(),
            [1.0, 0.0].into(),
        ];

        let mut trial = body_properties;
        trial.position += *steps
            .choose(rng)
            .expect("steps should have at least 1 element");
        trial
    }
}

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

impl CustomRandomWalk {
    /// Construct a new random walk simulation.
    fn new() -> anyhow::Result<CustomRandomWalk> {
        let n = 1000;
        let radius = 50.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let circle = Circle { radius };

        let microstate = Microstate::builder()
            .boundary(Closed(circle))
            .bodies(iter::repeat_n(Body::point(Cartesian::default()), n))
            .try_build()?;

        let translate_sweep = Sweep(Discrete);

        let hamiltonian = Zero;

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

impl Simulation for CustomRandomWalk {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();
        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<()> {
    let mut simulation = CustomRandomWalk::new()?;
    for _ in 0..1_000_000 {
        simulation.advance()?;
    }

    Ok(())
}

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

steps.choose(rng) chooses one of the elements of the array at random (with uniform probability). Arrays in Rust do not have the choose() method by default, it is provided by the IndexedRandom trait in the rand crate. choose() will return None when the input array is empty. That can never happen in this code, so you can safely unpack the option’s value with expect().

Important

Local trial moves MUST satisfy local detailed balance, as defined in Manousiouthakis & Deem.

The Simulation Model

The Random Walk tutorial used local variables in the main() function to store all the elements of the simulation model. You can certainly continue that practice in your applications. As you build more complex codes, however, you will need to move those elements to a struct so that the whole simulation model can be accessed in different modules.

The custom random walk model consists of the microstate, the Hamiltonian, the translation trial moves, and the temperature set point (in units of energy $kT$):

use rand::{Rng, seq::IndexedRandom};
use std::iter;

use hoomd_geometry::IsPointInside;
use hoomd_interaction::Zero;
use hoomd_mc::{LocalTrial, Sweep, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::AllPairs;
use hoomd_vector::{Cartesian, Metric};

/// Closed circular boundary condition.
struct Circle {
    radius: f64,
}

impl IsPointInside<Cartesian<2>> for Circle {
    fn is_point_inside(&self, point: &Cartesian<2>) -> bool {
        point.distance(&[0.0, 0.0].into()) < self.radius
    }
}

/// Take fixed steps left, right, down, or up.
struct Discrete;

impl LocalTrial<Point<Cartesian<2>>> for Discrete {
    fn propose<R: Rng>(
        &self,
        rng: &mut R,
        body_properties: Point<Cartesian<2>>,
    ) -> Point<Cartesian<2>> {
        let steps = [
            [0.0, -1.0].into(),
            [0.0, 1.0].into(),
            [-1.0, 0.0].into(),
            [1.0, 0.0].into(),
        ];

        let mut trial = body_properties;
        trial.position += *steps
            .choose(rng)
            .expect("steps should have at least 1 element");
        trial
    }
}

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

impl CustomRandomWalk {
    /// Construct a new random walk simulation.
    fn new() -> anyhow::Result<CustomRandomWalk> {
        let n = 1000;
        let radius = 50.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let circle = Circle { radius };

        let microstate = Microstate::builder()
            .boundary(Closed(circle))
            .bodies(iter::repeat_n(Body::point(Cartesian::default()), n))
            .try_build()?;

        let translate_sweep = Sweep(Discrete);

        let hamiltonian = Zero;

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

impl Simulation for CustomRandomWalk {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();
        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<()> {
    let mut simulation = CustomRandomWalk::new()?;
    for _ in 0..1_000_000 {
        simulation.advance()?;
    }

    Ok(())
}

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

Notice that this definition explicitly names the generic types of each of these fields. The Random Walk tutorial did not do this because the Rust compiler automatically determined which generic types were used. Structs can be made generic as well, but the details are beyond the scope of this tutorial. Consult the The Rust Programming Language to learn more.

Construct the Simulation Model

The new() method that constructs the CustomRandomWalk simulation model:

use rand::{Rng, seq::IndexedRandom};
use std::iter;

use hoomd_geometry::IsPointInside;
use hoomd_interaction::Zero;
use hoomd_mc::{LocalTrial, Sweep, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::AllPairs;
use hoomd_vector::{Cartesian, Metric};

/// Closed circular boundary condition.
struct Circle {
    radius: f64,
}

impl IsPointInside<Cartesian<2>> for Circle {
    fn is_point_inside(&self, point: &Cartesian<2>) -> bool {
        point.distance(&[0.0, 0.0].into()) < self.radius
    }
}

/// Take fixed steps left, right, down, or up.
struct Discrete;

impl LocalTrial<Point<Cartesian<2>>> for Discrete {
    fn propose<R: Rng>(
        &self,
        rng: &mut R,
        body_properties: Point<Cartesian<2>>,
    ) -> Point<Cartesian<2>> {
        let steps = [
            [0.0, -1.0].into(),
            [0.0, 1.0].into(),
            [-1.0, 0.0].into(),
            [1.0, 0.0].into(),
        ];

        let mut trial = body_properties;
        trial.position += *steps
            .choose(rng)
            .expect("steps should have at least 1 element");
        trial
    }
}

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

impl CustomRandomWalk {
    /// Construct a new random walk simulation.
    fn new() -> anyhow::Result<CustomRandomWalk> {
        let n = 1000;
        let radius = 50.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let circle = Circle { radius };

        let microstate = Microstate::builder()
            .boundary(Closed(circle))
            .bodies(iter::repeat_n(Body::point(Cartesian::default()), n))
            .try_build()?;

        let translate_sweep = Sweep(Discrete);

        let hamiltonian = Zero;

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

impl Simulation for CustomRandomWalk {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();
        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<()> {
    let mut simulation = CustomRandomWalk::new()?;
    for _ in 0..1_000_000 {
        simulation.advance()?;
    }

    Ok(())
}

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

One or more steps might fail, so return a Result<CustomRandomWalk>.

Parameters

Assign all the model parameters in one code block so that they are easy to modify:

use rand::{Rng, seq::IndexedRandom};
use std::iter;

use hoomd_geometry::IsPointInside;
use hoomd_interaction::Zero;
use hoomd_mc::{LocalTrial, Sweep, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::AllPairs;
use hoomd_vector::{Cartesian, Metric};

/// Closed circular boundary condition.
struct Circle {
    radius: f64,
}

impl IsPointInside<Cartesian<2>> for Circle {
    fn is_point_inside(&self, point: &Cartesian<2>) -> bool {
        point.distance(&[0.0, 0.0].into()) < self.radius
    }
}

/// Take fixed steps left, right, down, or up.
struct Discrete;

impl LocalTrial<Point<Cartesian<2>>> for Discrete {
    fn propose<R: Rng>(
        &self,
        rng: &mut R,
        body_properties: Point<Cartesian<2>>,
    ) -> Point<Cartesian<2>> {
        let steps = [
            [0.0, -1.0].into(),
            [0.0, 1.0].into(),
            [-1.0, 0.0].into(),
            [1.0, 0.0].into(),
        ];

        let mut trial = body_properties;
        trial.position += *steps
            .choose(rng)
            .expect("steps should have at least 1 element");
        trial
    }
}

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

impl CustomRandomWalk {
    /// Construct a new random walk simulation.
    fn new() -> anyhow::Result<CustomRandomWalk> {
        let n = 1000;
        let radius = 50.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let circle = Circle { radius };

        let microstate = Microstate::builder()
            .boundary(Closed(circle))
            .bodies(iter::repeat_n(Body::point(Cartesian::default()), n))
            .try_build()?;

        let translate_sweep = Sweep(Discrete);

        let hamiltonian = Zero;

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

impl Simulation for CustomRandomWalk {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();
        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<()> {
    let mut simulation = CustomRandomWalk::new()?;
    for _ in 0..1_000_000 {
        simulation.advance()?;
    }

    Ok(())
}

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

Microstate

Construct the Microstate with the circle boundary condition and place n bodies at the origin:

use rand::{Rng, seq::IndexedRandom};
use std::iter;

use hoomd_geometry::IsPointInside;
use hoomd_interaction::Zero;
use hoomd_mc::{LocalTrial, Sweep, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::AllPairs;
use hoomd_vector::{Cartesian, Metric};

/// Closed circular boundary condition.
struct Circle {
    radius: f64,
}

impl IsPointInside<Cartesian<2>> for Circle {
    fn is_point_inside(&self, point: &Cartesian<2>) -> bool {
        point.distance(&[0.0, 0.0].into()) < self.radius
    }
}

/// Take fixed steps left, right, down, or up.
struct Discrete;

impl LocalTrial<Point<Cartesian<2>>> for Discrete {
    fn propose<R: Rng>(
        &self,
        rng: &mut R,
        body_properties: Point<Cartesian<2>>,
    ) -> Point<Cartesian<2>> {
        let steps = [
            [0.0, -1.0].into(),
            [0.0, 1.0].into(),
            [-1.0, 0.0].into(),
            [1.0, 0.0].into(),
        ];

        let mut trial = body_properties;
        trial.position += *steps
            .choose(rng)
            .expect("steps should have at least 1 element");
        trial
    }
}

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

impl CustomRandomWalk {
    /// Construct a new random walk simulation.
    fn new() -> anyhow::Result<CustomRandomWalk> {
        let n = 1000;
        let radius = 50.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let circle = Circle { radius };

        let microstate = Microstate::builder()
            .boundary(Closed(circle))
            .bodies(iter::repeat_n(Body::point(Cartesian::default()), n))
            .try_build()?;

        let translate_sweep = Sweep(Discrete);

        let hamiltonian = Zero;

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

impl Simulation for CustomRandomWalk {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();
        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<()> {
    let mut simulation = CustomRandomWalk::new()?;
    for _ in 0..1_000_000 {
        simulation.advance()?;
    }

    Ok(())
}

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

The newtype Closed can wrap any type that implements IsPointInside (like the custom Circle) so that it can be used as a boundary condition.

Trial Moves

Apply the custom Discrete trial move to all bodies in the microstate:

use rand::{Rng, seq::IndexedRandom};
use std::iter;

use hoomd_geometry::IsPointInside;
use hoomd_interaction::Zero;
use hoomd_mc::{LocalTrial, Sweep, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::AllPairs;
use hoomd_vector::{Cartesian, Metric};

/// Closed circular boundary condition.
struct Circle {
    radius: f64,
}

impl IsPointInside<Cartesian<2>> for Circle {
    fn is_point_inside(&self, point: &Cartesian<2>) -> bool {
        point.distance(&[0.0, 0.0].into()) < self.radius
    }
}

/// Take fixed steps left, right, down, or up.
struct Discrete;

impl LocalTrial<Point<Cartesian<2>>> for Discrete {
    fn propose<R: Rng>(
        &self,
        rng: &mut R,
        body_properties: Point<Cartesian<2>>,
    ) -> Point<Cartesian<2>> {
        let steps = [
            [0.0, -1.0].into(),
            [0.0, 1.0].into(),
            [-1.0, 0.0].into(),
            [1.0, 0.0].into(),
        ];

        let mut trial = body_properties;
        trial.position += *steps
            .choose(rng)
            .expect("steps should have at least 1 element");
        trial
    }
}

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

impl CustomRandomWalk {
    /// Construct a new random walk simulation.
    fn new() -> anyhow::Result<CustomRandomWalk> {
        let n = 1000;
        let radius = 50.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let circle = Circle { radius };

        let microstate = Microstate::builder()
            .boundary(Closed(circle))
            .bodies(iter::repeat_n(Body::point(Cartesian::default()), n))
            .try_build()?;

        let translate_sweep = Sweep(Discrete);

        let hamiltonian = Zero;

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

impl Simulation for CustomRandomWalk {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();
        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<()> {
    let mut simulation = CustomRandomWalk::new()?;
    for _ in 0..1_000_000 {
        simulation.advance()?;
    }

    Ok(())
}

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

Sweep can wrap any type that implements LocalTrial.

Hamiltonian

As in the Random Walk tutorial, set $H=0$ so that bodies do not interact:

use rand::{Rng, seq::IndexedRandom};
use std::iter;

use hoomd_geometry::IsPointInside;
use hoomd_interaction::Zero;
use hoomd_mc::{LocalTrial, Sweep, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::AllPairs;
use hoomd_vector::{Cartesian, Metric};

/// Closed circular boundary condition.
struct Circle {
    radius: f64,
}

impl IsPointInside<Cartesian<2>> for Circle {
    fn is_point_inside(&self, point: &Cartesian<2>) -> bool {
        point.distance(&[0.0, 0.0].into()) < self.radius
    }
}

/// Take fixed steps left, right, down, or up.
struct Discrete;

impl LocalTrial<Point<Cartesian<2>>> for Discrete {
    fn propose<R: Rng>(
        &self,
        rng: &mut R,
        body_properties: Point<Cartesian<2>>,
    ) -> Point<Cartesian<2>> {
        let steps = [
            [0.0, -1.0].into(),
            [0.0, 1.0].into(),
            [-1.0, 0.0].into(),
            [1.0, 0.0].into(),
        ];

        let mut trial = body_properties;
        trial.position += *steps
            .choose(rng)
            .expect("steps should have at least 1 element");
        trial
    }
}

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

impl CustomRandomWalk {
    /// Construct a new random walk simulation.
    fn new() -> anyhow::Result<CustomRandomWalk> {
        let n = 1000;
        let radius = 50.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let circle = Circle { radius };

        let microstate = Microstate::builder()
            .boundary(Closed(circle))
            .bodies(iter::repeat_n(Body::point(Cartesian::default()), n))
            .try_build()?;

        let translate_sweep = Sweep(Discrete);

        let hamiltonian = Zero;

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

impl Simulation for CustomRandomWalk {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();
        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<()> {
    let mut simulation = CustomRandomWalk::new()?;
    for _ in 0..1_000_000 {
        simulation.advance()?;
    }

    Ok(())
}

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

Initialize the Struct

Now that all of the elements of the simulation model have been constructed, package them in a CustomRandomWalk struct:

use rand::{Rng, seq::IndexedRandom};
use std::iter;

use hoomd_geometry::IsPointInside;
use hoomd_interaction::Zero;
use hoomd_mc::{LocalTrial, Sweep, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::AllPairs;
use hoomd_vector::{Cartesian, Metric};

/// Closed circular boundary condition.
struct Circle {
    radius: f64,
}

impl IsPointInside<Cartesian<2>> for Circle {
    fn is_point_inside(&self, point: &Cartesian<2>) -> bool {
        point.distance(&[0.0, 0.0].into()) < self.radius
    }
}

/// Take fixed steps left, right, down, or up.
struct Discrete;

impl LocalTrial<Point<Cartesian<2>>> for Discrete {
    fn propose<R: Rng>(
        &self,
        rng: &mut R,
        body_properties: Point<Cartesian<2>>,
    ) -> Point<Cartesian<2>> {
        let steps = [
            [0.0, -1.0].into(),
            [0.0, 1.0].into(),
            [-1.0, 0.0].into(),
            [1.0, 0.0].into(),
        ];

        let mut trial = body_properties;
        trial.position += *steps
            .choose(rng)
            .expect("steps should have at least 1 element");
        trial
    }
}

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

impl CustomRandomWalk {
    /// Construct a new random walk simulation.
    fn new() -> anyhow::Result<CustomRandomWalk> {
        let n = 1000;
        let radius = 50.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let circle = Circle { radius };

        let microstate = Microstate::builder()
            .boundary(Closed(circle))
            .bodies(iter::repeat_n(Body::point(Cartesian::default()), n))
            .try_build()?;

        let translate_sweep = Sweep(Discrete);

        let hamiltonian = Zero;

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

impl Simulation for CustomRandomWalk {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();
        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<()> {
    let mut simulation = CustomRandomWalk::new()?;
    for _ in 0..1_000_000 {
        simulation.advance()?;
    }

    Ok(())
}

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

The struct is wrapped in Ok to indicate that the Result of this function is not an error.

Implement Simulation

The Simulation trait provides a common interface that all simulation models can follow:

use rand::{Rng, seq::IndexedRandom};
use std::iter;

use hoomd_geometry::IsPointInside;
use hoomd_interaction::Zero;
use hoomd_mc::{LocalTrial, Sweep, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::AllPairs;
use hoomd_vector::{Cartesian, Metric};

/// Closed circular boundary condition.
struct Circle {
    radius: f64,
}

impl IsPointInside<Cartesian<2>> for Circle {
    fn is_point_inside(&self, point: &Cartesian<2>) -> bool {
        point.distance(&[0.0, 0.0].into()) < self.radius
    }
}

/// Take fixed steps left, right, down, or up.
struct Discrete;

impl LocalTrial<Point<Cartesian<2>>> for Discrete {
    fn propose<R: Rng>(
        &self,
        rng: &mut R,
        body_properties: Point<Cartesian<2>>,
    ) -> Point<Cartesian<2>> {
        let steps = [
            [0.0, -1.0].into(),
            [0.0, 1.0].into(),
            [-1.0, 0.0].into(),
            [1.0, 0.0].into(),
        ];

        let mut trial = body_properties;
        trial.position += *steps
            .choose(rng)
            .expect("steps should have at least 1 element");
        trial
    }
}

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

impl CustomRandomWalk {
    /// Construct a new random walk simulation.
    fn new() -> anyhow::Result<CustomRandomWalk> {
        let n = 1000;
        let radius = 50.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let circle = Circle { radius };

        let microstate = Microstate::builder()
            .boundary(Closed(circle))
            .bodies(iter::repeat_n(Body::point(Cartesian::default()), n))
            .try_build()?;

        let translate_sweep = Sweep(Discrete);

        let hamiltonian = Zero;

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

impl Simulation for CustomRandomWalk {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();
        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<()> {
    let mut simulation = CustomRandomWalk::new()?;
    for _ in 0..1_000_000 {
        simulation.advance()?;
    }

    Ok(())
}

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

Advance the Simulation

The first method that all simulation models must have is an advance() method that moves the model forward one step. The implementation is identical to that in the Random Walk tutorial:

use rand::{Rng, seq::IndexedRandom};
use std::iter;

use hoomd_geometry::IsPointInside;
use hoomd_interaction::Zero;
use hoomd_mc::{LocalTrial, Sweep, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::AllPairs;
use hoomd_vector::{Cartesian, Metric};

/// Closed circular boundary condition.
struct Circle {
    radius: f64,
}

impl IsPointInside<Cartesian<2>> for Circle {
    fn is_point_inside(&self, point: &Cartesian<2>) -> bool {
        point.distance(&[0.0, 0.0].into()) < self.radius
    }
}

/// Take fixed steps left, right, down, or up.
struct Discrete;

impl LocalTrial<Point<Cartesian<2>>> for Discrete {
    fn propose<R: Rng>(
        &self,
        rng: &mut R,
        body_properties: Point<Cartesian<2>>,
    ) -> Point<Cartesian<2>> {
        let steps = [
            [0.0, -1.0].into(),
            [0.0, 1.0].into(),
            [-1.0, 0.0].into(),
            [1.0, 0.0].into(),
        ];

        let mut trial = body_properties;
        trial.position += *steps
            .choose(rng)
            .expect("steps should have at least 1 element");
        trial
    }
}

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

impl CustomRandomWalk {
    /// Construct a new random walk simulation.
    fn new() -> anyhow::Result<CustomRandomWalk> {
        let n = 1000;
        let radius = 50.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let circle = Circle { radius };

        let microstate = Microstate::builder()
            .boundary(Closed(circle))
            .bodies(iter::repeat_n(Body::point(Cartesian::default()), n))
            .try_build()?;

        let translate_sweep = Sweep(Discrete);

        let hamiltonian = Zero;

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

impl Simulation for CustomRandomWalk {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();
        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<()> {
    let mut simulation = CustomRandomWalk::new()?;
    for _ in 0..1_000_000 {
        simulation.advance()?;
    }

    Ok(())
}

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

All the complexity of the customizations is contained in the implementation of IsPointInside for Circle and LocalTrial for Discrete.

Get the Simulation Step

All simulation models must also implement a method that provides the current simulation step:

use rand::{Rng, seq::IndexedRandom};
use std::iter;

use hoomd_geometry::IsPointInside;
use hoomd_interaction::Zero;
use hoomd_mc::{LocalTrial, Sweep, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::AllPairs;
use hoomd_vector::{Cartesian, Metric};

/// Closed circular boundary condition.
struct Circle {
    radius: f64,
}

impl IsPointInside<Cartesian<2>> for Circle {
    fn is_point_inside(&self, point: &Cartesian<2>) -> bool {
        point.distance(&[0.0, 0.0].into()) < self.radius
    }
}

/// Take fixed steps left, right, down, or up.
struct Discrete;

impl LocalTrial<Point<Cartesian<2>>> for Discrete {
    fn propose<R: Rng>(
        &self,
        rng: &mut R,
        body_properties: Point<Cartesian<2>>,
    ) -> Point<Cartesian<2>> {
        let steps = [
            [0.0, -1.0].into(),
            [0.0, 1.0].into(),
            [-1.0, 0.0].into(),
            [1.0, 0.0].into(),
        ];

        let mut trial = body_properties;
        trial.position += *steps
            .choose(rng)
            .expect("steps should have at least 1 element");
        trial
    }
}

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

impl CustomRandomWalk {
    /// Construct a new random walk simulation.
    fn new() -> anyhow::Result<CustomRandomWalk> {
        let n = 1000;
        let radius = 50.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let circle = Circle { radius };

        let microstate = Microstate::builder()
            .boundary(Closed(circle))
            .bodies(iter::repeat_n(Body::point(Cartesian::default()), n))
            .try_build()?;

        let translate_sweep = Sweep(Discrete);

        let hamiltonian = Zero;

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

impl Simulation for CustomRandomWalk {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();
        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<()> {
    let mut simulation = CustomRandomWalk::new()?;
    for _ in 0..1_000_000 {
        simulation.advance()?;
    }

    Ok(())
}

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

Implement main()

To run the simulation, construct the CustomRandomWalk simulation model. Then call advance() many times:

use rand::{Rng, seq::IndexedRandom};
use std::iter;

use hoomd_geometry::IsPointInside;
use hoomd_interaction::Zero;
use hoomd_mc::{LocalTrial, Sweep, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::AllPairs;
use hoomd_vector::{Cartesian, Metric};

/// Closed circular boundary condition.
struct Circle {
    radius: f64,
}

impl IsPointInside<Cartesian<2>> for Circle {
    fn is_point_inside(&self, point: &Cartesian<2>) -> bool {
        point.distance(&[0.0, 0.0].into()) < self.radius
    }
}

/// Take fixed steps left, right, down, or up.
struct Discrete;

impl LocalTrial<Point<Cartesian<2>>> for Discrete {
    fn propose<R: Rng>(
        &self,
        rng: &mut R,
        body_properties: Point<Cartesian<2>>,
    ) -> Point<Cartesian<2>> {
        let steps = [
            [0.0, -1.0].into(),
            [0.0, 1.0].into(),
            [-1.0, 0.0].into(),
            [1.0, 0.0].into(),
        ];

        let mut trial = body_properties;
        trial.position += *steps
            .choose(rng)
            .expect("steps should have at least 1 element");
        trial
    }
}

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

impl CustomRandomWalk {
    /// Construct a new random walk simulation.
    fn new() -> anyhow::Result<CustomRandomWalk> {
        let n = 1000;
        let radius = 50.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let circle = Circle { radius };

        let microstate = Microstate::builder()
            .boundary(Closed(circle))
            .bodies(iter::repeat_n(Body::point(Cartesian::default()), n))
            .try_build()?;

        let translate_sweep = Sweep(Discrete);

        let hamiltonian = Zero;

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

impl Simulation for CustomRandomWalk {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();
        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<()> {
    let mut simulation = CustomRandomWalk::new()?;
    for _ in 0..1_000_000 {
        simulation.advance()?;
    }

    Ok(())
}

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

Later tutorials will show how you can write simulation trajectories and log properties of the simulation to files.

Note

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

Conclusion

This tutorial showed you how to customize the random walk simulation with new boundary conditions and apply your own trial moves to the points in it. Rust compiles your customizations into machine code and can inline them into the main simulation loop. This means that your custom simulations run just as fast as they do when using the built-in types.

Navigate to the top of the page to see the simulation in action. Notice that no points leave the boundary. Try pausing the simulation and advancing one step at a time. You should see that every body moves left, right, down, or up on every step.

You can also run the example in batch mode and then open the generated trajectory.gsd in Ovito or another visualization tool:

cargo run --release --example custom-random-walk

Complete Code

use rand::{Rng, seq::IndexedRandom};
use std::iter;

use hoomd_geometry::IsPointInside;
use hoomd_interaction::Zero;
use hoomd_mc::{LocalTrial, Sweep, Trial};
use hoomd_microstate::{
    Body, Microstate, SiteKey, boundary::Closed, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::AllPairs;
use hoomd_vector::{Cartesian, Metric};

/// Closed circular boundary condition.
struct Circle {
    radius: f64,
}

impl IsPointInside<Cartesian<2>> for Circle {
    fn is_point_inside(&self, point: &Cartesian<2>) -> bool {
        point.distance(&[0.0, 0.0].into()) < self.radius
    }
}

/// Take fixed steps left, right, down, or up.
struct Discrete;

impl LocalTrial<Point<Cartesian<2>>> for Discrete {
    fn propose<R: Rng>(
        &self,
        rng: &mut R,
        body_properties: Point<Cartesian<2>>,
    ) -> Point<Cartesian<2>> {
        let steps = [
            [0.0, -1.0].into(),
            [0.0, 1.0].into(),
            [-1.0, 0.0].into(),
            [1.0, 0.0].into(),
        ];

        let mut trial = body_properties;
        trial.position += *steps
            .choose(rng)
            .expect("steps should have at least 1 element");
        trial
    }
}

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

impl CustomRandomWalk {
    /// Construct a new random walk simulation.
    fn new() -> anyhow::Result<CustomRandomWalk> {
        let n = 1000;
        let radius = 50.0;
        let macrostate = Isothermal { temperature: 1.0 };

        let circle = Circle { radius };

        let microstate = Microstate::builder()
            .boundary(Closed(circle))
            .bodies(iter::repeat_n(Body::point(Cartesian::default()), n))
            .try_build()?;

        let translate_sweep = Sweep(Discrete);

        let hamiltonian = Zero;

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

impl Simulation for CustomRandomWalk {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        self.translate_sweep.apply(
            &mut self.microstate,
            &self.hamiltonian,
            &self.macrostate,
        );
        self.microstate.increment_step();
        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<()> {
    let mut simulation = CustomRandomWalk::new()?;
    for _ in 0..1_000_000 {
        simulation.advance()?;
    }

    Ok(())
}

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

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 to see the simulation in action. 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;

Tetronimoes

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

Overview

There are many ways you can model anisotropic bodies in hoomd-rs. This tutorial shows you how to place multiple sites in bodies that can translate and rotate. It uses the same site-site and site-field interactions as the Applying Interactions tutorial and the lattice moves of the Custom Random Walk tutorial. Combine those with tetronimo-shaped bodies, and you get something very interesting.

• Objective: Explain how to execute simulations with multi-site bodies.
• File: hoomd-rs/examples/mc-tutorial/tetronimoes.rs
• Run (interactively):

  cargo run --release --features "bevy" --example tetronimoes
  

• Run (in batch mode):

  cargo run --release --example tetronimoes
  

Use Declarations

use rand::{Rng, RngExt, seq::IndexedRandom};
use std::f64::consts::PI;

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

type PositionVector = Cartesian<2>;
type BodyProperties = OrientedPoint<PositionVector, Angle>;
type SiteProperties = Point<PositionVector>;

/// Take fixed steps left, right, down, up, rotate left, or rotate right.
struct DiscreteRotateOrTranslate;

impl LocalTrial<BodyProperties> for DiscreteRotateOrTranslate {
    fn propose<R: Rng>(
        &self,
        rng: &mut R,
        body_properties: BodyProperties,
    ) -> BodyProperties {
        let translate_steps = [
            [0.0, -1.0].into(),
            [0.0, 1.0].into(),
            [-1.0, 0.0].into(),
            [1.0, 0.0].into(),
        ];
        let rotate_steps = [-PI / 2.0, PI / 2.0];

        let mut trial = body_properties;
        if rng.random_bool(0.9) {
            trial.position += *translate_steps
                .choose(rng)
                .expect("translate_steps should have at least 1 element");
        } else {
            trial.orientation.theta += *rotate_steps
                .choose(rng)
                .expect("rotate_steps should have at least 1 element");
        }
        trial
    }
}

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

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

impl Tetronimoes {
    /// Construct a new tetronimo simulation.
    fn new() -> anyhow::Result<Tetronimoes> {
        let box_height = 30.0;
        let macrostate = Isothermal { temperature: 1.0 };
        let alpha = 1.0;
        let epsilon = 1000.0;
        let sigma = 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 vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let square = Rectangle::with_equal_edges(box_height.try_into()?);
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        let sweep = Sweep(DiscreteRotateOrTranslate);

        let template_sites = vec![
            // square
            vec![
                Point::new([-0.5, -0.5].into()),
                Point::new([0.5, -0.5].into()),
                Point::new([0.5, 0.5].into()),
                Point::new([-0.5, 0.5].into()),
            ],
            // line
            vec![
                Point::new([-1.5, 0.5].into()),
                Point::new([-0.5, 0.5].into()),
                Point::new([0.5, 0.5].into()),
                Point::new([1.5, 0.5].into()),
            ],
            // T
            vec![
                Point::new([-1.5, -0.5].into()),
                Point::new([-0.5, -0.5].into()),
                Point::new([0.5, -0.5].into()),
                Point::new([-0.5, 0.5].into()),
            ],
            // L1
            vec![
                Point::new([-1.5, -0.5].into()),
                Point::new([-0.5, -0.5].into()),
                Point::new([0.5, -0.5].into()),
                Point::new([0.5, 0.5].into()),
            ],
            // L2
            vec![
                Point::new([-1.5, 0.5].into()),
                Point::new([-0.5, 0.5].into()),
                Point::new([0.5, 0.5].into()),
                Point::new([0.5, -0.5].into()),
            ],
        ];

        Ok(Tetronimoes {
            microstate,
            hamiltonian,
            sweep,
            macrostate,
            template_sites,
        })
    }
}

impl Simulation for Tetronimoes {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        if self.microstate.step().is_multiple_of(100) {
            let mut rng = self.microstate.counter().make_rng();
            let sites = self
                .template_sites
                .choose(&mut rng)
                .expect("template_sites should have at least 1 element")
                .clone();

            let properties = OrientedPoint {
                position: [
                    0.0,
                    self.microstate.boundary().0.edge_lengths[1].get() / 2.0
                        - 2.0,
                ]
                .into(),
                orientation: Angle::from(0.0),
            };

            self.microstate.add_body(Body { sites, properties })?;
            self.microstate.increment_substep();
        }

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

        if self.hamiltonian.linear.total_energy(&self.microstate) > 100.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 = Tetronimoes::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("tetronimoes.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 tetronimoes_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use tetronimoes_interactive::main;

Type Aliases

Create type aliases for your model’s vector, body properties, and site properties types so that you don’t need to repeat the full generic type names throughout the code:

use rand::{Rng, RngExt, seq::IndexedRandom};
use std::f64::consts::PI;

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

type PositionVector = Cartesian<2>;
type BodyProperties = OrientedPoint<PositionVector, Angle>;
type SiteProperties = Point<PositionVector>;

/// Take fixed steps left, right, down, up, rotate left, or rotate right.
struct DiscreteRotateOrTranslate;

impl LocalTrial<BodyProperties> for DiscreteRotateOrTranslate {
    fn propose<R: Rng>(
        &self,
        rng: &mut R,
        body_properties: BodyProperties,
    ) -> BodyProperties {
        let translate_steps = [
            [0.0, -1.0].into(),
            [0.0, 1.0].into(),
            [-1.0, 0.0].into(),
            [1.0, 0.0].into(),
        ];
        let rotate_steps = [-PI / 2.0, PI / 2.0];

        let mut trial = body_properties;
        if rng.random_bool(0.9) {
            trial.position += *translate_steps
                .choose(rng)
                .expect("translate_steps should have at least 1 element");
        } else {
            trial.orientation.theta += *rotate_steps
                .choose(rng)
                .expect("rotate_steps should have at least 1 element");
        }
        trial
    }
}

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

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

impl Tetronimoes {
    /// Construct a new tetronimo simulation.
    fn new() -> anyhow::Result<Tetronimoes> {
        let box_height = 30.0;
        let macrostate = Isothermal { temperature: 1.0 };
        let alpha = 1.0;
        let epsilon = 1000.0;
        let sigma = 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 vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let square = Rectangle::with_equal_edges(box_height.try_into()?);
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        let sweep = Sweep(DiscreteRotateOrTranslate);

        let template_sites = vec![
            // square
            vec![
                Point::new([-0.5, -0.5].into()),
                Point::new([0.5, -0.5].into()),
                Point::new([0.5, 0.5].into()),
                Point::new([-0.5, 0.5].into()),
            ],
            // line
            vec![
                Point::new([-1.5, 0.5].into()),
                Point::new([-0.5, 0.5].into()),
                Point::new([0.5, 0.5].into()),
                Point::new([1.5, 0.5].into()),
            ],
            // T
            vec![
                Point::new([-1.5, -0.5].into()),
                Point::new([-0.5, -0.5].into()),
                Point::new([0.5, -0.5].into()),
                Point::new([-0.5, 0.5].into()),
            ],
            // L1
            vec![
                Point::new([-1.5, -0.5].into()),
                Point::new([-0.5, -0.5].into()),
                Point::new([0.5, -0.5].into()),
                Point::new([0.5, 0.5].into()),
            ],
            // L2
            vec![
                Point::new([-1.5, 0.5].into()),
                Point::new([-0.5, 0.5].into()),
                Point::new([0.5, 0.5].into()),
                Point::new([0.5, -0.5].into()),
            ],
        ];

        Ok(Tetronimoes {
            microstate,
            hamiltonian,
            sweep,
            macrostate,
            template_sites,
        })
    }
}

impl Simulation for Tetronimoes {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        if self.microstate.step().is_multiple_of(100) {
            let mut rng = self.microstate.counter().make_rng();
            let sites = self
                .template_sites
                .choose(&mut rng)
                .expect("template_sites should have at least 1 element")
                .clone();

            let properties = OrientedPoint {
                position: [
                    0.0,
                    self.microstate.boundary().0.edge_lengths[1].get() / 2.0
                        - 2.0,
                ]
                .into(),
                orientation: Angle::from(0.0),
            };

            self.microstate.add_body(Body { sites, properties })?;
            self.microstate.increment_substep();
        }

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

        if self.hamiltonian.linear.total_energy(&self.microstate) > 100.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 = Tetronimoes::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("tetronimoes.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 tetronimoes_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use tetronimoes_interactive::main;

The OrientedPoint type gives the tetronimo bodies in this tutorial both a position in space and an orientation that rotates about the origin of the body. These tetronimoes have disks at each site and therefore only need Point site properties.

Custom Trial Move

Implement a custom trial move to make the tetronimoes move in the way you might expect. Like in the Custom Random Walk, tetronimoes take discrete steps left, right, down, or up. Tetronimoes can also rotate by $\pm \pi /2$. The DiscreteRotateOrTranslate type implements this behavior:

use rand::{Rng, RngExt, seq::IndexedRandom};
use std::f64::consts::PI;

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

type PositionVector = Cartesian<2>;
type BodyProperties = OrientedPoint<PositionVector, Angle>;
type SiteProperties = Point<PositionVector>;

/// Take fixed steps left, right, down, up, rotate left, or rotate right.
struct DiscreteRotateOrTranslate;

impl LocalTrial<BodyProperties> for DiscreteRotateOrTranslate {
    fn propose<R: Rng>(
        &self,
        rng: &mut R,
        body_properties: BodyProperties,
    ) -> BodyProperties {
        let translate_steps = [
            [0.0, -1.0].into(),
            [0.0, 1.0].into(),
            [-1.0, 0.0].into(),
            [1.0, 0.0].into(),
        ];
        let rotate_steps = [-PI / 2.0, PI / 2.0];

        let mut trial = body_properties;
        if rng.random_bool(0.9) {
            trial.position += *translate_steps
                .choose(rng)
                .expect("translate_steps should have at least 1 element");
        } else {
            trial.orientation.theta += *rotate_steps
                .choose(rng)
                .expect("rotate_steps should have at least 1 element");
        }
        trial
    }
}

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

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

impl Tetronimoes {
    /// Construct a new tetronimo simulation.
    fn new() -> anyhow::Result<Tetronimoes> {
        let box_height = 30.0;
        let macrostate = Isothermal { temperature: 1.0 };
        let alpha = 1.0;
        let epsilon = 1000.0;
        let sigma = 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 vec_cell = VecCell::builder()
            .nominal_search_radius(
                hamiltonian.maximum_interaction_range().try_into()?,
            )
            .build();
        let square = Rectangle::with_equal_edges(box_height.try_into()?);
        let microstate = Microstate::builder()
            .spatial_data(vec_cell)
            .boundary(Closed(square))
            .try_build()?;

        let sweep = Sweep(DiscreteRotateOrTranslate);

        let template_sites = vec![
            // square
            vec![
                Point::new([-0.5, -0.5].into()),
                Point::new([0.5, -0.5].into()),
                Point::new([0.5, 0.5].into()),
                Point::new([-0.5, 0.5].into()),
            ],
            // line
            vec![
                Point::new([-1.5, 0.5].into()),
                Point::new([-0.5, 0.5].into()),
                Point::new([0.5, 0.5].into()),
                Point::new([1.5, 0.5].into()),
            ],
            // T
            vec![
                Point::new([-1.5, -0.5].into()),
                Point::new([-0.5, -0.5].into()),
                Point::new([0.5, -0.5].into()),
                Point::new([-0.5, 0.5].into()),
            ],
            // L1
            vec![
                Point::new([-1.5, -0.5].into()),
                Point::new([-0.5, -0.5].into()),
                Point::new([0.5, -0.5].into()),
                Point::new([0.5, 0.5].into()),
            ],
            // L2
            vec![
                Point::new([-1.5, 0.5].into()),
                Point::new([-0.5, 0.5].into()),
                Point::new([0.5, 0.5].into()),
                Point::new([0.5, -0.5].into()),
            ],
        ];

        Ok(Tetronimoes {
            microstate,
            hamiltonian,
            sweep,
            macrostate,
            template_sites,
        })
    }
}

impl Simulation for Tetronimoes {
    /// Advance the simulation forward one step.
    fn advance(&mut self) -> anyhow::Result<()> {
        if self.microstate.step().is_multiple_of(100) {
            let mut rng = self.microstate.counter().make_rng();
            let sites = self
                .template_sites
                .choose(&mut rng)
                .expect("template_sites should have at least 1 element")
                .clone();

            let properties = OrientedPoint {
                position: [
                    0.0,
                    self.microstate.boundary().0.edge_lengths[1].get() / 2.0
                        - 2.0,
                ]
                .into(),
                orientation: Angle::from(0.0),
            };

            self.microstate.add_body(Body { sites, properties })?;
            self.microstate.increment_substep();
        }

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

        if self.hamiltonian.linear.total_energy(&self.microstate) > 100.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 = Tetronimoes::new()?;
    let mut hoomd_gsd_file = HoomdGsdFile::create("tetronimoes.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 tetronimoes_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use tetronimoes_interactive::main;