Polydisperse Hard Disk Model
Click to focus. When focused, press m to open the options menu.
Refresh the page to restart the simulation.
Overview
In some models, every site is different from every other site. For example, you can model polydisperse hard disks where each site has its own radius. This tutorial shows you how to define a custom site type and compute a radius-dependent pairwise interactions.
- Objectives:
- Add a
radiusfield to a custom site properties type. - Show how to compute radius-dependent pairwise interactions.
- Execute a sample simulation where every site has a randomly chosen radius.
- Add a
- File:
hoomd-rs/examples/mc-tutorial/polydisperse-hard-disk-model.rs - Run (interactively):
cargo run --release --features "bevy" --example polydisperse-hard-disk-model - Run (in batch mode):
cargo run --release --example polydisperse-hard-disk-model
Use Declarations
use anyhow::{Context, anyhow};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use rand::{
SeedableRng,
distr::{Distribution, Uniform},
rngs::StdRng,
};
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
univariate::{Expanded, OverlapPenalty, UnivariateEnergy},
};
use hoomd_mc::{
BodyDistribution, QuickCompress, QuickInsert, Sweep, Translate, Trial,
Tune, TuneOptions,
};
use hoomd_microstate::{
AppendMicrostate, Body, Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_utility::valid::PositiveReal;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's radius.
radius: PositiveReal,
}
impl Transform<SiteProperties> for BodyProperties {
fn transform(&self, site_properties: &SiteProperties) -> SiteProperties {
SiteProperties {
position: self.position + site_properties.position,
..*site_properties
}
}
}
#[derive(MaximumInteractionRange)]
struct SitePairInteraction {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
if r < a.radius.get() + b.radius.get() {
f64::INFINITY
} else {
0.0
}
}
fn site_pair_energy_initial(
&self,
_a: &SiteProperties,
_b: &SiteProperties,
) -> f64 {
0.0
}
fn is_only_infinite_or_zero() -> bool {
true
}
}
#[derive(MaximumInteractionRange)]
struct SitePairOverlapPenalty {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairOverlapPenalty {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
let pair_interaction = Expanded {
delta: a.radius.get() + b.radius.get(),
f: OverlapPenalty::default(),
};
pair_interaction.energy(r)
}
}
struct PolydisperseBodyDistribution {
/// Radius of each disk to insert into the microstate.
radii: Vec<PositiveReal>,
/// Simulation boundary.
boundary: Periodic<Rectangle>,
}
impl BodyDistribution<Body<BodyProperties, SiteProperties>>
for PolydisperseBodyDistribution
{
fn sample<R: rand::Rng + ?Sized>(
&self,
index: usize,
rng: &mut R,
) -> Body<BodyProperties, SiteProperties> {
let properties = Point {
position: self.boundary.sample(rng),
};
let sites = vec![SiteProperties {
position: Cartesian::default(),
radius: self.radii[index],
}];
Body { properties, sites }
}
}
impl<X> AppendMicrostate<BodyProperties, SiteProperties, X, Periodic<Rectangle>>
for HoomdGsdFile
{
#[inline]
fn append_microstate(
&mut self,
microstate: &Microstate<
BodyProperties,
SiteProperties,
X,
Periodic<Rectangle>,
>,
) -> Result<hoomd_gsd::hoomd::Frame<'_>, hoomd_gsd::hoomd::AppendError>
{
self.append_frame(microstate.step())?
.configuration_box(microstate.boundary().shape().to_gsd_box())?
.configuration_dimensions(Dimensions::Two)?
.particles_position(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.position)
.map(|p| [p[0], p[1], 0.0].into()),
)?
.particles_diameter(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.radius.get() * 2.0),
)
}
}
impl PolydisperseHardDiskModel {
/// Construct a new hard disk self-assembly simulation.
fn new() -> anyhow::Result<PolydisperseHardDiskModel> {
let seed = 1;
let minimum_radius = 0.1;
let maximum_radius = 0.8;
let initial_packing_fraction = 0.6;
let target_packing_fraction = 0.72;
let n_disks = 64_usize.pow(2);
let maximum_distance = 0.07;
let macrostate = Isothermal { temperature: 1.0 };
let radius_distribution = Uniform::new(minimum_radius, maximum_radius)?;
let mut rng = StdRng::seed_from_u64(seed);
let mut radii = Vec::with_capacity(n_disks);
for r in radius_distribution.sample_iter(&mut rng).take(n_disks) {
radii.push(r.try_into()?);
}
let total_particle_area = radii.iter().fold(0.0, |total, r| {
let circle = Circle { radius: *r };
total + circle.volume()
});
let hamiltonian = PairwiseCutoff(SitePairInteraction {
maximum_interaction_range: maximum_radius * 2.0,
});
let overlap_penalty_hamiltonian =
PairwiseCutoff(SitePairOverlapPenalty {
maximum_interaction_range: maximum_radius * 2.0,
});
let initial_box_volume = total_particle_area / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.seed(seed as u32)
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let distribution = PolydisperseBodyDistribution {
boundary: microstate.boundary().clone(),
radii,
};
let quick_insert = QuickInsert::new(distribution, n_disks);
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume = total_particle_area / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(PolydisperseHardDiskModel {
microstate,
hamiltonian,
translate_sweep,
quick_insert,
quick_compress,
overlap_penalty_hamiltonian,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct PolydisperseHardDiskModel {
/// 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<SitePairInteraction>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<PolydisperseBodyDistribution>,
/// The current phase of the simulation.
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<SitePairOverlapPenalty>,
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for PolydisperseHardDiskModel {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl PolydisperseHardDiskModel {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = PolydisperseHardDiskModel::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("polydisperse-hard-disk-model.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 polydisperse_hard_disk_model_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use polydisperse_hard_disk_model_interactive::main;Type Aliases
Create type aliases for your model’s vector and body properties so that you don’t need to repeat the full generic type names throughout the code:
use anyhow::{Context, anyhow};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use rand::{
SeedableRng,
distr::{Distribution, Uniform},
rngs::StdRng,
};
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
univariate::{Expanded, OverlapPenalty, UnivariateEnergy},
};
use hoomd_mc::{
BodyDistribution, QuickCompress, QuickInsert, Sweep, Translate, Trial,
Tune, TuneOptions,
};
use hoomd_microstate::{
AppendMicrostate, Body, Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_utility::valid::PositiveReal;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's radius.
radius: PositiveReal,
}
impl Transform<SiteProperties> for BodyProperties {
fn transform(&self, site_properties: &SiteProperties) -> SiteProperties {
SiteProperties {
position: self.position + site_properties.position,
..*site_properties
}
}
}
#[derive(MaximumInteractionRange)]
struct SitePairInteraction {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
if r < a.radius.get() + b.radius.get() {
f64::INFINITY
} else {
0.0
}
}
fn site_pair_energy_initial(
&self,
_a: &SiteProperties,
_b: &SiteProperties,
) -> f64 {
0.0
}
fn is_only_infinite_or_zero() -> bool {
true
}
}
#[derive(MaximumInteractionRange)]
struct SitePairOverlapPenalty {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairOverlapPenalty {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
let pair_interaction = Expanded {
delta: a.radius.get() + b.radius.get(),
f: OverlapPenalty::default(),
};
pair_interaction.energy(r)
}
}
struct PolydisperseBodyDistribution {
/// Radius of each disk to insert into the microstate.
radii: Vec<PositiveReal>,
/// Simulation boundary.
boundary: Periodic<Rectangle>,
}
impl BodyDistribution<Body<BodyProperties, SiteProperties>>
for PolydisperseBodyDistribution
{
fn sample<R: rand::Rng + ?Sized>(
&self,
index: usize,
rng: &mut R,
) -> Body<BodyProperties, SiteProperties> {
let properties = Point {
position: self.boundary.sample(rng),
};
let sites = vec![SiteProperties {
position: Cartesian::default(),
radius: self.radii[index],
}];
Body { properties, sites }
}
}
impl<X> AppendMicrostate<BodyProperties, SiteProperties, X, Periodic<Rectangle>>
for HoomdGsdFile
{
#[inline]
fn append_microstate(
&mut self,
microstate: &Microstate<
BodyProperties,
SiteProperties,
X,
Periodic<Rectangle>,
>,
) -> Result<hoomd_gsd::hoomd::Frame<'_>, hoomd_gsd::hoomd::AppendError>
{
self.append_frame(microstate.step())?
.configuration_box(microstate.boundary().shape().to_gsd_box())?
.configuration_dimensions(Dimensions::Two)?
.particles_position(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.position)
.map(|p| [p[0], p[1], 0.0].into()),
)?
.particles_diameter(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.radius.get() * 2.0),
)
}
}
impl PolydisperseHardDiskModel {
/// Construct a new hard disk self-assembly simulation.
fn new() -> anyhow::Result<PolydisperseHardDiskModel> {
let seed = 1;
let minimum_radius = 0.1;
let maximum_radius = 0.8;
let initial_packing_fraction = 0.6;
let target_packing_fraction = 0.72;
let n_disks = 64_usize.pow(2);
let maximum_distance = 0.07;
let macrostate = Isothermal { temperature: 1.0 };
let radius_distribution = Uniform::new(minimum_radius, maximum_radius)?;
let mut rng = StdRng::seed_from_u64(seed);
let mut radii = Vec::with_capacity(n_disks);
for r in radius_distribution.sample_iter(&mut rng).take(n_disks) {
radii.push(r.try_into()?);
}
let total_particle_area = radii.iter().fold(0.0, |total, r| {
let circle = Circle { radius: *r };
total + circle.volume()
});
let hamiltonian = PairwiseCutoff(SitePairInteraction {
maximum_interaction_range: maximum_radius * 2.0,
});
let overlap_penalty_hamiltonian =
PairwiseCutoff(SitePairOverlapPenalty {
maximum_interaction_range: maximum_radius * 2.0,
});
let initial_box_volume = total_particle_area / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.seed(seed as u32)
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let distribution = PolydisperseBodyDistribution {
boundary: microstate.boundary().clone(),
radii,
};
let quick_insert = QuickInsert::new(distribution, n_disks);
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume = total_particle_area / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(PolydisperseHardDiskModel {
microstate,
hamiltonian,
translate_sweep,
quick_insert,
quick_compress,
overlap_penalty_hamiltonian,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct PolydisperseHardDiskModel {
/// 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<SitePairInteraction>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<PolydisperseBodyDistribution>,
/// The current phase of the simulation.
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<SitePairOverlapPenalty>,
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for PolydisperseHardDiskModel {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl PolydisperseHardDiskModel {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = PolydisperseHardDiskModel::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("polydisperse-hard-disk-model.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 polydisperse_hard_disk_model_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use polydisperse_hard_disk_model_interactive::main;The sites are in this tutorial are placed at points in space and interact via
an isotropic interaction. Therefore, use Point for the body properties.
Site Properties
Define SiteProperties
No built-in hoomd-rs site type has a radius field. Define a custom type to hold it along with the site’s position:
use anyhow::{Context, anyhow};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use rand::{
SeedableRng,
distr::{Distribution, Uniform},
rngs::StdRng,
};
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
univariate::{Expanded, OverlapPenalty, UnivariateEnergy},
};
use hoomd_mc::{
BodyDistribution, QuickCompress, QuickInsert, Sweep, Translate, Trial,
Tune, TuneOptions,
};
use hoomd_microstate::{
AppendMicrostate, Body, Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_utility::valid::PositiveReal;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's radius.
radius: PositiveReal,
}
impl Transform<SiteProperties> for BodyProperties {
fn transform(&self, site_properties: &SiteProperties) -> SiteProperties {
SiteProperties {
position: self.position + site_properties.position,
..*site_properties
}
}
}
#[derive(MaximumInteractionRange)]
struct SitePairInteraction {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
if r < a.radius.get() + b.radius.get() {
f64::INFINITY
} else {
0.0
}
}
fn site_pair_energy_initial(
&self,
_a: &SiteProperties,
_b: &SiteProperties,
) -> f64 {
0.0
}
fn is_only_infinite_or_zero() -> bool {
true
}
}
#[derive(MaximumInteractionRange)]
struct SitePairOverlapPenalty {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairOverlapPenalty {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
let pair_interaction = Expanded {
delta: a.radius.get() + b.radius.get(),
f: OverlapPenalty::default(),
};
pair_interaction.energy(r)
}
}
struct PolydisperseBodyDistribution {
/// Radius of each disk to insert into the microstate.
radii: Vec<PositiveReal>,
/// Simulation boundary.
boundary: Periodic<Rectangle>,
}
impl BodyDistribution<Body<BodyProperties, SiteProperties>>
for PolydisperseBodyDistribution
{
fn sample<R: rand::Rng + ?Sized>(
&self,
index: usize,
rng: &mut R,
) -> Body<BodyProperties, SiteProperties> {
let properties = Point {
position: self.boundary.sample(rng),
};
let sites = vec![SiteProperties {
position: Cartesian::default(),
radius: self.radii[index],
}];
Body { properties, sites }
}
}
impl<X> AppendMicrostate<BodyProperties, SiteProperties, X, Periodic<Rectangle>>
for HoomdGsdFile
{
#[inline]
fn append_microstate(
&mut self,
microstate: &Microstate<
BodyProperties,
SiteProperties,
X,
Periodic<Rectangle>,
>,
) -> Result<hoomd_gsd::hoomd::Frame<'_>, hoomd_gsd::hoomd::AppendError>
{
self.append_frame(microstate.step())?
.configuration_box(microstate.boundary().shape().to_gsd_box())?
.configuration_dimensions(Dimensions::Two)?
.particles_position(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.position)
.map(|p| [p[0], p[1], 0.0].into()),
)?
.particles_diameter(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.radius.get() * 2.0),
)
}
}
impl PolydisperseHardDiskModel {
/// Construct a new hard disk self-assembly simulation.
fn new() -> anyhow::Result<PolydisperseHardDiskModel> {
let seed = 1;
let minimum_radius = 0.1;
let maximum_radius = 0.8;
let initial_packing_fraction = 0.6;
let target_packing_fraction = 0.72;
let n_disks = 64_usize.pow(2);
let maximum_distance = 0.07;
let macrostate = Isothermal { temperature: 1.0 };
let radius_distribution = Uniform::new(minimum_radius, maximum_radius)?;
let mut rng = StdRng::seed_from_u64(seed);
let mut radii = Vec::with_capacity(n_disks);
for r in radius_distribution.sample_iter(&mut rng).take(n_disks) {
radii.push(r.try_into()?);
}
let total_particle_area = radii.iter().fold(0.0, |total, r| {
let circle = Circle { radius: *r };
total + circle.volume()
});
let hamiltonian = PairwiseCutoff(SitePairInteraction {
maximum_interaction_range: maximum_radius * 2.0,
});
let overlap_penalty_hamiltonian =
PairwiseCutoff(SitePairOverlapPenalty {
maximum_interaction_range: maximum_radius * 2.0,
});
let initial_box_volume = total_particle_area / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.seed(seed as u32)
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let distribution = PolydisperseBodyDistribution {
boundary: microstate.boundary().clone(),
radii,
};
let quick_insert = QuickInsert::new(distribution, n_disks);
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume = total_particle_area / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(PolydisperseHardDiskModel {
microstate,
hamiltonian,
translate_sweep,
quick_insert,
quick_compress,
overlap_penalty_hamiltonian,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct PolydisperseHardDiskModel {
/// 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<SitePairInteraction>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<PolydisperseBodyDistribution>,
/// The current phase of the simulation.
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<SitePairOverlapPenalty>,
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for PolydisperseHardDiskModel {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl PolydisperseHardDiskModel {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = PolydisperseHardDiskModel::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("polydisperse-hard-disk-model.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 polydisperse_hard_disk_model_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use polydisperse_hard_disk_model_interactive::main;As in Type-dependent Interactions, #[derive] the traits that are needed.
Transforming Sites
Implement Transform<SiteProperties> for the body properties type:
use anyhow::{Context, anyhow};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use rand::{
SeedableRng,
distr::{Distribution, Uniform},
rngs::StdRng,
};
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
univariate::{Expanded, OverlapPenalty, UnivariateEnergy},
};
use hoomd_mc::{
BodyDistribution, QuickCompress, QuickInsert, Sweep, Translate, Trial,
Tune, TuneOptions,
};
use hoomd_microstate::{
AppendMicrostate, Body, Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_utility::valid::PositiveReal;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's radius.
radius: PositiveReal,
}
impl Transform<SiteProperties> for BodyProperties {
fn transform(&self, site_properties: &SiteProperties) -> SiteProperties {
SiteProperties {
position: self.position + site_properties.position,
..*site_properties
}
}
}
#[derive(MaximumInteractionRange)]
struct SitePairInteraction {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
if r < a.radius.get() + b.radius.get() {
f64::INFINITY
} else {
0.0
}
}
fn site_pair_energy_initial(
&self,
_a: &SiteProperties,
_b: &SiteProperties,
) -> f64 {
0.0
}
fn is_only_infinite_or_zero() -> bool {
true
}
}
#[derive(MaximumInteractionRange)]
struct SitePairOverlapPenalty {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairOverlapPenalty {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
let pair_interaction = Expanded {
delta: a.radius.get() + b.radius.get(),
f: OverlapPenalty::default(),
};
pair_interaction.energy(r)
}
}
struct PolydisperseBodyDistribution {
/// Radius of each disk to insert into the microstate.
radii: Vec<PositiveReal>,
/// Simulation boundary.
boundary: Periodic<Rectangle>,
}
impl BodyDistribution<Body<BodyProperties, SiteProperties>>
for PolydisperseBodyDistribution
{
fn sample<R: rand::Rng + ?Sized>(
&self,
index: usize,
rng: &mut R,
) -> Body<BodyProperties, SiteProperties> {
let properties = Point {
position: self.boundary.sample(rng),
};
let sites = vec![SiteProperties {
position: Cartesian::default(),
radius: self.radii[index],
}];
Body { properties, sites }
}
}
impl<X> AppendMicrostate<BodyProperties, SiteProperties, X, Periodic<Rectangle>>
for HoomdGsdFile
{
#[inline]
fn append_microstate(
&mut self,
microstate: &Microstate<
BodyProperties,
SiteProperties,
X,
Periodic<Rectangle>,
>,
) -> Result<hoomd_gsd::hoomd::Frame<'_>, hoomd_gsd::hoomd::AppendError>
{
self.append_frame(microstate.step())?
.configuration_box(microstate.boundary().shape().to_gsd_box())?
.configuration_dimensions(Dimensions::Two)?
.particles_position(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.position)
.map(|p| [p[0], p[1], 0.0].into()),
)?
.particles_diameter(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.radius.get() * 2.0),
)
}
}
impl PolydisperseHardDiskModel {
/// Construct a new hard disk self-assembly simulation.
fn new() -> anyhow::Result<PolydisperseHardDiskModel> {
let seed = 1;
let minimum_radius = 0.1;
let maximum_radius = 0.8;
let initial_packing_fraction = 0.6;
let target_packing_fraction = 0.72;
let n_disks = 64_usize.pow(2);
let maximum_distance = 0.07;
let macrostate = Isothermal { temperature: 1.0 };
let radius_distribution = Uniform::new(minimum_radius, maximum_radius)?;
let mut rng = StdRng::seed_from_u64(seed);
let mut radii = Vec::with_capacity(n_disks);
for r in radius_distribution.sample_iter(&mut rng).take(n_disks) {
radii.push(r.try_into()?);
}
let total_particle_area = radii.iter().fold(0.0, |total, r| {
let circle = Circle { radius: *r };
total + circle.volume()
});
let hamiltonian = PairwiseCutoff(SitePairInteraction {
maximum_interaction_range: maximum_radius * 2.0,
});
let overlap_penalty_hamiltonian =
PairwiseCutoff(SitePairOverlapPenalty {
maximum_interaction_range: maximum_radius * 2.0,
});
let initial_box_volume = total_particle_area / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.seed(seed as u32)
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let distribution = PolydisperseBodyDistribution {
boundary: microstate.boundary().clone(),
radii,
};
let quick_insert = QuickInsert::new(distribution, n_disks);
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume = total_particle_area / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(PolydisperseHardDiskModel {
microstate,
hamiltonian,
translate_sweep,
quick_insert,
quick_compress,
overlap_penalty_hamiltonian,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct PolydisperseHardDiskModel {
/// 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<SitePairInteraction>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<PolydisperseBodyDistribution>,
/// The current phase of the simulation.
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<SitePairOverlapPenalty>,
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for PolydisperseHardDiskModel {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl PolydisperseHardDiskModel {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = PolydisperseHardDiskModel::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("polydisperse-hard-disk-model.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 polydisperse_hard_disk_model_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use polydisperse_hard_disk_model_interactive::main;Site-site Interactions
The site_pair_energy method is a function only of two sites. It does not
reference any particular microstate. At the same time, the interaction type
must be aware of the maximum radius it can expect so that it can be used
with PairwiseCutoff. You can meet these requirements by storing the maximum
interaction range in a field of the site pair interaction type.
Define the Hard Shape SitePairInteraction
Define the SitePairInteraction type:
use anyhow::{Context, anyhow};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use rand::{
SeedableRng,
distr::{Distribution, Uniform},
rngs::StdRng,
};
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
univariate::{Expanded, OverlapPenalty, UnivariateEnergy},
};
use hoomd_mc::{
BodyDistribution, QuickCompress, QuickInsert, Sweep, Translate, Trial,
Tune, TuneOptions,
};
use hoomd_microstate::{
AppendMicrostate, Body, Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_utility::valid::PositiveReal;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's radius.
radius: PositiveReal,
}
impl Transform<SiteProperties> for BodyProperties {
fn transform(&self, site_properties: &SiteProperties) -> SiteProperties {
SiteProperties {
position: self.position + site_properties.position,
..*site_properties
}
}
}
#[derive(MaximumInteractionRange)]
struct SitePairInteraction {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
if r < a.radius.get() + b.radius.get() {
f64::INFINITY
} else {
0.0
}
}
fn site_pair_energy_initial(
&self,
_a: &SiteProperties,
_b: &SiteProperties,
) -> f64 {
0.0
}
fn is_only_infinite_or_zero() -> bool {
true
}
}
#[derive(MaximumInteractionRange)]
struct SitePairOverlapPenalty {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairOverlapPenalty {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
let pair_interaction = Expanded {
delta: a.radius.get() + b.radius.get(),
f: OverlapPenalty::default(),
};
pair_interaction.energy(r)
}
}
struct PolydisperseBodyDistribution {
/// Radius of each disk to insert into the microstate.
radii: Vec<PositiveReal>,
/// Simulation boundary.
boundary: Periodic<Rectangle>,
}
impl BodyDistribution<Body<BodyProperties, SiteProperties>>
for PolydisperseBodyDistribution
{
fn sample<R: rand::Rng + ?Sized>(
&self,
index: usize,
rng: &mut R,
) -> Body<BodyProperties, SiteProperties> {
let properties = Point {
position: self.boundary.sample(rng),
};
let sites = vec![SiteProperties {
position: Cartesian::default(),
radius: self.radii[index],
}];
Body { properties, sites }
}
}
impl<X> AppendMicrostate<BodyProperties, SiteProperties, X, Periodic<Rectangle>>
for HoomdGsdFile
{
#[inline]
fn append_microstate(
&mut self,
microstate: &Microstate<
BodyProperties,
SiteProperties,
X,
Periodic<Rectangle>,
>,
) -> Result<hoomd_gsd::hoomd::Frame<'_>, hoomd_gsd::hoomd::AppendError>
{
self.append_frame(microstate.step())?
.configuration_box(microstate.boundary().shape().to_gsd_box())?
.configuration_dimensions(Dimensions::Two)?
.particles_position(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.position)
.map(|p| [p[0], p[1], 0.0].into()),
)?
.particles_diameter(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.radius.get() * 2.0),
)
}
}
impl PolydisperseHardDiskModel {
/// Construct a new hard disk self-assembly simulation.
fn new() -> anyhow::Result<PolydisperseHardDiskModel> {
let seed = 1;
let minimum_radius = 0.1;
let maximum_radius = 0.8;
let initial_packing_fraction = 0.6;
let target_packing_fraction = 0.72;
let n_disks = 64_usize.pow(2);
let maximum_distance = 0.07;
let macrostate = Isothermal { temperature: 1.0 };
let radius_distribution = Uniform::new(minimum_radius, maximum_radius)?;
let mut rng = StdRng::seed_from_u64(seed);
let mut radii = Vec::with_capacity(n_disks);
for r in radius_distribution.sample_iter(&mut rng).take(n_disks) {
radii.push(r.try_into()?);
}
let total_particle_area = radii.iter().fold(0.0, |total, r| {
let circle = Circle { radius: *r };
total + circle.volume()
});
let hamiltonian = PairwiseCutoff(SitePairInteraction {
maximum_interaction_range: maximum_radius * 2.0,
});
let overlap_penalty_hamiltonian =
PairwiseCutoff(SitePairOverlapPenalty {
maximum_interaction_range: maximum_radius * 2.0,
});
let initial_box_volume = total_particle_area / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.seed(seed as u32)
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let distribution = PolydisperseBodyDistribution {
boundary: microstate.boundary().clone(),
radii,
};
let quick_insert = QuickInsert::new(distribution, n_disks);
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume = total_particle_area / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(PolydisperseHardDiskModel {
microstate,
hamiltonian,
translate_sweep,
quick_insert,
quick_compress,
overlap_penalty_hamiltonian,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct PolydisperseHardDiskModel {
/// 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<SitePairInteraction>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<PolydisperseBodyDistribution>,
/// The current phase of the simulation.
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<SitePairOverlapPenalty>,
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for PolydisperseHardDiskModel {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl PolydisperseHardDiskModel {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = PolydisperseHardDiskModel::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("polydisperse-hard-disk-model.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 polydisperse_hard_disk_model_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use polydisperse_hard_disk_model_interactive::main;Every site-site interaction type must implement MaximumInteractionRange
which sets the distance above which the interactions go to zero. You can implement
this trait directly, or add the maximum_interaction_range field and
#derive[(MaximumInteractionRange)] as shown here.
Implement SitePairEnergy
PairwiseCutoff uses the SitePairEnergy trait to calculate the interaction
energy that each pair of sites contributes to the total. Implement the trait
and compute the radius-dependent interactions:
use anyhow::{Context, anyhow};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use rand::{
SeedableRng,
distr::{Distribution, Uniform},
rngs::StdRng,
};
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
univariate::{Expanded, OverlapPenalty, UnivariateEnergy},
};
use hoomd_mc::{
BodyDistribution, QuickCompress, QuickInsert, Sweep, Translate, Trial,
Tune, TuneOptions,
};
use hoomd_microstate::{
AppendMicrostate, Body, Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_utility::valid::PositiveReal;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's radius.
radius: PositiveReal,
}
impl Transform<SiteProperties> for BodyProperties {
fn transform(&self, site_properties: &SiteProperties) -> SiteProperties {
SiteProperties {
position: self.position + site_properties.position,
..*site_properties
}
}
}
#[derive(MaximumInteractionRange)]
struct SitePairInteraction {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
if r < a.radius.get() + b.radius.get() {
f64::INFINITY
} else {
0.0
}
}
fn site_pair_energy_initial(
&self,
_a: &SiteProperties,
_b: &SiteProperties,
) -> f64 {
0.0
}
fn is_only_infinite_or_zero() -> bool {
true
}
}
#[derive(MaximumInteractionRange)]
struct SitePairOverlapPenalty {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairOverlapPenalty {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
let pair_interaction = Expanded {
delta: a.radius.get() + b.radius.get(),
f: OverlapPenalty::default(),
};
pair_interaction.energy(r)
}
}
struct PolydisperseBodyDistribution {
/// Radius of each disk to insert into the microstate.
radii: Vec<PositiveReal>,
/// Simulation boundary.
boundary: Periodic<Rectangle>,
}
impl BodyDistribution<Body<BodyProperties, SiteProperties>>
for PolydisperseBodyDistribution
{
fn sample<R: rand::Rng + ?Sized>(
&self,
index: usize,
rng: &mut R,
) -> Body<BodyProperties, SiteProperties> {
let properties = Point {
position: self.boundary.sample(rng),
};
let sites = vec![SiteProperties {
position: Cartesian::default(),
radius: self.radii[index],
}];
Body { properties, sites }
}
}
impl<X> AppendMicrostate<BodyProperties, SiteProperties, X, Periodic<Rectangle>>
for HoomdGsdFile
{
#[inline]
fn append_microstate(
&mut self,
microstate: &Microstate<
BodyProperties,
SiteProperties,
X,
Periodic<Rectangle>,
>,
) -> Result<hoomd_gsd::hoomd::Frame<'_>, hoomd_gsd::hoomd::AppendError>
{
self.append_frame(microstate.step())?
.configuration_box(microstate.boundary().shape().to_gsd_box())?
.configuration_dimensions(Dimensions::Two)?
.particles_position(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.position)
.map(|p| [p[0], p[1], 0.0].into()),
)?
.particles_diameter(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.radius.get() * 2.0),
)
}
}
impl PolydisperseHardDiskModel {
/// Construct a new hard disk self-assembly simulation.
fn new() -> anyhow::Result<PolydisperseHardDiskModel> {
let seed = 1;
let minimum_radius = 0.1;
let maximum_radius = 0.8;
let initial_packing_fraction = 0.6;
let target_packing_fraction = 0.72;
let n_disks = 64_usize.pow(2);
let maximum_distance = 0.07;
let macrostate = Isothermal { temperature: 1.0 };
let radius_distribution = Uniform::new(minimum_radius, maximum_radius)?;
let mut rng = StdRng::seed_from_u64(seed);
let mut radii = Vec::with_capacity(n_disks);
for r in radius_distribution.sample_iter(&mut rng).take(n_disks) {
radii.push(r.try_into()?);
}
let total_particle_area = radii.iter().fold(0.0, |total, r| {
let circle = Circle { radius: *r };
total + circle.volume()
});
let hamiltonian = PairwiseCutoff(SitePairInteraction {
maximum_interaction_range: maximum_radius * 2.0,
});
let overlap_penalty_hamiltonian =
PairwiseCutoff(SitePairOverlapPenalty {
maximum_interaction_range: maximum_radius * 2.0,
});
let initial_box_volume = total_particle_area / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.seed(seed as u32)
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let distribution = PolydisperseBodyDistribution {
boundary: microstate.boundary().clone(),
radii,
};
let quick_insert = QuickInsert::new(distribution, n_disks);
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume = total_particle_area / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(PolydisperseHardDiskModel {
microstate,
hamiltonian,
translate_sweep,
quick_insert,
quick_compress,
overlap_penalty_hamiltonian,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct PolydisperseHardDiskModel {
/// 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<SitePairInteraction>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<PolydisperseBodyDistribution>,
/// The current phase of the simulation.
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<SitePairOverlapPenalty>,
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for PolydisperseHardDiskModel {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl PolydisperseHardDiskModel {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = PolydisperseHardDiskModel::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("polydisperse-hard-disk-model.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 polydisperse_hard_disk_model_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use polydisperse_hard_disk_model_interactive::main;First, compute the distance between the two sites. Return infinity when the disks overlap and 0 when they do not.
Special Methods for Hard Shapes
Hard potentials (potentials that always result of a value of 0 or infinity)
should implement two additional methods. site_pair_energy_initial should
return 0 and is_only_infinite_or_zero should return true:
use anyhow::{Context, anyhow};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use rand::{
SeedableRng,
distr::{Distribution, Uniform},
rngs::StdRng,
};
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
univariate::{Expanded, OverlapPenalty, UnivariateEnergy},
};
use hoomd_mc::{
BodyDistribution, QuickCompress, QuickInsert, Sweep, Translate, Trial,
Tune, TuneOptions,
};
use hoomd_microstate::{
AppendMicrostate, Body, Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_utility::valid::PositiveReal;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's radius.
radius: PositiveReal,
}
impl Transform<SiteProperties> for BodyProperties {
fn transform(&self, site_properties: &SiteProperties) -> SiteProperties {
SiteProperties {
position: self.position + site_properties.position,
..*site_properties
}
}
}
#[derive(MaximumInteractionRange)]
struct SitePairInteraction {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
if r < a.radius.get() + b.radius.get() {
f64::INFINITY
} else {
0.0
}
}
fn site_pair_energy_initial(
&self,
_a: &SiteProperties,
_b: &SiteProperties,
) -> f64 {
0.0
}
fn is_only_infinite_or_zero() -> bool {
true
}
}
#[derive(MaximumInteractionRange)]
struct SitePairOverlapPenalty {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairOverlapPenalty {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
let pair_interaction = Expanded {
delta: a.radius.get() + b.radius.get(),
f: OverlapPenalty::default(),
};
pair_interaction.energy(r)
}
}
struct PolydisperseBodyDistribution {
/// Radius of each disk to insert into the microstate.
radii: Vec<PositiveReal>,
/// Simulation boundary.
boundary: Periodic<Rectangle>,
}
impl BodyDistribution<Body<BodyProperties, SiteProperties>>
for PolydisperseBodyDistribution
{
fn sample<R: rand::Rng + ?Sized>(
&self,
index: usize,
rng: &mut R,
) -> Body<BodyProperties, SiteProperties> {
let properties = Point {
position: self.boundary.sample(rng),
};
let sites = vec![SiteProperties {
position: Cartesian::default(),
radius: self.radii[index],
}];
Body { properties, sites }
}
}
impl<X> AppendMicrostate<BodyProperties, SiteProperties, X, Periodic<Rectangle>>
for HoomdGsdFile
{
#[inline]
fn append_microstate(
&mut self,
microstate: &Microstate<
BodyProperties,
SiteProperties,
X,
Periodic<Rectangle>,
>,
) -> Result<hoomd_gsd::hoomd::Frame<'_>, hoomd_gsd::hoomd::AppendError>
{
self.append_frame(microstate.step())?
.configuration_box(microstate.boundary().shape().to_gsd_box())?
.configuration_dimensions(Dimensions::Two)?
.particles_position(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.position)
.map(|p| [p[0], p[1], 0.0].into()),
)?
.particles_diameter(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.radius.get() * 2.0),
)
}
}
impl PolydisperseHardDiskModel {
/// Construct a new hard disk self-assembly simulation.
fn new() -> anyhow::Result<PolydisperseHardDiskModel> {
let seed = 1;
let minimum_radius = 0.1;
let maximum_radius = 0.8;
let initial_packing_fraction = 0.6;
let target_packing_fraction = 0.72;
let n_disks = 64_usize.pow(2);
let maximum_distance = 0.07;
let macrostate = Isothermal { temperature: 1.0 };
let radius_distribution = Uniform::new(minimum_radius, maximum_radius)?;
let mut rng = StdRng::seed_from_u64(seed);
let mut radii = Vec::with_capacity(n_disks);
for r in radius_distribution.sample_iter(&mut rng).take(n_disks) {
radii.push(r.try_into()?);
}
let total_particle_area = radii.iter().fold(0.0, |total, r| {
let circle = Circle { radius: *r };
total + circle.volume()
});
let hamiltonian = PairwiseCutoff(SitePairInteraction {
maximum_interaction_range: maximum_radius * 2.0,
});
let overlap_penalty_hamiltonian =
PairwiseCutoff(SitePairOverlapPenalty {
maximum_interaction_range: maximum_radius * 2.0,
});
let initial_box_volume = total_particle_area / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.seed(seed as u32)
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let distribution = PolydisperseBodyDistribution {
boundary: microstate.boundary().clone(),
radii,
};
let quick_insert = QuickInsert::new(distribution, n_disks);
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume = total_particle_area / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(PolydisperseHardDiskModel {
microstate,
hamiltonian,
translate_sweep,
quick_insert,
quick_compress,
overlap_penalty_hamiltonian,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct PolydisperseHardDiskModel {
/// 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<SitePairInteraction>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<PolydisperseBodyDistribution>,
/// The current phase of the simulation.
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<SitePairOverlapPenalty>,
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for PolydisperseHardDiskModel {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl PolydisperseHardDiskModel {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = PolydisperseHardDiskModel::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("polydisperse-hard-disk-model.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 polydisperse_hard_disk_model_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use polydisperse_hard_disk_model_interactive::main;Note
Your simulations will run correctly without these methods (the default
site_pair_energy_initialcallssite_pair_energyand the defaultis_only_infinite_or_zeroreturnsfalse). However, these implementations will make hard shape simulations execute faster. There is no need to compute the energy of the initial state before a trial move as it will always be zero for a hard potential. Even if the initial energy were infinity, a trial move with an initial infinite energy and a finite final energy would always be accepted.
Site-site Overlap Penalty
Define SitePairOverlapPenalty
To use QuickInsert and QuickCompress we need to define another site-site
interaction potential. This one uses OverlapPenalty to allow partially
overlapping sites during initialization:
use anyhow::{Context, anyhow};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use rand::{
SeedableRng,
distr::{Distribution, Uniform},
rngs::StdRng,
};
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
univariate::{Expanded, OverlapPenalty, UnivariateEnergy},
};
use hoomd_mc::{
BodyDistribution, QuickCompress, QuickInsert, Sweep, Translate, Trial,
Tune, TuneOptions,
};
use hoomd_microstate::{
AppendMicrostate, Body, Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_utility::valid::PositiveReal;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's radius.
radius: PositiveReal,
}
impl Transform<SiteProperties> for BodyProperties {
fn transform(&self, site_properties: &SiteProperties) -> SiteProperties {
SiteProperties {
position: self.position + site_properties.position,
..*site_properties
}
}
}
#[derive(MaximumInteractionRange)]
struct SitePairInteraction {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
if r < a.radius.get() + b.radius.get() {
f64::INFINITY
} else {
0.0
}
}
fn site_pair_energy_initial(
&self,
_a: &SiteProperties,
_b: &SiteProperties,
) -> f64 {
0.0
}
fn is_only_infinite_or_zero() -> bool {
true
}
}
#[derive(MaximumInteractionRange)]
struct SitePairOverlapPenalty {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairOverlapPenalty {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
let pair_interaction = Expanded {
delta: a.radius.get() + b.radius.get(),
f: OverlapPenalty::default(),
};
pair_interaction.energy(r)
}
}
struct PolydisperseBodyDistribution {
/// Radius of each disk to insert into the microstate.
radii: Vec<PositiveReal>,
/// Simulation boundary.
boundary: Periodic<Rectangle>,
}
impl BodyDistribution<Body<BodyProperties, SiteProperties>>
for PolydisperseBodyDistribution
{
fn sample<R: rand::Rng + ?Sized>(
&self,
index: usize,
rng: &mut R,
) -> Body<BodyProperties, SiteProperties> {
let properties = Point {
position: self.boundary.sample(rng),
};
let sites = vec![SiteProperties {
position: Cartesian::default(),
radius: self.radii[index],
}];
Body { properties, sites }
}
}
impl<X> AppendMicrostate<BodyProperties, SiteProperties, X, Periodic<Rectangle>>
for HoomdGsdFile
{
#[inline]
fn append_microstate(
&mut self,
microstate: &Microstate<
BodyProperties,
SiteProperties,
X,
Periodic<Rectangle>,
>,
) -> Result<hoomd_gsd::hoomd::Frame<'_>, hoomd_gsd::hoomd::AppendError>
{
self.append_frame(microstate.step())?
.configuration_box(microstate.boundary().shape().to_gsd_box())?
.configuration_dimensions(Dimensions::Two)?
.particles_position(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.position)
.map(|p| [p[0], p[1], 0.0].into()),
)?
.particles_diameter(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.radius.get() * 2.0),
)
}
}
impl PolydisperseHardDiskModel {
/// Construct a new hard disk self-assembly simulation.
fn new() -> anyhow::Result<PolydisperseHardDiskModel> {
let seed = 1;
let minimum_radius = 0.1;
let maximum_radius = 0.8;
let initial_packing_fraction = 0.6;
let target_packing_fraction = 0.72;
let n_disks = 64_usize.pow(2);
let maximum_distance = 0.07;
let macrostate = Isothermal { temperature: 1.0 };
let radius_distribution = Uniform::new(minimum_radius, maximum_radius)?;
let mut rng = StdRng::seed_from_u64(seed);
let mut radii = Vec::with_capacity(n_disks);
for r in radius_distribution.sample_iter(&mut rng).take(n_disks) {
radii.push(r.try_into()?);
}
let total_particle_area = radii.iter().fold(0.0, |total, r| {
let circle = Circle { radius: *r };
total + circle.volume()
});
let hamiltonian = PairwiseCutoff(SitePairInteraction {
maximum_interaction_range: maximum_radius * 2.0,
});
let overlap_penalty_hamiltonian =
PairwiseCutoff(SitePairOverlapPenalty {
maximum_interaction_range: maximum_radius * 2.0,
});
let initial_box_volume = total_particle_area / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.seed(seed as u32)
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let distribution = PolydisperseBodyDistribution {
boundary: microstate.boundary().clone(),
radii,
};
let quick_insert = QuickInsert::new(distribution, n_disks);
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume = total_particle_area / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(PolydisperseHardDiskModel {
microstate,
hamiltonian,
translate_sweep,
quick_insert,
quick_compress,
overlap_penalty_hamiltonian,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct PolydisperseHardDiskModel {
/// 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<SitePairInteraction>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<PolydisperseBodyDistribution>,
/// The current phase of the simulation.
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<SitePairOverlapPenalty>,
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for PolydisperseHardDiskModel {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl PolydisperseHardDiskModel {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = PolydisperseHardDiskModel::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("polydisperse-hard-disk-model.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 polydisperse_hard_disk_model_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use polydisperse_hard_disk_model_interactive::main;As with SitePairInteraction, this type also needs to have a maximum interaction range.
Implement the Overlap Penalty Computation
As in Hard Disk Self-Assembly, you can use an Expanded<OverlapPenalty> to compute
the overlap penalty shifted to the surface of the disks. When the disks are
polydisperse, you need to set delta to the sum of the radii unique to each
pair interaction:
use anyhow::{Context, anyhow};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use rand::{
SeedableRng,
distr::{Distribution, Uniform},
rngs::StdRng,
};
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
univariate::{Expanded, OverlapPenalty, UnivariateEnergy},
};
use hoomd_mc::{
BodyDistribution, QuickCompress, QuickInsert, Sweep, Translate, Trial,
Tune, TuneOptions,
};
use hoomd_microstate::{
AppendMicrostate, Body, Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_utility::valid::PositiveReal;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's radius.
radius: PositiveReal,
}
impl Transform<SiteProperties> for BodyProperties {
fn transform(&self, site_properties: &SiteProperties) -> SiteProperties {
SiteProperties {
position: self.position + site_properties.position,
..*site_properties
}
}
}
#[derive(MaximumInteractionRange)]
struct SitePairInteraction {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
if r < a.radius.get() + b.radius.get() {
f64::INFINITY
} else {
0.0
}
}
fn site_pair_energy_initial(
&self,
_a: &SiteProperties,
_b: &SiteProperties,
) -> f64 {
0.0
}
fn is_only_infinite_or_zero() -> bool {
true
}
}
#[derive(MaximumInteractionRange)]
struct SitePairOverlapPenalty {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairOverlapPenalty {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
let pair_interaction = Expanded {
delta: a.radius.get() + b.radius.get(),
f: OverlapPenalty::default(),
};
pair_interaction.energy(r)
}
}
struct PolydisperseBodyDistribution {
/// Radius of each disk to insert into the microstate.
radii: Vec<PositiveReal>,
/// Simulation boundary.
boundary: Periodic<Rectangle>,
}
impl BodyDistribution<Body<BodyProperties, SiteProperties>>
for PolydisperseBodyDistribution
{
fn sample<R: rand::Rng + ?Sized>(
&self,
index: usize,
rng: &mut R,
) -> Body<BodyProperties, SiteProperties> {
let properties = Point {
position: self.boundary.sample(rng),
};
let sites = vec![SiteProperties {
position: Cartesian::default(),
radius: self.radii[index],
}];
Body { properties, sites }
}
}
impl<X> AppendMicrostate<BodyProperties, SiteProperties, X, Periodic<Rectangle>>
for HoomdGsdFile
{
#[inline]
fn append_microstate(
&mut self,
microstate: &Microstate<
BodyProperties,
SiteProperties,
X,
Periodic<Rectangle>,
>,
) -> Result<hoomd_gsd::hoomd::Frame<'_>, hoomd_gsd::hoomd::AppendError>
{
self.append_frame(microstate.step())?
.configuration_box(microstate.boundary().shape().to_gsd_box())?
.configuration_dimensions(Dimensions::Two)?
.particles_position(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.position)
.map(|p| [p[0], p[1], 0.0].into()),
)?
.particles_diameter(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.radius.get() * 2.0),
)
}
}
impl PolydisperseHardDiskModel {
/// Construct a new hard disk self-assembly simulation.
fn new() -> anyhow::Result<PolydisperseHardDiskModel> {
let seed = 1;
let minimum_radius = 0.1;
let maximum_radius = 0.8;
let initial_packing_fraction = 0.6;
let target_packing_fraction = 0.72;
let n_disks = 64_usize.pow(2);
let maximum_distance = 0.07;
let macrostate = Isothermal { temperature: 1.0 };
let radius_distribution = Uniform::new(minimum_radius, maximum_radius)?;
let mut rng = StdRng::seed_from_u64(seed);
let mut radii = Vec::with_capacity(n_disks);
for r in radius_distribution.sample_iter(&mut rng).take(n_disks) {
radii.push(r.try_into()?);
}
let total_particle_area = radii.iter().fold(0.0, |total, r| {
let circle = Circle { radius: *r };
total + circle.volume()
});
let hamiltonian = PairwiseCutoff(SitePairInteraction {
maximum_interaction_range: maximum_radius * 2.0,
});
let overlap_penalty_hamiltonian =
PairwiseCutoff(SitePairOverlapPenalty {
maximum_interaction_range: maximum_radius * 2.0,
});
let initial_box_volume = total_particle_area / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.seed(seed as u32)
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let distribution = PolydisperseBodyDistribution {
boundary: microstate.boundary().clone(),
radii,
};
let quick_insert = QuickInsert::new(distribution, n_disks);
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume = total_particle_area / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(PolydisperseHardDiskModel {
microstate,
hamiltonian,
translate_sweep,
quick_insert,
quick_compress,
overlap_penalty_hamiltonian,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct PolydisperseHardDiskModel {
/// 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<SitePairInteraction>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<PolydisperseBodyDistribution>,
/// The current phase of the simulation.
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<SitePairOverlapPenalty>,
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for PolydisperseHardDiskModel {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl PolydisperseHardDiskModel {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = PolydisperseHardDiskModel::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("polydisperse-hard-disk-model.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 polydisperse_hard_disk_model_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use polydisperse_hard_disk_model_interactive::main;Body Distribution
QuickInsert randomly draws bodies from a distribution and places them in
the microstate. Hard Ellipse Self-Assembly and similar tutorials use the
UniformIn distribution which randomizes the body’s position (and orientation
when present) but keeps all the other body and site properties fixed.
Define PolydisperseBodyDistribution
Define a custom distribution type that samples bodies whose sites have random radii:
use anyhow::{Context, anyhow};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use rand::{
SeedableRng,
distr::{Distribution, Uniform},
rngs::StdRng,
};
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
univariate::{Expanded, OverlapPenalty, UnivariateEnergy},
};
use hoomd_mc::{
BodyDistribution, QuickCompress, QuickInsert, Sweep, Translate, Trial,
Tune, TuneOptions,
};
use hoomd_microstate::{
AppendMicrostate, Body, Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_utility::valid::PositiveReal;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's radius.
radius: PositiveReal,
}
impl Transform<SiteProperties> for BodyProperties {
fn transform(&self, site_properties: &SiteProperties) -> SiteProperties {
SiteProperties {
position: self.position + site_properties.position,
..*site_properties
}
}
}
#[derive(MaximumInteractionRange)]
struct SitePairInteraction {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
if r < a.radius.get() + b.radius.get() {
f64::INFINITY
} else {
0.0
}
}
fn site_pair_energy_initial(
&self,
_a: &SiteProperties,
_b: &SiteProperties,
) -> f64 {
0.0
}
fn is_only_infinite_or_zero() -> bool {
true
}
}
#[derive(MaximumInteractionRange)]
struct SitePairOverlapPenalty {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairOverlapPenalty {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
let pair_interaction = Expanded {
delta: a.radius.get() + b.radius.get(),
f: OverlapPenalty::default(),
};
pair_interaction.energy(r)
}
}
struct PolydisperseBodyDistribution {
/// Radius of each disk to insert into the microstate.
radii: Vec<PositiveReal>,
/// Simulation boundary.
boundary: Periodic<Rectangle>,
}
impl BodyDistribution<Body<BodyProperties, SiteProperties>>
for PolydisperseBodyDistribution
{
fn sample<R: rand::Rng + ?Sized>(
&self,
index: usize,
rng: &mut R,
) -> Body<BodyProperties, SiteProperties> {
let properties = Point {
position: self.boundary.sample(rng),
};
let sites = vec![SiteProperties {
position: Cartesian::default(),
radius: self.radii[index],
}];
Body { properties, sites }
}
}
impl<X> AppendMicrostate<BodyProperties, SiteProperties, X, Periodic<Rectangle>>
for HoomdGsdFile
{
#[inline]
fn append_microstate(
&mut self,
microstate: &Microstate<
BodyProperties,
SiteProperties,
X,
Periodic<Rectangle>,
>,
) -> Result<hoomd_gsd::hoomd::Frame<'_>, hoomd_gsd::hoomd::AppendError>
{
self.append_frame(microstate.step())?
.configuration_box(microstate.boundary().shape().to_gsd_box())?
.configuration_dimensions(Dimensions::Two)?
.particles_position(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.position)
.map(|p| [p[0], p[1], 0.0].into()),
)?
.particles_diameter(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.radius.get() * 2.0),
)
}
}
impl PolydisperseHardDiskModel {
/// Construct a new hard disk self-assembly simulation.
fn new() -> anyhow::Result<PolydisperseHardDiskModel> {
let seed = 1;
let minimum_radius = 0.1;
let maximum_radius = 0.8;
let initial_packing_fraction = 0.6;
let target_packing_fraction = 0.72;
let n_disks = 64_usize.pow(2);
let maximum_distance = 0.07;
let macrostate = Isothermal { temperature: 1.0 };
let radius_distribution = Uniform::new(minimum_radius, maximum_radius)?;
let mut rng = StdRng::seed_from_u64(seed);
let mut radii = Vec::with_capacity(n_disks);
for r in radius_distribution.sample_iter(&mut rng).take(n_disks) {
radii.push(r.try_into()?);
}
let total_particle_area = radii.iter().fold(0.0, |total, r| {
let circle = Circle { radius: *r };
total + circle.volume()
});
let hamiltonian = PairwiseCutoff(SitePairInteraction {
maximum_interaction_range: maximum_radius * 2.0,
});
let overlap_penalty_hamiltonian =
PairwiseCutoff(SitePairOverlapPenalty {
maximum_interaction_range: maximum_radius * 2.0,
});
let initial_box_volume = total_particle_area / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.seed(seed as u32)
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let distribution = PolydisperseBodyDistribution {
boundary: microstate.boundary().clone(),
radii,
};
let quick_insert = QuickInsert::new(distribution, n_disks);
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume = total_particle_area / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(PolydisperseHardDiskModel {
microstate,
hamiltonian,
translate_sweep,
quick_insert,
quick_compress,
overlap_penalty_hamiltonian,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct PolydisperseHardDiskModel {
/// 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<SitePairInteraction>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<PolydisperseBodyDistribution>,
/// The current phase of the simulation.
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<SitePairOverlapPenalty>,
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for PolydisperseHardDiskModel {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl PolydisperseHardDiskModel {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = PolydisperseHardDiskModel::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("polydisperse-hard-disk-model.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 polydisperse_hard_disk_model_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use polydisperse_hard_disk_model_interactive::main;Warning
You cannot simply draw a new random radius for each sample. If you did, the larger sites would be more likely to overlap with existing ones and
QuickInsertwould strongly bias toward placing smaller sites.
To ensure that the inserted sites have an unbiased distribution of radii, precompute the radii and store them in the distribution type.
Implement BodyDistribution
Implement the BodyDistribution trait:
use anyhow::{Context, anyhow};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use rand::{
SeedableRng,
distr::{Distribution, Uniform},
rngs::StdRng,
};
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
univariate::{Expanded, OverlapPenalty, UnivariateEnergy},
};
use hoomd_mc::{
BodyDistribution, QuickCompress, QuickInsert, Sweep, Translate, Trial,
Tune, TuneOptions,
};
use hoomd_microstate::{
AppendMicrostate, Body, Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_utility::valid::PositiveReal;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's radius.
radius: PositiveReal,
}
impl Transform<SiteProperties> for BodyProperties {
fn transform(&self, site_properties: &SiteProperties) -> SiteProperties {
SiteProperties {
position: self.position + site_properties.position,
..*site_properties
}
}
}
#[derive(MaximumInteractionRange)]
struct SitePairInteraction {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
if r < a.radius.get() + b.radius.get() {
f64::INFINITY
} else {
0.0
}
}
fn site_pair_energy_initial(
&self,
_a: &SiteProperties,
_b: &SiteProperties,
) -> f64 {
0.0
}
fn is_only_infinite_or_zero() -> bool {
true
}
}
#[derive(MaximumInteractionRange)]
struct SitePairOverlapPenalty {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairOverlapPenalty {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
let pair_interaction = Expanded {
delta: a.radius.get() + b.radius.get(),
f: OverlapPenalty::default(),
};
pair_interaction.energy(r)
}
}
struct PolydisperseBodyDistribution {
/// Radius of each disk to insert into the microstate.
radii: Vec<PositiveReal>,
/// Simulation boundary.
boundary: Periodic<Rectangle>,
}
impl BodyDistribution<Body<BodyProperties, SiteProperties>>
for PolydisperseBodyDistribution
{
fn sample<R: rand::Rng + ?Sized>(
&self,
index: usize,
rng: &mut R,
) -> Body<BodyProperties, SiteProperties> {
let properties = Point {
position: self.boundary.sample(rng),
};
let sites = vec![SiteProperties {
position: Cartesian::default(),
radius: self.radii[index],
}];
Body { properties, sites }
}
}
impl<X> AppendMicrostate<BodyProperties, SiteProperties, X, Periodic<Rectangle>>
for HoomdGsdFile
{
#[inline]
fn append_microstate(
&mut self,
microstate: &Microstate<
BodyProperties,
SiteProperties,
X,
Periodic<Rectangle>,
>,
) -> Result<hoomd_gsd::hoomd::Frame<'_>, hoomd_gsd::hoomd::AppendError>
{
self.append_frame(microstate.step())?
.configuration_box(microstate.boundary().shape().to_gsd_box())?
.configuration_dimensions(Dimensions::Two)?
.particles_position(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.position)
.map(|p| [p[0], p[1], 0.0].into()),
)?
.particles_diameter(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.radius.get() * 2.0),
)
}
}
impl PolydisperseHardDiskModel {
/// Construct a new hard disk self-assembly simulation.
fn new() -> anyhow::Result<PolydisperseHardDiskModel> {
let seed = 1;
let minimum_radius = 0.1;
let maximum_radius = 0.8;
let initial_packing_fraction = 0.6;
let target_packing_fraction = 0.72;
let n_disks = 64_usize.pow(2);
let maximum_distance = 0.07;
let macrostate = Isothermal { temperature: 1.0 };
let radius_distribution = Uniform::new(minimum_radius, maximum_radius)?;
let mut rng = StdRng::seed_from_u64(seed);
let mut radii = Vec::with_capacity(n_disks);
for r in radius_distribution.sample_iter(&mut rng).take(n_disks) {
radii.push(r.try_into()?);
}
let total_particle_area = radii.iter().fold(0.0, |total, r| {
let circle = Circle { radius: *r };
total + circle.volume()
});
let hamiltonian = PairwiseCutoff(SitePairInteraction {
maximum_interaction_range: maximum_radius * 2.0,
});
let overlap_penalty_hamiltonian =
PairwiseCutoff(SitePairOverlapPenalty {
maximum_interaction_range: maximum_radius * 2.0,
});
let initial_box_volume = total_particle_area / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.seed(seed as u32)
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let distribution = PolydisperseBodyDistribution {
boundary: microstate.boundary().clone(),
radii,
};
let quick_insert = QuickInsert::new(distribution, n_disks);
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume = total_particle_area / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(PolydisperseHardDiskModel {
microstate,
hamiltonian,
translate_sweep,
quick_insert,
quick_compress,
overlap_penalty_hamiltonian,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct PolydisperseHardDiskModel {
/// 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<SitePairInteraction>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<PolydisperseBodyDistribution>,
/// The current phase of the simulation.
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<SitePairOverlapPenalty>,
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for PolydisperseHardDiskModel {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl PolydisperseHardDiskModel {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = PolydisperseHardDiskModel::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("polydisperse-hard-disk-model.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 polydisperse_hard_disk_model_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use polydisperse_hard_disk_model_interactive::main;QuickInsert calls this sample method each time it attempts to place a new body.
It increments the index argument after each successfully placed site.
This implementation of sample chooses a position randomly in the simulation boundary
and chooses a unique (but predetermined) radius based on index.
Implement AppendMicrostate
As in Type-dependent Interactions, implement the AppendMicrostate trait
to write the custom SiteProperties type. In this case, use the
particles_diameter method to write the unique diameter of each site:
use anyhow::{Context, anyhow};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use rand::{
SeedableRng,
distr::{Distribution, Uniform},
rngs::StdRng,
};
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
univariate::{Expanded, OverlapPenalty, UnivariateEnergy},
};
use hoomd_mc::{
BodyDistribution, QuickCompress, QuickInsert, Sweep, Translate, Trial,
Tune, TuneOptions,
};
use hoomd_microstate::{
AppendMicrostate, Body, Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_utility::valid::PositiveReal;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's radius.
radius: PositiveReal,
}
impl Transform<SiteProperties> for BodyProperties {
fn transform(&self, site_properties: &SiteProperties) -> SiteProperties {
SiteProperties {
position: self.position + site_properties.position,
..*site_properties
}
}
}
#[derive(MaximumInteractionRange)]
struct SitePairInteraction {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
if r < a.radius.get() + b.radius.get() {
f64::INFINITY
} else {
0.0
}
}
fn site_pair_energy_initial(
&self,
_a: &SiteProperties,
_b: &SiteProperties,
) -> f64 {
0.0
}
fn is_only_infinite_or_zero() -> bool {
true
}
}
#[derive(MaximumInteractionRange)]
struct SitePairOverlapPenalty {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairOverlapPenalty {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
let pair_interaction = Expanded {
delta: a.radius.get() + b.radius.get(),
f: OverlapPenalty::default(),
};
pair_interaction.energy(r)
}
}
struct PolydisperseBodyDistribution {
/// Radius of each disk to insert into the microstate.
radii: Vec<PositiveReal>,
/// Simulation boundary.
boundary: Periodic<Rectangle>,
}
impl BodyDistribution<Body<BodyProperties, SiteProperties>>
for PolydisperseBodyDistribution
{
fn sample<R: rand::Rng + ?Sized>(
&self,
index: usize,
rng: &mut R,
) -> Body<BodyProperties, SiteProperties> {
let properties = Point {
position: self.boundary.sample(rng),
};
let sites = vec![SiteProperties {
position: Cartesian::default(),
radius: self.radii[index],
}];
Body { properties, sites }
}
}
impl<X> AppendMicrostate<BodyProperties, SiteProperties, X, Periodic<Rectangle>>
for HoomdGsdFile
{
#[inline]
fn append_microstate(
&mut self,
microstate: &Microstate<
BodyProperties,
SiteProperties,
X,
Periodic<Rectangle>,
>,
) -> Result<hoomd_gsd::hoomd::Frame<'_>, hoomd_gsd::hoomd::AppendError>
{
self.append_frame(microstate.step())?
.configuration_box(microstate.boundary().shape().to_gsd_box())?
.configuration_dimensions(Dimensions::Two)?
.particles_position(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.position)
.map(|p| [p[0], p[1], 0.0].into()),
)?
.particles_diameter(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.radius.get() * 2.0),
)
}
}
impl PolydisperseHardDiskModel {
/// Construct a new hard disk self-assembly simulation.
fn new() -> anyhow::Result<PolydisperseHardDiskModel> {
let seed = 1;
let minimum_radius = 0.1;
let maximum_radius = 0.8;
let initial_packing_fraction = 0.6;
let target_packing_fraction = 0.72;
let n_disks = 64_usize.pow(2);
let maximum_distance = 0.07;
let macrostate = Isothermal { temperature: 1.0 };
let radius_distribution = Uniform::new(minimum_radius, maximum_radius)?;
let mut rng = StdRng::seed_from_u64(seed);
let mut radii = Vec::with_capacity(n_disks);
for r in radius_distribution.sample_iter(&mut rng).take(n_disks) {
radii.push(r.try_into()?);
}
let total_particle_area = radii.iter().fold(0.0, |total, r| {
let circle = Circle { radius: *r };
total + circle.volume()
});
let hamiltonian = PairwiseCutoff(SitePairInteraction {
maximum_interaction_range: maximum_radius * 2.0,
});
let overlap_penalty_hamiltonian =
PairwiseCutoff(SitePairOverlapPenalty {
maximum_interaction_range: maximum_radius * 2.0,
});
let initial_box_volume = total_particle_area / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.seed(seed as u32)
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let distribution = PolydisperseBodyDistribution {
boundary: microstate.boundary().clone(),
radii,
};
let quick_insert = QuickInsert::new(distribution, n_disks);
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume = total_particle_area / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(PolydisperseHardDiskModel {
microstate,
hamiltonian,
translate_sweep,
quick_insert,
quick_compress,
overlap_penalty_hamiltonian,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct PolydisperseHardDiskModel {
/// 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<SitePairInteraction>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<PolydisperseBodyDistribution>,
/// The current phase of the simulation.
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<SitePairOverlapPenalty>,
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for PolydisperseHardDiskModel {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl PolydisperseHardDiskModel {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = PolydisperseHardDiskModel::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("polydisperse-hard-disk-model.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 polydisperse_hard_disk_model_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use polydisperse_hard_disk_model_interactive::main;Construct the Simulation Model
The new() method constructs a new simulation model:
use anyhow::{Context, anyhow};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use rand::{
SeedableRng,
distr::{Distribution, Uniform},
rngs::StdRng,
};
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
univariate::{Expanded, OverlapPenalty, UnivariateEnergy},
};
use hoomd_mc::{
BodyDistribution, QuickCompress, QuickInsert, Sweep, Translate, Trial,
Tune, TuneOptions,
};
use hoomd_microstate::{
AppendMicrostate, Body, Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_utility::valid::PositiveReal;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's radius.
radius: PositiveReal,
}
impl Transform<SiteProperties> for BodyProperties {
fn transform(&self, site_properties: &SiteProperties) -> SiteProperties {
SiteProperties {
position: self.position + site_properties.position,
..*site_properties
}
}
}
#[derive(MaximumInteractionRange)]
struct SitePairInteraction {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
if r < a.radius.get() + b.radius.get() {
f64::INFINITY
} else {
0.0
}
}
fn site_pair_energy_initial(
&self,
_a: &SiteProperties,
_b: &SiteProperties,
) -> f64 {
0.0
}
fn is_only_infinite_or_zero() -> bool {
true
}
}
#[derive(MaximumInteractionRange)]
struct SitePairOverlapPenalty {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairOverlapPenalty {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
let pair_interaction = Expanded {
delta: a.radius.get() + b.radius.get(),
f: OverlapPenalty::default(),
};
pair_interaction.energy(r)
}
}
struct PolydisperseBodyDistribution {
/// Radius of each disk to insert into the microstate.
radii: Vec<PositiveReal>,
/// Simulation boundary.
boundary: Periodic<Rectangle>,
}
impl BodyDistribution<Body<BodyProperties, SiteProperties>>
for PolydisperseBodyDistribution
{
fn sample<R: rand::Rng + ?Sized>(
&self,
index: usize,
rng: &mut R,
) -> Body<BodyProperties, SiteProperties> {
let properties = Point {
position: self.boundary.sample(rng),
};
let sites = vec![SiteProperties {
position: Cartesian::default(),
radius: self.radii[index],
}];
Body { properties, sites }
}
}
impl<X> AppendMicrostate<BodyProperties, SiteProperties, X, Periodic<Rectangle>>
for HoomdGsdFile
{
#[inline]
fn append_microstate(
&mut self,
microstate: &Microstate<
BodyProperties,
SiteProperties,
X,
Periodic<Rectangle>,
>,
) -> Result<hoomd_gsd::hoomd::Frame<'_>, hoomd_gsd::hoomd::AppendError>
{
self.append_frame(microstate.step())?
.configuration_box(microstate.boundary().shape().to_gsd_box())?
.configuration_dimensions(Dimensions::Two)?
.particles_position(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.position)
.map(|p| [p[0], p[1], 0.0].into()),
)?
.particles_diameter(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.radius.get() * 2.0),
)
}
}
impl PolydisperseHardDiskModel {
/// Construct a new hard disk self-assembly simulation.
fn new() -> anyhow::Result<PolydisperseHardDiskModel> {
let seed = 1;
let minimum_radius = 0.1;
let maximum_radius = 0.8;
let initial_packing_fraction = 0.6;
let target_packing_fraction = 0.72;
let n_disks = 64_usize.pow(2);
let maximum_distance = 0.07;
let macrostate = Isothermal { temperature: 1.0 };
let radius_distribution = Uniform::new(minimum_radius, maximum_radius)?;
let mut rng = StdRng::seed_from_u64(seed);
let mut radii = Vec::with_capacity(n_disks);
for r in radius_distribution.sample_iter(&mut rng).take(n_disks) {
radii.push(r.try_into()?);
}
let total_particle_area = radii.iter().fold(0.0, |total, r| {
let circle = Circle { radius: *r };
total + circle.volume()
});
let hamiltonian = PairwiseCutoff(SitePairInteraction {
maximum_interaction_range: maximum_radius * 2.0,
});
let overlap_penalty_hamiltonian =
PairwiseCutoff(SitePairOverlapPenalty {
maximum_interaction_range: maximum_radius * 2.0,
});
let initial_box_volume = total_particle_area / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.seed(seed as u32)
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let distribution = PolydisperseBodyDistribution {
boundary: microstate.boundary().clone(),
radii,
};
let quick_insert = QuickInsert::new(distribution, n_disks);
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume = total_particle_area / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(PolydisperseHardDiskModel {
microstate,
hamiltonian,
translate_sweep,
quick_insert,
quick_compress,
overlap_penalty_hamiltonian,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct PolydisperseHardDiskModel {
/// 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<SitePairInteraction>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<PolydisperseBodyDistribution>,
/// The current phase of the simulation.
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<SitePairOverlapPenalty>,
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for PolydisperseHardDiskModel {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl PolydisperseHardDiskModel {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = PolydisperseHardDiskModel::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("polydisperse-hard-disk-model.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 polydisperse_hard_disk_model_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use polydisperse_hard_disk_model_interactive::main;Parameters
Assign all the model parameters in one code block:
use anyhow::{Context, anyhow};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use rand::{
SeedableRng,
distr::{Distribution, Uniform},
rngs::StdRng,
};
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
univariate::{Expanded, OverlapPenalty, UnivariateEnergy},
};
use hoomd_mc::{
BodyDistribution, QuickCompress, QuickInsert, Sweep, Translate, Trial,
Tune, TuneOptions,
};
use hoomd_microstate::{
AppendMicrostate, Body, Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_utility::valid::PositiveReal;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's radius.
radius: PositiveReal,
}
impl Transform<SiteProperties> for BodyProperties {
fn transform(&self, site_properties: &SiteProperties) -> SiteProperties {
SiteProperties {
position: self.position + site_properties.position,
..*site_properties
}
}
}
#[derive(MaximumInteractionRange)]
struct SitePairInteraction {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
if r < a.radius.get() + b.radius.get() {
f64::INFINITY
} else {
0.0
}
}
fn site_pair_energy_initial(
&self,
_a: &SiteProperties,
_b: &SiteProperties,
) -> f64 {
0.0
}
fn is_only_infinite_or_zero() -> bool {
true
}
}
#[derive(MaximumInteractionRange)]
struct SitePairOverlapPenalty {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairOverlapPenalty {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
let pair_interaction = Expanded {
delta: a.radius.get() + b.radius.get(),
f: OverlapPenalty::default(),
};
pair_interaction.energy(r)
}
}
struct PolydisperseBodyDistribution {
/// Radius of each disk to insert into the microstate.
radii: Vec<PositiveReal>,
/// Simulation boundary.
boundary: Periodic<Rectangle>,
}
impl BodyDistribution<Body<BodyProperties, SiteProperties>>
for PolydisperseBodyDistribution
{
fn sample<R: rand::Rng + ?Sized>(
&self,
index: usize,
rng: &mut R,
) -> Body<BodyProperties, SiteProperties> {
let properties = Point {
position: self.boundary.sample(rng),
};
let sites = vec![SiteProperties {
position: Cartesian::default(),
radius: self.radii[index],
}];
Body { properties, sites }
}
}
impl<X> AppendMicrostate<BodyProperties, SiteProperties, X, Periodic<Rectangle>>
for HoomdGsdFile
{
#[inline]
fn append_microstate(
&mut self,
microstate: &Microstate<
BodyProperties,
SiteProperties,
X,
Periodic<Rectangle>,
>,
) -> Result<hoomd_gsd::hoomd::Frame<'_>, hoomd_gsd::hoomd::AppendError>
{
self.append_frame(microstate.step())?
.configuration_box(microstate.boundary().shape().to_gsd_box())?
.configuration_dimensions(Dimensions::Two)?
.particles_position(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.position)
.map(|p| [p[0], p[1], 0.0].into()),
)?
.particles_diameter(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.radius.get() * 2.0),
)
}
}
impl PolydisperseHardDiskModel {
/// Construct a new hard disk self-assembly simulation.
fn new() -> anyhow::Result<PolydisperseHardDiskModel> {
let seed = 1;
let minimum_radius = 0.1;
let maximum_radius = 0.8;
let initial_packing_fraction = 0.6;
let target_packing_fraction = 0.72;
let n_disks = 64_usize.pow(2);
let maximum_distance = 0.07;
let macrostate = Isothermal { temperature: 1.0 };
let radius_distribution = Uniform::new(minimum_radius, maximum_radius)?;
let mut rng = StdRng::seed_from_u64(seed);
let mut radii = Vec::with_capacity(n_disks);
for r in radius_distribution.sample_iter(&mut rng).take(n_disks) {
radii.push(r.try_into()?);
}
let total_particle_area = radii.iter().fold(0.0, |total, r| {
let circle = Circle { radius: *r };
total + circle.volume()
});
let hamiltonian = PairwiseCutoff(SitePairInteraction {
maximum_interaction_range: maximum_radius * 2.0,
});
let overlap_penalty_hamiltonian =
PairwiseCutoff(SitePairOverlapPenalty {
maximum_interaction_range: maximum_radius * 2.0,
});
let initial_box_volume = total_particle_area / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.seed(seed as u32)
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let distribution = PolydisperseBodyDistribution {
boundary: microstate.boundary().clone(),
radii,
};
let quick_insert = QuickInsert::new(distribution, n_disks);
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume = total_particle_area / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(PolydisperseHardDiskModel {
microstate,
hamiltonian,
translate_sweep,
quick_insert,
quick_compress,
overlap_penalty_hamiltonian,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct PolydisperseHardDiskModel {
/// 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<SitePairInteraction>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<PolydisperseBodyDistribution>,
/// The current phase of the simulation.
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<SitePairOverlapPenalty>,
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for PolydisperseHardDiskModel {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl PolydisperseHardDiskModel {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = PolydisperseHardDiskModel::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("polydisperse-hard-disk-model.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 polydisperse_hard_disk_model_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use polydisperse_hard_disk_model_interactive::main;seed is the random number seed. Change this value to select different random
distribution of radii. minimum_radius and maximum_radius set the limits of
the radii distribution. initial_packing_fraction is the area of the disks
divided by the area of the simulation boundary in the initial state. Choose
this value so that disks can be placed easily in the microstate. During the
Initialize phase, the microstate will be compressed until it reaches the
packing fraction target_packing_fraction. n_disks is the number of disks
to add, maximum_distance is the largest distance a translation trial move can
take (initially), and macrostate holds the current temperature set point (in
units of energy).
Precompute the Radii
Use the rand crate to randomly sample radii from a uniform distribution:
use anyhow::{Context, anyhow};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use rand::{
SeedableRng,
distr::{Distribution, Uniform},
rngs::StdRng,
};
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
univariate::{Expanded, OverlapPenalty, UnivariateEnergy},
};
use hoomd_mc::{
BodyDistribution, QuickCompress, QuickInsert, Sweep, Translate, Trial,
Tune, TuneOptions,
};
use hoomd_microstate::{
AppendMicrostate, Body, Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_utility::valid::PositiveReal;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's radius.
radius: PositiveReal,
}
impl Transform<SiteProperties> for BodyProperties {
fn transform(&self, site_properties: &SiteProperties) -> SiteProperties {
SiteProperties {
position: self.position + site_properties.position,
..*site_properties
}
}
}
#[derive(MaximumInteractionRange)]
struct SitePairInteraction {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
if r < a.radius.get() + b.radius.get() {
f64::INFINITY
} else {
0.0
}
}
fn site_pair_energy_initial(
&self,
_a: &SiteProperties,
_b: &SiteProperties,
) -> f64 {
0.0
}
fn is_only_infinite_or_zero() -> bool {
true
}
}
#[derive(MaximumInteractionRange)]
struct SitePairOverlapPenalty {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairOverlapPenalty {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
let pair_interaction = Expanded {
delta: a.radius.get() + b.radius.get(),
f: OverlapPenalty::default(),
};
pair_interaction.energy(r)
}
}
struct PolydisperseBodyDistribution {
/// Radius of each disk to insert into the microstate.
radii: Vec<PositiveReal>,
/// Simulation boundary.
boundary: Periodic<Rectangle>,
}
impl BodyDistribution<Body<BodyProperties, SiteProperties>>
for PolydisperseBodyDistribution
{
fn sample<R: rand::Rng + ?Sized>(
&self,
index: usize,
rng: &mut R,
) -> Body<BodyProperties, SiteProperties> {
let properties = Point {
position: self.boundary.sample(rng),
};
let sites = vec![SiteProperties {
position: Cartesian::default(),
radius: self.radii[index],
}];
Body { properties, sites }
}
}
impl<X> AppendMicrostate<BodyProperties, SiteProperties, X, Periodic<Rectangle>>
for HoomdGsdFile
{
#[inline]
fn append_microstate(
&mut self,
microstate: &Microstate<
BodyProperties,
SiteProperties,
X,
Periodic<Rectangle>,
>,
) -> Result<hoomd_gsd::hoomd::Frame<'_>, hoomd_gsd::hoomd::AppendError>
{
self.append_frame(microstate.step())?
.configuration_box(microstate.boundary().shape().to_gsd_box())?
.configuration_dimensions(Dimensions::Two)?
.particles_position(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.position)
.map(|p| [p[0], p[1], 0.0].into()),
)?
.particles_diameter(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.radius.get() * 2.0),
)
}
}
impl PolydisperseHardDiskModel {
/// Construct a new hard disk self-assembly simulation.
fn new() -> anyhow::Result<PolydisperseHardDiskModel> {
let seed = 1;
let minimum_radius = 0.1;
let maximum_radius = 0.8;
let initial_packing_fraction = 0.6;
let target_packing_fraction = 0.72;
let n_disks = 64_usize.pow(2);
let maximum_distance = 0.07;
let macrostate = Isothermal { temperature: 1.0 };
let radius_distribution = Uniform::new(minimum_radius, maximum_radius)?;
let mut rng = StdRng::seed_from_u64(seed);
let mut radii = Vec::with_capacity(n_disks);
for r in radius_distribution.sample_iter(&mut rng).take(n_disks) {
radii.push(r.try_into()?);
}
let total_particle_area = radii.iter().fold(0.0, |total, r| {
let circle = Circle { radius: *r };
total + circle.volume()
});
let hamiltonian = PairwiseCutoff(SitePairInteraction {
maximum_interaction_range: maximum_radius * 2.0,
});
let overlap_penalty_hamiltonian =
PairwiseCutoff(SitePairOverlapPenalty {
maximum_interaction_range: maximum_radius * 2.0,
});
let initial_box_volume = total_particle_area / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.seed(seed as u32)
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let distribution = PolydisperseBodyDistribution {
boundary: microstate.boundary().clone(),
radii,
};
let quick_insert = QuickInsert::new(distribution, n_disks);
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume = total_particle_area / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(PolydisperseHardDiskModel {
microstate,
hamiltonian,
translate_sweep,
quick_insert,
quick_compress,
overlap_penalty_hamiltonian,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct PolydisperseHardDiskModel {
/// 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<SitePairInteraction>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<PolydisperseBodyDistribution>,
/// The current phase of the simulation.
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<SitePairOverlapPenalty>,
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for PolydisperseHardDiskModel {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl PolydisperseHardDiskModel {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = PolydisperseHardDiskModel::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("polydisperse-hard-disk-model.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 polydisperse_hard_disk_model_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use polydisperse_hard_disk_model_interactive::main;Place the sampled radii in a Vec.
Tip
This example uses a uniform distribution for simplicity. A normal distribution would likely be more appropriate for production use.
Particle Area
Sum the area of all sites:
use anyhow::{Context, anyhow};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use rand::{
SeedableRng,
distr::{Distribution, Uniform},
rngs::StdRng,
};
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
univariate::{Expanded, OverlapPenalty, UnivariateEnergy},
};
use hoomd_mc::{
BodyDistribution, QuickCompress, QuickInsert, Sweep, Translate, Trial,
Tune, TuneOptions,
};
use hoomd_microstate::{
AppendMicrostate, Body, Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_utility::valid::PositiveReal;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's radius.
radius: PositiveReal,
}
impl Transform<SiteProperties> for BodyProperties {
fn transform(&self, site_properties: &SiteProperties) -> SiteProperties {
SiteProperties {
position: self.position + site_properties.position,
..*site_properties
}
}
}
#[derive(MaximumInteractionRange)]
struct SitePairInteraction {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
if r < a.radius.get() + b.radius.get() {
f64::INFINITY
} else {
0.0
}
}
fn site_pair_energy_initial(
&self,
_a: &SiteProperties,
_b: &SiteProperties,
) -> f64 {
0.0
}
fn is_only_infinite_or_zero() -> bool {
true
}
}
#[derive(MaximumInteractionRange)]
struct SitePairOverlapPenalty {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairOverlapPenalty {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
let pair_interaction = Expanded {
delta: a.radius.get() + b.radius.get(),
f: OverlapPenalty::default(),
};
pair_interaction.energy(r)
}
}
struct PolydisperseBodyDistribution {
/// Radius of each disk to insert into the microstate.
radii: Vec<PositiveReal>,
/// Simulation boundary.
boundary: Periodic<Rectangle>,
}
impl BodyDistribution<Body<BodyProperties, SiteProperties>>
for PolydisperseBodyDistribution
{
fn sample<R: rand::Rng + ?Sized>(
&self,
index: usize,
rng: &mut R,
) -> Body<BodyProperties, SiteProperties> {
let properties = Point {
position: self.boundary.sample(rng),
};
let sites = vec![SiteProperties {
position: Cartesian::default(),
radius: self.radii[index],
}];
Body { properties, sites }
}
}
impl<X> AppendMicrostate<BodyProperties, SiteProperties, X, Periodic<Rectangle>>
for HoomdGsdFile
{
#[inline]
fn append_microstate(
&mut self,
microstate: &Microstate<
BodyProperties,
SiteProperties,
X,
Periodic<Rectangle>,
>,
) -> Result<hoomd_gsd::hoomd::Frame<'_>, hoomd_gsd::hoomd::AppendError>
{
self.append_frame(microstate.step())?
.configuration_box(microstate.boundary().shape().to_gsd_box())?
.configuration_dimensions(Dimensions::Two)?
.particles_position(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.position)
.map(|p| [p[0], p[1], 0.0].into()),
)?
.particles_diameter(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.radius.get() * 2.0),
)
}
}
impl PolydisperseHardDiskModel {
/// Construct a new hard disk self-assembly simulation.
fn new() -> anyhow::Result<PolydisperseHardDiskModel> {
let seed = 1;
let minimum_radius = 0.1;
let maximum_radius = 0.8;
let initial_packing_fraction = 0.6;
let target_packing_fraction = 0.72;
let n_disks = 64_usize.pow(2);
let maximum_distance = 0.07;
let macrostate = Isothermal { temperature: 1.0 };
let radius_distribution = Uniform::new(minimum_radius, maximum_radius)?;
let mut rng = StdRng::seed_from_u64(seed);
let mut radii = Vec::with_capacity(n_disks);
for r in radius_distribution.sample_iter(&mut rng).take(n_disks) {
radii.push(r.try_into()?);
}
let total_particle_area = radii.iter().fold(0.0, |total, r| {
let circle = Circle { radius: *r };
total + circle.volume()
});
let hamiltonian = PairwiseCutoff(SitePairInteraction {
maximum_interaction_range: maximum_radius * 2.0,
});
let overlap_penalty_hamiltonian =
PairwiseCutoff(SitePairOverlapPenalty {
maximum_interaction_range: maximum_radius * 2.0,
});
let initial_box_volume = total_particle_area / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.seed(seed as u32)
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let distribution = PolydisperseBodyDistribution {
boundary: microstate.boundary().clone(),
radii,
};
let quick_insert = QuickInsert::new(distribution, n_disks);
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume = total_particle_area / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(PolydisperseHardDiskModel {
microstate,
hamiltonian,
translate_sweep,
quick_insert,
quick_compress,
overlap_penalty_hamiltonian,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct PolydisperseHardDiskModel {
/// 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<SitePairInteraction>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<PolydisperseBodyDistribution>,
/// The current phase of the simulation.
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<SitePairOverlapPenalty>,
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for PolydisperseHardDiskModel {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl PolydisperseHardDiskModel {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = PolydisperseHardDiskModel::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("polydisperse-hard-disk-model.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 polydisperse_hard_disk_model_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use polydisperse_hard_disk_model_interactive::main;This value will later be used to compute the total area of the simulation boundary from the packing fraction.
Hamiltonian
Construct the SitePairInteraction type and use it with PairwiseCutoff
to form the system’s Hamiltonian. Similarly construct the overlap penalty
Hamiltonian using the SitePairOverlapPenalty implemented above:
use anyhow::{Context, anyhow};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use rand::{
SeedableRng,
distr::{Distribution, Uniform},
rngs::StdRng,
};
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
univariate::{Expanded, OverlapPenalty, UnivariateEnergy},
};
use hoomd_mc::{
BodyDistribution, QuickCompress, QuickInsert, Sweep, Translate, Trial,
Tune, TuneOptions,
};
use hoomd_microstate::{
AppendMicrostate, Body, Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_utility::valid::PositiveReal;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's radius.
radius: PositiveReal,
}
impl Transform<SiteProperties> for BodyProperties {
fn transform(&self, site_properties: &SiteProperties) -> SiteProperties {
SiteProperties {
position: self.position + site_properties.position,
..*site_properties
}
}
}
#[derive(MaximumInteractionRange)]
struct SitePairInteraction {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
if r < a.radius.get() + b.radius.get() {
f64::INFINITY
} else {
0.0
}
}
fn site_pair_energy_initial(
&self,
_a: &SiteProperties,
_b: &SiteProperties,
) -> f64 {
0.0
}
fn is_only_infinite_or_zero() -> bool {
true
}
}
#[derive(MaximumInteractionRange)]
struct SitePairOverlapPenalty {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairOverlapPenalty {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
let pair_interaction = Expanded {
delta: a.radius.get() + b.radius.get(),
f: OverlapPenalty::default(),
};
pair_interaction.energy(r)
}
}
struct PolydisperseBodyDistribution {
/// Radius of each disk to insert into the microstate.
radii: Vec<PositiveReal>,
/// Simulation boundary.
boundary: Periodic<Rectangle>,
}
impl BodyDistribution<Body<BodyProperties, SiteProperties>>
for PolydisperseBodyDistribution
{
fn sample<R: rand::Rng + ?Sized>(
&self,
index: usize,
rng: &mut R,
) -> Body<BodyProperties, SiteProperties> {
let properties = Point {
position: self.boundary.sample(rng),
};
let sites = vec![SiteProperties {
position: Cartesian::default(),
radius: self.radii[index],
}];
Body { properties, sites }
}
}
impl<X> AppendMicrostate<BodyProperties, SiteProperties, X, Periodic<Rectangle>>
for HoomdGsdFile
{
#[inline]
fn append_microstate(
&mut self,
microstate: &Microstate<
BodyProperties,
SiteProperties,
X,
Periodic<Rectangle>,
>,
) -> Result<hoomd_gsd::hoomd::Frame<'_>, hoomd_gsd::hoomd::AppendError>
{
self.append_frame(microstate.step())?
.configuration_box(microstate.boundary().shape().to_gsd_box())?
.configuration_dimensions(Dimensions::Two)?
.particles_position(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.position)
.map(|p| [p[0], p[1], 0.0].into()),
)?
.particles_diameter(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.radius.get() * 2.0),
)
}
}
impl PolydisperseHardDiskModel {
/// Construct a new hard disk self-assembly simulation.
fn new() -> anyhow::Result<PolydisperseHardDiskModel> {
let seed = 1;
let minimum_radius = 0.1;
let maximum_radius = 0.8;
let initial_packing_fraction = 0.6;
let target_packing_fraction = 0.72;
let n_disks = 64_usize.pow(2);
let maximum_distance = 0.07;
let macrostate = Isothermal { temperature: 1.0 };
let radius_distribution = Uniform::new(minimum_radius, maximum_radius)?;
let mut rng = StdRng::seed_from_u64(seed);
let mut radii = Vec::with_capacity(n_disks);
for r in radius_distribution.sample_iter(&mut rng).take(n_disks) {
radii.push(r.try_into()?);
}
let total_particle_area = radii.iter().fold(0.0, |total, r| {
let circle = Circle { radius: *r };
total + circle.volume()
});
let hamiltonian = PairwiseCutoff(SitePairInteraction {
maximum_interaction_range: maximum_radius * 2.0,
});
let overlap_penalty_hamiltonian =
PairwiseCutoff(SitePairOverlapPenalty {
maximum_interaction_range: maximum_radius * 2.0,
});
let initial_box_volume = total_particle_area / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.seed(seed as u32)
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let distribution = PolydisperseBodyDistribution {
boundary: microstate.boundary().clone(),
radii,
};
let quick_insert = QuickInsert::new(distribution, n_disks);
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume = total_particle_area / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(PolydisperseHardDiskModel {
microstate,
hamiltonian,
translate_sweep,
quick_insert,
quick_compress,
overlap_penalty_hamiltonian,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct PolydisperseHardDiskModel {
/// 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<SitePairInteraction>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<PolydisperseBodyDistribution>,
/// The current phase of the simulation.
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<SitePairOverlapPenalty>,
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for PolydisperseHardDiskModel {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl PolydisperseHardDiskModel {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = PolydisperseHardDiskModel::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("polydisperse-hard-disk-model.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 polydisperse_hard_disk_model_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use polydisperse_hard_disk_model_interactive::main;With a uniform distribution, the maximum interaction range between any two sites is twice the maximum radius of any individual site.
Tip
If you use a normal distribution, you need to compute this maximum.
Warning
If you set the maximum interaction range too small, hoomd-rs will miss some interactions that should be computed. hoomd-rs has no idea what your custom
site_pair_energycomputes, so it cannot validate or provide an error when you choose amaximum_interaction_rangethat is incommensurate with your implementation.
Construct the Microstate
use anyhow::{Context, anyhow};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use rand::{
SeedableRng,
distr::{Distribution, Uniform},
rngs::StdRng,
};
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
univariate::{Expanded, OverlapPenalty, UnivariateEnergy},
};
use hoomd_mc::{
BodyDistribution, QuickCompress, QuickInsert, Sweep, Translate, Trial,
Tune, TuneOptions,
};
use hoomd_microstate::{
AppendMicrostate, Body, Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_utility::valid::PositiveReal;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's radius.
radius: PositiveReal,
}
impl Transform<SiteProperties> for BodyProperties {
fn transform(&self, site_properties: &SiteProperties) -> SiteProperties {
SiteProperties {
position: self.position + site_properties.position,
..*site_properties
}
}
}
#[derive(MaximumInteractionRange)]
struct SitePairInteraction {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
if r < a.radius.get() + b.radius.get() {
f64::INFINITY
} else {
0.0
}
}
fn site_pair_energy_initial(
&self,
_a: &SiteProperties,
_b: &SiteProperties,
) -> f64 {
0.0
}
fn is_only_infinite_or_zero() -> bool {
true
}
}
#[derive(MaximumInteractionRange)]
struct SitePairOverlapPenalty {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairOverlapPenalty {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
let pair_interaction = Expanded {
delta: a.radius.get() + b.radius.get(),
f: OverlapPenalty::default(),
};
pair_interaction.energy(r)
}
}
struct PolydisperseBodyDistribution {
/// Radius of each disk to insert into the microstate.
radii: Vec<PositiveReal>,
/// Simulation boundary.
boundary: Periodic<Rectangle>,
}
impl BodyDistribution<Body<BodyProperties, SiteProperties>>
for PolydisperseBodyDistribution
{
fn sample<R: rand::Rng + ?Sized>(
&self,
index: usize,
rng: &mut R,
) -> Body<BodyProperties, SiteProperties> {
let properties = Point {
position: self.boundary.sample(rng),
};
let sites = vec![SiteProperties {
position: Cartesian::default(),
radius: self.radii[index],
}];
Body { properties, sites }
}
}
impl<X> AppendMicrostate<BodyProperties, SiteProperties, X, Periodic<Rectangle>>
for HoomdGsdFile
{
#[inline]
fn append_microstate(
&mut self,
microstate: &Microstate<
BodyProperties,
SiteProperties,
X,
Periodic<Rectangle>,
>,
) -> Result<hoomd_gsd::hoomd::Frame<'_>, hoomd_gsd::hoomd::AppendError>
{
self.append_frame(microstate.step())?
.configuration_box(microstate.boundary().shape().to_gsd_box())?
.configuration_dimensions(Dimensions::Two)?
.particles_position(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.position)
.map(|p| [p[0], p[1], 0.0].into()),
)?
.particles_diameter(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.radius.get() * 2.0),
)
}
}
impl PolydisperseHardDiskModel {
/// Construct a new hard disk self-assembly simulation.
fn new() -> anyhow::Result<PolydisperseHardDiskModel> {
let seed = 1;
let minimum_radius = 0.1;
let maximum_radius = 0.8;
let initial_packing_fraction = 0.6;
let target_packing_fraction = 0.72;
let n_disks = 64_usize.pow(2);
let maximum_distance = 0.07;
let macrostate = Isothermal { temperature: 1.0 };
let radius_distribution = Uniform::new(minimum_radius, maximum_radius)?;
let mut rng = StdRng::seed_from_u64(seed);
let mut radii = Vec::with_capacity(n_disks);
for r in radius_distribution.sample_iter(&mut rng).take(n_disks) {
radii.push(r.try_into()?);
}
let total_particle_area = radii.iter().fold(0.0, |total, r| {
let circle = Circle { radius: *r };
total + circle.volume()
});
let hamiltonian = PairwiseCutoff(SitePairInteraction {
maximum_interaction_range: maximum_radius * 2.0,
});
let overlap_penalty_hamiltonian =
PairwiseCutoff(SitePairOverlapPenalty {
maximum_interaction_range: maximum_radius * 2.0,
});
let initial_box_volume = total_particle_area / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.seed(seed as u32)
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let distribution = PolydisperseBodyDistribution {
boundary: microstate.boundary().clone(),
radii,
};
let quick_insert = QuickInsert::new(distribution, n_disks);
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume = total_particle_area / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(PolydisperseHardDiskModel {
microstate,
hamiltonian,
translate_sweep,
quick_insert,
quick_compress,
overlap_penalty_hamiltonian,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct PolydisperseHardDiskModel {
/// 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<SitePairInteraction>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<PolydisperseBodyDistribution>,
/// The current phase of the simulation.
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<SitePairOverlapPenalty>,
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for PolydisperseHardDiskModel {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl PolydisperseHardDiskModel {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = PolydisperseHardDiskModel::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("polydisperse-hard-disk-model.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 polydisperse_hard_disk_model_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use polydisperse_hard_disk_model_interactive::main;Place Polydisperse Bodies with QuickInsert
Construct the custom PolydisperseBodyDistribution, pass it the predetermined
radii. Then make a new QuickInsert with this distribution:
use anyhow::{Context, anyhow};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use rand::{
SeedableRng,
distr::{Distribution, Uniform},
rngs::StdRng,
};
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
univariate::{Expanded, OverlapPenalty, UnivariateEnergy},
};
use hoomd_mc::{
BodyDistribution, QuickCompress, QuickInsert, Sweep, Translate, Trial,
Tune, TuneOptions,
};
use hoomd_microstate::{
AppendMicrostate, Body, Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_utility::valid::PositiveReal;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's radius.
radius: PositiveReal,
}
impl Transform<SiteProperties> for BodyProperties {
fn transform(&self, site_properties: &SiteProperties) -> SiteProperties {
SiteProperties {
position: self.position + site_properties.position,
..*site_properties
}
}
}
#[derive(MaximumInteractionRange)]
struct SitePairInteraction {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
if r < a.radius.get() + b.radius.get() {
f64::INFINITY
} else {
0.0
}
}
fn site_pair_energy_initial(
&self,
_a: &SiteProperties,
_b: &SiteProperties,
) -> f64 {
0.0
}
fn is_only_infinite_or_zero() -> bool {
true
}
}
#[derive(MaximumInteractionRange)]
struct SitePairOverlapPenalty {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairOverlapPenalty {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
let pair_interaction = Expanded {
delta: a.radius.get() + b.radius.get(),
f: OverlapPenalty::default(),
};
pair_interaction.energy(r)
}
}
struct PolydisperseBodyDistribution {
/// Radius of each disk to insert into the microstate.
radii: Vec<PositiveReal>,
/// Simulation boundary.
boundary: Periodic<Rectangle>,
}
impl BodyDistribution<Body<BodyProperties, SiteProperties>>
for PolydisperseBodyDistribution
{
fn sample<R: rand::Rng + ?Sized>(
&self,
index: usize,
rng: &mut R,
) -> Body<BodyProperties, SiteProperties> {
let properties = Point {
position: self.boundary.sample(rng),
};
let sites = vec![SiteProperties {
position: Cartesian::default(),
radius: self.radii[index],
}];
Body { properties, sites }
}
}
impl<X> AppendMicrostate<BodyProperties, SiteProperties, X, Periodic<Rectangle>>
for HoomdGsdFile
{
#[inline]
fn append_microstate(
&mut self,
microstate: &Microstate<
BodyProperties,
SiteProperties,
X,
Periodic<Rectangle>,
>,
) -> Result<hoomd_gsd::hoomd::Frame<'_>, hoomd_gsd::hoomd::AppendError>
{
self.append_frame(microstate.step())?
.configuration_box(microstate.boundary().shape().to_gsd_box())?
.configuration_dimensions(Dimensions::Two)?
.particles_position(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.position)
.map(|p| [p[0], p[1], 0.0].into()),
)?
.particles_diameter(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.radius.get() * 2.0),
)
}
}
impl PolydisperseHardDiskModel {
/// Construct a new hard disk self-assembly simulation.
fn new() -> anyhow::Result<PolydisperseHardDiskModel> {
let seed = 1;
let minimum_radius = 0.1;
let maximum_radius = 0.8;
let initial_packing_fraction = 0.6;
let target_packing_fraction = 0.72;
let n_disks = 64_usize.pow(2);
let maximum_distance = 0.07;
let macrostate = Isothermal { temperature: 1.0 };
let radius_distribution = Uniform::new(minimum_radius, maximum_radius)?;
let mut rng = StdRng::seed_from_u64(seed);
let mut radii = Vec::with_capacity(n_disks);
for r in radius_distribution.sample_iter(&mut rng).take(n_disks) {
radii.push(r.try_into()?);
}
let total_particle_area = radii.iter().fold(0.0, |total, r| {
let circle = Circle { radius: *r };
total + circle.volume()
});
let hamiltonian = PairwiseCutoff(SitePairInteraction {
maximum_interaction_range: maximum_radius * 2.0,
});
let overlap_penalty_hamiltonian =
PairwiseCutoff(SitePairOverlapPenalty {
maximum_interaction_range: maximum_radius * 2.0,
});
let initial_box_volume = total_particle_area / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.seed(seed as u32)
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let distribution = PolydisperseBodyDistribution {
boundary: microstate.boundary().clone(),
radii,
};
let quick_insert = QuickInsert::new(distribution, n_disks);
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume = total_particle_area / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(PolydisperseHardDiskModel {
microstate,
hamiltonian,
translate_sweep,
quick_insert,
quick_compress,
overlap_penalty_hamiltonian,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct PolydisperseHardDiskModel {
/// 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<SitePairInteraction>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<PolydisperseBodyDistribution>,
/// The current phase of the simulation.
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<SitePairOverlapPenalty>,
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for PolydisperseHardDiskModel {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl PolydisperseHardDiskModel {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = PolydisperseHardDiskModel::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("polydisperse-hard-disk-model.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 polydisperse_hard_disk_model_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use polydisperse_hard_disk_model_interactive::main;Initialization and Simulation
The remaining initialization and simulation code is very similar to that in the Hard Ellipse Self-Assembly tutorial. The only difference is that rotation moves are not present here (see the complete code below).
Execute the Simulation in Batch Mode
In a terminal, execute the following command to run the simulation in batch mode:
cargo run --release --example polydisperse-hard-disk-model
Visualize the Simulation Results
Open the generated polydisperse-hard-disk-model.gsd in Ovito or another
visualization tool to see the simulation results. By default. Ovito displays
the sites with the diameter set in the GSD file.
Render with Tachyon, and you should see something like:
Note
The diameters in the GSD file take precedence in OVITO. When they are set, there is no way to override the display radius. For example, you can no longer set the radius smaller to more easily see self-assembled crystal planes.
Conclusion
This tutorial showed you how to add a site radius field, how to use that when computing custom site-site interactions, and how to write each site’s diameter to a GSD file.
Navigate to the top of the page to see the simulation in action. Notice how
the sites have different sizes and are randomly distributed in the simulation
by QuickInsert.
Complete Code
use anyhow::{Context, anyhow};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use rand::{
SeedableRng,
distr::{Distribution, Uniform},
rngs::StdRng,
};
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
univariate::{Expanded, OverlapPenalty, UnivariateEnergy},
};
use hoomd_mc::{
BodyDistribution, QuickCompress, QuickInsert, Sweep, Translate, Trial,
Tune, TuneOptions,
};
use hoomd_microstate::{
AppendMicrostate, Body, Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_utility::valid::PositiveReal;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's radius.
radius: PositiveReal,
}
impl Transform<SiteProperties> for BodyProperties {
fn transform(&self, site_properties: &SiteProperties) -> SiteProperties {
SiteProperties {
position: self.position + site_properties.position,
..*site_properties
}
}
}
#[derive(MaximumInteractionRange)]
struct SitePairInteraction {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
if r < a.radius.get() + b.radius.get() {
f64::INFINITY
} else {
0.0
}
}
fn site_pair_energy_initial(
&self,
_a: &SiteProperties,
_b: &SiteProperties,
) -> f64 {
0.0
}
fn is_only_infinite_or_zero() -> bool {
true
}
}
#[derive(MaximumInteractionRange)]
struct SitePairOverlapPenalty {
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairOverlapPenalty {
fn site_pair_energy(&self, a: &SiteProperties, b: &SiteProperties) -> f64 {
let r = a.position().distance(b.position());
let pair_interaction = Expanded {
delta: a.radius.get() + b.radius.get(),
f: OverlapPenalty::default(),
};
pair_interaction.energy(r)
}
}
struct PolydisperseBodyDistribution {
/// Radius of each disk to insert into the microstate.
radii: Vec<PositiveReal>,
/// Simulation boundary.
boundary: Periodic<Rectangle>,
}
impl BodyDistribution<Body<BodyProperties, SiteProperties>>
for PolydisperseBodyDistribution
{
fn sample<R: rand::Rng + ?Sized>(
&self,
index: usize,
rng: &mut R,
) -> Body<BodyProperties, SiteProperties> {
let properties = Point {
position: self.boundary.sample(rng),
};
let sites = vec![SiteProperties {
position: Cartesian::default(),
radius: self.radii[index],
}];
Body { properties, sites }
}
}
impl<X> AppendMicrostate<BodyProperties, SiteProperties, X, Periodic<Rectangle>>
for HoomdGsdFile
{
#[inline]
fn append_microstate(
&mut self,
microstate: &Microstate<
BodyProperties,
SiteProperties,
X,
Periodic<Rectangle>,
>,
) -> Result<hoomd_gsd::hoomd::Frame<'_>, hoomd_gsd::hoomd::AppendError>
{
self.append_frame(microstate.step())?
.configuration_box(microstate.boundary().shape().to_gsd_box())?
.configuration_dimensions(Dimensions::Two)?
.particles_position(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.position)
.map(|p| [p[0], p[1], 0.0].into()),
)?
.particles_diameter(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.radius.get() * 2.0),
)
}
}
impl PolydisperseHardDiskModel {
/// Construct a new hard disk self-assembly simulation.
fn new() -> anyhow::Result<PolydisperseHardDiskModel> {
let seed = 1;
let minimum_radius = 0.1;
let maximum_radius = 0.8;
let initial_packing_fraction = 0.6;
let target_packing_fraction = 0.72;
let n_disks = 64_usize.pow(2);
let maximum_distance = 0.07;
let macrostate = Isothermal { temperature: 1.0 };
let radius_distribution = Uniform::new(minimum_radius, maximum_radius)?;
let mut rng = StdRng::seed_from_u64(seed);
let mut radii = Vec::with_capacity(n_disks);
for r in radius_distribution.sample_iter(&mut rng).take(n_disks) {
radii.push(r.try_into()?);
}
let total_particle_area = radii.iter().fold(0.0, |total, r| {
let circle = Circle { radius: *r };
total + circle.volume()
});
let hamiltonian = PairwiseCutoff(SitePairInteraction {
maximum_interaction_range: maximum_radius * 2.0,
});
let overlap_penalty_hamiltonian =
PairwiseCutoff(SitePairOverlapPenalty {
maximum_interaction_range: maximum_radius * 2.0,
});
let initial_box_volume = total_particle_area / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.seed(seed as u32)
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let distribution = PolydisperseBodyDistribution {
boundary: microstate.boundary().clone(),
radii,
};
let quick_insert = QuickInsert::new(distribution, n_disks);
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume = total_particle_area / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(PolydisperseHardDiskModel {
microstate,
hamiltonian,
translate_sweep,
quick_insert,
quick_compress,
overlap_penalty_hamiltonian,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct PolydisperseHardDiskModel {
/// 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<SitePairInteraction>,
/// Trial moves to apply.
translate_sweep: Sweep<Translate<PositionVector>>,
/// Temperature set point.
macrostate: Isothermal,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// Quick insert algorithm.
quick_insert: QuickInsert<PolydisperseBodyDistribution>,
/// The current phase of the simulation.
/// How sites interact when inserted and compressed.
overlap_penalty_hamiltonian: PairwiseCutoff<SitePairOverlapPenalty>,
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for PolydisperseHardDiskModel {
/// Advance the simulation forward one step.
fn advance(&mut self) -> anyhow::Result<()> {
match self.phase {
Phase::Initialize => {
self.initialize().context("failed to initialize")?
}
Phase::Equilibrate => self.equilibrate(),
}
self.microstate.increment_step();
Ok(())
}
/// Get the current simulation step.
fn step(&self) -> u64 {
self.microstate.step()
}
}
impl PolydisperseHardDiskModel {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert.is_complete() {
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
}
self.translate_sweep.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
&Isothermal { temperature: 1.0 },
);
if self.quick_compress.is_complete() {
self.translate_sweep.tune_with_options(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
&TuneOptions::default(),
);
self.phase = Phase::Equilibrate;
println!(
"Initialization complete at step {}.",
self.microstate.step()
);
}
if self.step() >= 20_000 {
let n = self.microstate.bodies().len();
let target_n = self.quick_insert.target();
let volume = self.microstate.boundary().volume();
let target_volume = self.quick_compress.target_volume();
return Err(anyhow!(
"inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}"
));
}
Ok(())
}
fn equilibrate(&mut self) {
self.translate_sweep.apply(
&mut self.microstate,
&self.hamiltonian,
&self.macrostate,
);
}
}
// 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 = PolydisperseHardDiskModel::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("polydisperse-hard-disk-model.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 polydisperse_hard_disk_model_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use polydisperse_hard_disk_model_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.
