Type-dependent Interactions
Click to focus. When focused, press m to open the options menu.
Refresh the page to restart the simulation.
Overview
Some models label each site with a site type and apply interactions as a function of the site type. For example, you can model coarse-grained phase separation with attractive A-A and B-B interactions and purely repulsive A-B interactions. This tutorial shows you how to define site types and compute type-dependent pairwise interactions.
- Objectives:
- Show how to use an
enumto name the possible site types. - Define a custom site properties struct that includes the site type.
- Show how to compute type-dependent pairwise interactions.
- Demonstrate phase separation of A and B site types.
- Show how to use an
- File:
hoomd-rs/examples/mc-tutorial/type-dependent-interactions.rs - Run (interactively):
cargo run --release --features "bevy" --example type-dependent-interactions - Run (in batch mode):
cargo run --release --example type-dependent-interactions
Use Declarations
use anyhow::{Context, anyhow};
use strum::VariantNames;
use strum_macros::VariantNames;
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
pairwise::Isotropic,
univariate::{
Expanded, LennardJones, OverlapPenalty, UnivariateEnergy,
WeeksChandlerAnderson,
},
};
use hoomd_mc::{
QuickCompress, QuickInsert, Sweep, Translate, Trial, Tune, UniformIn,
};
use hoomd_microstate::{AppendMicrostate, Site};
use hoomd_microstate::{
Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, PartialEq, VariantNames)]
enum SiteType {
#[default]
A,
B,
}
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's type.
site_type: SiteType,
}
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 {
lj_aa: LennardJones<12, 6>,
wca_ab: WeeksChandlerAnderson,
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(
&self,
site_properties_i: &SiteProperties,
site_properties_j: &SiteProperties,
) -> f64 {
let r = site_properties_i
.position
.distance(&site_properties_j.position);
match (site_properties_i.site_type, site_properties_j.site_type) {
(SiteType::A, SiteType::A) => self.lj_aa.energy(r),
(SiteType::A, SiteType::B) | (SiteType::B, SiteType::A) => {
self.wca_ab.energy(r)
}
(SiteType::B, SiteType::B) => {
1.0 / r.powi(12) - f64::exp(-1.0 / 2.0 * r.powi(2))
}
}
}
}
impl TypeDependentInteractions {
/// Construct a new type-dependent interactions simulation.
fn new() -> anyhow::Result<TypeDependentInteractions> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.5;
let n_disks = 512;
let maximum_distance = 0.07;
let sigma = 1.0;
let macrostate = Isothermal { temperature: 1.0 };
let lj_aa = LennardJones {
epsilon: 2.0,
sigma: 1.0,
};
let wca_ab = WeeksChandlerAnderson {
epsilon: 1.0,
sigma: 1.0,
};
let hamiltonian = PairwiseCutoff(SitePairInteraction {
lj_aa,
wca_ab,
maximum_interaction_range: 2.5,
});
let overlap_penalty = Isotropic {
interaction: Expanded {
delta: sigma,
f: OverlapPenalty::default(),
},
r_cut: sigma,
};
let overlap_penalty_hamiltonian = PairwiseCutoff(overlap_penalty);
let circle = Circle {
radius: (sigma / 2.0).try_into()?,
};
let initial_box_volume =
n_disks as f64 * circle.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::A,
..SiteProperties::default()
}],
};
let quick_insert_a = QuickInsert::new(distribution, n_disks / 2);
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::B,
..SiteProperties::default()
}],
};
let quick_insert_b =
QuickInsert::new(distribution, n_disks - quick_insert_a.target());
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume =
n_disks as f64 * circle.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(TypeDependentInteractions {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
quick_insert_a,
quick_insert_b,
quick_compress,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct TypeDependentInteractions {
/// 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 insert algorithm for A bodies.
quick_insert_a: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick insert algorithm for B bodies.
quick_insert_b: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// How sites interact during compression.
overlap_penalty_hamiltonian:
PairwiseCutoff<Isotropic<Expanded<OverlapPenalty>>>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for TypeDependentInteractions {
/// 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 TypeDependentInteractions {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert_a.is_complete()
&& self.quick_insert_b.is_complete()
{
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert_a
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
self.quick_insert_b
.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_default(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
);
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_a.target() + self.quick_insert_b.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,
);
}
}
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_type_id(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.site_type as u32),
)?
.particles_types(SiteType::VARIANTS.iter().copied())
}
}
// 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_interaction::TotalEnergy;
let mut simulation = TypeDependentInteractions::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("type-dependent-interactions.gsd")?;
for _ in 0..100_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.log_scalar(
"potential_energy",
simulation.hamiltonian.total_energy(&simulation.microstate),
)?
.log_scalars(
"particles/site_energy",
simulation
.microstate
.iter_sites_tag_order()
.map(|s| site_energy(&simulation, s)),
)?
.end()?;
}
}
Ok(())
}
#[allow(dead_code, reason = "site_energy is used in the batch mode main()")]
fn site_energy(
simulation: &TypeDependentInteractions,
site: &Site<SiteProperties>,
) -> f64 {
use hoomd_interaction::DeltaEnergyRemove;
let body_index = simulation.microstate.body_indices()[site.body_tag]
.expect("site's parent body should be in the microstate");
-simulation
.hamiltonian
.delta_energy_remove(&simulation.microstate, body_index)
}
#[cfg(feature = "bevy")]
mod type_dependent_interactions_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use type_dependent_interactions_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 strum::VariantNames;
use strum_macros::VariantNames;
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
pairwise::Isotropic,
univariate::{
Expanded, LennardJones, OverlapPenalty, UnivariateEnergy,
WeeksChandlerAnderson,
},
};
use hoomd_mc::{
QuickCompress, QuickInsert, Sweep, Translate, Trial, Tune, UniformIn,
};
use hoomd_microstate::{AppendMicrostate, Site};
use hoomd_microstate::{
Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, PartialEq, VariantNames)]
enum SiteType {
#[default]
A,
B,
}
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's type.
site_type: SiteType,
}
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 {
lj_aa: LennardJones<12, 6>,
wca_ab: WeeksChandlerAnderson,
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(
&self,
site_properties_i: &SiteProperties,
site_properties_j: &SiteProperties,
) -> f64 {
let r = site_properties_i
.position
.distance(&site_properties_j.position);
match (site_properties_i.site_type, site_properties_j.site_type) {
(SiteType::A, SiteType::A) => self.lj_aa.energy(r),
(SiteType::A, SiteType::B) | (SiteType::B, SiteType::A) => {
self.wca_ab.energy(r)
}
(SiteType::B, SiteType::B) => {
1.0 / r.powi(12) - f64::exp(-1.0 / 2.0 * r.powi(2))
}
}
}
}
impl TypeDependentInteractions {
/// Construct a new type-dependent interactions simulation.
fn new() -> anyhow::Result<TypeDependentInteractions> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.5;
let n_disks = 512;
let maximum_distance = 0.07;
let sigma = 1.0;
let macrostate = Isothermal { temperature: 1.0 };
let lj_aa = LennardJones {
epsilon: 2.0,
sigma: 1.0,
};
let wca_ab = WeeksChandlerAnderson {
epsilon: 1.0,
sigma: 1.0,
};
let hamiltonian = PairwiseCutoff(SitePairInteraction {
lj_aa,
wca_ab,
maximum_interaction_range: 2.5,
});
let overlap_penalty = Isotropic {
interaction: Expanded {
delta: sigma,
f: OverlapPenalty::default(),
},
r_cut: sigma,
};
let overlap_penalty_hamiltonian = PairwiseCutoff(overlap_penalty);
let circle = Circle {
radius: (sigma / 2.0).try_into()?,
};
let initial_box_volume =
n_disks as f64 * circle.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::A,
..SiteProperties::default()
}],
};
let quick_insert_a = QuickInsert::new(distribution, n_disks / 2);
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::B,
..SiteProperties::default()
}],
};
let quick_insert_b =
QuickInsert::new(distribution, n_disks - quick_insert_a.target());
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume =
n_disks as f64 * circle.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(TypeDependentInteractions {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
quick_insert_a,
quick_insert_b,
quick_compress,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct TypeDependentInteractions {
/// 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 insert algorithm for A bodies.
quick_insert_a: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick insert algorithm for B bodies.
quick_insert_b: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// How sites interact during compression.
overlap_penalty_hamiltonian:
PairwiseCutoff<Isotropic<Expanded<OverlapPenalty>>>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for TypeDependentInteractions {
/// 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 TypeDependentInteractions {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert_a.is_complete()
&& self.quick_insert_b.is_complete()
{
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert_a
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
self.quick_insert_b
.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_default(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
);
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_a.target() + self.quick_insert_b.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,
);
}
}
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_type_id(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.site_type as u32),
)?
.particles_types(SiteType::VARIANTS.iter().copied())
}
}
// 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_interaction::TotalEnergy;
let mut simulation = TypeDependentInteractions::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("type-dependent-interactions.gsd")?;
for _ in 0..100_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.log_scalar(
"potential_energy",
simulation.hamiltonian.total_energy(&simulation.microstate),
)?
.log_scalars(
"particles/site_energy",
simulation
.microstate
.iter_sites_tag_order()
.map(|s| site_energy(&simulation, s)),
)?
.end()?;
}
}
Ok(())
}
#[allow(dead_code, reason = "site_energy is used in the batch mode main()")]
fn site_energy(
simulation: &TypeDependentInteractions,
site: &Site<SiteProperties>,
) -> f64 {
use hoomd_interaction::DeltaEnergyRemove;
let body_index = simulation.microstate.body_indices()[site.body_tag]
.expect("site's parent body should be in the microstate");
-simulation
.hamiltonian
.delta_energy_remove(&simulation.microstate, body_index)
}
#[cfg(feature = "bevy")]
mod type_dependent_interactions_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use type_dependent_interactions_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
You might be familiar with simulation tools where you assign types to sites based on a numerical index or string name. In hoomd-rs, you can assign the site type (and any other site-specific parameter(s)) using any Rust datatype.
SiteType enum
Using an enum ensures that every site always has a well-defined site type.
If you fail to assign a type or set an invalid one, Rust will issue an error at
compile time. In this example, SiteType enumerates all the site types:
use anyhow::{Context, anyhow};
use strum::VariantNames;
use strum_macros::VariantNames;
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
pairwise::Isotropic,
univariate::{
Expanded, LennardJones, OverlapPenalty, UnivariateEnergy,
WeeksChandlerAnderson,
},
};
use hoomd_mc::{
QuickCompress, QuickInsert, Sweep, Translate, Trial, Tune, UniformIn,
};
use hoomd_microstate::{AppendMicrostate, Site};
use hoomd_microstate::{
Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, PartialEq, VariantNames)]
enum SiteType {
#[default]
A,
B,
}
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's type.
site_type: SiteType,
}
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 {
lj_aa: LennardJones<12, 6>,
wca_ab: WeeksChandlerAnderson,
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(
&self,
site_properties_i: &SiteProperties,
site_properties_j: &SiteProperties,
) -> f64 {
let r = site_properties_i
.position
.distance(&site_properties_j.position);
match (site_properties_i.site_type, site_properties_j.site_type) {
(SiteType::A, SiteType::A) => self.lj_aa.energy(r),
(SiteType::A, SiteType::B) | (SiteType::B, SiteType::A) => {
self.wca_ab.energy(r)
}
(SiteType::B, SiteType::B) => {
1.0 / r.powi(12) - f64::exp(-1.0 / 2.0 * r.powi(2))
}
}
}
}
impl TypeDependentInteractions {
/// Construct a new type-dependent interactions simulation.
fn new() -> anyhow::Result<TypeDependentInteractions> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.5;
let n_disks = 512;
let maximum_distance = 0.07;
let sigma = 1.0;
let macrostate = Isothermal { temperature: 1.0 };
let lj_aa = LennardJones {
epsilon: 2.0,
sigma: 1.0,
};
let wca_ab = WeeksChandlerAnderson {
epsilon: 1.0,
sigma: 1.0,
};
let hamiltonian = PairwiseCutoff(SitePairInteraction {
lj_aa,
wca_ab,
maximum_interaction_range: 2.5,
});
let overlap_penalty = Isotropic {
interaction: Expanded {
delta: sigma,
f: OverlapPenalty::default(),
},
r_cut: sigma,
};
let overlap_penalty_hamiltonian = PairwiseCutoff(overlap_penalty);
let circle = Circle {
radius: (sigma / 2.0).try_into()?,
};
let initial_box_volume =
n_disks as f64 * circle.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::A,
..SiteProperties::default()
}],
};
let quick_insert_a = QuickInsert::new(distribution, n_disks / 2);
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::B,
..SiteProperties::default()
}],
};
let quick_insert_b =
QuickInsert::new(distribution, n_disks - quick_insert_a.target());
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume =
n_disks as f64 * circle.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(TypeDependentInteractions {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
quick_insert_a,
quick_insert_b,
quick_compress,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct TypeDependentInteractions {
/// 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 insert algorithm for A bodies.
quick_insert_a: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick insert algorithm for B bodies.
quick_insert_b: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// How sites interact during compression.
overlap_penalty_hamiltonian:
PairwiseCutoff<Isotropic<Expanded<OverlapPenalty>>>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for TypeDependentInteractions {
/// 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 TypeDependentInteractions {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert_a.is_complete()
&& self.quick_insert_b.is_complete()
{
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert_a
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
self.quick_insert_b
.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_default(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
);
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_a.target() + self.quick_insert_b.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,
);
}
}
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_type_id(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.site_type as u32),
)?
.particles_types(SiteType::VARIANTS.iter().copied())
}
}
// 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_interaction::TotalEnergy;
let mut simulation = TypeDependentInteractions::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("type-dependent-interactions.gsd")?;
for _ in 0..100_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.log_scalar(
"potential_energy",
simulation.hamiltonian.total_energy(&simulation.microstate),
)?
.log_scalars(
"particles/site_energy",
simulation
.microstate
.iter_sites_tag_order()
.map(|s| site_energy(&simulation, s)),
)?
.end()?;
}
}
Ok(())
}
#[allow(dead_code, reason = "site_energy is used in the batch mode main()")]
fn site_energy(
simulation: &TypeDependentInteractions,
site: &Site<SiteProperties>,
) -> f64 {
use hoomd_interaction::DeltaEnergyRemove;
let body_index = simulation.microstate.body_indices()[site.body_tag]
.expect("site's parent body should be in the microstate");
-simulation
.hamiltonian
.delta_energy_remove(&simulation.microstate, body_index)
}
#[cfg(feature = "bevy")]
mod type_dependent_interactions_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use type_dependent_interactions_interactive::main;Define SiteProperties
Previous tutorials used one of the built-in structures (Point or
OrientedPoint) to represent the site properties. These types are limited
as they represent only a site’s position (or position and orientation). This
tutorial defines a new SiteProperties struct that gives each site a position
in space and a given site type:
use anyhow::{Context, anyhow};
use strum::VariantNames;
use strum_macros::VariantNames;
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
pairwise::Isotropic,
univariate::{
Expanded, LennardJones, OverlapPenalty, UnivariateEnergy,
WeeksChandlerAnderson,
},
};
use hoomd_mc::{
QuickCompress, QuickInsert, Sweep, Translate, Trial, Tune, UniformIn,
};
use hoomd_microstate::{AppendMicrostate, Site};
use hoomd_microstate::{
Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, PartialEq, VariantNames)]
enum SiteType {
#[default]
A,
B,
}
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's type.
site_type: SiteType,
}
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 {
lj_aa: LennardJones<12, 6>,
wca_ab: WeeksChandlerAnderson,
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(
&self,
site_properties_i: &SiteProperties,
site_properties_j: &SiteProperties,
) -> f64 {
let r = site_properties_i
.position
.distance(&site_properties_j.position);
match (site_properties_i.site_type, site_properties_j.site_type) {
(SiteType::A, SiteType::A) => self.lj_aa.energy(r),
(SiteType::A, SiteType::B) | (SiteType::B, SiteType::A) => {
self.wca_ab.energy(r)
}
(SiteType::B, SiteType::B) => {
1.0 / r.powi(12) - f64::exp(-1.0 / 2.0 * r.powi(2))
}
}
}
}
impl TypeDependentInteractions {
/// Construct a new type-dependent interactions simulation.
fn new() -> anyhow::Result<TypeDependentInteractions> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.5;
let n_disks = 512;
let maximum_distance = 0.07;
let sigma = 1.0;
let macrostate = Isothermal { temperature: 1.0 };
let lj_aa = LennardJones {
epsilon: 2.0,
sigma: 1.0,
};
let wca_ab = WeeksChandlerAnderson {
epsilon: 1.0,
sigma: 1.0,
};
let hamiltonian = PairwiseCutoff(SitePairInteraction {
lj_aa,
wca_ab,
maximum_interaction_range: 2.5,
});
let overlap_penalty = Isotropic {
interaction: Expanded {
delta: sigma,
f: OverlapPenalty::default(),
},
r_cut: sigma,
};
let overlap_penalty_hamiltonian = PairwiseCutoff(overlap_penalty);
let circle = Circle {
radius: (sigma / 2.0).try_into()?,
};
let initial_box_volume =
n_disks as f64 * circle.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::A,
..SiteProperties::default()
}],
};
let quick_insert_a = QuickInsert::new(distribution, n_disks / 2);
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::B,
..SiteProperties::default()
}],
};
let quick_insert_b =
QuickInsert::new(distribution, n_disks - quick_insert_a.target());
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume =
n_disks as f64 * circle.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(TypeDependentInteractions {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
quick_insert_a,
quick_insert_b,
quick_compress,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct TypeDependentInteractions {
/// 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 insert algorithm for A bodies.
quick_insert_a: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick insert algorithm for B bodies.
quick_insert_b: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// How sites interact during compression.
overlap_penalty_hamiltonian:
PairwiseCutoff<Isotropic<Expanded<OverlapPenalty>>>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for TypeDependentInteractions {
/// 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 TypeDependentInteractions {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert_a.is_complete()
&& self.quick_insert_b.is_complete()
{
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert_a
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
self.quick_insert_b
.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_default(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
);
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_a.target() + self.quick_insert_b.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,
);
}
}
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_type_id(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.site_type as u32),
)?
.particles_types(SiteType::VARIANTS.iter().copied())
}
}
// 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_interaction::TotalEnergy;
let mut simulation = TypeDependentInteractions::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("type-dependent-interactions.gsd")?;
for _ in 0..100_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.log_scalar(
"potential_energy",
simulation.hamiltonian.total_energy(&simulation.microstate),
)?
.log_scalars(
"particles/site_energy",
simulation
.microstate
.iter_sites_tag_order()
.map(|s| site_energy(&simulation, s)),
)?
.end()?;
}
}
Ok(())
}
#[allow(dead_code, reason = "site_energy is used in the batch mode main()")]
fn site_energy(
simulation: &TypeDependentInteractions,
site: &Site<SiteProperties>,
) -> f64 {
use hoomd_interaction::DeltaEnergyRemove;
let body_index = simulation.microstate.body_indices()[site.body_tag]
.expect("site's parent body should be in the microstate");
-simulation
.hamiltonian
.delta_energy_remove(&simulation.microstate, body_index)
}
#[cfg(feature = "bevy")]
mod type_dependent_interactions_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use type_dependent_interactions_interactive::main;Initialization, MC, and MD methods operate on a generic site properties type
S with trait bounds set as needed. The derive macro implements the
listed traits for SiteProperties automatically. All the types listed in the
above code block are required by at least one method in this example. Rust
provides the Clone, Copy, Default traits (and their derive macros).
hoomd_microstate defines Position and Orientation along with their
corresponding derive macros.
Note
You can add any fields to your
SitePropertiestype and use those fields when computing interactions.positionis the only required field.
Transforming Sites
Body stores each of its sites in the body frame of reference. The
Transform trait (implemented for the body properties type) transforms
site properties from the body frame to the system frame. You must implement
Transform<SiteProperties> for the body properties type to use Microstate:
use anyhow::{Context, anyhow};
use strum::VariantNames;
use strum_macros::VariantNames;
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
pairwise::Isotropic,
univariate::{
Expanded, LennardJones, OverlapPenalty, UnivariateEnergy,
WeeksChandlerAnderson,
},
};
use hoomd_mc::{
QuickCompress, QuickInsert, Sweep, Translate, Trial, Tune, UniformIn,
};
use hoomd_microstate::{AppendMicrostate, Site};
use hoomd_microstate::{
Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, PartialEq, VariantNames)]
enum SiteType {
#[default]
A,
B,
}
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's type.
site_type: SiteType,
}
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 {
lj_aa: LennardJones<12, 6>,
wca_ab: WeeksChandlerAnderson,
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(
&self,
site_properties_i: &SiteProperties,
site_properties_j: &SiteProperties,
) -> f64 {
let r = site_properties_i
.position
.distance(&site_properties_j.position);
match (site_properties_i.site_type, site_properties_j.site_type) {
(SiteType::A, SiteType::A) => self.lj_aa.energy(r),
(SiteType::A, SiteType::B) | (SiteType::B, SiteType::A) => {
self.wca_ab.energy(r)
}
(SiteType::B, SiteType::B) => {
1.0 / r.powi(12) - f64::exp(-1.0 / 2.0 * r.powi(2))
}
}
}
}
impl TypeDependentInteractions {
/// Construct a new type-dependent interactions simulation.
fn new() -> anyhow::Result<TypeDependentInteractions> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.5;
let n_disks = 512;
let maximum_distance = 0.07;
let sigma = 1.0;
let macrostate = Isothermal { temperature: 1.0 };
let lj_aa = LennardJones {
epsilon: 2.0,
sigma: 1.0,
};
let wca_ab = WeeksChandlerAnderson {
epsilon: 1.0,
sigma: 1.0,
};
let hamiltonian = PairwiseCutoff(SitePairInteraction {
lj_aa,
wca_ab,
maximum_interaction_range: 2.5,
});
let overlap_penalty = Isotropic {
interaction: Expanded {
delta: sigma,
f: OverlapPenalty::default(),
},
r_cut: sigma,
};
let overlap_penalty_hamiltonian = PairwiseCutoff(overlap_penalty);
let circle = Circle {
radius: (sigma / 2.0).try_into()?,
};
let initial_box_volume =
n_disks as f64 * circle.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::A,
..SiteProperties::default()
}],
};
let quick_insert_a = QuickInsert::new(distribution, n_disks / 2);
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::B,
..SiteProperties::default()
}],
};
let quick_insert_b =
QuickInsert::new(distribution, n_disks - quick_insert_a.target());
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume =
n_disks as f64 * circle.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(TypeDependentInteractions {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
quick_insert_a,
quick_insert_b,
quick_compress,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct TypeDependentInteractions {
/// 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 insert algorithm for A bodies.
quick_insert_a: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick insert algorithm for B bodies.
quick_insert_b: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// How sites interact during compression.
overlap_penalty_hamiltonian:
PairwiseCutoff<Isotropic<Expanded<OverlapPenalty>>>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for TypeDependentInteractions {
/// 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 TypeDependentInteractions {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert_a.is_complete()
&& self.quick_insert_b.is_complete()
{
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert_a
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
self.quick_insert_b
.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_default(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
);
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_a.target() + self.quick_insert_b.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,
);
}
}
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_type_id(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.site_type as u32),
)?
.particles_types(SiteType::VARIANTS.iter().copied())
}
}
// 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_interaction::TotalEnergy;
let mut simulation = TypeDependentInteractions::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("type-dependent-interactions.gsd")?;
for _ in 0..100_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.log_scalar(
"potential_energy",
simulation.hamiltonian.total_energy(&simulation.microstate),
)?
.log_scalars(
"particles/site_energy",
simulation
.microstate
.iter_sites_tag_order()
.map(|s| site_energy(&simulation, s)),
)?
.end()?;
}
}
Ok(())
}
#[allow(dead_code, reason = "site_energy is used in the batch mode main()")]
fn site_energy(
simulation: &TypeDependentInteractions,
site: &Site<SiteProperties>,
) -> f64 {
use hoomd_interaction::DeltaEnergyRemove;
let body_index = simulation.microstate.body_indices()[site.body_tag]
.expect("site's parent body should be in the microstate");
-simulation
.hamiltonian
.delta_energy_remove(&simulation.microstate, body_index)
}
#[cfg(feature = "bevy")]
mod type_dependent_interactions_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use type_dependent_interactions_interactive::main;This implementation is suitable for point bodies with point sites in Euclidean
space as it transforms from one frame to the other by vector addition and copies
all other fields. See the hoomd_microstate::property module documentation for
code that works with oriented bodies and/or sites.
Site-site Interactions
For demonstration purposes, this tutorial shows you how to implement a coarse-grained model where A-A interactions attract via the Lennard-Jones potential, B-B interact via the sum of a power law and a Gaussian, and A-B interact with the Weeks-Chandler-Anderson potential.
Define SitePairInteraction
Define a new struct to hold the interaction parameters. Use the provided
types for the LennardJones and WeeksChandlerAnderson interactions:
use anyhow::{Context, anyhow};
use strum::VariantNames;
use strum_macros::VariantNames;
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
pairwise::Isotropic,
univariate::{
Expanded, LennardJones, OverlapPenalty, UnivariateEnergy,
WeeksChandlerAnderson,
},
};
use hoomd_mc::{
QuickCompress, QuickInsert, Sweep, Translate, Trial, Tune, UniformIn,
};
use hoomd_microstate::{AppendMicrostate, Site};
use hoomd_microstate::{
Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, PartialEq, VariantNames)]
enum SiteType {
#[default]
A,
B,
}
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's type.
site_type: SiteType,
}
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 {
lj_aa: LennardJones<12, 6>,
wca_ab: WeeksChandlerAnderson,
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(
&self,
site_properties_i: &SiteProperties,
site_properties_j: &SiteProperties,
) -> f64 {
let r = site_properties_i
.position
.distance(&site_properties_j.position);
match (site_properties_i.site_type, site_properties_j.site_type) {
(SiteType::A, SiteType::A) => self.lj_aa.energy(r),
(SiteType::A, SiteType::B) | (SiteType::B, SiteType::A) => {
self.wca_ab.energy(r)
}
(SiteType::B, SiteType::B) => {
1.0 / r.powi(12) - f64::exp(-1.0 / 2.0 * r.powi(2))
}
}
}
}
impl TypeDependentInteractions {
/// Construct a new type-dependent interactions simulation.
fn new() -> anyhow::Result<TypeDependentInteractions> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.5;
let n_disks = 512;
let maximum_distance = 0.07;
let sigma = 1.0;
let macrostate = Isothermal { temperature: 1.0 };
let lj_aa = LennardJones {
epsilon: 2.0,
sigma: 1.0,
};
let wca_ab = WeeksChandlerAnderson {
epsilon: 1.0,
sigma: 1.0,
};
let hamiltonian = PairwiseCutoff(SitePairInteraction {
lj_aa,
wca_ab,
maximum_interaction_range: 2.5,
});
let overlap_penalty = Isotropic {
interaction: Expanded {
delta: sigma,
f: OverlapPenalty::default(),
},
r_cut: sigma,
};
let overlap_penalty_hamiltonian = PairwiseCutoff(overlap_penalty);
let circle = Circle {
radius: (sigma / 2.0).try_into()?,
};
let initial_box_volume =
n_disks as f64 * circle.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::A,
..SiteProperties::default()
}],
};
let quick_insert_a = QuickInsert::new(distribution, n_disks / 2);
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::B,
..SiteProperties::default()
}],
};
let quick_insert_b =
QuickInsert::new(distribution, n_disks - quick_insert_a.target());
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume =
n_disks as f64 * circle.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(TypeDependentInteractions {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
quick_insert_a,
quick_insert_b,
quick_compress,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct TypeDependentInteractions {
/// 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 insert algorithm for A bodies.
quick_insert_a: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick insert algorithm for B bodies.
quick_insert_b: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// How sites interact during compression.
overlap_penalty_hamiltonian:
PairwiseCutoff<Isotropic<Expanded<OverlapPenalty>>>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for TypeDependentInteractions {
/// 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 TypeDependentInteractions {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert_a.is_complete()
&& self.quick_insert_b.is_complete()
{
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert_a
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
self.quick_insert_b
.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_default(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
);
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_a.target() + self.quick_insert_b.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,
);
}
}
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_type_id(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.site_type as u32),
)?
.particles_types(SiteType::VARIANTS.iter().copied())
}
}
// 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_interaction::TotalEnergy;
let mut simulation = TypeDependentInteractions::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("type-dependent-interactions.gsd")?;
for _ in 0..100_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.log_scalar(
"potential_energy",
simulation.hamiltonian.total_energy(&simulation.microstate),
)?
.log_scalars(
"particles/site_energy",
simulation
.microstate
.iter_sites_tag_order()
.map(|s| site_energy(&simulation, s)),
)?
.end()?;
}
}
Ok(())
}
#[allow(dead_code, reason = "site_energy is used in the batch mode main()")]
fn site_energy(
simulation: &TypeDependentInteractions,
site: &Site<SiteProperties>,
) -> f64 {
use hoomd_interaction::DeltaEnergyRemove;
let body_index = simulation.microstate.body_indices()[site.body_tag]
.expect("site's parent body should be in the microstate");
-simulation
.hamiltonian
.delta_energy_remove(&simulation.microstate, body_index)
}
#[cfg(feature = "bevy")]
mod type_dependent_interactions_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use type_dependent_interactions_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 type-dependent interactions:
use anyhow::{Context, anyhow};
use strum::VariantNames;
use strum_macros::VariantNames;
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
pairwise::Isotropic,
univariate::{
Expanded, LennardJones, OverlapPenalty, UnivariateEnergy,
WeeksChandlerAnderson,
},
};
use hoomd_mc::{
QuickCompress, QuickInsert, Sweep, Translate, Trial, Tune, UniformIn,
};
use hoomd_microstate::{AppendMicrostate, Site};
use hoomd_microstate::{
Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, PartialEq, VariantNames)]
enum SiteType {
#[default]
A,
B,
}
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's type.
site_type: SiteType,
}
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 {
lj_aa: LennardJones<12, 6>,
wca_ab: WeeksChandlerAnderson,
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(
&self,
site_properties_i: &SiteProperties,
site_properties_j: &SiteProperties,
) -> f64 {
let r = site_properties_i
.position
.distance(&site_properties_j.position);
match (site_properties_i.site_type, site_properties_j.site_type) {
(SiteType::A, SiteType::A) => self.lj_aa.energy(r),
(SiteType::A, SiteType::B) | (SiteType::B, SiteType::A) => {
self.wca_ab.energy(r)
}
(SiteType::B, SiteType::B) => {
1.0 / r.powi(12) - f64::exp(-1.0 / 2.0 * r.powi(2))
}
}
}
}
impl TypeDependentInteractions {
/// Construct a new type-dependent interactions simulation.
fn new() -> anyhow::Result<TypeDependentInteractions> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.5;
let n_disks = 512;
let maximum_distance = 0.07;
let sigma = 1.0;
let macrostate = Isothermal { temperature: 1.0 };
let lj_aa = LennardJones {
epsilon: 2.0,
sigma: 1.0,
};
let wca_ab = WeeksChandlerAnderson {
epsilon: 1.0,
sigma: 1.0,
};
let hamiltonian = PairwiseCutoff(SitePairInteraction {
lj_aa,
wca_ab,
maximum_interaction_range: 2.5,
});
let overlap_penalty = Isotropic {
interaction: Expanded {
delta: sigma,
f: OverlapPenalty::default(),
},
r_cut: sigma,
};
let overlap_penalty_hamiltonian = PairwiseCutoff(overlap_penalty);
let circle = Circle {
radius: (sigma / 2.0).try_into()?,
};
let initial_box_volume =
n_disks as f64 * circle.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::A,
..SiteProperties::default()
}],
};
let quick_insert_a = QuickInsert::new(distribution, n_disks / 2);
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::B,
..SiteProperties::default()
}],
};
let quick_insert_b =
QuickInsert::new(distribution, n_disks - quick_insert_a.target());
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume =
n_disks as f64 * circle.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(TypeDependentInteractions {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
quick_insert_a,
quick_insert_b,
quick_compress,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct TypeDependentInteractions {
/// 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 insert algorithm for A bodies.
quick_insert_a: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick insert algorithm for B bodies.
quick_insert_b: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// How sites interact during compression.
overlap_penalty_hamiltonian:
PairwiseCutoff<Isotropic<Expanded<OverlapPenalty>>>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for TypeDependentInteractions {
/// 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 TypeDependentInteractions {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert_a.is_complete()
&& self.quick_insert_b.is_complete()
{
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert_a
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
self.quick_insert_b
.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_default(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
);
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_a.target() + self.quick_insert_b.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,
);
}
}
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_type_id(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.site_type as u32),
)?
.particles_types(SiteType::VARIANTS.iter().copied())
}
}
// 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_interaction::TotalEnergy;
let mut simulation = TypeDependentInteractions::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("type-dependent-interactions.gsd")?;
for _ in 0..100_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.log_scalar(
"potential_energy",
simulation.hamiltonian.total_energy(&simulation.microstate),
)?
.log_scalars(
"particles/site_energy",
simulation
.microstate
.iter_sites_tag_order()
.map(|s| site_energy(&simulation, s)),
)?
.end()?;
}
}
Ok(())
}
#[allow(dead_code, reason = "site_energy is used in the batch mode main()")]
fn site_energy(
simulation: &TypeDependentInteractions,
site: &Site<SiteProperties>,
) -> f64 {
use hoomd_interaction::DeltaEnergyRemove;
let body_index = simulation.microstate.body_indices()[site.body_tag]
.expect("site's parent body should be in the microstate");
-simulation
.hamiltonian
.delta_energy_remove(&simulation.microstate, body_index)
}
#[cfg(feature = "bevy")]
mod type_dependent_interactions_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use type_dependent_interactions_interactive::main;First, compute the distance between the two sites. Then use Rust’s match
expression to compute the desired potential as a function of the two site
types. Rust will produce a compile error should you miss one or more cases.
You can compute interactions that are not included in hoomd-interaction by
writing the expression directly, as demonstrated in the B-B case.
Important
Ensure that
site_pair_energy(i,j)computes the same energy assite_pair_energy(j,i). Rust does not enforce this symmetry and hoomd-rs cannot detect when there is a problem.
Construct the Simulation Model
The new() method constructs a new simulation model:
use anyhow::{Context, anyhow};
use strum::VariantNames;
use strum_macros::VariantNames;
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
pairwise::Isotropic,
univariate::{
Expanded, LennardJones, OverlapPenalty, UnivariateEnergy,
WeeksChandlerAnderson,
},
};
use hoomd_mc::{
QuickCompress, QuickInsert, Sweep, Translate, Trial, Tune, UniformIn,
};
use hoomd_microstate::{AppendMicrostate, Site};
use hoomd_microstate::{
Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, PartialEq, VariantNames)]
enum SiteType {
#[default]
A,
B,
}
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's type.
site_type: SiteType,
}
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 {
lj_aa: LennardJones<12, 6>,
wca_ab: WeeksChandlerAnderson,
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(
&self,
site_properties_i: &SiteProperties,
site_properties_j: &SiteProperties,
) -> f64 {
let r = site_properties_i
.position
.distance(&site_properties_j.position);
match (site_properties_i.site_type, site_properties_j.site_type) {
(SiteType::A, SiteType::A) => self.lj_aa.energy(r),
(SiteType::A, SiteType::B) | (SiteType::B, SiteType::A) => {
self.wca_ab.energy(r)
}
(SiteType::B, SiteType::B) => {
1.0 / r.powi(12) - f64::exp(-1.0 / 2.0 * r.powi(2))
}
}
}
}
impl TypeDependentInteractions {
/// Construct a new type-dependent interactions simulation.
fn new() -> anyhow::Result<TypeDependentInteractions> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.5;
let n_disks = 512;
let maximum_distance = 0.07;
let sigma = 1.0;
let macrostate = Isothermal { temperature: 1.0 };
let lj_aa = LennardJones {
epsilon: 2.0,
sigma: 1.0,
};
let wca_ab = WeeksChandlerAnderson {
epsilon: 1.0,
sigma: 1.0,
};
let hamiltonian = PairwiseCutoff(SitePairInteraction {
lj_aa,
wca_ab,
maximum_interaction_range: 2.5,
});
let overlap_penalty = Isotropic {
interaction: Expanded {
delta: sigma,
f: OverlapPenalty::default(),
},
r_cut: sigma,
};
let overlap_penalty_hamiltonian = PairwiseCutoff(overlap_penalty);
let circle = Circle {
radius: (sigma / 2.0).try_into()?,
};
let initial_box_volume =
n_disks as f64 * circle.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::A,
..SiteProperties::default()
}],
};
let quick_insert_a = QuickInsert::new(distribution, n_disks / 2);
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::B,
..SiteProperties::default()
}],
};
let quick_insert_b =
QuickInsert::new(distribution, n_disks - quick_insert_a.target());
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume =
n_disks as f64 * circle.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(TypeDependentInteractions {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
quick_insert_a,
quick_insert_b,
quick_compress,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct TypeDependentInteractions {
/// 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 insert algorithm for A bodies.
quick_insert_a: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick insert algorithm for B bodies.
quick_insert_b: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// How sites interact during compression.
overlap_penalty_hamiltonian:
PairwiseCutoff<Isotropic<Expanded<OverlapPenalty>>>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for TypeDependentInteractions {
/// 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 TypeDependentInteractions {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert_a.is_complete()
&& self.quick_insert_b.is_complete()
{
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert_a
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
self.quick_insert_b
.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_default(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
);
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_a.target() + self.quick_insert_b.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,
);
}
}
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_type_id(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.site_type as u32),
)?
.particles_types(SiteType::VARIANTS.iter().copied())
}
}
// 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_interaction::TotalEnergy;
let mut simulation = TypeDependentInteractions::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("type-dependent-interactions.gsd")?;
for _ in 0..100_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.log_scalar(
"potential_energy",
simulation.hamiltonian.total_energy(&simulation.microstate),
)?
.log_scalars(
"particles/site_energy",
simulation
.microstate
.iter_sites_tag_order()
.map(|s| site_energy(&simulation, s)),
)?
.end()?;
}
}
Ok(())
}
#[allow(dead_code, reason = "site_energy is used in the batch mode main()")]
fn site_energy(
simulation: &TypeDependentInteractions,
site: &Site<SiteProperties>,
) -> f64 {
use hoomd_interaction::DeltaEnergyRemove;
let body_index = simulation.microstate.body_indices()[site.body_tag]
.expect("site's parent body should be in the microstate");
-simulation
.hamiltonian
.delta_energy_remove(&simulation.microstate, body_index)
}
#[cfg(feature = "bevy")]
mod type_dependent_interactions_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use type_dependent_interactions_interactive::main;Parameters
Assign all the model parameters in one code block:
use anyhow::{Context, anyhow};
use strum::VariantNames;
use strum_macros::VariantNames;
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
pairwise::Isotropic,
univariate::{
Expanded, LennardJones, OverlapPenalty, UnivariateEnergy,
WeeksChandlerAnderson,
},
};
use hoomd_mc::{
QuickCompress, QuickInsert, Sweep, Translate, Trial, Tune, UniformIn,
};
use hoomd_microstate::{AppendMicrostate, Site};
use hoomd_microstate::{
Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, PartialEq, VariantNames)]
enum SiteType {
#[default]
A,
B,
}
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's type.
site_type: SiteType,
}
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 {
lj_aa: LennardJones<12, 6>,
wca_ab: WeeksChandlerAnderson,
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(
&self,
site_properties_i: &SiteProperties,
site_properties_j: &SiteProperties,
) -> f64 {
let r = site_properties_i
.position
.distance(&site_properties_j.position);
match (site_properties_i.site_type, site_properties_j.site_type) {
(SiteType::A, SiteType::A) => self.lj_aa.energy(r),
(SiteType::A, SiteType::B) | (SiteType::B, SiteType::A) => {
self.wca_ab.energy(r)
}
(SiteType::B, SiteType::B) => {
1.0 / r.powi(12) - f64::exp(-1.0 / 2.0 * r.powi(2))
}
}
}
}
impl TypeDependentInteractions {
/// Construct a new type-dependent interactions simulation.
fn new() -> anyhow::Result<TypeDependentInteractions> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.5;
let n_disks = 512;
let maximum_distance = 0.07;
let sigma = 1.0;
let macrostate = Isothermal { temperature: 1.0 };
let lj_aa = LennardJones {
epsilon: 2.0,
sigma: 1.0,
};
let wca_ab = WeeksChandlerAnderson {
epsilon: 1.0,
sigma: 1.0,
};
let hamiltonian = PairwiseCutoff(SitePairInteraction {
lj_aa,
wca_ab,
maximum_interaction_range: 2.5,
});
let overlap_penalty = Isotropic {
interaction: Expanded {
delta: sigma,
f: OverlapPenalty::default(),
},
r_cut: sigma,
};
let overlap_penalty_hamiltonian = PairwiseCutoff(overlap_penalty);
let circle = Circle {
radius: (sigma / 2.0).try_into()?,
};
let initial_box_volume =
n_disks as f64 * circle.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::A,
..SiteProperties::default()
}],
};
let quick_insert_a = QuickInsert::new(distribution, n_disks / 2);
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::B,
..SiteProperties::default()
}],
};
let quick_insert_b =
QuickInsert::new(distribution, n_disks - quick_insert_a.target());
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume =
n_disks as f64 * circle.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(TypeDependentInteractions {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
quick_insert_a,
quick_insert_b,
quick_compress,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct TypeDependentInteractions {
/// 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 insert algorithm for A bodies.
quick_insert_a: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick insert algorithm for B bodies.
quick_insert_b: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// How sites interact during compression.
overlap_penalty_hamiltonian:
PairwiseCutoff<Isotropic<Expanded<OverlapPenalty>>>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for TypeDependentInteractions {
/// 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 TypeDependentInteractions {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert_a.is_complete()
&& self.quick_insert_b.is_complete()
{
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert_a
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
self.quick_insert_b
.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_default(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
);
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_a.target() + self.quick_insert_b.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,
);
}
}
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_type_id(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.site_type as u32),
)?
.particles_types(SiteType::VARIANTS.iter().copied())
}
}
// 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_interaction::TotalEnergy;
let mut simulation = TypeDependentInteractions::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("type-dependent-interactions.gsd")?;
for _ in 0..100_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.log_scalar(
"potential_energy",
simulation.hamiltonian.total_energy(&simulation.microstate),
)?
.log_scalars(
"particles/site_energy",
simulation
.microstate
.iter_sites_tag_order()
.map(|s| site_energy(&simulation, s)),
)?
.end()?;
}
}
Ok(())
}
#[allow(dead_code, reason = "site_energy is used in the batch mode main()")]
fn site_energy(
simulation: &TypeDependentInteractions,
site: &Site<SiteProperties>,
) -> f64 {
use hoomd_interaction::DeltaEnergyRemove;
let body_index = simulation.microstate.body_indices()[site.body_tag]
.expect("site's parent body should be in the microstate");
-simulation
.hamiltonian
.delta_energy_remove(&simulation.microstate, body_index)
}
#[cfg(feature = "bevy")]
mod type_dependent_interactions_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use type_dependent_interactions_interactive::main;initial_packing_fraction is the volume of the disks divided by the volume
of the simulation boundary in the initial state. Choose this value so
that 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), sigma is the disk diameter, and macrostate holds the current
temperature set point (in units of energy).
Hamiltonian
Construct the SitePairInteraction type and use it with PairwiseCutoff to
form the system’s Hamiltonian:
use anyhow::{Context, anyhow};
use strum::VariantNames;
use strum_macros::VariantNames;
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
pairwise::Isotropic,
univariate::{
Expanded, LennardJones, OverlapPenalty, UnivariateEnergy,
WeeksChandlerAnderson,
},
};
use hoomd_mc::{
QuickCompress, QuickInsert, Sweep, Translate, Trial, Tune, UniformIn,
};
use hoomd_microstate::{AppendMicrostate, Site};
use hoomd_microstate::{
Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, PartialEq, VariantNames)]
enum SiteType {
#[default]
A,
B,
}
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's type.
site_type: SiteType,
}
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 {
lj_aa: LennardJones<12, 6>,
wca_ab: WeeksChandlerAnderson,
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(
&self,
site_properties_i: &SiteProperties,
site_properties_j: &SiteProperties,
) -> f64 {
let r = site_properties_i
.position
.distance(&site_properties_j.position);
match (site_properties_i.site_type, site_properties_j.site_type) {
(SiteType::A, SiteType::A) => self.lj_aa.energy(r),
(SiteType::A, SiteType::B) | (SiteType::B, SiteType::A) => {
self.wca_ab.energy(r)
}
(SiteType::B, SiteType::B) => {
1.0 / r.powi(12) - f64::exp(-1.0 / 2.0 * r.powi(2))
}
}
}
}
impl TypeDependentInteractions {
/// Construct a new type-dependent interactions simulation.
fn new() -> anyhow::Result<TypeDependentInteractions> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.5;
let n_disks = 512;
let maximum_distance = 0.07;
let sigma = 1.0;
let macrostate = Isothermal { temperature: 1.0 };
let lj_aa = LennardJones {
epsilon: 2.0,
sigma: 1.0,
};
let wca_ab = WeeksChandlerAnderson {
epsilon: 1.0,
sigma: 1.0,
};
let hamiltonian = PairwiseCutoff(SitePairInteraction {
lj_aa,
wca_ab,
maximum_interaction_range: 2.5,
});
let overlap_penalty = Isotropic {
interaction: Expanded {
delta: sigma,
f: OverlapPenalty::default(),
},
r_cut: sigma,
};
let overlap_penalty_hamiltonian = PairwiseCutoff(overlap_penalty);
let circle = Circle {
radius: (sigma / 2.0).try_into()?,
};
let initial_box_volume =
n_disks as f64 * circle.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::A,
..SiteProperties::default()
}],
};
let quick_insert_a = QuickInsert::new(distribution, n_disks / 2);
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::B,
..SiteProperties::default()
}],
};
let quick_insert_b =
QuickInsert::new(distribution, n_disks - quick_insert_a.target());
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume =
n_disks as f64 * circle.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(TypeDependentInteractions {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
quick_insert_a,
quick_insert_b,
quick_compress,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct TypeDependentInteractions {
/// 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 insert algorithm for A bodies.
quick_insert_a: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick insert algorithm for B bodies.
quick_insert_b: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// How sites interact during compression.
overlap_penalty_hamiltonian:
PairwiseCutoff<Isotropic<Expanded<OverlapPenalty>>>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for TypeDependentInteractions {
/// 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 TypeDependentInteractions {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert_a.is_complete()
&& self.quick_insert_b.is_complete()
{
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert_a
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
self.quick_insert_b
.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_default(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
);
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_a.target() + self.quick_insert_b.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,
);
}
}
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_type_id(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.site_type as u32),
)?
.particles_types(SiteType::VARIANTS.iter().copied())
}
}
// 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_interaction::TotalEnergy;
let mut simulation = TypeDependentInteractions::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("type-dependent-interactions.gsd")?;
for _ in 0..100_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.log_scalar(
"potential_energy",
simulation.hamiltonian.total_energy(&simulation.microstate),
)?
.log_scalars(
"particles/site_energy",
simulation
.microstate
.iter_sites_tag_order()
.map(|s| site_energy(&simulation, s)),
)?
.end()?;
}
}
Ok(())
}
#[allow(dead_code, reason = "site_energy is used in the batch mode main()")]
fn site_energy(
simulation: &TypeDependentInteractions,
site: &Site<SiteProperties>,
) -> f64 {
use hoomd_interaction::DeltaEnergyRemove;
let body_index = simulation.microstate.body_indices()[site.body_tag]
.expect("site's parent body should be in the microstate");
-simulation
.hamiltonian
.delta_energy_remove(&simulation.microstate, body_index)
}
#[cfg(feature = "bevy")]
mod type_dependent_interactions_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use type_dependent_interactions_interactive::main;Overlap Penalty Hamiltonian
Even though SitePairInteraction treats all sites as points with well-defined
potentials at all distances , it is helpful to place sites where
as it can take many steps to relax high energy states.
Use the hard disk OverlapPenalty as a Hamiltonian when initializing:
use anyhow::{Context, anyhow};
use strum::VariantNames;
use strum_macros::VariantNames;
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
pairwise::Isotropic,
univariate::{
Expanded, LennardJones, OverlapPenalty, UnivariateEnergy,
WeeksChandlerAnderson,
},
};
use hoomd_mc::{
QuickCompress, QuickInsert, Sweep, Translate, Trial, Tune, UniformIn,
};
use hoomd_microstate::{AppendMicrostate, Site};
use hoomd_microstate::{
Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, PartialEq, VariantNames)]
enum SiteType {
#[default]
A,
B,
}
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's type.
site_type: SiteType,
}
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 {
lj_aa: LennardJones<12, 6>,
wca_ab: WeeksChandlerAnderson,
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(
&self,
site_properties_i: &SiteProperties,
site_properties_j: &SiteProperties,
) -> f64 {
let r = site_properties_i
.position
.distance(&site_properties_j.position);
match (site_properties_i.site_type, site_properties_j.site_type) {
(SiteType::A, SiteType::A) => self.lj_aa.energy(r),
(SiteType::A, SiteType::B) | (SiteType::B, SiteType::A) => {
self.wca_ab.energy(r)
}
(SiteType::B, SiteType::B) => {
1.0 / r.powi(12) - f64::exp(-1.0 / 2.0 * r.powi(2))
}
}
}
}
impl TypeDependentInteractions {
/// Construct a new type-dependent interactions simulation.
fn new() -> anyhow::Result<TypeDependentInteractions> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.5;
let n_disks = 512;
let maximum_distance = 0.07;
let sigma = 1.0;
let macrostate = Isothermal { temperature: 1.0 };
let lj_aa = LennardJones {
epsilon: 2.0,
sigma: 1.0,
};
let wca_ab = WeeksChandlerAnderson {
epsilon: 1.0,
sigma: 1.0,
};
let hamiltonian = PairwiseCutoff(SitePairInteraction {
lj_aa,
wca_ab,
maximum_interaction_range: 2.5,
});
let overlap_penalty = Isotropic {
interaction: Expanded {
delta: sigma,
f: OverlapPenalty::default(),
},
r_cut: sigma,
};
let overlap_penalty_hamiltonian = PairwiseCutoff(overlap_penalty);
let circle = Circle {
radius: (sigma / 2.0).try_into()?,
};
let initial_box_volume =
n_disks as f64 * circle.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::A,
..SiteProperties::default()
}],
};
let quick_insert_a = QuickInsert::new(distribution, n_disks / 2);
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::B,
..SiteProperties::default()
}],
};
let quick_insert_b =
QuickInsert::new(distribution, n_disks - quick_insert_a.target());
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume =
n_disks as f64 * circle.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(TypeDependentInteractions {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
quick_insert_a,
quick_insert_b,
quick_compress,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct TypeDependentInteractions {
/// 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 insert algorithm for A bodies.
quick_insert_a: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick insert algorithm for B bodies.
quick_insert_b: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// How sites interact during compression.
overlap_penalty_hamiltonian:
PairwiseCutoff<Isotropic<Expanded<OverlapPenalty>>>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for TypeDependentInteractions {
/// 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 TypeDependentInteractions {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert_a.is_complete()
&& self.quick_insert_b.is_complete()
{
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert_a
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
self.quick_insert_b
.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_default(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
);
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_a.target() + self.quick_insert_b.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,
);
}
}
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_type_id(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.site_type as u32),
)?
.particles_types(SiteType::VARIANTS.iter().copied())
}
}
// 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_interaction::TotalEnergy;
let mut simulation = TypeDependentInteractions::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("type-dependent-interactions.gsd")?;
for _ in 0..100_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.log_scalar(
"potential_energy",
simulation.hamiltonian.total_energy(&simulation.microstate),
)?
.log_scalars(
"particles/site_energy",
simulation
.microstate
.iter_sites_tag_order()
.map(|s| site_energy(&simulation, s)),
)?
.end()?;
}
}
Ok(())
}
#[allow(dead_code, reason = "site_energy is used in the batch mode main()")]
fn site_energy(
simulation: &TypeDependentInteractions,
site: &Site<SiteProperties>,
) -> f64 {
use hoomd_interaction::DeltaEnergyRemove;
let body_index = simulation.microstate.body_indices()[site.body_tag]
.expect("site's parent body should be in the microstate");
-simulation
.hamiltonian
.delta_energy_remove(&simulation.microstate, body_index)
}
#[cfg(feature = "bevy")]
mod type_dependent_interactions_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use type_dependent_interactions_interactive::main;Construct the Boundary
Place all bodies in periodic square boundary conditions at the chosen packing fraction:
use anyhow::{Context, anyhow};
use strum::VariantNames;
use strum_macros::VariantNames;
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
pairwise::Isotropic,
univariate::{
Expanded, LennardJones, OverlapPenalty, UnivariateEnergy,
WeeksChandlerAnderson,
},
};
use hoomd_mc::{
QuickCompress, QuickInsert, Sweep, Translate, Trial, Tune, UniformIn,
};
use hoomd_microstate::{AppendMicrostate, Site};
use hoomd_microstate::{
Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, PartialEq, VariantNames)]
enum SiteType {
#[default]
A,
B,
}
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's type.
site_type: SiteType,
}
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 {
lj_aa: LennardJones<12, 6>,
wca_ab: WeeksChandlerAnderson,
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(
&self,
site_properties_i: &SiteProperties,
site_properties_j: &SiteProperties,
) -> f64 {
let r = site_properties_i
.position
.distance(&site_properties_j.position);
match (site_properties_i.site_type, site_properties_j.site_type) {
(SiteType::A, SiteType::A) => self.lj_aa.energy(r),
(SiteType::A, SiteType::B) | (SiteType::B, SiteType::A) => {
self.wca_ab.energy(r)
}
(SiteType::B, SiteType::B) => {
1.0 / r.powi(12) - f64::exp(-1.0 / 2.0 * r.powi(2))
}
}
}
}
impl TypeDependentInteractions {
/// Construct a new type-dependent interactions simulation.
fn new() -> anyhow::Result<TypeDependentInteractions> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.5;
let n_disks = 512;
let maximum_distance = 0.07;
let sigma = 1.0;
let macrostate = Isothermal { temperature: 1.0 };
let lj_aa = LennardJones {
epsilon: 2.0,
sigma: 1.0,
};
let wca_ab = WeeksChandlerAnderson {
epsilon: 1.0,
sigma: 1.0,
};
let hamiltonian = PairwiseCutoff(SitePairInteraction {
lj_aa,
wca_ab,
maximum_interaction_range: 2.5,
});
let overlap_penalty = Isotropic {
interaction: Expanded {
delta: sigma,
f: OverlapPenalty::default(),
},
r_cut: sigma,
};
let overlap_penalty_hamiltonian = PairwiseCutoff(overlap_penalty);
let circle = Circle {
radius: (sigma / 2.0).try_into()?,
};
let initial_box_volume =
n_disks as f64 * circle.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::A,
..SiteProperties::default()
}],
};
let quick_insert_a = QuickInsert::new(distribution, n_disks / 2);
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::B,
..SiteProperties::default()
}],
};
let quick_insert_b =
QuickInsert::new(distribution, n_disks - quick_insert_a.target());
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume =
n_disks as f64 * circle.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(TypeDependentInteractions {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
quick_insert_a,
quick_insert_b,
quick_compress,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct TypeDependentInteractions {
/// 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 insert algorithm for A bodies.
quick_insert_a: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick insert algorithm for B bodies.
quick_insert_b: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// How sites interact during compression.
overlap_penalty_hamiltonian:
PairwiseCutoff<Isotropic<Expanded<OverlapPenalty>>>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for TypeDependentInteractions {
/// 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 TypeDependentInteractions {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert_a.is_complete()
&& self.quick_insert_b.is_complete()
{
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert_a
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
self.quick_insert_b
.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_default(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
);
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_a.target() + self.quick_insert_b.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,
);
}
}
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_type_id(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.site_type as u32),
)?
.particles_types(SiteType::VARIANTS.iter().copied())
}
}
// 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_interaction::TotalEnergy;
let mut simulation = TypeDependentInteractions::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("type-dependent-interactions.gsd")?;
for _ in 0..100_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.log_scalar(
"potential_energy",
simulation.hamiltonian.total_energy(&simulation.microstate),
)?
.log_scalars(
"particles/site_energy",
simulation
.microstate
.iter_sites_tag_order()
.map(|s| site_energy(&simulation, s)),
)?
.end()?;
}
}
Ok(())
}
#[allow(dead_code, reason = "site_energy is used in the batch mode main()")]
fn site_energy(
simulation: &TypeDependentInteractions,
site: &Site<SiteProperties>,
) -> f64 {
use hoomd_interaction::DeltaEnergyRemove;
let body_index = simulation.microstate.body_indices()[site.body_tag]
.expect("site's parent body should be in the microstate");
-simulation
.hamiltonian
.delta_energy_remove(&simulation.microstate, body_index)
}
#[cfg(feature = "bevy")]
mod type_dependent_interactions_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use type_dependent_interactions_interactive::main;Place A and B bodies with QuickInsert
Use one QuickInsert to place half of the bodies with site type A and a
second to place the remainder with site type B:
use anyhow::{Context, anyhow};
use strum::VariantNames;
use strum_macros::VariantNames;
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
pairwise::Isotropic,
univariate::{
Expanded, LennardJones, OverlapPenalty, UnivariateEnergy,
WeeksChandlerAnderson,
},
};
use hoomd_mc::{
QuickCompress, QuickInsert, Sweep, Translate, Trial, Tune, UniformIn,
};
use hoomd_microstate::{AppendMicrostate, Site};
use hoomd_microstate::{
Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, PartialEq, VariantNames)]
enum SiteType {
#[default]
A,
B,
}
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's type.
site_type: SiteType,
}
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 {
lj_aa: LennardJones<12, 6>,
wca_ab: WeeksChandlerAnderson,
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(
&self,
site_properties_i: &SiteProperties,
site_properties_j: &SiteProperties,
) -> f64 {
let r = site_properties_i
.position
.distance(&site_properties_j.position);
match (site_properties_i.site_type, site_properties_j.site_type) {
(SiteType::A, SiteType::A) => self.lj_aa.energy(r),
(SiteType::A, SiteType::B) | (SiteType::B, SiteType::A) => {
self.wca_ab.energy(r)
}
(SiteType::B, SiteType::B) => {
1.0 / r.powi(12) - f64::exp(-1.0 / 2.0 * r.powi(2))
}
}
}
}
impl TypeDependentInteractions {
/// Construct a new type-dependent interactions simulation.
fn new() -> anyhow::Result<TypeDependentInteractions> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.5;
let n_disks = 512;
let maximum_distance = 0.07;
let sigma = 1.0;
let macrostate = Isothermal { temperature: 1.0 };
let lj_aa = LennardJones {
epsilon: 2.0,
sigma: 1.0,
};
let wca_ab = WeeksChandlerAnderson {
epsilon: 1.0,
sigma: 1.0,
};
let hamiltonian = PairwiseCutoff(SitePairInteraction {
lj_aa,
wca_ab,
maximum_interaction_range: 2.5,
});
let overlap_penalty = Isotropic {
interaction: Expanded {
delta: sigma,
f: OverlapPenalty::default(),
},
r_cut: sigma,
};
let overlap_penalty_hamiltonian = PairwiseCutoff(overlap_penalty);
let circle = Circle {
radius: (sigma / 2.0).try_into()?,
};
let initial_box_volume =
n_disks as f64 * circle.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::A,
..SiteProperties::default()
}],
};
let quick_insert_a = QuickInsert::new(distribution, n_disks / 2);
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::B,
..SiteProperties::default()
}],
};
let quick_insert_b =
QuickInsert::new(distribution, n_disks - quick_insert_a.target());
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume =
n_disks as f64 * circle.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(TypeDependentInteractions {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
quick_insert_a,
quick_insert_b,
quick_compress,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct TypeDependentInteractions {
/// 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 insert algorithm for A bodies.
quick_insert_a: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick insert algorithm for B bodies.
quick_insert_b: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// How sites interact during compression.
overlap_penalty_hamiltonian:
PairwiseCutoff<Isotropic<Expanded<OverlapPenalty>>>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for TypeDependentInteractions {
/// 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 TypeDependentInteractions {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert_a.is_complete()
&& self.quick_insert_b.is_complete()
{
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert_a
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
self.quick_insert_b
.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_default(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
);
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_a.target() + self.quick_insert_b.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,
);
}
}
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_type_id(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.site_type as u32),
)?
.particles_types(SiteType::VARIANTS.iter().copied())
}
}
// 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_interaction::TotalEnergy;
let mut simulation = TypeDependentInteractions::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("type-dependent-interactions.gsd")?;
for _ in 0..100_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.log_scalar(
"potential_energy",
simulation.hamiltonian.total_energy(&simulation.microstate),
)?
.log_scalars(
"particles/site_energy",
simulation
.microstate
.iter_sites_tag_order()
.map(|s| site_energy(&simulation, s)),
)?
.end()?;
}
}
Ok(())
}
#[allow(dead_code, reason = "site_energy is used in the batch mode main()")]
fn site_energy(
simulation: &TypeDependentInteractions,
site: &Site<SiteProperties>,
) -> f64 {
use hoomd_interaction::DeltaEnergyRemove;
let body_index = simulation.microstate.body_indices()[site.body_tag]
.expect("site's parent body should be in the microstate");
-simulation
.hamiltonian
.delta_energy_remove(&simulation.microstate, body_index)
}
#[cfg(feature = "bevy")]
mod type_dependent_interactions_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use type_dependent_interactions_interactive::main;Initialization and Simulation
The remaining initialization and simulation code is very similar to that in the
Hard Ellipse Self-Assembly tutorial. The differences are that rotation moves
are not present here and there are two quick_insert methods to apply instead
of one (see also the complete code below).
Execute the Simulation in Batch Mode
Implement AppendMicrostate
Patchy Particle Self-assembly and previous tutorials could call:
hoomd_gsd_file.append_microstate(&simulation.microstate) and it just worked.
hoomd-rs has implemented the AppendMicrostate trait for typical combinations
of Point and OrientedPoint site properties (in 2D and 3D) and various
boundary conditions.
When you customize the site properties and/or boundary types, you also need
to implement AppendMicrostate to write GSD files. Your implementation must
to append a frame at the current step, transform your boundary to a GSD box
definition, write particle positions projected to Cartesian 3-vectors, and can
also write any other data chunks you would like in the frame:
use anyhow::{Context, anyhow};
use strum::VariantNames;
use strum_macros::VariantNames;
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
pairwise::Isotropic,
univariate::{
Expanded, LennardJones, OverlapPenalty, UnivariateEnergy,
WeeksChandlerAnderson,
},
};
use hoomd_mc::{
QuickCompress, QuickInsert, Sweep, Translate, Trial, Tune, UniformIn,
};
use hoomd_microstate::{AppendMicrostate, Site};
use hoomd_microstate::{
Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, PartialEq, VariantNames)]
enum SiteType {
#[default]
A,
B,
}
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's type.
site_type: SiteType,
}
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 {
lj_aa: LennardJones<12, 6>,
wca_ab: WeeksChandlerAnderson,
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(
&self,
site_properties_i: &SiteProperties,
site_properties_j: &SiteProperties,
) -> f64 {
let r = site_properties_i
.position
.distance(&site_properties_j.position);
match (site_properties_i.site_type, site_properties_j.site_type) {
(SiteType::A, SiteType::A) => self.lj_aa.energy(r),
(SiteType::A, SiteType::B) | (SiteType::B, SiteType::A) => {
self.wca_ab.energy(r)
}
(SiteType::B, SiteType::B) => {
1.0 / r.powi(12) - f64::exp(-1.0 / 2.0 * r.powi(2))
}
}
}
}
impl TypeDependentInteractions {
/// Construct a new type-dependent interactions simulation.
fn new() -> anyhow::Result<TypeDependentInteractions> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.5;
let n_disks = 512;
let maximum_distance = 0.07;
let sigma = 1.0;
let macrostate = Isothermal { temperature: 1.0 };
let lj_aa = LennardJones {
epsilon: 2.0,
sigma: 1.0,
};
let wca_ab = WeeksChandlerAnderson {
epsilon: 1.0,
sigma: 1.0,
};
let hamiltonian = PairwiseCutoff(SitePairInteraction {
lj_aa,
wca_ab,
maximum_interaction_range: 2.5,
});
let overlap_penalty = Isotropic {
interaction: Expanded {
delta: sigma,
f: OverlapPenalty::default(),
},
r_cut: sigma,
};
let overlap_penalty_hamiltonian = PairwiseCutoff(overlap_penalty);
let circle = Circle {
radius: (sigma / 2.0).try_into()?,
};
let initial_box_volume =
n_disks as f64 * circle.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::A,
..SiteProperties::default()
}],
};
let quick_insert_a = QuickInsert::new(distribution, n_disks / 2);
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::B,
..SiteProperties::default()
}],
};
let quick_insert_b =
QuickInsert::new(distribution, n_disks - quick_insert_a.target());
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume =
n_disks as f64 * circle.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(TypeDependentInteractions {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
quick_insert_a,
quick_insert_b,
quick_compress,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct TypeDependentInteractions {
/// 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 insert algorithm for A bodies.
quick_insert_a: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick insert algorithm for B bodies.
quick_insert_b: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// How sites interact during compression.
overlap_penalty_hamiltonian:
PairwiseCutoff<Isotropic<Expanded<OverlapPenalty>>>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for TypeDependentInteractions {
/// 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 TypeDependentInteractions {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert_a.is_complete()
&& self.quick_insert_b.is_complete()
{
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert_a
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
self.quick_insert_b
.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_default(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
);
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_a.target() + self.quick_insert_b.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,
);
}
}
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_type_id(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.site_type as u32),
)?
.particles_types(SiteType::VARIANTS.iter().copied())
}
}
// 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_interaction::TotalEnergy;
let mut simulation = TypeDependentInteractions::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("type-dependent-interactions.gsd")?;
for _ in 0..100_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.log_scalar(
"potential_energy",
simulation.hamiltonian.total_energy(&simulation.microstate),
)?
.log_scalars(
"particles/site_energy",
simulation
.microstate
.iter_sites_tag_order()
.map(|s| site_energy(&simulation, s)),
)?
.end()?;
}
}
Ok(())
}
#[allow(dead_code, reason = "site_energy is used in the batch mode main()")]
fn site_energy(
simulation: &TypeDependentInteractions,
site: &Site<SiteProperties>,
) -> f64 {
use hoomd_interaction::DeltaEnergyRemove;
let body_index = simulation.microstate.body_indices()[site.body_tag]
.expect("site's parent body should be in the microstate");
-simulation
.hamiltonian
.delta_energy_remove(&simulation.microstate, body_index)
}
#[cfg(feature = "bevy")]
mod type_dependent_interactions_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use type_dependent_interactions_interactive::main;append_frame returns a Frame type. You write data to the current frame by
adding methods to the call chain. Notice how the code above calls
.configuration_box()?
.configuration_dimensions()?
.particles_position()?...`)
The frame ends once the call chain completes and the Frame is dropped.
The arguments of each method pass the data to write. For example,
configuration_dimensions takes the number of dimensions (2 or 3). The
particles_position method takes an IntoIterator that produces
Cartesian<3> vectors. Use Microstate’s iter_sites_tag_order to ensure that
all sites are written to the file in increasing tag order.
Enumerations in Rust internally number the variants starting at 0 (by default).
This conveniently allows us to convert a site’s type to a 0-based type id
with the cast: s.properties.site_type as u32. We also need to write the
string type names to the file. Rust itself doesn’t provide any mechanism
to do so. This example uses the strum crate’s #[derive(VariantNames)]
macro on the SiteType type to create the SiteType::VARIANTS array
of string type names.
Tip
Use a built-in implementation in
hoomd-microstate/src/append.rsas a starting point when you need to implementAppendMicrostatefor your custom site and/or boundary types.
Implement main()
To run the simulation, construct the TypeDependentInteractions simulation model.
Then call advance() many times:
use anyhow::{Context, anyhow};
use strum::VariantNames;
use strum_macros::VariantNames;
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
pairwise::Isotropic,
univariate::{
Expanded, LennardJones, OverlapPenalty, UnivariateEnergy,
WeeksChandlerAnderson,
},
};
use hoomd_mc::{
QuickCompress, QuickInsert, Sweep, Translate, Trial, Tune, UniformIn,
};
use hoomd_microstate::{AppendMicrostate, Site};
use hoomd_microstate::{
Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, PartialEq, VariantNames)]
enum SiteType {
#[default]
A,
B,
}
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's type.
site_type: SiteType,
}
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 {
lj_aa: LennardJones<12, 6>,
wca_ab: WeeksChandlerAnderson,
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(
&self,
site_properties_i: &SiteProperties,
site_properties_j: &SiteProperties,
) -> f64 {
let r = site_properties_i
.position
.distance(&site_properties_j.position);
match (site_properties_i.site_type, site_properties_j.site_type) {
(SiteType::A, SiteType::A) => self.lj_aa.energy(r),
(SiteType::A, SiteType::B) | (SiteType::B, SiteType::A) => {
self.wca_ab.energy(r)
}
(SiteType::B, SiteType::B) => {
1.0 / r.powi(12) - f64::exp(-1.0 / 2.0 * r.powi(2))
}
}
}
}
impl TypeDependentInteractions {
/// Construct a new type-dependent interactions simulation.
fn new() -> anyhow::Result<TypeDependentInteractions> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.5;
let n_disks = 512;
let maximum_distance = 0.07;
let sigma = 1.0;
let macrostate = Isothermal { temperature: 1.0 };
let lj_aa = LennardJones {
epsilon: 2.0,
sigma: 1.0,
};
let wca_ab = WeeksChandlerAnderson {
epsilon: 1.0,
sigma: 1.0,
};
let hamiltonian = PairwiseCutoff(SitePairInteraction {
lj_aa,
wca_ab,
maximum_interaction_range: 2.5,
});
let overlap_penalty = Isotropic {
interaction: Expanded {
delta: sigma,
f: OverlapPenalty::default(),
},
r_cut: sigma,
};
let overlap_penalty_hamiltonian = PairwiseCutoff(overlap_penalty);
let circle = Circle {
radius: (sigma / 2.0).try_into()?,
};
let initial_box_volume =
n_disks as f64 * circle.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::A,
..SiteProperties::default()
}],
};
let quick_insert_a = QuickInsert::new(distribution, n_disks / 2);
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::B,
..SiteProperties::default()
}],
};
let quick_insert_b =
QuickInsert::new(distribution, n_disks - quick_insert_a.target());
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume =
n_disks as f64 * circle.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(TypeDependentInteractions {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
quick_insert_a,
quick_insert_b,
quick_compress,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct TypeDependentInteractions {
/// 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 insert algorithm for A bodies.
quick_insert_a: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick insert algorithm for B bodies.
quick_insert_b: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// How sites interact during compression.
overlap_penalty_hamiltonian:
PairwiseCutoff<Isotropic<Expanded<OverlapPenalty>>>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for TypeDependentInteractions {
/// 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 TypeDependentInteractions {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert_a.is_complete()
&& self.quick_insert_b.is_complete()
{
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert_a
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
self.quick_insert_b
.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_default(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
);
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_a.target() + self.quick_insert_b.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,
);
}
}
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_type_id(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.site_type as u32),
)?
.particles_types(SiteType::VARIANTS.iter().copied())
}
}
// 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_interaction::TotalEnergy;
let mut simulation = TypeDependentInteractions::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("type-dependent-interactions.gsd")?;
for _ in 0..100_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.log_scalar(
"potential_energy",
simulation.hamiltonian.total_energy(&simulation.microstate),
)?
.log_scalars(
"particles/site_energy",
simulation
.microstate
.iter_sites_tag_order()
.map(|s| site_energy(&simulation, s)),
)?
.end()?;
}
}
Ok(())
}
#[allow(dead_code, reason = "site_energy is used in the batch mode main()")]
fn site_energy(
simulation: &TypeDependentInteractions,
site: &Site<SiteProperties>,
) -> f64 {
use hoomd_interaction::DeltaEnergyRemove;
let body_index = simulation.microstate.body_indices()[site.body_tag]
.expect("site's parent body should be in the microstate");
-simulation
.hamiltonian
.delta_energy_remove(&simulation.microstate, body_index)
}
#[cfg(feature = "bevy")]
mod type_dependent_interactions_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use type_dependent_interactions_interactive::main;See Applying Interactions for a step-by-step explanation of this code.
Note
This
main()function runs in batch mode. There is a differentmain()(not shown here) used in the interactive example. The interactive example does not write the GSD file.
Log to the GSD File
To log custom quantities to the GSD file, add .log_scalar, .log_scalars and/or
.log_arrays to the call chain after .append_microstate:
use anyhow::{Context, anyhow};
use strum::VariantNames;
use strum_macros::VariantNames;
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
pairwise::Isotropic,
univariate::{
Expanded, LennardJones, OverlapPenalty, UnivariateEnergy,
WeeksChandlerAnderson,
},
};
use hoomd_mc::{
QuickCompress, QuickInsert, Sweep, Translate, Trial, Tune, UniformIn,
};
use hoomd_microstate::{AppendMicrostate, Site};
use hoomd_microstate::{
Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, PartialEq, VariantNames)]
enum SiteType {
#[default]
A,
B,
}
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's type.
site_type: SiteType,
}
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 {
lj_aa: LennardJones<12, 6>,
wca_ab: WeeksChandlerAnderson,
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(
&self,
site_properties_i: &SiteProperties,
site_properties_j: &SiteProperties,
) -> f64 {
let r = site_properties_i
.position
.distance(&site_properties_j.position);
match (site_properties_i.site_type, site_properties_j.site_type) {
(SiteType::A, SiteType::A) => self.lj_aa.energy(r),
(SiteType::A, SiteType::B) | (SiteType::B, SiteType::A) => {
self.wca_ab.energy(r)
}
(SiteType::B, SiteType::B) => {
1.0 / r.powi(12) - f64::exp(-1.0 / 2.0 * r.powi(2))
}
}
}
}
impl TypeDependentInteractions {
/// Construct a new type-dependent interactions simulation.
fn new() -> anyhow::Result<TypeDependentInteractions> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.5;
let n_disks = 512;
let maximum_distance = 0.07;
let sigma = 1.0;
let macrostate = Isothermal { temperature: 1.0 };
let lj_aa = LennardJones {
epsilon: 2.0,
sigma: 1.0,
};
let wca_ab = WeeksChandlerAnderson {
epsilon: 1.0,
sigma: 1.0,
};
let hamiltonian = PairwiseCutoff(SitePairInteraction {
lj_aa,
wca_ab,
maximum_interaction_range: 2.5,
});
let overlap_penalty = Isotropic {
interaction: Expanded {
delta: sigma,
f: OverlapPenalty::default(),
},
r_cut: sigma,
};
let overlap_penalty_hamiltonian = PairwiseCutoff(overlap_penalty);
let circle = Circle {
radius: (sigma / 2.0).try_into()?,
};
let initial_box_volume =
n_disks as f64 * circle.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::A,
..SiteProperties::default()
}],
};
let quick_insert_a = QuickInsert::new(distribution, n_disks / 2);
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::B,
..SiteProperties::default()
}],
};
let quick_insert_b =
QuickInsert::new(distribution, n_disks - quick_insert_a.target());
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume =
n_disks as f64 * circle.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(TypeDependentInteractions {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
quick_insert_a,
quick_insert_b,
quick_compress,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct TypeDependentInteractions {
/// 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 insert algorithm for A bodies.
quick_insert_a: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick insert algorithm for B bodies.
quick_insert_b: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// How sites interact during compression.
overlap_penalty_hamiltonian:
PairwiseCutoff<Isotropic<Expanded<OverlapPenalty>>>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for TypeDependentInteractions {
/// 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 TypeDependentInteractions {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert_a.is_complete()
&& self.quick_insert_b.is_complete()
{
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert_a
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
self.quick_insert_b
.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_default(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
);
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_a.target() + self.quick_insert_b.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,
);
}
}
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_type_id(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.site_type as u32),
)?
.particles_types(SiteType::VARIANTS.iter().copied())
}
}
// 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_interaction::TotalEnergy;
let mut simulation = TypeDependentInteractions::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("type-dependent-interactions.gsd")?;
for _ in 0..100_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.log_scalar(
"potential_energy",
simulation.hamiltonian.total_energy(&simulation.microstate),
)?
.log_scalars(
"particles/site_energy",
simulation
.microstate
.iter_sites_tag_order()
.map(|s| site_energy(&simulation, s)),
)?
.end()?;
}
}
Ok(())
}
#[allow(dead_code, reason = "site_energy is used in the batch mode main()")]
fn site_energy(
simulation: &TypeDependentInteractions,
site: &Site<SiteProperties>,
) -> f64 {
use hoomd_interaction::DeltaEnergyRemove;
let body_index = simulation.microstate.body_indices()[site.body_tag]
.expect("site's parent body should be in the microstate");
-simulation
.hamiltonian
.delta_energy_remove(&simulation.microstate, body_index)
}
#[cfg(feature = "bevy")]
mod type_dependent_interactions_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use type_dependent_interactions_interactive::main;The given quantities are stored in the frame alongside the particle properties.
Log names that start with “particles/” store per-particle values. This tutorial
logs the energy each site contributes to the total computed by the site_energy
helper method:
use anyhow::{Context, anyhow};
use strum::VariantNames;
use strum_macros::VariantNames;
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
pairwise::Isotropic,
univariate::{
Expanded, LennardJones, OverlapPenalty, UnivariateEnergy,
WeeksChandlerAnderson,
},
};
use hoomd_mc::{
QuickCompress, QuickInsert, Sweep, Translate, Trial, Tune, UniformIn,
};
use hoomd_microstate::{AppendMicrostate, Site};
use hoomd_microstate::{
Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, PartialEq, VariantNames)]
enum SiteType {
#[default]
A,
B,
}
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's type.
site_type: SiteType,
}
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 {
lj_aa: LennardJones<12, 6>,
wca_ab: WeeksChandlerAnderson,
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(
&self,
site_properties_i: &SiteProperties,
site_properties_j: &SiteProperties,
) -> f64 {
let r = site_properties_i
.position
.distance(&site_properties_j.position);
match (site_properties_i.site_type, site_properties_j.site_type) {
(SiteType::A, SiteType::A) => self.lj_aa.energy(r),
(SiteType::A, SiteType::B) | (SiteType::B, SiteType::A) => {
self.wca_ab.energy(r)
}
(SiteType::B, SiteType::B) => {
1.0 / r.powi(12) - f64::exp(-1.0 / 2.0 * r.powi(2))
}
}
}
}
impl TypeDependentInteractions {
/// Construct a new type-dependent interactions simulation.
fn new() -> anyhow::Result<TypeDependentInteractions> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.5;
let n_disks = 512;
let maximum_distance = 0.07;
let sigma = 1.0;
let macrostate = Isothermal { temperature: 1.0 };
let lj_aa = LennardJones {
epsilon: 2.0,
sigma: 1.0,
};
let wca_ab = WeeksChandlerAnderson {
epsilon: 1.0,
sigma: 1.0,
};
let hamiltonian = PairwiseCutoff(SitePairInteraction {
lj_aa,
wca_ab,
maximum_interaction_range: 2.5,
});
let overlap_penalty = Isotropic {
interaction: Expanded {
delta: sigma,
f: OverlapPenalty::default(),
},
r_cut: sigma,
};
let overlap_penalty_hamiltonian = PairwiseCutoff(overlap_penalty);
let circle = Circle {
radius: (sigma / 2.0).try_into()?,
};
let initial_box_volume =
n_disks as f64 * circle.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::A,
..SiteProperties::default()
}],
};
let quick_insert_a = QuickInsert::new(distribution, n_disks / 2);
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::B,
..SiteProperties::default()
}],
};
let quick_insert_b =
QuickInsert::new(distribution, n_disks - quick_insert_a.target());
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume =
n_disks as f64 * circle.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(TypeDependentInteractions {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
quick_insert_a,
quick_insert_b,
quick_compress,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct TypeDependentInteractions {
/// 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 insert algorithm for A bodies.
quick_insert_a: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick insert algorithm for B bodies.
quick_insert_b: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// How sites interact during compression.
overlap_penalty_hamiltonian:
PairwiseCutoff<Isotropic<Expanded<OverlapPenalty>>>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for TypeDependentInteractions {
/// 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 TypeDependentInteractions {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert_a.is_complete()
&& self.quick_insert_b.is_complete()
{
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert_a
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
self.quick_insert_b
.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_default(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
);
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_a.target() + self.quick_insert_b.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,
);
}
}
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_type_id(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.site_type as u32),
)?
.particles_types(SiteType::VARIANTS.iter().copied())
}
}
// 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_interaction::TotalEnergy;
let mut simulation = TypeDependentInteractions::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("type-dependent-interactions.gsd")?;
for _ in 0..100_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.log_scalar(
"potential_energy",
simulation.hamiltonian.total_energy(&simulation.microstate),
)?
.log_scalars(
"particles/site_energy",
simulation
.microstate
.iter_sites_tag_order()
.map(|s| site_energy(&simulation, s)),
)?
.end()?;
}
}
Ok(())
}
#[allow(dead_code, reason = "site_energy is used in the batch mode main()")]
fn site_energy(
simulation: &TypeDependentInteractions,
site: &Site<SiteProperties>,
) -> f64 {
use hoomd_interaction::DeltaEnergyRemove;
let body_index = simulation.microstate.body_indices()[site.body_tag]
.expect("site's parent body should be in the microstate");
-simulation
.hamiltonian
.delta_energy_remove(&simulation.microstate, body_index)
}
#[cfg(feature = "bevy")]
mod type_dependent_interactions_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use type_dependent_interactions_interactive::main;The method computes the negative of the energy delta that results when removing the site’s parent body. The bodies in this tutorial are all single-site bodies, so this is equivalent to the site’s contribution. You can use any of the APIs in hoomd-rs to compute and log any site-specific quantity you are interested in.
Note
Use GSD to log per-particle quantities. Parquet and other column-oriented formats cannot easily store thousands of scalar columns.
Prefer parquet when logging a few scalar quantities over large numbers of frames. GSD uses 32 bytes of overhead per data chunk per frame, so a GSD file with only scalar log quantities will be about 5 times larger than an equivalent Parquet file.
Run the Simulation
In a terminal, execute the following command to run the simulation in batch mode:
cargo run --release --example type-dependent-interactions
Visualize the Simulation Results
Open the generated type-dependent-interactions.gsd in Ovito or another
visualization tool to see the simulation results. By default. Ovito displays
the A particles with an appropriate size, but the B particles are too large.
To fix this and color the sites by the logged site_energy quantity, import
this modifier snippet (use the copy button to copy the entire text to your
clipboard):
{"description": "OVITO Modifier Snippet: Edit types (Particle Type) | Color coding (site_energy)", "payload": "AAAHinjafRTLThNR9LS0pS2FykOQRwxqTNyIIrJiAQKNYUFAwC7YkMvMZTownWnmAcGY6ELjwo3xB/wA+QCWbtzAzi/AjQtd6AeY6DkzZ9pbGJjkzL33vN/lo6+fi39PNgGCOQDII/TBClRhCTbwHIdVcMFBbAdCGaKvBFMwAZMwjf+HjEsjpBAy0ZH+wO/bEQ5GEAbBAxN8kLCFYCO4YMChIh2aCS8ZBqJkybPZf4/fGbJUxXuO3l/uvX91/7RwChe/To4kFcsivGFdL/ksEX7VdRrS9cmDCWiAQH989FADC33z4AFG344bx6wcIqdEiUcXqB6M3/oY5uxyqfgrkG+xg0VCpObufprRN8eO8NFFXqaA/Uyn6NrdxPXEsp1cFJRtab7GiIEXZxGil1OZDYuLYSuG+sN3S3ggTPVZ8329aX4QYTg1+yfCz/44fvbr7c3j+ETUEDkU1h7m8X8DYYYtT3HSKYTCqnB9U7OkF3bFVbkaJjVXWBzhXogMjrIxcmGsSZlsxUZhpH8zR4/nBK4mlQ6IeCijIvCdJ/pu4PlrwjakSuvek7KxLi2p+aZjt1E0x3Lcp67QTWn7oLQfnXlp61VhBW3Kip6PuWhDU/93Sd30pb5x2AhzFNY8x7RL3Ka4h7SarJuasCqWrKMLW1teTejOATdRyFOqC89LIgwjqiGfe3JZerUFiiSJ607INS+0vR2hyYXAskzbqNhi25J6Ev9ozTRqFoK/ToIV3ZCJ1rtDvUmU3n39YA1TGiTK9biXksrnsqFmKkNZUBGFphUVm3PbUFSCskwONtdsgCRKyRZ1mSjiHNjSPV9MomRDZXh5HWMypERlKdhBXbqmtqSre69o6hituWNKV93M09xHNCT9YX0XHB2Lt4x/5s0zC23sEYWlarqmbnpKZ2eZk+azt4LNGraqoqrADLRNSvG8E1PkeqIniv4is9By6uL6sXQnk2iPFReFL1a2d3EWFSFaCn2LOKGG8NsDjL2iBZlX8LEk7cnCmtzZEK4h/dY8whzPa4ozuslZTfP7hBumg9/fePAzPOr7XIifvJmI/J3XPsV0k3j/A200N3Q="}
Alternately, you can:
- Add an Edit types modification to the Ovito pipeline.
- Select type B.
- Set the Display radius to 0.5.
- Add a Color coding modification to the pipeline.
- Set the
site_energyInput property and choose theViridisColor gradient
Render with Tachyon, and you should see something like:
Conclusion
This tutorial showed you how to implement custom site properties and site-site interaction types. Specifically, it demonstrated the addition of a site type field and showed how you can use that when computing the interaction energies.
Navigate to the top of the page and refresh to see the simulation in action again. Notice that the sites quickly phase separate. Wait long enough and the system should form two stripes. Due to the differently shaped potentials, the A and B domains are distinct. The A sites form a hexagonal solid and B form a lower density fluid.
Complete Code
use anyhow::{Context, anyhow};
use strum::VariantNames;
use strum_macros::VariantNames;
use hoomd_geometry::{
Volume,
shape::{Circle, Rectangle},
};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use hoomd_interaction::{
MaximumInteractionRange, PairwiseCutoff, SitePairEnergy,
pairwise::Isotropic,
univariate::{
Expanded, LennardJones, OverlapPenalty, UnivariateEnergy,
WeeksChandlerAnderson,
},
};
use hoomd_mc::{
QuickCompress, QuickInsert, Sweep, Translate, Trial, Tune, UniformIn,
};
use hoomd_microstate::{AppendMicrostate, Site};
use hoomd_microstate::{
Microstate, SiteKey, Transform,
boundary::Periodic,
property::{Point, Position},
};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{Cartesian, Metric};
type PositionVector = Cartesian<2>;
type BodyProperties = Point<PositionVector>;
#[derive(Clone, Copy, Default, PartialEq, VariantNames)]
enum SiteType {
#[default]
A,
B,
}
#[derive(Clone, Copy, Default, Position)]
struct SiteProperties {
/// The site's position.
position: PositionVector,
/// The site's type.
site_type: SiteType,
}
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 {
lj_aa: LennardJones<12, 6>,
wca_ab: WeeksChandlerAnderson,
maximum_interaction_range: f64,
}
impl SitePairEnergy<SiteProperties> for SitePairInteraction {
fn site_pair_energy(
&self,
site_properties_i: &SiteProperties,
site_properties_j: &SiteProperties,
) -> f64 {
let r = site_properties_i
.position
.distance(&site_properties_j.position);
match (site_properties_i.site_type, site_properties_j.site_type) {
(SiteType::A, SiteType::A) => self.lj_aa.energy(r),
(SiteType::A, SiteType::B) | (SiteType::B, SiteType::A) => {
self.wca_ab.energy(r)
}
(SiteType::B, SiteType::B) => {
1.0 / r.powi(12) - f64::exp(-1.0 / 2.0 * r.powi(2))
}
}
}
}
impl TypeDependentInteractions {
/// Construct a new type-dependent interactions simulation.
fn new() -> anyhow::Result<TypeDependentInteractions> {
let initial_packing_fraction = 0.3;
let target_packing_fraction = 0.5;
let n_disks = 512;
let maximum_distance = 0.07;
let sigma = 1.0;
let macrostate = Isothermal { temperature: 1.0 };
let lj_aa = LennardJones {
epsilon: 2.0,
sigma: 1.0,
};
let wca_ab = WeeksChandlerAnderson {
epsilon: 1.0,
sigma: 1.0,
};
let hamiltonian = PairwiseCutoff(SitePairInteraction {
lj_aa,
wca_ab,
maximum_interaction_range: 2.5,
});
let overlap_penalty = Isotropic {
interaction: Expanded {
delta: sigma,
f: OverlapPenalty::default(),
},
r_cut: sigma,
};
let overlap_penalty_hamiltonian = PairwiseCutoff(overlap_penalty);
let circle = Circle {
radius: (sigma / 2.0).try_into()?,
};
let initial_box_volume =
n_disks as f64 * circle.volume() / initial_packing_fraction;
let initial_box_edge_length = initial_box_volume.sqrt();
let square =
Rectangle::with_equal_edges(initial_box_edge_length.try_into()?);
let periodic_square =
Periodic::new(hamiltonian.maximum_interaction_range(), square)?;
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::A,
..SiteProperties::default()
}],
};
let quick_insert_a = QuickInsert::new(distribution, n_disks / 2);
let distribution = UniformIn {
boundary: periodic_square.clone(),
template_sites: vec![SiteProperties {
site_type: SiteType::B,
..SiteProperties::default()
}],
};
let quick_insert_b =
QuickInsert::new(distribution, n_disks - quick_insert_a.target());
let vec_cell = VecCell::builder()
.nominal_search_radius(
hamiltonian.maximum_interaction_range().try_into()?,
)
.build();
let microstate = Microstate::builder()
.boundary(periodic_square)
.spatial_data(vec_cell)
.try_build()?;
let translate =
Translate::with_maximum_distance(maximum_distance.try_into()?);
let translate_sweep = Sweep(translate);
let target_box_volume =
n_disks as f64 * circle.volume() / target_packing_fraction;
let quick_compress =
QuickCompress::with_target_volume(target_box_volume.try_into()?);
Ok(TypeDependentInteractions {
microstate,
overlap_penalty_hamiltonian,
hamiltonian,
translate_sweep,
quick_insert_a,
quick_insert_b,
quick_compress,
macrostate,
phase: Phase::Initialize,
})
}
}
#[cfg_attr(feature = "bevy", derive(Resource))]
struct TypeDependentInteractions {
/// 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 insert algorithm for A bodies.
quick_insert_a: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick insert algorithm for B bodies.
quick_insert_b: QuickInsert<UniformIn<SiteProperties, Periodic<Rectangle>>>,
/// Quick compress algorithm
quick_compress: QuickCompress<Periodic<Rectangle>>,
/// How sites interact during compression.
overlap_penalty_hamiltonian:
PairwiseCutoff<Isotropic<Expanded<OverlapPenalty>>>,
/// The current phase of the simulation.
phase: Phase,
}
enum Phase {
Initialize,
Equilibrate,
}
impl Simulation for TypeDependentInteractions {
/// 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 TypeDependentInteractions {
fn initialize(&mut self) -> anyhow::Result<()> {
if self.quick_insert_a.is_complete()
&& self.quick_insert_b.is_complete()
{
self.quick_compress.apply(
&mut self.microstate,
&self.overlap_penalty_hamiltonian,
|_| true,
);
} else {
self.quick_insert_a
.apply(&mut self.microstate, &self.overlap_penalty_hamiltonian);
self.quick_insert_b
.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_default(
&self.microstate,
&self.hamiltonian,
&self.macrostate,
);
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_a.target() + self.quick_insert_b.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,
);
}
}
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_type_id(
microstate
.iter_sites_tag_order()
.map(|s| s.properties.site_type as u32),
)?
.particles_types(SiteType::VARIANTS.iter().copied())
}
}
// 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_interaction::TotalEnergy;
let mut simulation = TypeDependentInteractions::new()?;
let mut hoomd_gsd_file =
HoomdGsdFile::create("type-dependent-interactions.gsd")?;
for _ in 0..100_000 {
simulation.advance()?;
if simulation.step().is_multiple_of(10_000) {
hoomd_gsd_file
.append_microstate(&simulation.microstate)?
.log_scalar(
"potential_energy",
simulation.hamiltonian.total_energy(&simulation.microstate),
)?
.log_scalars(
"particles/site_energy",
simulation
.microstate
.iter_sites_tag_order()
.map(|s| site_energy(&simulation, s)),
)?
.end()?;
}
}
Ok(())
}
#[allow(dead_code, reason = "site_energy is used in the batch mode main()")]
fn site_energy(
simulation: &TypeDependentInteractions,
site: &Site<SiteProperties>,
) -> f64 {
use hoomd_interaction::DeltaEnergyRemove;
let body_index = simulation.microstate.body_indices()[site.body_tag]
.expect("site's parent body should be in the microstate");
-simulation
.hamiltonian
.delta_energy_remove(&simulation.microstate, body_index)
}
#[cfg(feature = "bevy")]
mod type_dependent_interactions_interactive;
#[cfg(feature = "bevy")]
use bevy::prelude::Resource;
#[cfg(feature = "bevy")]
use type_dependent_interactions_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.
