hoomd_md/
update_net_force.rs

1// Copyright (c) 2024-2026 The Regents of the University of Michigan.
2// Part of hoomd-rs, released under the BSD 3-Clause License.
3
4//! Define `UpdateNetForce`
5
6use rayon::prelude::*;
7
8use hoomd_interaction::{NetBodyForceAndVirial, NetBodyForceVirialAndTorque};
9use hoomd_microstate::{
10    Microstate,
11    property::{NetForce, NetTorque, NetVirial},
12};
13use hoomd_vector::{Outer, Vector, Wedge};
14
15/// Compute the net force and virial given by an interaction model and apply it
16/// to each body in the microstate.
17///
18/// Given an interaction model that implements [`NetBodyForceAndVirial`],
19/// [`UpdateNetForceAndVirial`] sets the [`NetForce`] and [`NetVirial`]
20/// properties of each body in the microstate to the one computed by the
21/// interaction model.
22///
23/// [`NetBodyForceAndVirial`]: hoomd_interaction::NetBodyForceAndVirial
24/// [`NetForce`]: hoomd_microstate::property::NetForce
25/// [`NetVirial`]: hoomd_microstate::property::NetVirial
26///
27/// # Example
28/// ```
29/// use hoomd_interaction::{
30///     PairwiseCutoff, Rigid, pairwise::Isotropic, univariate::LennardJones,
31/// };
32/// use hoomd_md::UpdateNetForceAndVirial;
33/// use hoomd_microstate::{
34///     Body, Microstate,
35///     property::{DynamicPoint, Point},
36/// };
37/// use hoomd_vector::Cartesian;
38///
39/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
40/// let mut microstate = Microstate::builder()
41///     .bodies([
42///         Body::single_site(DynamicPoint::default(), Point::default()),
43///         Body::single_site(
44///             DynamicPoint {
45///                 position: Cartesian::<2>::from([2.0, 0.0]),
46///                 ..Default::default()
47///             },
48///             Point::default(),
49///         ),
50///     ])
51///     .try_build()?;
52///
53/// let lennard_jones = LennardJones::<12, 6>::default();
54/// let pairwise_cutoff = PairwiseCutoff(Isotropic {
55///     interaction: lennard_jones,
56///     r_cut: 2.5,
57/// });
58/// let rigid = Rigid(pairwise_cutoff);
59///
60/// microstate.update_net_force_and_virial(&rigid);
61/// #   Ok(())
62/// # }
63/// ```
64pub trait UpdateNetForceAndVirial<E> {
65    /// Compute and set the net force and virial on each body.
66    fn update_net_force_and_virial(&mut self, interaction_model: &E);
67}
68
69/// Compute the net force, virial, and torque given by an interaction model and
70/// apply them to each body in the microstate.
71///
72/// Given an interaction model that implements [`NetBodyForceVirialAndTorque`],
73/// [`UpdateNetForceVirialAndTorque`] sets the [`NetForce`], [`NetVirial`] and
74/// [`NetTorque`] properties of each body in the microstate to the ones computed
75/// by the interaction model.
76///
77/// [`NetBodyForceVirialAndTorque`]: hoomd_interaction::NetBodyForceVirialAndTorque
78/// [`NetForce`]: hoomd_microstate::property::NetForce
79/// [`NetVirial`]: hoomd_microstate::property::NetVirial
80/// [`NetTorque`]: hoomd_microstate::property::NetTorque
81///
82/// # Example
83///
84/// ```
85/// use hoomd_interaction::{
86///     PairwiseCutoff, Rigid, pairwise::Isotropic, univariate::LennardJones,
87/// };
88/// use hoomd_md::UpdateNetForceVirialAndTorque;
89/// use hoomd_microstate::{
90///     Body, Microstate,
91///     property::{DynamicOrientedPoint, Point},
92/// };
93/// use hoomd_vector::{Angle, Cartesian};
94///
95/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
96/// let mut microstate: Microstate<
97///     DynamicOrientedPoint<Cartesian<2>, Angle>,
98///     Point<Cartesian<2>>,
99///     _,
100///     _,
101/// > = Microstate::builder()
102///     .bodies([
103///         Body::single_site(
104///             DynamicOrientedPoint {
105///                 position: Cartesian::<2>::from([0.0, -1.0]),
106///                 ..Default::default()
107///             },
108///             Point::new([0.0, 1.0].into()),
109///         ),
110///         Body::single_site(
111///             DynamicOrientedPoint {
112///                 position: Cartesian::<2>::from([2.0, -2.0]),
113///                 ..Default::default()
114///             },
115///             Point::new([0.0, 2.0].into()),
116///         ),
117///     ])
118///     .try_build()?;
119///
120/// let lennard_jones = LennardJones::<12, 6>::default();
121/// let pairwise_cutoff = PairwiseCutoff(Isotropic {
122///     interaction: lennard_jones,
123///     r_cut: 2.5,
124/// });
125/// let rigid = Rigid(pairwise_cutoff);
126///
127/// microstate.update_net_force_virial_and_torque(&rigid);
128/// #   Ok(())
129/// # }
130/// ```
131pub trait UpdateNetForceVirialAndTorque<E> {
132    /// Compute and set the net force, virial, and torque on each body.
133    fn update_net_force_virial_and_torque(&mut self, interaction_model: &E);
134}
135
136impl<V, B, S, X, C, E> UpdateNetForceAndVirial<E> for Microstate<B, S, X, C>
137where
138    V: Default + Vector + Outer + Send,
139    V::Tensor: Copy + Send,
140    B: NetForce<NetForce = V> + NetVirial<NetVirial = V::Tensor> + Sync,
141    S: Sync,
142    X: Sync,
143    C: Sync,
144    E: NetBodyForceAndVirial<B, S, X, C, Force = V> + Sync,
145{
146    #[inline]
147    fn update_net_force_and_virial(&mut self, interaction_model: &E) {
148        let mut net_force_and_virial_tmp = Vec::new();
149
150        (0..self.bodies().len())
151            .into_par_iter()
152            .map(|body_index| interaction_model.net_body_force_and_virial(self, body_index))
153            .collect_into_vec(&mut net_force_and_virial_tmp);
154
155        for (body_index, (net_force, net_virial)) in net_force_and_virial_tmp.iter().enumerate() {
156            self.set_body_net_force(body_index, *net_force);
157            self.set_body_net_virial(body_index, *net_virial);
158        }
159    }
160}
161
162impl<V, B, S, X, C, E> UpdateNetForceVirialAndTorque<E> for Microstate<B, S, X, C>
163where
164    V: Default + Vector + Wedge + Outer + Send,
165    V::Tensor: Copy + Send,
166    V::Bivector: Copy + Send,
167    B: NetForce<NetForce = V>
168        + NetVirial<NetVirial = V::Tensor>
169        + NetTorque<NetTorque = V::Bivector>
170        + Sync,
171    S: Sync,
172    X: Sync,
173    C: Sync,
174    E: NetBodyForceVirialAndTorque<B, S, X, C, Force = V> + Sync,
175{
176    #[inline]
177    fn update_net_force_virial_and_torque(&mut self, interaction_model: &E) {
178        let mut net_force_virial_and_torque_tmp = Vec::new();
179
180        (0..self.bodies().len())
181            .into_par_iter()
182            .map(|body_index| interaction_model.net_body_force_virial_and_torque(self, body_index))
183            .collect_into_vec(&mut net_force_virial_and_torque_tmp);
184
185        for (body_index, (net_force, net_virial, net_torque)) in
186            net_force_virial_and_torque_tmp.iter().enumerate()
187        {
188            self.set_body_net_force(body_index, *net_force);
189            self.set_body_net_virial(body_index, *net_virial);
190            self.set_body_net_torque(body_index, *net_torque);
191        }
192    }
193}
194
195#[cfg(test)]
196mod test {
197    use super::*;
198    use approxim::assert_relative_eq;
199    use assert2::check;
200
201    use hoomd_interaction::{PairwiseCutoff, Rigid, pairwise::Isotropic, univariate::LennardJones};
202    use hoomd_microstate::{
203        Body,
204        property::{DynamicOrientedPoint, DynamicPoint, Point},
205    };
206    use hoomd_vector::{Angle, Cartesian, Versor};
207
208    // TODO: add virial tests
209
210    #[test]
211    fn net_force_2d() -> anyhow::Result<()> {
212        let mut microstate = Microstate::builder()
213            .bodies([
214                Body::single_site(DynamicPoint::default(), Point::default()),
215                Body::single_site(
216                    DynamicPoint {
217                        position: Cartesian::<2>::from([2.0, 0.0]),
218                        ..Default::default()
219                    },
220                    Point::default(),
221                ),
222            ])
223            .try_build()?;
224
225        let lennard_jones = LennardJones::<12, 6>::default();
226        let pairwise_cutoff = PairwiseCutoff(Isotropic {
227            interaction: lennard_jones,
228            r_cut: 2.5,
229        });
230        let rigid = Rigid(pairwise_cutoff);
231
232        check!(microstate.bodies()[0].item.properties.net_force == [0.0, 0.0].into());
233        check!(microstate.bodies()[1].item.properties.net_force == [0.0, 0.0].into());
234
235        microstate.update_net_force_and_virial(&rigid);
236
237        assert_relative_eq!(
238            microstate.bodies()[0].item.properties.net_force,
239            [93.0 / 512.0, 0.0].into()
240        );
241        assert_relative_eq!(
242            microstate.bodies()[1].item.properties.net_force,
243            [-93.0 / 512.0, 0.0].into()
244        );
245
246        Ok(())
247    }
248
249    #[test]
250    fn net_force_3d() -> anyhow::Result<()> {
251        let mut microstate = Microstate::builder()
252            .bodies([
253                Body::single_site(DynamicPoint::default(), Point::default()),
254                Body::single_site(
255                    DynamicPoint {
256                        position: Cartesian::<3>::from([2.0, 0.0, 0.0]),
257                        ..Default::default()
258                    },
259                    Point::default(),
260                ),
261            ])
262            .try_build()?;
263
264        let lennard_jones = LennardJones::<12, 6>::default();
265        let pairwise_cutoff = PairwiseCutoff(Isotropic {
266            interaction: lennard_jones,
267            r_cut: 2.5,
268        });
269        let rigid = Rigid(pairwise_cutoff);
270
271        check!(microstate.bodies()[0].item.properties.net_force == [0.0, 0.0, 0.0].into());
272        check!(microstate.bodies()[1].item.properties.net_force == [0.0, 0.0, 0.0].into());
273
274        microstate.update_net_force_and_virial(&rigid);
275
276        assert_relative_eq!(
277            microstate.bodies()[0].item.properties.net_force,
278            [93.0 / 512.0, 0.0, 0.0].into()
279        );
280        assert_relative_eq!(
281            microstate.bodies()[1].item.properties.net_force,
282            [-93.0 / 512.0, 0.0, 0.0].into()
283        );
284
285        Ok(())
286    }
287
288    #[test]
289    fn net_force_and_torque_2d() -> anyhow::Result<()> {
290        let mut microstate: Microstate<
291            DynamicOrientedPoint<Cartesian<2>, Angle>,
292            Point<Cartesian<2>>,
293            _,
294            _,
295        > = Microstate::builder()
296            .bodies([
297                Body::single_site(
298                    DynamicOrientedPoint {
299                        position: Cartesian::<2>::from([0.0, -1.0]),
300                        ..Default::default()
301                    },
302                    Point::new([0.0, 1.0].into()),
303                ),
304                Body::single_site(
305                    DynamicOrientedPoint {
306                        position: Cartesian::<2>::from([2.0, -2.0]),
307                        ..Default::default()
308                    },
309                    Point::new([0.0, 2.0].into()),
310                ),
311            ])
312            .try_build()?;
313
314        let lennard_jones = LennardJones::<12, 6>::default();
315        let pairwise_cutoff = PairwiseCutoff(Isotropic {
316            interaction: lennard_jones,
317            r_cut: 2.5,
318        });
319        let rigid = Rigid(pairwise_cutoff);
320
321        check!(microstate.bodies()[0].item.properties.net_force == [0.0, 0.0].into());
322        check!(microstate.bodies()[1].item.properties.net_force == [0.0, 0.0].into());
323        check!(microstate.bodies()[0].item.properties.net_torque == 0.0);
324        check!(microstate.bodies()[1].item.properties.net_torque == 0.0);
325
326        microstate.update_net_force_virial_and_torque(&rigid);
327
328        assert_relative_eq!(
329            microstate.bodies()[0].item.properties.net_force,
330            [93.0 / 512.0, 0.0].into()
331        );
332        assert_relative_eq!(
333            microstate.bodies()[1].item.properties.net_force,
334            [-93.0 / 512.0, 0.0].into()
335        );
336
337        assert_relative_eq!(
338            microstate.bodies()[0].item.properties.net_torque,
339            -93.0 / 512.0
340        );
341        assert_relative_eq!(
342            microstate.bodies()[1].item.properties.net_torque,
343            2.0 * 93.0 / 512.0
344        );
345
346        Ok(())
347    }
348
349    #[test]
350    fn net_force_and_torque_3d() -> anyhow::Result<()> {
351        let mut microstate: Microstate<
352            DynamicOrientedPoint<Cartesian<3>, Versor>,
353            Point<Cartesian<3>>,
354            _,
355            _,
356        > = Microstate::builder()
357            .bodies([
358                Body::single_site(
359                    DynamicOrientedPoint {
360                        position: Cartesian::<3>::from([0.0, -1.0, 0.0]),
361                        ..Default::default()
362                    },
363                    Point::new([0.0, 1.0, 0.0].into()),
364                ),
365                Body::single_site(
366                    DynamicOrientedPoint {
367                        position: Cartesian::<3>::from([2.0, -2.0, 0.0]),
368                        ..Default::default()
369                    },
370                    Point::new([0.0, 2.0, 0.0].into()),
371                ),
372            ])
373            .try_build()?;
374
375        let lennard_jones = LennardJones::<12, 6>::default();
376        let pairwise_cutoff = PairwiseCutoff(Isotropic {
377            interaction: lennard_jones,
378            r_cut: 2.5,
379        });
380        let rigid = Rigid(pairwise_cutoff);
381
382        check!(microstate.bodies()[0].item.properties.net_force == [0.0, 0.0, 0.0].into());
383        check!(microstate.bodies()[1].item.properties.net_force == [0.0, 0.0, 0.0].into());
384        check!(microstate.bodies()[0].item.properties.net_torque == [0.0, 0.0, 0.0].into());
385        check!(microstate.bodies()[1].item.properties.net_torque == [0.0, 0.0, 0.0].into());
386
387        microstate.update_net_force_virial_and_torque(&rigid);
388
389        assert_relative_eq!(
390            microstate.bodies()[0].item.properties.net_force,
391            [93.0 / 512.0, 0.0, 0.0].into()
392        );
393        assert_relative_eq!(
394            microstate.bodies()[1].item.properties.net_force,
395            [-93.0 / 512.0, 0.0, 0.0].into()
396        );
397
398        assert_relative_eq!(
399            microstate.bodies()[0].item.properties.net_torque,
400            [0.0, 0.0, -93.0 / 512.0].into()
401        );
402        assert_relative_eq!(
403            microstate.bodies()[1].item.properties.net_torque,
404            [0.0, 0.0, 2.0 * 93.0 / 512.0].into()
405        );
406
407        Ok(())
408    }
409}
hoomd_md/update_net_force.rs

hoomd_md/
update_net_force.rs
