Hard Tetrahedron Self-Assembly
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 represents sites with hard convex polyhedra. When compressed to a sufficiently high packing fraction, systems of hard tetrahedra self-assemble into a quasicrystal: 10.1038/nature08641.
- Objectives:
- Explain how to model system of hard convex polytopes.
- Demonstrate the self-assembly of hard tetrahedra.
- File:
hoomd-rs/examples/mc-tutorial/hard-tetrahedron-self-assembly.rs - Run (interactively):
cargo run --release --features "bevy" --example hard-tetrahedron-self-assembly - Run (in batch mode):
cargo run --release --example hard-tetrahedron-self-assembly
Use Declarations
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Convex, Volume,
shape::{ConvexPolyhedron, Cuboid},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff,
pairwise::{Anisotropic, ApproximateShapeOverlap, HardShape},
univariate::OverlapPenalty,
};
use hoomd_mc::{
QuickCompress, QuickInsert, Rotate, Sweep, Translate, Trial, Tune,
TuneOptions, UniformIn,
};
use hoomd_microstate::{
Microstate, SiteKey, boundary::Periodic, property::OrientedPoint,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{self, Cartesian, Versor};
type PositionVector = Cartesian<3>;
type Orientation = Versor;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
impl HardTetrahedronSelfAssembly {
/// Construct a new hard tetrahedron self-assembly simulation.
fn new() -> anyhow::Result<HardTetrahedronSelfAssembly> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.50;
let n_bodies = 256;
let maximum_distance = 0.04;
let maximum_rotation = 0.04;
let macrostate = Isothermal { temperature: 1.0 };
let a = 1.0_f64;
let h = 6.0_f64.sqrt() / 3.0 * a;
let tetrahedron_volume = 1.0 / 12.0 * 2.0_f64.sqrt() * a.powi(3);
let tetrahedron = ConvexPolyhedron::with_vertices(vec![
[3.0_f64.sqrt() / 3.0 * a, 0.0, -h / 4.0].into(),
[-3.0_f64.sqrt() / 6.0 * a, 0.5 * a, -h / 4.0].into(),
[-3.0_f64.sqrt() / 6.0 * a, -0.5 * a, -h / 4.0].into(),
[0.0, 0.0, 3.0 * h / 4.0].into(),
])?;
let hamiltonian =
PairwiseCutoff(HardShape(Convex(tetrahedron.clone())));
let initial_box_volume =
n_bodies as f64 * tetrahedron_volume / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.cbrt();
let cube =
Cuboid::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_cube =
Periodic::new(hamiltonian.maximum_interaction_range(), cube)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_cube)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let rotate =
Rotate::with_maximum_rotation(maximum_rotation.try_into()?);
let rotate_sweep = Sweep(rotate);
let distribution = UniformIn {
boundary: microstate.boundary().clone(),
template_sites: vec![SiteProperties::default()],
};
let quick_insert = QuickInsert::new(distribution, n_bodies);
let target_box_volume =
n_bodies as f64 * tetrahedron_volume / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
Convex(tetrahedron),
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: hamiltonian.maximum_interaction_range(),
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardTetrahedronSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardTetrahedronSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 3>,
Periodic<Cuboid>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Convex<ConvexPolyhedron>>>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Trial moves to apply.
rotate_sweep: Sweep<Rotate<Orientation>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm.
quick_compress: QuickCompress<Periodic<Cuboid>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Cuboid>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<
ApproximateShapeOverlap<OverlapPenalty, Convex<ConvexPolyhedron>>,
>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for HardTetrahedronSelfAssembly {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl HardTetrahedronSelfAssembly {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
self.rotate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.rotate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
self.rotate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = HardTetrahedronSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-tetrahedron-self-assembly.gsd")?;
for _ in 0..40_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.end()?;
}
}
Ok(())
}
#[cfg(feature = "bevy")]
mod hard_tetrahedron_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_tetrahedron_self_assembly_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 anyhow::{Context, anyhow};
use hoomd_geometry::{
Convex, Volume,
shape::{ConvexPolyhedron, Cuboid},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff,
pairwise::{Anisotropic, ApproximateShapeOverlap, HardShape},
univariate::OverlapPenalty,
};
use hoomd_mc::{
QuickCompress, QuickInsert, Rotate, Sweep, Translate, Trial, Tune,
TuneOptions, UniformIn,
};
use hoomd_microstate::{
Microstate, SiteKey, boundary::Periodic, property::OrientedPoint,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{self, Cartesian, Versor};
type PositionVector = Cartesian<3>;
type Orientation = Versor;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
impl HardTetrahedronSelfAssembly {
/// Construct a new hard tetrahedron self-assembly simulation.
fn new() -> anyhow::Result<HardTetrahedronSelfAssembly> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.50;
let n_bodies = 256;
let maximum_distance = 0.04;
let maximum_rotation = 0.04;
let macrostate = Isothermal { temperature: 1.0 };
let a = 1.0_f64;
let h = 6.0_f64.sqrt() / 3.0 * a;
let tetrahedron_volume = 1.0 / 12.0 * 2.0_f64.sqrt() * a.powi(3);
let tetrahedron = ConvexPolyhedron::with_vertices(vec![
[3.0_f64.sqrt() / 3.0 * a, 0.0, -h / 4.0].into(),
[-3.0_f64.sqrt() / 6.0 * a, 0.5 * a, -h / 4.0].into(),
[-3.0_f64.sqrt() / 6.0 * a, -0.5 * a, -h / 4.0].into(),
[0.0, 0.0, 3.0 * h / 4.0].into(),
])?;
let hamiltonian =
PairwiseCutoff(HardShape(Convex(tetrahedron.clone())));
let initial_box_volume =
n_bodies as f64 * tetrahedron_volume / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.cbrt();
let cube =
Cuboid::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_cube =
Periodic::new(hamiltonian.maximum_interaction_range(), cube)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_cube)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let rotate =
Rotate::with_maximum_rotation(maximum_rotation.try_into()?);
let rotate_sweep = Sweep(rotate);
let distribution = UniformIn {
boundary: microstate.boundary().clone(),
template_sites: vec![SiteProperties::default()],
};
let quick_insert = QuickInsert::new(distribution, n_bodies);
let target_box_volume =
n_bodies as f64 * tetrahedron_volume / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
Convex(tetrahedron),
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: hamiltonian.maximum_interaction_range(),
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardTetrahedronSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardTetrahedronSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 3>,
Periodic<Cuboid>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Convex<ConvexPolyhedron>>>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Trial moves to apply.
rotate_sweep: Sweep<Rotate<Orientation>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm.
quick_compress: QuickCompress<Periodic<Cuboid>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Cuboid>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<
ApproximateShapeOverlap<OverlapPenalty, Convex<ConvexPolyhedron>>,
>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for HardTetrahedronSelfAssembly {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl HardTetrahedronSelfAssembly {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
self.rotate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.rotate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
self.rotate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = HardTetrahedronSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-tetrahedron-self-assembly.gsd")?;
for _ in 0..40_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.end()?;
}
}
Ok(())
}
#[cfg(feature = "bevy")]
mod hard_tetrahedron_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_tetrahedron_self_assembly_interactive::main;The sites are in this tutorial are represented by tetrahedra with both
position and orientation. Therefore, use OrientedPoint for both the body
and site properties. Use Versor to represent rotations in 3D.
Construct the Simulation Model
The new() method constructs a new simulation model:
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Convex, Volume,
shape::{ConvexPolyhedron, Cuboid},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff,
pairwise::{Anisotropic, ApproximateShapeOverlap, HardShape},
univariate::OverlapPenalty,
};
use hoomd_mc::{
QuickCompress, QuickInsert, Rotate, Sweep, Translate, Trial, Tune,
TuneOptions, UniformIn,
};
use hoomd_microstate::{
Microstate, SiteKey, boundary::Periodic, property::OrientedPoint,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{self, Cartesian, Versor};
type PositionVector = Cartesian<3>;
type Orientation = Versor;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
impl HardTetrahedronSelfAssembly {
/// Construct a new hard tetrahedron self-assembly simulation.
fn new() -> anyhow::Result<HardTetrahedronSelfAssembly> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.50;
let n_bodies = 256;
let maximum_distance = 0.04;
let maximum_rotation = 0.04;
let macrostate = Isothermal { temperature: 1.0 };
let a = 1.0_f64;
let h = 6.0_f64.sqrt() / 3.0 * a;
let tetrahedron_volume = 1.0 / 12.0 * 2.0_f64.sqrt() * a.powi(3);
let tetrahedron = ConvexPolyhedron::with_vertices(vec![
[3.0_f64.sqrt() / 3.0 * a, 0.0, -h / 4.0].into(),
[-3.0_f64.sqrt() / 6.0 * a, 0.5 * a, -h / 4.0].into(),
[-3.0_f64.sqrt() / 6.0 * a, -0.5 * a, -h / 4.0].into(),
[0.0, 0.0, 3.0 * h / 4.0].into(),
])?;
let hamiltonian =
PairwiseCutoff(HardShape(Convex(tetrahedron.clone())));
let initial_box_volume =
n_bodies as f64 * tetrahedron_volume / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.cbrt();
let cube =
Cuboid::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_cube =
Periodic::new(hamiltonian.maximum_interaction_range(), cube)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_cube)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let rotate =
Rotate::with_maximum_rotation(maximum_rotation.try_into()?);
let rotate_sweep = Sweep(rotate);
let distribution = UniformIn {
boundary: microstate.boundary().clone(),
template_sites: vec![SiteProperties::default()],
};
let quick_insert = QuickInsert::new(distribution, n_bodies);
let target_box_volume =
n_bodies as f64 * tetrahedron_volume / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
Convex(tetrahedron),
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: hamiltonian.maximum_interaction_range(),
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardTetrahedronSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardTetrahedronSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 3>,
Periodic<Cuboid>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Convex<ConvexPolyhedron>>>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Trial moves to apply.
rotate_sweep: Sweep<Rotate<Orientation>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm.
quick_compress: QuickCompress<Periodic<Cuboid>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Cuboid>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<
ApproximateShapeOverlap<OverlapPenalty, Convex<ConvexPolyhedron>>,
>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for HardTetrahedronSelfAssembly {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl HardTetrahedronSelfAssembly {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
self.rotate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.rotate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
self.rotate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = HardTetrahedronSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-tetrahedron-self-assembly.gsd")?;
for _ in 0..40_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.end()?;
}
}
Ok(())
}
#[cfg(feature = "bevy")]
mod hard_tetrahedron_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_tetrahedron_self_assembly_interactive::main;Parameters
Assign all the model parameters in one code block:
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Convex, Volume,
shape::{ConvexPolyhedron, Cuboid},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff,
pairwise::{Anisotropic, ApproximateShapeOverlap, HardShape},
univariate::OverlapPenalty,
};
use hoomd_mc::{
QuickCompress, QuickInsert, Rotate, Sweep, Translate, Trial, Tune,
TuneOptions, UniformIn,
};
use hoomd_microstate::{
Microstate, SiteKey, boundary::Periodic, property::OrientedPoint,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{self, Cartesian, Versor};
type PositionVector = Cartesian<3>;
type Orientation = Versor;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
impl HardTetrahedronSelfAssembly {
/// Construct a new hard tetrahedron self-assembly simulation.
fn new() -> anyhow::Result<HardTetrahedronSelfAssembly> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.50;
let n_bodies = 256;
let maximum_distance = 0.04;
let maximum_rotation = 0.04;
let macrostate = Isothermal { temperature: 1.0 };
let a = 1.0_f64;
let h = 6.0_f64.sqrt() / 3.0 * a;
let tetrahedron_volume = 1.0 / 12.0 * 2.0_f64.sqrt() * a.powi(3);
let tetrahedron = ConvexPolyhedron::with_vertices(vec![
[3.0_f64.sqrt() / 3.0 * a, 0.0, -h / 4.0].into(),
[-3.0_f64.sqrt() / 6.0 * a, 0.5 * a, -h / 4.0].into(),
[-3.0_f64.sqrt() / 6.0 * a, -0.5 * a, -h / 4.0].into(),
[0.0, 0.0, 3.0 * h / 4.0].into(),
])?;
let hamiltonian =
PairwiseCutoff(HardShape(Convex(tetrahedron.clone())));
let initial_box_volume =
n_bodies as f64 * tetrahedron_volume / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.cbrt();
let cube =
Cuboid::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_cube =
Periodic::new(hamiltonian.maximum_interaction_range(), cube)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_cube)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let rotate =
Rotate::with_maximum_rotation(maximum_rotation.try_into()?);
let rotate_sweep = Sweep(rotate);
let distribution = UniformIn {
boundary: microstate.boundary().clone(),
template_sites: vec![SiteProperties::default()],
};
let quick_insert = QuickInsert::new(distribution, n_bodies);
let target_box_volume =
n_bodies as f64 * tetrahedron_volume / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
Convex(tetrahedron),
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: hamiltonian.maximum_interaction_range(),
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardTetrahedronSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardTetrahedronSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 3>,
Periodic<Cuboid>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Convex<ConvexPolyhedron>>>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Trial moves to apply.
rotate_sweep: Sweep<Rotate<Orientation>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm.
quick_compress: QuickCompress<Periodic<Cuboid>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Cuboid>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<
ApproximateShapeOverlap<OverlapPenalty, Convex<ConvexPolyhedron>>,
>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for HardTetrahedronSelfAssembly {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl HardTetrahedronSelfAssembly {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
self.rotate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.rotate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
self.rotate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = HardTetrahedronSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-tetrahedron-self-assembly.gsd")?;
for _ in 0..40_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.end()?;
}
}
Ok(())
}
#[cfg(feature = "bevy")]
mod hard_tetrahedron_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_tetrahedron_self_assembly_interactive::main;initial_packing_fraction is the volume of the tetrahedra divided by the
volume of the simulation boundary in the initial state. Choose this value so
that tetrahedra can be placed easily in the microstate. During the Initialize
phase, the microstate will be compressed until it reaches the packing fraction
target_packing_fraction. n_bodies is the number of tetrahedra to add,
maximum_distance is the largest distance a translation trial move can take
(initially), maximum_rotationcontrols the size of the rotation trial moves
(initially), and macrostate holds the temperature set point (in units of
energy).
Hamiltonian
A ConvexPolyhedron is the shape given by the convex hull of the given vertices.
Wrap it in the Convex newtype for use with the HardShape Hamiltonian:
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Convex, Volume,
shape::{ConvexPolyhedron, Cuboid},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff,
pairwise::{Anisotropic, ApproximateShapeOverlap, HardShape},
univariate::OverlapPenalty,
};
use hoomd_mc::{
QuickCompress, QuickInsert, Rotate, Sweep, Translate, Trial, Tune,
TuneOptions, UniformIn,
};
use hoomd_microstate::{
Microstate, SiteKey, boundary::Periodic, property::OrientedPoint,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{self, Cartesian, Versor};
type PositionVector = Cartesian<3>;
type Orientation = Versor;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
impl HardTetrahedronSelfAssembly {
/// Construct a new hard tetrahedron self-assembly simulation.
fn new() -> anyhow::Result<HardTetrahedronSelfAssembly> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.50;
let n_bodies = 256;
let maximum_distance = 0.04;
let maximum_rotation = 0.04;
let macrostate = Isothermal { temperature: 1.0 };
let a = 1.0_f64;
let h = 6.0_f64.sqrt() / 3.0 * a;
let tetrahedron_volume = 1.0 / 12.0 * 2.0_f64.sqrt() * a.powi(3);
let tetrahedron = ConvexPolyhedron::with_vertices(vec![
[3.0_f64.sqrt() / 3.0 * a, 0.0, -h / 4.0].into(),
[-3.0_f64.sqrt() / 6.0 * a, 0.5 * a, -h / 4.0].into(),
[-3.0_f64.sqrt() / 6.0 * a, -0.5 * a, -h / 4.0].into(),
[0.0, 0.0, 3.0 * h / 4.0].into(),
])?;
let hamiltonian =
PairwiseCutoff(HardShape(Convex(tetrahedron.clone())));
let initial_box_volume =
n_bodies as f64 * tetrahedron_volume / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.cbrt();
let cube =
Cuboid::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_cube =
Periodic::new(hamiltonian.maximum_interaction_range(), cube)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_cube)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let rotate =
Rotate::with_maximum_rotation(maximum_rotation.try_into()?);
let rotate_sweep = Sweep(rotate);
let distribution = UniformIn {
boundary: microstate.boundary().clone(),
template_sites: vec![SiteProperties::default()],
};
let quick_insert = QuickInsert::new(distribution, n_bodies);
let target_box_volume =
n_bodies as f64 * tetrahedron_volume / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
Convex(tetrahedron),
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: hamiltonian.maximum_interaction_range(),
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardTetrahedronSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardTetrahedronSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 3>,
Periodic<Cuboid>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Convex<ConvexPolyhedron>>>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Trial moves to apply.
rotate_sweep: Sweep<Rotate<Orientation>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm.
quick_compress: QuickCompress<Periodic<Cuboid>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Cuboid>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<
ApproximateShapeOverlap<OverlapPenalty, Convex<ConvexPolyhedron>>,
>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for HardTetrahedronSelfAssembly {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl HardTetrahedronSelfAssembly {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
self.rotate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.rotate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
self.rotate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = HardTetrahedronSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-tetrahedron-self-assembly.gsd")?;
for _ in 0..40_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.end()?;
}
}
Ok(())
}
#[cfg(feature = "bevy")]
mod hard_tetrahedron_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_tetrahedron_self_assembly_interactive::main;Coxeter is one of many tools you can use to construct the vertices of a convex polyhedron.
Initialization and Simulation
See the Hard Ellipse Self-Assembly tutorial for a complete explanation of remaining initialization and simulation code (see also the complete code below).
Execute the Simulation in Batch Mode
Implement main()
To run the simulation, construct the HardTetrahedronSelfAssembly 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 anyhow::{Context, anyhow};
use hoomd_geometry::{
Convex, Volume,
shape::{ConvexPolyhedron, Cuboid},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff,
pairwise::{Anisotropic, ApproximateShapeOverlap, HardShape},
univariate::OverlapPenalty,
};
use hoomd_mc::{
QuickCompress, QuickInsert, Rotate, Sweep, Translate, Trial, Tune,
TuneOptions, UniformIn,
};
use hoomd_microstate::{
Microstate, SiteKey, boundary::Periodic, property::OrientedPoint,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{self, Cartesian, Versor};
type PositionVector = Cartesian<3>;
type Orientation = Versor;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
impl HardTetrahedronSelfAssembly {
/// Construct a new hard tetrahedron self-assembly simulation.
fn new() -> anyhow::Result<HardTetrahedronSelfAssembly> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.50;
let n_bodies = 256;
let maximum_distance = 0.04;
let maximum_rotation = 0.04;
let macrostate = Isothermal { temperature: 1.0 };
let a = 1.0_f64;
let h = 6.0_f64.sqrt() / 3.0 * a;
let tetrahedron_volume = 1.0 / 12.0 * 2.0_f64.sqrt() * a.powi(3);
let tetrahedron = ConvexPolyhedron::with_vertices(vec![
[3.0_f64.sqrt() / 3.0 * a, 0.0, -h / 4.0].into(),
[-3.0_f64.sqrt() / 6.0 * a, 0.5 * a, -h / 4.0].into(),
[-3.0_f64.sqrt() / 6.0 * a, -0.5 * a, -h / 4.0].into(),
[0.0, 0.0, 3.0 * h / 4.0].into(),
])?;
let hamiltonian =
PairwiseCutoff(HardShape(Convex(tetrahedron.clone())));
let initial_box_volume =
n_bodies as f64 * tetrahedron_volume / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.cbrt();
let cube =
Cuboid::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_cube =
Periodic::new(hamiltonian.maximum_interaction_range(), cube)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_cube)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let rotate =
Rotate::with_maximum_rotation(maximum_rotation.try_into()?);
let rotate_sweep = Sweep(rotate);
let distribution = UniformIn {
boundary: microstate.boundary().clone(),
template_sites: vec![SiteProperties::default()],
};
let quick_insert = QuickInsert::new(distribution, n_bodies);
let target_box_volume =
n_bodies as f64 * tetrahedron_volume / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
Convex(tetrahedron),
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: hamiltonian.maximum_interaction_range(),
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardTetrahedronSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardTetrahedronSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 3>,
Periodic<Cuboid>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Convex<ConvexPolyhedron>>>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Trial moves to apply.
rotate_sweep: Sweep<Rotate<Orientation>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm.
quick_compress: QuickCompress<Periodic<Cuboid>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Cuboid>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<
ApproximateShapeOverlap<OverlapPenalty, Convex<ConvexPolyhedron>>,
>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for HardTetrahedronSelfAssembly {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl HardTetrahedronSelfAssembly {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
self.rotate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.rotate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
self.rotate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = HardTetrahedronSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-tetrahedron-self-assembly.gsd")?;
for _ in 0..40_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.end()?;
}
}
Ok(())
}
#[cfg(feature = "bevy")]
mod hard_tetrahedron_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_tetrahedron_self_assembly_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 differentmain()(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 hard-tetrahedron-self-assembly
Visualize the Simulation Results
Open the generated hard-tetrahedron-self-assembly.gsd in Ovito or another
visualization tool to see the simulation results. To render the tetrahedra in
Ovito, import this modifier snippet (use the copy button to copy the entire
text to your clipboard):
{"description": "OVITO Modifier Snippet: Edit types (Particle Type)", "payload": "AAAHKXjafVU9TBRBFH573HG/cAL+oYVnaAwFosTEBMyJeAUFgeBJYUOG3eF2de/2srsHucTCwspEY6hsDIWJjRY2diRUhoTERkysTIwWFhYYO4mJzpt9cyzrwiTvZu5972/ee/O2+Prdq9yf7TsAresAkBHUD7OwANNQFXsJ5sAFR3C7BBUhWAUYgxG4BFfE7yjxtIBw6+qn/xnaE0QPab9PewFl5lynyV2/Lc4j0AQm/PlggQ42cPDgoojgIK8kImsLSS40Lv+HelA6vybjPlxLLbxTSqM/SUFZLYSmMErJ6JYnPKY7p4zSRVY2qpuLMvKYCo3unCTfPXRG+71os/yzo1DseDomaFAr7+auNV98ykV3AfZhFEGeJ2X9AMbJ8hhlGj1l55jrW7rNPfHnzJEJGhB0+giPx1VSpEO1TpBTXCdVeaEui2KG5DDUXZL9rrK/+aH98m2+eKH8JWBsbpc2rg73/y5/DfgTb0D7/HxrR/GV/GZEfmP97940e7Q+8RTgx8et+fK3Z4+Hzz0ZhyNXkgoJoY5VK03VShAW5WsRW4oPB+1gLRM6yee5YfncqLabshqyzZKE9XpOy9V56F10TJ3STV63dGZXbF7nDX9x0TOZ4axSsoNHVWeeFwcMClaT3/b4DPfMKcd23DipISl1g+n3lpnOp1q2bTVqlQZbsrkRJ3/WtGqmLci/hYoVo8ZjvfdIu3FI34qxOs8MqxWr1+seCg3EuBbsBwrPSpd4W0ysFipjSiLhzHYHbsL9UOTxt0asR4/LHyKFBqvzWBVntcHdaFURSUljKnLkJNFIeA7lVyyPai5iTK5p9AARy1mGYFvLFnfDXfeL2moQU1wR7SabbcYxlGSaBPAxF9RsQKEgiH2s6lqsUbMpk+hTYjgk8hQUqSkIR1DuJvPZ7NJdrvvBuJQIjrU0IguWFwzOhJpxmVBsShxHW3aeL1eZW+M+3a4T2KTuWys87EO+u9H9z5FcC6FPEK496gR6pIlp+hqYNLfwIjs01THP7xH8BwoEH1E="}
Alternately, you can:
- Download tetrahedron.obj.
- Add an Edit types modification to the Ovito pipeline.
- Choose the
Mesh/user-definedshape. - Set the Display radius to 1.
- Click Load geometry file… and choose
tetrahedron.obj. To save time when visualizing many similar systems, you can save the session state or use modifier templates.
Render with Tachyon, and you should see something like:
Tip
Use Coxeter, Blender, OnShape, or any modeling tool you prefer to create
.objfiles that best represent your model’s sites.
Conclusion
This tutorial showed you how to perform hard tetrahedron self-assembly simulations using a shape overlap potential.
Navigate to the top of the page to see the simulation in action. Notice that tetrahedra are first added in a large batch. Once all the overlaps are removed, another batch appears. After all all tetrahedra are in the microstate and not overlapping, the simulation compresses to a higher packing fraction. After that, it speeds up as it begins using the more efficient hard particle overlap Hamiltonian. Watch the simulation long enough and you should see dimers and pentamers form. These motifs organize to form a quasicrystal, but only after very long simulation times with at least 4096 tetrahedra: 10.1038/nature08641.
Complete Code
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Convex, Volume,
shape::{ConvexPolyhedron, Cuboid},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff,
pairwise::{Anisotropic, ApproximateShapeOverlap, HardShape},
univariate::OverlapPenalty,
};
use hoomd_mc::{
QuickCompress, QuickInsert, Rotate, Sweep, Translate, Trial, Tune,
TuneOptions, UniformIn,
};
use hoomd_microstate::{
Microstate, SiteKey, boundary::Periodic, property::OrientedPoint,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{self, Cartesian, Versor};
type PositionVector = Cartesian<3>;
type Orientation = Versor;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
impl HardTetrahedronSelfAssembly {
/// Construct a new hard tetrahedron self-assembly simulation.
fn new() -> anyhow::Result<HardTetrahedronSelfAssembly> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.50;
let n_bodies = 256;
let maximum_distance = 0.04;
let maximum_rotation = 0.04;
let macrostate = Isothermal { temperature: 1.0 };
let a = 1.0_f64;
let h = 6.0_f64.sqrt() / 3.0 * a;
let tetrahedron_volume = 1.0 / 12.0 * 2.0_f64.sqrt() * a.powi(3);
let tetrahedron = ConvexPolyhedron::with_vertices(vec![
[3.0_f64.sqrt() / 3.0 * a, 0.0, -h / 4.0].into(),
[-3.0_f64.sqrt() / 6.0 * a, 0.5 * a, -h / 4.0].into(),
[-3.0_f64.sqrt() / 6.0 * a, -0.5 * a, -h / 4.0].into(),
[0.0, 0.0, 3.0 * h / 4.0].into(),
])?;
let hamiltonian =
PairwiseCutoff(HardShape(Convex(tetrahedron.clone())));
let initial_box_volume =
n_bodies as f64 * tetrahedron_volume / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.cbrt();
let cube =
Cuboid::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_cube =
Periodic::new(hamiltonian.maximum_interaction_range(), cube)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_cube)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let rotate =
Rotate::with_maximum_rotation(maximum_rotation.try_into()?);
let rotate_sweep = Sweep(rotate);
let distribution = UniformIn {
boundary: microstate.boundary().clone(),
template_sites: vec![SiteProperties::default()],
};
let quick_insert = QuickInsert::new(distribution, n_bodies);
let target_box_volume =
n_bodies as f64 * tetrahedron_volume / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
Convex(tetrahedron),
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: hamiltonian.maximum_interaction_range(),
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardTetrahedronSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardTetrahedronSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 3>,
Periodic<Cuboid>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Convex<ConvexPolyhedron>>>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Trial moves to apply.
rotate_sweep: Sweep<Rotate<Orientation>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm.
quick_compress: QuickCompress<Periodic<Cuboid>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Cuboid>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<
ApproximateShapeOverlap<OverlapPenalty, Convex<ConvexPolyhedron>>,
>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for HardTetrahedronSelfAssembly {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl HardTetrahedronSelfAssembly {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
self.rotate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.rotate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
self.rotate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = HardTetrahedronSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-tetrahedron-self-assembly.gsd")?;
for _ in 0..40_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.end()?;
}
}
Ok(())
}
#[cfg(feature = "bevy")]
mod hard_tetrahedron_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_tetrahedron_self_assembly_interactive::main;Development of hoomd-rs is led by the Glotzer Group at the University of Michigan.
Copyright © 2024-2026 The Regents of the University of Michigan.
