Hard Ellipse 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 ellipses. You can apply the same techniques to any hard shape. When compressed to a sufficiently high packing fraction, systems of hard ellipses self-assemble into the nematic phase.
- Objectives:
- Show how to quickly insert bodies into the microstate.
- Demonstrate the self-assembly of hard ellipses into the nematic phase.
- File:
hoomd-rs/examples/mc-tutorial/hard-ellipse-self-assembly.rs - Run (interactively):
cargo run --release --features "bevy" --example hard-ellipse-self-assembly - Run (in batch mode):
cargo run --release --example hard-ellipse-self-assembly
Use Declarations
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_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::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_self_assembly_interactive::main;The sites are in this tutorial are represented by ellipses with both
position and orientation. Therefore, use OrientedPoint for both the body
and site properties.
The Simulation Model
Here is the type that holds the simulation model:
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_self_assembly_interactive::main;The simulation will consist of two phases:
- Initialize: Add new ellipses to the microstate and compress to the target packing fraction.
- Equilibrate: Perform hard particle Monte Carlo to self-assembly the nematic phase.
The phase field tracks the current phase of the simulation. It stores an enum:
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_self_assembly_interactive::main;Construct the Simulation Model
The new() method constructs a new simulation model:
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_self_assembly_interactive::main;Parameters
Assign all the model parameters in one code block:
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_self_assembly_interactive::main;initial_packing_fraction is the volume of the ellipses divided by the volume
of the simulation boundary in the initial state. Choose this value so that
ellipses 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 ellipses to add,
maximum_distance is the largest distance a translation trial move can take
(initially), maximum_rotation is the largest angle possible in a rotation
trial move (initially), sigma is the major axis of the ellipse, aspect is
the ellipse aspect ratio and macrostate holds the temperature set point (in
units of energy).
To ensure that sigma is the major axis, aspect must be greater than or equal
to 1.0.
Hamiltonian
HardShape represents each site with the given shape. The site pair
energy is infinite when the two sites overlap and 0 when they do
not. Use PairwiseCutoff with the HardShape interaction as the Hamiltonian:
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_self_assembly_interactive::main;Periodic Boundary Conditions
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_self_assembly_interactive::main;For the case of hard ellipses, the largest distance between the
centers of two potentially overlapping ellipses is sigma — when two
ellipses a distance sigma apart rotated so their their long axes just touch.
HardShape computes this in the maximum_interaction_range() method.
Microstate
Construct a microstate with the periodic boundary conditions and the VecCell
spatial data structure:
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_self_assembly_interactive::main;Start with no bodies in the microstate.
Trial Moves
Apply both Translate and Rotate trial moves to the bodies:
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_self_assembly_interactive::main;In 2D simulations, Rotate uniformly selects a random angle between
-maximum_rotation and maximum_rotation and rotates the body by that angle.
Quickly Insert Bodies
The Hard Disk Self-Assembly tutorial placed disks on a square lattice. While
the same can be done for hard ellipses, the process is not as simple. Instead,
this tutorial uses QuickInsert to randomly place the ellipses.
QuickInsert will add up to n_bodies new bodies to the microstate
drawn randomly from the given distribution. UniformIn generates bodies
with positions uniformly distributed in the given boundary and orientations
uniformly distributed among all possible orientations:
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_self_assembly_interactive::main;UniformIn clones template_sites for each new body. In this case, a body
is represented by a single site at the body’s origin and has a default
orientation.
Quickly Compress the Microstate
QuickCompress will irreversibly scale the microstate toward the target boundary
volume:
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff,
pairwise::{HardSphere, Isotropic},
univariate::{Expanded, OverlapPenalty},
};
use hoomd_mc::{QuickCompress, Sweep, Translate, Trial};
use hoomd_microstate::{
Body, Microstate, SiteKey, boundary::Periodic, property::Point,
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::Cartesian;
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
type SiteProperties = Point<PositionVector>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardDiskSelfAssembly {
/// Positions of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardSphere>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// How sites interact during compression.
overlap_penalty_hamiltonian:
PairwiseCutoff<Isotropic<Expanded<OverlapPenalty>>>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Compress,
Equilibrate,
}
impl HardDiskSelfAssembly {
/// Construct a new hard disk self-assembly simulation.
fn new() -> anyhow::Result<HardDiskSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.76;
let n_disks = 64_usize.pow(2);
let maximum_distance = 0.07;
let sigma = 1.0;
let macrostate = Isothermal { temperature: 1.0 };
let hamiltonian = PairwiseCutoff(HardSphere { diameter: sigma });
let circle = Circle {
radius: (sigma / 2.0).try_into()?,
};
let initial_box_volume =
n_disks as f64 * circle.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let mut microstate = Microstate::builder()
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let n_on_side_f64 = (n_disks as f64).sqrt().ceil();
let a = initial_box_edge_length / n_on_side_f64;
let n_on_side = n_on_side_f64 as usize;
for j in 0..n_on_side {
let y = -initial_box_edge_length / 2.0 + j as f64 * a;
for i in 0..n_on_side {
let x = -initial_box_edge_length / 2.0 + i as f64 * a;
if microstate.bodies().len() < n_disks {
microstate
.add_body(Body::point(Cartesian::from([x, y])))?;
}
}
}
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume =
n_disks as f64 * circle.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let overlap_penalty = Isotropic {
interaction: Expanded {
delta: sigma,
f: OverlapPenalty::default(),
},
r_cut: sigma,
};
let overlap_penalty_hamiltonian = PairwiseCutoff(overlap_penalty);
Ok(HardDiskSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
quick_compress,
macrostate,
phase: Phase::Compress,
})
}
}
impl Simulation for HardDiskSelfAssembly {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Compress => self.compress().context("failed to compress")?,
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl HardDiskSelfAssembly {
fn compress(&mut self) -> anyhow::Result<()> {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.phase = Phase::Equilibrate;
println!(
"Compression complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 10_000 {
let current = self.microstate.boundary().volume();
let target = self.quick_compress.target_volume();
let step = self.microstate.step();
return Err(anyhow!(
"Achieved volume {current} after {step} steps. The target was {target}."
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_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 = HardDiskSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-disk-self-assembly.gsd")?;
for _ in 0..100_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_disk_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_disk_self_assembly_interactive::main;Overlap Penalty Hamiltonian
QuickInsert will only add a body when the change in energy due to the addition
is finite. You could use the hard particle hamiltonian with QuickInsert
and ensure that no ellipses in the microstate overlap. However, random body
insertions do not pack densely.
One way around this problem is to allow bodies to overlap a little when
inserted and allow later translate and rotate trial moves to remove that
overlap. The OverlapPenalty potential consists of an infinite energy core
followed by a harmonic potential added to a step function. The infinite core
prevents inserted bodies from overlapping too much, the harmonic potential
encourages the trial moves to separate bodies, and the step function prevents
the trial moves from moving non-overlapping sites into overlapping configurations.
Express this Hamiltonian using PairwiseCutoff with an Anisotropic interaction:
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_self_assembly_interactive::main;ApproximateShapeOverlap computes the approximate amount of overlap between
a pair of shapes, OverlapPenalty applies the potential described above,
and the Anisotropic PairwiseCutoff computes this potential on pairs of
sites.
Important
Use
ApproximateShapeOverlaponly to remove overlaps during initialization. It does not compute the exact amount of overlap and is therefore not appropriate for use in production sampling.
Initialize the Struct
Package all these values into a struct to represent the simulation:
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_self_assembly_interactive::main;Begin the simulation in the Initialize phase.
Implement Simulation
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_self_assembly_interactive::main;Advance the Simulation
advance calls self.initialize() to advance the simulation when in the
initialization phase and self.equilibrate() when in the equilibration phase:
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_self_assembly_interactive::main;Each method advances the simulation one step and potentially changes the
phase. The initialize method might return an error (see below). The anyhow
method context adds additional information to the error message.
Get the simulation step
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_self_assembly_interactive::main;Implement HardEllipseSelfAssembly
Place the model-specific methods in the inherent implementation for
HardEllipseSelfAssembly:
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_self_assembly_interactive::main;Initialize
Implement the initialize phase:
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_self_assembly_interactive::main;Add New Bodies and Compress
The quick_insert.apply method adds new randomly placed bodies to the microstate and
quick_compress.apply irreversibly compresses the microstate toward the
target boundary volume. Insert all n_bodies ellipses first, then compress
to the target volume:
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_self_assembly_interactive::main;Separate Overlapping Bodies
QuickInsert is quick because it only inserts bodies (it never
removes them) that may overlap with others (as determined by the
overlap_penalty_hamiltonian). Apply translation and rotation trial moves
to separate overlapping pairs and make free space available for more body
insertions:
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_self_assembly_interactive::main;Use overlap_penalty_hamiltonian for trial moves during the initialize step.
The harmonic part of OverlapPenalty allows overlaps to be removed over many
simulation steps. Pass a fixed temperature=1.0 because the energy scale in
OverlapPenalty has no relation to that in hamiltonian.
Warning
If you use
hamiltonianhere, then a trial move would need to remove an overlap in one step. That might not be possible depending on the amount of overlap, the trial move size, and the density of the system.
Transition to the Equilibrate State
After many steps, QuickInsert should add all the requested bodies,
QuickCompress should achieve the target boundary volume, and all
overlaps will be removed (the total energy of overlap_penalty_hamiltonian is 0).
When all those are true quick_compress.is_complete() will return true. Before
proceeding to the equilibrate phase, call tune_with_options to adjust the
trial move sizes and achieve a 20% move acceptance rate:
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_self_assembly_interactive::main;Detect Failures
It might happen that after a long time, QuickInsert fails to add the target
number of bodies or QuickCompress is unable to achieve the target volume.
Instead of running the simulation for an infinitely long time, detect this
condition and report an error:
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_self_assembly_interactive::main;You can use the anyhow crate to report errors with context. Set
initial_packing_fraction to 1.0, run the example, and you should get an
error similar to:
Error: failed to initialize
Caused by:
inserted 391/512 bodies and compressed to 80.4247719318987 / 114.89253133128388
Equilibrate
The equilibration phase of the simulation applies the translate and rotate
trial moves with the hard overlap Hamiltonian (hamiltonian):
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_self_assembly_interactive::main;Equilibration never ends in this tutorial. In your own simulations, you might transition to a production phase after a certain number of steps and eventually complete the simulation.
Execute the Simulation in Batch Mode
Implement main()
To run the simulation, construct the HardEllipseSelfAssembly 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::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_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-ellipse-self-assembly
Visualize the Simulation Results
Open the generated hard-ellipse-self-assembly.gsd in Ovito or another
visualization tool to see the simulation results. To render the ellipses in
Ovito, import this modifier snippet (use the copy button to copy the entire
text to your clipboard):
{"description": "OVITO Modifier Snippet: Compute property (Aspherical Shape)", "payload": "AAAC6njahVLNSsNAEP4a21p/aP0BFS8WEfQgpUW8K+rBg1q0FO0ttts00GbDZgP6BD6BV/GuB1/Ciw/iUwg6m0yaiIoD3+zm29mZLzNbeX59mv546wDhHoASYQFnaOMYLVqraEJBEjtBqCC2WeyghgZ2ydeZy8UwC174e5uQJ2wQqthHAB8DCMroogsbQ2IviLGJFxRjEZa4loXLyF9FvsPsOq9Fqluj+smu8W2XZ1hIrUCYM7c/M5Ye54zeO9ZdlqH2Q91U0hdK3yLzi+WuHPnSE54+tUciyB7NiBtfiSBwpTfmjdZSTwyFY2vzh4+5jLRFT7jO4Fqqox8XLZZkmfg1wvIB1Q21SDSdyJ7bd4WKa0RhZhJbTVtptzsUwR/xh6mWAt9bNTNn3vWcTOYiR6yY3v2SYZLPN43I/wvmOdyMoZQpk6SZJ0ydi37LVo7QaQPplYzfV2Rt7lAy33ceQpHnbPb3ceKomQ/GfQEXHnuD"}
Alternately, you can:
- Add a Compute Property modification to the pipeline.
- Set Property Name to Aspherical Shape.
- Set X=0.5, Y=0.1, Z=0.1. 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:
Conclusion
This tutorial showed you how to perform hard-ellipse self-assembly simulations
using a shape overlap potential, periodic boundary conditions, and QuickInsert
to add bodies.
Navigate to the top of the page to see the simulation in action. Notice that ellipses are first added in a large batch. Once all the overlaps are removed, another batch appears. After all all ellipses 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 domains form where all the ellipses point in roughly the same direction while at the same time there is no translational order. This is the nematic phase.
Complete Code
use anyhow::{Context, anyhow};
use hoomd_geometry::{
Volume,
shape::{Ellipse, Rectangle},
};
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, Angle, Cartesian};
type PositionVector = Cartesian<2>;
type Orientation = Angle;
type BodyProperties = OrientedPoint<PositionVector, Orientation>;
type SiteProperties = OrientedPoint<PositionVector, Orientation>;
#[cfg_attr(feature = "bevy", derive(Resource))]
struct HardEllipseSelfAssembly {
/// Positions and orientations of all the bodies in the simulation.
microstate: Microstate<
BodyProperties,
SiteProperties,
VecCell<SiteKey, 2>,
Periodic<Rectangle>,
>,
/// How sites interact with other sites and fields.
hamiltonian: PairwiseCutoff<HardShape<Ellipse>>,
/// 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<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<
Anisotropic<ApproximateShapeOverlap<OverlapPenalty, Ellipse>>,
>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl HardEllipseSelfAssembly {
/// Construct a new hard ellipse self-assembly simulation.
fn new() -> anyhow::Result<HardEllipseSelfAssembly> {
let initial_packing_fraction = 0.4;
let target_packing_fraction = 0.7;
let n_bodies = 512;
let maximum_distance = 0.07;
let maximum_rotation = 0.3;
let sigma = 1.0;
let aspect = 5.0;
let macrostate = Isothermal { temperature: 1.0 };
assert!(aspect >= 1.0);
let ellipse = Ellipse::with_semi_axes([
(sigma / 2.0).try_into()?,
(sigma / aspect / 2.0).try_into()?,
]);
let hamiltonian = PairwiseCutoff(HardShape(ellipse.clone()));
let initial_box_volume =
n_bodies as f64 * ellipse.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.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 * ellipse.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
let approximate_shape_overlap = Anisotropic {
interaction: ApproximateShapeOverlap::new(
ellipse,
OverlapPenalty::default(),
0.01.try_into()?,
),
r_cut: sigma,
};
let overlap_penalty_hamiltonian =
PairwiseCutoff(approximate_shape_overlap);
Ok(HardEllipseSelfAssembly {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
rotate_sweep,
quick_compress,
quick_insert,
macrostate,
phase: Phase::Initialize,
})
}
}
impl Simulation for HardEllipseSelfAssembly {
/// 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 HardEllipseSelfAssembly {
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 = HardEllipseSelfAssembly::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("hard-ellipse-self-assembly.gsd")?;
for _ in 0..100_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_ellipse_self_assembly_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use hard_ellipse_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.
