hoomd_md/compute/
rotational_kinetic_energy.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//! Implement `RotationalKineticEnergy`
5
6use super::RotationalKineticEnergy;
7use hoomd_microstate::{Body, Microstate, Tagged, property::DynamicOrientedPoint};
8use hoomd_vector::{Angle, Outer, Versor, Wedge};
9
10impl<P, S, X, C> RotationalKineticEnergy<DynamicOrientedPoint<P, Angle>, S>
11    for Microstate<DynamicOrientedPoint<P, Angle>, S, X, C>
12where
13    P: Wedge + Outer,
14{
15    #[inline]
16    fn rotational_kinetic_energy_with_filter<
17        F: Fn(&Tagged<Body<DynamicOrientedPoint<P, Angle>, S>>) -> bool,
18    >(
19        &self,
20        should_sum_body: F,
21    ) -> (f64, usize) {
22        self.bodies()
23            .iter()
24            .filter(|&body| should_sum_body(body))
25            .fold((0.0, 0), |(total, count), body| {
26                let moment_of_inertia = body.item.properties.moment_of_inertia;
27                let angular_momentum = body.item.properties.angular_momentum;
28
29                if moment_of_inertia == 0.0 {
30                    (total, count)
31                } else {
32                    (
33                        total + angular_momentum.powi(2) / (2.0 * moment_of_inertia),
34                        count + 1,
35                    )
36                }
37            })
38    }
39}
40
41impl<P, S, X, C> RotationalKineticEnergy<DynamicOrientedPoint<P, Versor>, S>
42    for Microstate<DynamicOrientedPoint<P, Versor>, S, X, C>
43where
44    P: Wedge + Outer,
45{
46    #[inline]
47    fn rotational_kinetic_energy_with_filter<
48        F: Fn(&Tagged<Body<DynamicOrientedPoint<P, Versor>, S>>) -> bool,
49    >(
50        &self,
51        should_sum_body: F,
52    ) -> (f64, usize) {
53        self.bodies()
54            .iter()
55            .filter(|&body| should_sum_body(body))
56            .fold((0.0, 0), |(mut total, mut count), body| {
57                let moment_of_inertia = body.item.properties.moment_of_inertia;
58                let angular_momentum = body.item.properties.angular_momentum;
59
60                for (momentum, inertia) in
61                    angular_momentum.coordinates.iter().zip(moment_of_inertia)
62                {
63                    if inertia != 0.0 {
64                        total += momentum.powi(2) / (2.0 * inertia);
65                        count += 1;
66                    }
67                }
68
69                (total, count)
70            })
71    }
72}
73
74#[cfg(test)]
75mod test {
76    use super::*;
77    use approxim::assert_relative_eq;
78    use assert2::check;
79
80    use hoomd_microstate::{
81        Body,
82        property::{DynamicOrientedPoint, Point},
83    };
84    use hoomd_vector::{Angle, Cartesian, Versor};
85
86    #[test]
87    fn kinetic_energy_2d() -> anyhow::Result<()> {
88        let microstate: Microstate<DynamicOrientedPoint<Cartesian<2>, Angle>, _, _, _> =
89            Microstate::builder()
90                .bodies([
91                    Body::single_site(DynamicOrientedPoint::default(), Point::default()),
92                    Body::single_site(
93                        DynamicOrientedPoint {
94                            moment_of_inertia: 0.0,
95                            ..Default::default()
96                        },
97                        Point::default(),
98                    ),
99                    Body::single_site(
100                        DynamicOrientedPoint {
101                            moment_of_inertia: 2.0,
102                            angular_momentum: 8.0,
103                            ..Default::default()
104                        },
105                        Point::default(),
106                    ),
107                    Body::single_site(
108                        DynamicOrientedPoint {
109                            moment_of_inertia: 4.0,
110                            angular_momentum: 3.0,
111                            ..Default::default()
112                        },
113                        Point::default(),
114                    ),
115                    Body::single_site(
116                        DynamicOrientedPoint {
117                            moment_of_inertia: 3.0,
118                            angular_momentum: 2.0,
119                            ..Default::default()
120                        },
121                        Point::default(),
122                    ),
123                ])
124                .try_build()?;
125
126        let (total_kinetic_energy, total_degrees_of_freedom) =
127            microstate.rotational_kinetic_energy();
128        check!(total_degrees_of_freedom == 4);
129        assert_relative_eq!(total_kinetic_energy, 64.0 / 4.0 + 9.0 / 8.0 + 4.0 / 6.0);
130
131        let (filtered_kinetic_energy, filtered_degrees_of_freedom) =
132            microstate.rotational_kinetic_energy_with_filter(|b| b.tag <= 2);
133        check!(filtered_degrees_of_freedom == 2);
134        assert_relative_eq!(filtered_kinetic_energy, 64.0 / 4.0);
135
136        Ok(())
137    }
138
139    #[test]
140    fn kinetic_energy_3d() -> anyhow::Result<()> {
141        let microstate: Microstate<DynamicOrientedPoint<Cartesian<3>, Versor>, _, _, _> =
142            Microstate::builder()
143                .bodies([
144                    Body::single_site(DynamicOrientedPoint::default(), Point::default()),
145                    Body::single_site(
146                        DynamicOrientedPoint {
147                            moment_of_inertia: [0.0, 0.0, 0.0],
148                            angular_momentum: [1.0, 1.0, 1.0].into(),
149                            ..Default::default()
150                        },
151                        Point::default(),
152                    ),
153                    Body::single_site(
154                        DynamicOrientedPoint {
155                            moment_of_inertia: [2.0, 0.0, 0.0],
156                            angular_momentum: [8.0, 1.0, 1.0].into(),
157                            ..Default::default()
158                        },
159                        Point::default(),
160                    ),
161                    Body::single_site(
162                        DynamicOrientedPoint {
163                            moment_of_inertia: [0.0, 6.0, 0.0],
164                            angular_momentum: [1.0, 3.0, 1.0].into(),
165                            ..Default::default()
166                        },
167                        Point::default(),
168                    ),
169                    Body::single_site(
170                        DynamicOrientedPoint {
171                            moment_of_inertia: [0.0, 0.0, 3.0],
172                            angular_momentum: [1.0, 1.0, -4.0].into(),
173                            ..Default::default()
174                        },
175                        Point::default(),
176                    ),
177                    Body::single_site(
178                        DynamicOrientedPoint {
179                            moment_of_inertia: [2.0, 4.0, 6.0],
180                            angular_momentum: [3.0, 2.0, -4.0].into(),
181                            ..Default::default()
182                        },
183                        Point::default(),
184                    ),
185                ])
186                .try_build()?;
187
188        let (total_kinetic_energy, total_degrees_of_freedom) =
189            microstate.rotational_kinetic_energy();
190        check!(total_degrees_of_freedom == 9);
191        assert_relative_eq!(
192            total_kinetic_energy,
193            64.0 / 4.0 + 9.0 / 12.0 + 16.0 / 6.0 + 9.0 / 4.0 + 4.0 / 8.0 + 16.0 / 12.0
194        );
195
196        let (filtered_kinetic_energy, filtered_degrees_of_freedom) =
197            microstate.rotational_kinetic_energy_with_filter(|b| b.tag <= 2);
198        check!(filtered_degrees_of_freedom == 4);
199        assert_relative_eq!(filtered_kinetic_energy, 64.0 / 4.0);
200
201        Ok(())
202    }
203}
hoomd_md/compute/rotational_kinetic_energy.rs

hoomd_md/compute/
rotational_kinetic_energy.rs
