diff --git a/Cargo.lock b/Cargo.lock index b364297c..6c6f212d 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -3061,6 +3061,7 @@ name = "examples" version = "1.2.0" dependencies = [ "anyhow", + "arrayvec", "bevy", "bevy_egui", "hoomd-bevy", diff --git a/doc/src/SUMMARY.md b/doc/src/SUMMARY.md index 0a2ec649..1ed9a417 100644 --- a/doc/src/SUMMARY.md +++ b/doc/src/SUMMARY.md @@ -39,6 +39,7 @@ - [Grandcanonical ensemble simulation of ellipses (2D)](mc-examples/ellipse-gcmc-2d.md) - [Melt an ideal hexagonal structure of hexagons](mc-examples/crystal-stability.md) - [Multi-site patchy shapes](mc-examples/multi-site-patchy-shapes.md) + - [3D shapes confined to a 2D plane](mc-examples/capsules-on-a-plane.md) # Reference diff --git a/doc/src/mc-examples/capsules-on-a-plane.md b/doc/src/mc-examples/capsules-on-a-plane.md new file mode 100644 index 00000000..b9316197 --- /dev/null +++ b/doc/src/mc-examples/capsules-on-a-plane.md @@ -0,0 +1,44 @@ +# 3D shapes confined to a 2D plane + + +{{#include ../../scripts/canvas.html}} + +## Overview + +* Objective: Self-assemble 3D capsules that are allowed to freely rotate while their + center is confined to a 2D plane. +* File: `hoomd-rs/examples/mc-examples/capsules-on-a-plane.rs` +* Run (interactively): + ```shell + cargo run --release --features "bevy" --example capsules-on-a-plane + ``` +* Run (in batch mode): + ```shell + cargo run --release --example capsules-on-a-plane + ``` + +## Strategy + +1. Represent bodies with `Cartesian<2>` positions and `Versor` orientations. +2. Define a custom `SiteProperties` with `Cartesian<3>` positions and `Versor` + orientations. Transform sites from the body frame to the simulation frame by first + lifting the body position into 3D: `(body_x, body_y, 0)` and then transforming + as normal. +3. Define a custom `Boundary` type that wraps `Periodic`. + Implement `Wrap`, `Volume`, `MapPoint`, `Scale`, and `Distribution` + for the custom boundary by calling the same methods on the inner type. + Implement `Wrap` and `GenerateGhosts` by projecting + the 3D site position into 2D, calling the same methods on the wrapped type, + then lift the result back into 3D. + +See the example code for details. It implements all of these steps in a general fashion. +You can copy and paste this code and use it for Monte Carlo simulations with any 3D +interaction model (including multi-site rigid bodies) where the body centers should +be confined to the _xy_ plane. + +## Complete Code +```rust,ignore +{{#rustdoc_include ../../../examples/mc-examples/capsules-on-a-plane.rs:all}} diff --git a/doc/src/release-notes.md b/doc/src/release-notes.md index 79968362..6042f8b1 100644 --- a/doc/src/release-notes.md +++ b/doc/src/release-notes.md @@ -4,13 +4,18 @@ *Added:* +* `[examples]`: Add "Melt an ideal hexagonal structure of hexagons" example (#393). +* `[examples]`: Add "Multi-site patchy shapes" example (#395). +* `[examples]`: Add "3D shapes confined to a 2D plane" example (#397). +* `[examples]`: Use the mouse wheel to zoom examples with 3D cameras (#397). +* `[examples]`: Click and drag to orbit the camera in examples with 3D cameras (#397). * `[hoomd-geometry]`: Add `Hyperparallelepiped`, `Triclinic`, and `Rhomboid` shapes (#88). * `[hoomd-linear-algebra]` Add `matrix::qr` module with methods to compute the QR factorization of matrices (#88). * `[hoomd-linear-algebra]` Add methods `iter_column_slice`, `iter_column_slice_mut`, `iter_submatrix`, to `Matrix` (#88). * `[hoomd-macrostate]`: Add `Fugacity` trait that access the fugacity of a macrostate (#354). * `[hoomd-macrostate]`: Add `IsothermoalIsofugacity` type that stores the system's temperature and fugacity (#354). * `[hoomd-mc]`: Add `GrandCanonical` trial move type. It inserts and removes bodies and samples from the constant chemical potential, constant volume, constant temperature ensemble (#354). -* `[hoomd-mc`]: Add `InsertReMoveCount` type that counts the accepted and rejected insertion and removal moves (#354). +* `[hoomd-mc`]: Add `InsertRemoveCount` type that counts the accepted and rejected insertion and removal moves (#354). * `[hoomd-microstate`]: Implement all traits necessary to use `Hyperparallelepiped`, `Triclinic`, and `Rhomboid` as simulation boundaries (#88). * `[hoomd-microstate`]: Implement `AppendMicrostate` for microstates with `Triclinic` and `Rhomboid` boundaries (#88). * `[hoomd-microstate`]: Add `Microstate::iter_bodies_tag_order` method that iterates over all bodies in tag order (#393). @@ -19,12 +24,14 @@ * `[hoomd-vector`]: Implement `Mul> for f64`, `Mul for Cartesian`, `Mul> for PositiveReal`, `MulAssign for Cartesian`, and `MulAssign for Cartesian` (#88). * `[hoomd-vector`]: Add `Cartesian::basis` method that constructs Cartesian basis vectors (#393). * `[hoomd-workspace`]: Add the `hoomd-workspace` crate. Use it to initialize and read state points in a signac workspace (#359). -* *examples*: Add "Melt an ideal hexagonal structure of hexagons" example (#393). -* *examples*: Add "Multi-site patchy shapes" example (#395). *Changed:* * Build the documentation with mdBook 0.5.4 and KaTeX 0.18.1 (#376). +* `[hoomd-microstate]`: Relax trait bounds to allow bodie and site positions to + have different types (#397). +* `[hoomd-mc]`: Relax trait bounds to allow bodie and site positions to + have different types (#397). *Deprecated:* @@ -32,10 +39,10 @@ *Fixed:* +* `[examples]`: Properly auto size egui parameter windows (#354). +* `[examples]`: Reduce the zoom speed (#393). * `[hoomd-interaction]`: Fix typos in documentation (#358). * `[hoomd-microstate]`: Fix typos in documentation (#358). -* *examples*: Properly auto size egui parameter windows (#354). -* *examples*: Reduce the zoom speed (#393). ## 1.2.0 (2026-07-07) diff --git a/examples/Cargo.toml b/examples/Cargo.toml index 0d3684a3..43182db3 100644 --- a/examples/Cargo.toml +++ b/examples/Cargo.toml @@ -39,6 +39,7 @@ hoomd-spatial.workspace = true hoomd-utility.workspace = true hoomd-vector.workspace = true +arrayvec.workspace = true anyhow.workspace = true bevy = { workspace = true, optional = true } bevy_egui = { workspace = true, optional = true } @@ -169,3 +170,10 @@ path = "mc-examples/multi-site-patchy-shapes.rs" [package.metadata.example.multi-site-patchy-shapes] path = "mc-examples" + +[[example]] +name = "capsules-on-a-plane" +path = "mc-examples/capsules-on-a-plane.rs" + +[package.metadata.example.capsules-on-a-plane] +path = "mc-examples" diff --git a/examples/mc-examples/capsules-on-a-plane.rs b/examples/mc-examples/capsules-on-a-plane.rs new file mode 100644 index 00000000..19f5e633 --- /dev/null +++ b/examples/mc-examples/capsules-on-a-plane.rs @@ -0,0 +1,418 @@ +// ANCHOR: all +use anyhow::{Context, anyhow}; +use arrayvec::ArrayVec; + +use hoomd_geometry::{ + Convex, MapPoint, Scale, Volume, + shape::{Capsule, Rectangle}, +}; +use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile}; +use hoomd_interaction::{ + MaximumInteractionRange, PairwiseCutoff, + pairwise::{Anisotropic, ApproximateShapeOverlap, HardShape}, + univariate::OverlapPenalty, +}; +use hoomd_mc::{ + QuickCompress, QuickInsert, Rotate, Sweep, Translate, Trial, Tune, + TuneOptions, UniformIn, +}; +use hoomd_microstate::{ + AppendMicrostate, Microstate, SiteKey, Transform, + boundary::{GenerateGhosts, MAX_GHOSTS, Periodic, Wrap}, + property::{Orientation, OrientedPoint, Point, Position}, +}; +use hoomd_simulation::{Simulation, macrostate::Isothermal}; +use hoomd_spatial::VecCell; +use hoomd_vector::{self, Cartesian, Rotate as _, Rotation, Versor}; + +type BodyProperties = OrientedPoint, Versor>; + +#[derive(Clone, Copy, Default, Position, Orientation)] +struct SiteProperties { + position: Cartesian<3>, + orientation: Versor, +} + +#[derive(Clone, PartialEq)] +struct Boundary(Periodic); + +impl Transform for OrientedPoint, Versor> { + fn transform(&self, site_properties: &SiteProperties) -> SiteProperties { + let lifted_body_position = + Cartesian::from([self.position[0], self.position[1], 0.0]); + + SiteProperties { + position: lifted_body_position + + self.orientation.rotate(&site_properties.position), + orientation: self.orientation.combine(&site_properties.orientation), + } + } +} + +impl Wrap for Boundary { + fn wrap( + &self, + properties: SiteProperties, + ) -> Result { + let wrapped_projection = self.0.wrap(Point { + position: Cartesian::from([ + properties.position[0], + properties.position[1], + ]), + })?; + Ok(SiteProperties { + position: Cartesian::from([ + wrapped_projection.position[0], + wrapped_projection.position[1], + properties.position[2], + ]), + ..properties + }) + } +} + +impl Wrap for Boundary { + fn wrap( + &self, + properties: BodyProperties, + ) -> Result { + self.0.wrap(properties) + } +} + +impl GenerateGhosts for Boundary { + fn maximum_interaction_range(&self) -> f64 { + self.0.maximum_interaction_range() + } + + fn generate_ghosts( + &self, + site_properties: &SiteProperties, + ) -> ArrayVec { + let projected_site = Point { + position: Cartesian::from([ + site_properties.position[0], + site_properties.position[1], + ]), + }; + let projected_ghosts = self.0.generate_ghosts(&projected_site); + projected_ghosts + .iter() + .map(|s| SiteProperties { + position: Cartesian::from([ + s.position[0], + s.position[1], + site_properties.position[2], + ]), + ..*site_properties + }) + .collect() + } +} + +impl Volume for Boundary { + fn volume(&self) -> f64 { + self.0.volume() + } +} + +impl MapPoint> for Boundary { + fn map_point( + &self, + point: Cartesian<2>, + other: &Self, + ) -> Result, hoomd_geometry::Error> { + self.0.map_point(point, &other.0) + } +} + +impl Scale for Boundary { + fn scale_length(&self, v: hoomd_utility::valid::PositiveReal) -> Self { + Boundary(self.0.scale_length(v)) + } + + fn scale_volume(&self, v: hoomd_utility::valid::PositiveReal) -> Self { + Boundary(self.0.scale_volume(v)) + } +} + +impl Distribution> for Boundary { + fn sample(&self, rng: &mut R) -> Cartesian<2> { + self.0.sample(rng) + } +} + +impl Quasi2dCapsuleSelfAssembly { + /// Construct a new hard tetrahedron self-assembly simulation. + fn new() -> anyhow::Result { + let initial_number_density = 0.12; + let target_number_density = 0.22; + let n_bodies = 256; + let maximum_distance = 0.04; + let maximum_rotation = 0.04; + let macrostate = Isothermal { temperature: 1.0 }; + + let capsule = Capsule { + radius: 1.0.try_into()?, + height: 5.0.try_into()?, + }; + let hamiltonian = PairwiseCutoff(HardShape(capsule.clone())); + + let initial_box_volume = n_bodies as f64 / initial_number_density; + let initial_box_edge_length = initial_box_volume.sqrt(); + let rectangle = + Rectangle::with_equal_edges(initial_box_edge_length.try_into()?); + let periodic_rectangle = Boundary(Periodic::new( + hamiltonian.maximum_interaction_range(), + rectangle, + )?); + + let vec_cell = VecCell::builder() + .nominal_search_radius( + hamiltonian.maximum_interaction_range().try_into()?, + ) + .build(); + let microstate = Microstate::builder() + .boundary(periodic_rectangle) + .spatial_data(vec_cell) + .try_build()?; + + let translate = + Translate::with_maximum_distance(maximum_distance.try_into()?); + let translate_sweep = Sweep(translate); + + let rotate = + Rotate::with_maximum_rotation(maximum_rotation.try_into()?); + let rotate_sweep = Sweep(rotate); + + let distribution = UniformIn { + boundary: microstate.boundary().clone(), + template_sites: vec![SiteProperties::default()], + }; + let quick_insert = QuickInsert::new(distribution, n_bodies); + + let target_box_volume = n_bodies as f64 / target_number_density; + let quick_compress = + QuickCompress::with_target_volume(target_box_volume.try_into()?); + + let approximate_shape_overlap = Anisotropic { + interaction: ApproximateShapeOverlap::new( + Convex(capsule), + OverlapPenalty::default(), + 0.01.try_into()?, + ), + r_cut: hamiltonian.maximum_interaction_range(), + }; + + let overlap_penalty_hamiltonian = + PairwiseCutoff(approximate_shape_overlap); + + Ok(Quasi2dCapsuleSelfAssembly { + microstate, + overlap_penalty_hamiltonian, + hamiltonian, + translate_sweep, + rotate_sweep, + quick_compress, + quick_insert, + macrostate, + phase: Phase::Initialize, + }) + } +} + +#[cfg_attr(feature = "bevy", derive(Resource))] +struct Quasi2dCapsuleSelfAssembly { + /// Positions and orientations of all the bodies in the simulation. + microstate: Microstate< + BodyProperties, + SiteProperties, + VecCell, + Boundary, + >, + /// How sites interact with other sites and fields. + hamiltonian: PairwiseCutoff>>, + /// Trial moves to apply. + translate_sweep: Sweep>>, + /// Trial moves to apply. + rotate_sweep: Sweep>, + /// Temperature set point. + macrostate: Isothermal, + /// Quick compress algorithm. + quick_compress: QuickCompress, + /// Quick insert algorithm. + quick_insert: QuickInsert>, + /// How sites interact when inserted and compressed. + overlap_penalty_hamiltonian: PairwiseCutoff< + Anisotropic< + ApproximateShapeOverlap>>, + >, + >, + /// The current phase of the simulation. + phase: Phase, +} + +enum Phase { + Initialize, + Equilibrate, +} + +impl Simulation for Quasi2dCapsuleSelfAssembly { + /// 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 Quasi2dCapsuleSelfAssembly { + fn initialize(&mut self) -> anyhow::Result<()> { + if self.quick_insert.is_complete() { + self.quick_compress.apply( + &mut self.microstate, + &self.overlap_penalty_hamiltonian, + |_| true, + ); + } else { + self.quick_insert + .apply(&mut self.microstate, &self.overlap_penalty_hamiltonian); + } + + self.translate_sweep.apply( + &mut self.microstate, + &self.overlap_penalty_hamiltonian, + &Isothermal { temperature: 1.0 }, + ); + + self.rotate_sweep.apply( + &mut self.microstate, + &self.overlap_penalty_hamiltonian, + &Isothermal { temperature: 1.0 }, + ); + + if self.quick_compress.is_complete() { + self.translate_sweep.tune_with_options( + &self.microstate, + &self.hamiltonian, + &self.macrostate, + &TuneOptions::default(), + ); + self.rotate_sweep.tune_with_options( + &self.microstate, + &self.hamiltonian, + &self.macrostate, + &TuneOptions::default(), + ); + + self.phase = Phase::Equilibrate; + println!( + "Initialization complete at step {}.", + self.microstate.step() + ); + } + + if self.step() >= 20_000 { + let n = self.microstate.bodies().len(); + let target_n = self.quick_insert.target(); + let volume = self.microstate.boundary().0.volume(); + let target_volume = self.quick_compress.target_volume(); + return Err(anyhow!( + "inserted {n}/{target_n} bodies and compressed to {volume} / {target_volume}" + )); + } + + Ok(()) + } + + fn equilibrate(&mut self) { + self.translate_sweep.apply( + &mut self.microstate, + &self.hamiltonian, + &self.macrostate, + ); + + self.rotate_sweep.apply( + &mut self.microstate, + &self.hamiltonian, + &self.macrostate, + ); + } +} + +impl AppendMicrostate + for HoomdGsdFile +{ + #[inline] + fn append_microstate( + &mut self, + microstate: &Microstate, + ) -> Result, hoomd_gsd::hoomd::AppendError> + { + let edge_lengths = microstate.boundary().0.shape().edge_lengths; + + self.append_frame(microstate.step())? + .configuration_box([ + edge_lengths[0].get(), + edge_lengths[1].get(), + 5.0, + 0.0, + 0.0, + 0.0, + ])? + .configuration_dimensions(Dimensions::Three)? + .particles_orientation( + microstate + .iter_sites_tag_order() + .map(|s| s.properties.orientation), + )? + .particles_position( + microstate + .iter_sites_tag_order() + .map(|s| s.properties.position), + ) + } +} + +// Remove the cfg(not(...)) line when using this code outside the hoomd-rs/examples directory. +#[cfg(not(feature = "bevy"))] +fn main() -> anyhow::Result<()> { + use hoomd_gsd::hoomd::HoomdGsdFile; + use hoomd_microstate::AppendMicrostate; + + let mut simulation = Quasi2dCapsuleSelfAssembly::new()?; + let mut hoomd_gsd_file = HoomdGsdFile::create("capsules-on-a-plane.gsd")?; + + for _ in 0..40_000 { + simulation.advance()?; + + if simulation.step().is_multiple_of(10_000) { + hoomd_gsd_file + .append_microstate(&simulation.microstate)? + .end()?; + } + } + + Ok(()) +} +// ANCHOR_END: all + +#[cfg(feature = "bevy")] +mod capsules_on_a_plane_interactive; +#[cfg(feature = "bevy")] +use bevy::prelude::Resource; +#[cfg(feature = "bevy")] +use capsules_on_a_plane_interactive::main; +use rand::{Rng, distr::Distribution}; diff --git a/examples/mc-examples/capsules_on_a_plane_interactive.rs b/examples/mc-examples/capsules_on_a_plane_interactive.rs new file mode 100644 index 00000000..dde1a8d3 --- /dev/null +++ b/examples/mc-examples/capsules_on_a_plane_interactive.rs @@ -0,0 +1,106 @@ +use std::f64::consts::PI; + +use hoomd_bevy::{ + AdvanceSet, HoomdBevyPlugin, InitialCamera, PRIMARY_COLOR_3D, Settings, + representation::surface_mesh, +}; + +use anyhow::Context; +use bevy::prelude::*; +use bevy_egui::EguiPlugin; +use hoomd_vector::{Rotation, Versor}; + +use super::Quasi2dCapsuleSelfAssembly; + +/// Mark the tetrahedron representation type. +struct A; + +pub(crate) fn main() -> anyhow::Result<()> { + let simulation = Quasi2dCapsuleSelfAssembly::new() + .context("failed to setup simulation")?; + + let l = + simulation.microstate.boundary().0.shape().edge_lengths[1].get() as f32; + let hoomd_bevy_plugin = HoomdBevyPlugin { + initial_settings: Settings { + camera: InitialCamera::Orthographic3d(l + 1.0), + ..default() + }, + simulation, + }; + + let mut app = App::new(); + hoomd_bevy::add_default_plugins(&mut app); + app.add_plugins(EguiPlugin::default()); + hoomd_bevy_plugin.build(&mut app); + + let capsule = Capsule3d { + radius: 1.0, + half_length: 2.5, + }; + let capsule_material = StandardMaterial { + base_color: PRIMARY_COLOR_3D, + perceptual_roughness: 0.2, + ..default() + }; + + app.add_systems( + Startup, + (move || (capsule.mesh().build(), capsule_material.clone())) + .pipe(surface_mesh::SurfaceMesh::::setup), + ); + + app.add_systems( + Update, + (sync_sites,) + .run_if(resource_changed::) + .after(AdvanceSet), + ); + + app.run(); + + Ok(()) +} + +/// Copy the current positions of simulation sites to bevy entities. +fn sync_sites( + mut commands: Commands, + site_representation: Res>, + site_query: Query< + (Entity, &mut Transform), + With>, + >, + simulation: Res, +) { + let sites = simulation.microstate.sites(); + + surface_mesh::SurfaceMesh::sync( + &mut commands, + site_representation, + site_query, + sites.iter().map(|site| { + let orientation = site.properties.orientation; + let rotation = Versor::from_axis_angle( + [1.0, 0.0, 0.0] + .try_into() + .expect("hard-coded vector should be non-zero length"), + PI / 2.0, + ); + let orientation = orientation.combine(&rotation); + + ( + Vec3::new( + site.properties.position[0] as f32, + site.properties.position[1] as f32, + site.properties.position[2] as f32, + ), + Quat::from_xyzw( + orientation.get().vector[0] as f32, + orientation.get().vector[1] as f32, + orientation.get().vector[2] as f32, + orientation.get().scalar as f32, + ), + ) + }), + ); +} diff --git a/examples/mc-tutorial/hard_tetrahedron_self_assembly_interactive.rs b/examples/mc-tutorial/hard_tetrahedron_self_assembly_interactive.rs index a18a4b00..b79a6e0d 100644 --- a/examples/mc-tutorial/hard_tetrahedron_self_assembly_interactive.rs +++ b/examples/mc-tutorial/hard_tetrahedron_self_assembly_interactive.rs @@ -56,7 +56,11 @@ pub(crate) fn main() -> anyhow::Result<()> { ), ], }; - let tetrahedron_material = StandardMaterial::from(PRIMARY_COLOR_3D); + let tetrahedron_material = StandardMaterial { + base_color: PRIMARY_COLOR_3D, + perceptual_roughness: 0.2, + ..default() + }; app.add_systems( Startup, diff --git a/hoomd-bevy/src/lib.rs b/hoomd-bevy/src/lib.rs index 559043c2..9a618a5e 100644 --- a/hoomd-bevy/src/lib.rs +++ b/hoomd-bevy/src/lib.rs @@ -205,7 +205,8 @@ pub enum InitialCamera { /// automatically based on the window dimensions. /// /// Controls: - /// * TODO + /// * Left click and drag to rotate. + /// * Scroll to zoom. Orthographic3d(f32), } @@ -254,6 +255,22 @@ pub struct CameraControl2d { dragging: bool, } +/// Settings used by the 3d camera controls. +#[derive(Debug, Default, Resource)] +pub struct CameraControl3d { + /// Coordinates clicked in the window. + last_click_position: Vec2, + + /// Initial camera transform when clicked + last_camera_transform: Transform, + + /// Initial light transform when clicked + last_light_transform: Transform, + + /// Track whether the user is dragging the view. + dragging: bool, +} + /// The overlay UI root node. #[derive(Component)] struct OverlayRoot; @@ -575,10 +592,6 @@ where } /// Left click and drag to pan the 2D camera. - /// - /// # Panics - /// - /// Panics when the 2D camera viewport is invalid. fn camera_mouse_pan_control_2d( camera: Single< (&Camera, &GlobalTransform, &mut Transform, &mut Projection), @@ -679,6 +692,97 @@ where } } + /// Zoom the 3d camera using the mouse wheel or trackpad scroll gesture. + fn camera_mouse_zoom_control_3d( + time: Res