Hard Disk Self-Assembly
Click to focus. When focused, press m to open the options menu.
Refresh the page to restart the simulation.
Overview
When compressed to a sufficiently high packing fraction, systems of hard particles self-assemble into ordered structures. Use periodic boundary conditions to model the behavior of the bulk material.
- Objectives:
- Explain how to execute simulations with periodic boundary conditions.
- Show how to quickly compress the microstate to a target packing fraction.
- Demonstrate the self-assembly of hard disks into the hexagonal phase.
- File:
hoomd-rs/examples/mc-tutorial/hard-disk-self-assembly.rs - Run (interactively):
cargo run --release --features "bevy" --example hard-disk-self-assembly - Run (in batch mode):
cargo run --release --example hard-disk-self-assembly
Use Declarations
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;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::{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;The sites are in this tutorial are represented by isotropic disks.
Therefore, use Point 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::{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;The simulation will consist of two phases:
- Compress: Decrease the volume of the microstate until it reaches a target.
- Equilibrate: Perform hard particle Monte Carlo to self-assembly the hexagonal phase.
The phase field tracks the current phase of the simulation. It stores an enum:
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;Construct the Simulation Model
The new() method constructs a new simulation model:
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;Parameters
Assign all the model parameters in one code block:
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;initial_packing_fraction is the volume of the disks divided by the volume of
the simulation boundary in the initial state. Choose this value so that the
placed disks do not overlap. During the Compress phase, the microstate will
be compressed until it reaches the packing fraction target_packing_fraction.
n_disks is the number of disks to place in the microstate, maximum_distance
is the largest distance a translation trial move can take, sigma is the
diameter of each disk, and macrostate holds the temperature set point (in
units of energy).
Hamiltonian
HardSphere represents each site with a hard sphere of the given diameter. The hard
sphere site pair energy is infinite when the two sites overlap and 0 when
they do not. Use PairwiseCutoff with the HardSphere interaction as the Hamiltonian:
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;Periodic Boundary Conditions
Use periodic boundary conditions via the Periodic type to represent an
infinitely repeating system. Provide the underlying shape and the maximum
interaction range between sites to construct Periodic:
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;Periodic uses this distance to generate ghost sites outside the boundary that
are periodic images of sites inside. Methods like PairwiseCutoff will compute
interactions between sites inside the boundary and all other sites (whether
they are ghosts or not). All pairs separated by a distance larger than the maximum
interaction range are assumed to be non-overlapping. You must choose this value
appropriately for your shape(s). For the case of hard disks, the largest distance
between the centers of two potentially overlapping shapes is sigma. HardSphere
provides this via the maximum_interaction_range() method.
Important
In hoomd-rs, it is YOUR responsibility to determine the appropriate
maximum_interaction_range. You might be used to other simulation codes, HOOMD-blue for example, that automatically determine this maximum for you. That is not possible in hoomd-rs as your model’s interactions and/or any analysis methods could be any arbitrary code.
Warning
If you set
maximum_interaction_rangetoo small,PairwiseCutoff(and similar methods) will miss interactions that should be computed.
Microstate
Construct a microstate with the periodic boundary conditions and the VecCell
spatial data structure:
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;Place Disks
Place n_disks disks in the microstate on a square lattice:
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;Trial Moves
Apply both Translate trial moves to the bodies:
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;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
QuickCompress will introduce strain in the system each time it compresses
the microstate. Translation trial moves can relieve that strain over many steps.
It is possible to use QuickCompress with pure hard shape interactions, though
you must choose the translation trial move distance very carefully.
The OverlapPenalty potential works around this problem. It consists of an infinite energy core
followed by a harmonic potential added to a step function. The infinite core
prevents bodies from overlapping too much on compression, 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 Isotropic, Expanded interaction:
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;OverlapPenalty applies the potential described above (centered on )
which Expanded shifts to the surface of the sphere at .
Initialize the Struct
Package all these values into a struct to represent the simulation:
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;Begin the simulation in the Compress phase.
Implement Simulation
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;Advance the Simulation
advance calls self.apply() to advance the simulation when in the
compress phase and self.equilibrate() when in the equilibrate phase:
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;Each method advances the simulation one step and potentially changes the
phase. The compress 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::{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;Implement HardDiskSelfAssembly
Place the model-specific methods in the inherent implementation for
HardDiskSelfAssembly:
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;Compress
Implement the compress phase:
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;Compress the Microstate
The quick_compress.apply method irreversibly compresses 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;To avoid jamming the system, QuickCompress waits until the total energy of
the given Hamiltonian (overlap_penalty_hamiltonian in this case) is 0 before
compressing.
Separate Overlapping Bodies
QuickCompress is quick because it only scales the system toward the
target (never away from) and may cause some pairs of sites to overlap. Apply
translation and rotation trial moves to separate overlapping pairs and make free
space available for further compression:
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;Use overlap_penalty_hamiltonian for trial moves during the compress 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 packing fraction of the system.
Transition to the Equilibrate State
After many steps, QuickCompress should achieve the target boundary volume and all
overlaps will be removed (the total energy of overlap_penalty_hamiltonian is 0).
When both those are true quick_compress.is_complete() will return true and the
simulation can proceed to the equilibrate phase:
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;Detect Failures
It might happen QuickCompress fails to reach the target boundary volume, even after
many steps. 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::{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;You can use the anyhow crate to report errors with context. Set
target_packing_fraction to 1.0, run the example, and you should get an error
similar to:
Error: failed at step: 10000
Caused by:
0: failed to compress
1: Achieved volume 3971.206835460344 after 10000 steps. The target was 3216.990877275948.
Equilibrate
The equilibration phase of the simulation applies the translate
trial moves with the hard overlap Hamiltonian (hamiltonian):
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;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 HardDiskSelfAssembly 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::{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;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.
Run the Simulation
In a terminal, execute the following command to run the simulation in batch mode:
cargo run --release --example hard-particle-self-assembly
Visualize the Simulation Results
Open the generated tetronimoes.gsd in Ovito or another visualization
tool to see the simulation results. Ovito will render the disks as spheres
with the expected diameter of 1 by default.
Render with Tachyon, and you should see something like:
Conclusion
This tutorial showed you how to perform hard disk self-assembly simulations
using periodic boundary conditions and QuickCompress to achieve a target
packing fraction.
Navigate to the top of the page to see the simulation in action. Notice that the disks start on an evenly spaced grid and are quickly scaled to a higher packing fraction. The many grain boundaries are caused by the quick compression. Over time, local trial moves will heal these grain boundaries leaving a single crystal.
Complete Code
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;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.
