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;

Enumerate Possible Moves

First, enumerate the possible moves:

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;

Choose and Propose a Move

Then, choose a random move and mutates the body properties accordingly:

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;

Using random_bool(), this code proposes translate moves more often than rotate moves because the result is more visually interesting.

The Simulation Model

Here is the type that holds the simulation model along with the corresponding Hamiltonian:

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;

Construct the Simulation Model

The new() method constructs a new simulation model:

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;

Parameters

Assign all the model parameters in one code block:

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;

Hamiltonian

Use the same Hamiltonian as the Applying Interactions tutorial:

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;

These interactions are applied to the sites in the microstate.

Microstate

Use the VecCell spatial data structure, confine the bodies and sites inside of a closed square, and start with no bodies in the microstate:

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;

Trial Moves

Apply sweeps of the custom DiscreteRotateOrTranslate trial move:

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;

Prepare the Tetronimo Shapes

new() also prepares a list of tetronimo shapes for later use. There are five types of tetronimoes that are each represented by a vector of four points:

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;

Initialize the Struct

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;

Implement Simulation

The Simulation implementation closely follows that in Applying Interactions.

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;

Advance the Simulation

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;

Add New Tetronimoes

The code that adds tetronimoes is more complex than that for disks:

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;

It first chooses a random tetronimo from the template_sites, then it adds the body near the top of the boundary with a default orientation of $\theta =0$ and clone of the chosen sites. Each body has four sites in this example.

hoomd-rs uses a counter based random number generator. Whenever you need to use random numbers in your code, you can get a Rng to generate them by calling microstate.counter().make_rng().

Important

Whenever you use counter.make_rng, You MUST indicate that your substep is complete by calling microstate.increment_substep() so that the next substep will use a different set of random numbers.

Apply Trial Moves

Apply the custom trial move to each body in the microstate:

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;

Reset the Simulation

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;

Get the Simulation Step

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;

Execute the Simulation in Batch Mode

Implement main()

To run the simulation, construct the Tetronimoes simulation model. Then call advance() many times and write the sites to a GSD file periodically so that you can inspect the results of the simulation:

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;

See Applying Interactions for a step-by-step explanation of this code.

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 tetronimoes

Visualize the Simulation Results

Open the generated tetronimoes.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 add bodies with multiple sites and how they can be translated and rotated by trial moves.

Navigate to the top of the page to see the simulation in action. Notice how the randomly generated tetronimoes fall to the bottom while randomly rotating.

Complete 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;