hoomd_md/modify/
zero_center_angular_momentum.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 `ZeroCenterAngularMomentum`
5
6use super::ZeroCenterAngularMomentum;
7use hoomd_linear_algebra::{GeneralMatrix, MatMul, matrix::Matrix};
8use hoomd_microstate::{
9    Body, Microstate, SiteKey, Tagged, Transform,
10    boundary::{GenerateGhosts, Wrap},
11    property::{
12        DynamicOrientedPoint, DynamicPoint, Mass, Momentum, Position, RotationalMotionTypes,
13    },
14};
15use hoomd_spatial::PointUpdate;
16use hoomd_vector::{Cartesian, InnerProduct, Outer, Wedge};
17
18/// Zero a 3D microstate's angular momentum.
19#[inline]
20fn zero_angular_momentum_3d<B, S, X, C, F>(
21    microstate: &mut Microstate<B, S, X, C>,
22    should_zero_body: F,
23) where
24    B: Position<Position = Cartesian<3>>
25        + Mass
26        + Momentum<Momentum = Cartesian<3>>
27        + Transform<S>
28        + Clone,
29    S: Position<Position = Cartesian<3>> + Default,
30    X: PointUpdate<Cartesian<3>, SiteKey>,
31    C: Wrap<B> + Wrap<S> + GenerateGhosts<S>,
32    F: Fn(&Tagged<Body<B, S>>) -> bool,
33{
34    let mut center_of_mass = Cartesian::default();
35    let mut total_mass = 0.0;
36
37    for body in microstate.bodies() {
38        if !should_zero_body(body) {
39            continue;
40        }
41
42        let position = body.item.properties.position();
43        let mass = body.item.properties.mass();
44
45        center_of_mass += *position * mass;
46        total_mass += mass;
47    }
48    center_of_mass /= total_mass;
49
50    let mut angular_momentum_center = Cartesian::default();
51    let mut moment_of_inertia_center = Matrix::<3, 3>::zeros();
52    for body in microstate.bodies() {
53        if !should_zero_body(body) {
54            continue;
55        }
56
57        let position = body.item.properties.position();
58        let momentum = body.item.properties.momentum();
59        let mass = body.item.properties.mass();
60
61        let r = *position - center_of_mass;
62        angular_momentum_center += r.wedge(momentum);
63
64        moment_of_inertia_center +=
65            (Matrix::with_diagonal([r.norm_squared(); 3]) - r.outer(&r)) * mass;
66    }
67
68    let center_angular_momentum_matrix = angular_momentum_center.to_row_matrix();
69    let (u, s, vt) = moment_of_inertia_center.svd();
70
71    // If the system do not rotate w. r. t. the principle axis (I_principal=0),
72    // set the omega component to 0 by setting the corresponding s^-1 to 0.
73    let mut s_inv_dense = Matrix::<3, 3>::zeros();
74    if s[0] > 0.0 {
75        s_inv_dense.rows[0][0] = 1.0 / s[0];
76    }
77    if s[1] > 0.0 {
78        s_inv_dense.rows[1][1] = 1.0 / s[1];
79    }
80    if s[2] > 0.0 {
81        s_inv_dense.rows[2][2] = 1.0 / s[2];
82    }
83
84    // omega = L * v * s^-1 * u^t (omega and L are row matrices)
85    let omega = center_angular_momentum_matrix
86        .matmul(&vt.transpose())
87        .matmul(&s_inv_dense)
88        .matmul(&u.transpose());
89    let center_angular_velocity = Cartesian::from(omega.rows[0]);
90
91    for body_index in 0..microstate.bodies().len() {
92        let body = &microstate.bodies()[body_index];
93        if !should_zero_body(body) {
94            continue;
95        }
96
97        let mut body_properties = body.item.properties.clone();
98
99        let position = body_properties.position();
100        let mass = body_properties.mass();
101
102        let r = *position - center_of_mass;
103
104        *body_properties.momentum_mut() -= center_angular_velocity.wedge(&r) * mass;
105
106        microstate
107            .update_body_properties(body_index, body_properties)
108            .expect("Bodies and sites should remain in simulation boundary.");
109    }
110}
111
112/// Zero a 2D microstate's angular momentum.
113#[inline]
114fn zero_angular_momentum_2d<B, S, X, C, F>(
115    microstate: &mut Microstate<B, S, X, C>,
116    should_zero_body: F,
117) where
118    B: Position<Position = Cartesian<2>>
119        + Mass
120        + Momentum<Momentum = Cartesian<2>>
121        + Transform<S>
122        + Clone,
123    S: Position<Position = Cartesian<2>> + Default,
124    X: PointUpdate<Cartesian<2>, SiteKey>,
125    C: Wrap<B> + Wrap<S> + GenerateGhosts<S>,
126    F: Fn(&Tagged<Body<B, S>>) -> bool,
127{
128    let mut center_of_mass = Cartesian::default();
129    let mut total_mass = 0.0;
130
131    for body in microstate.bodies() {
132        if !should_zero_body(body) {
133            continue;
134        }
135
136        let position = body.item.properties.position();
137        let mass = body.item.properties.mass();
138
139        center_of_mass += *position * mass;
140        total_mass += mass;
141    }
142    center_of_mass /= total_mass;
143
144    let mut angular_momentum_center = 0.0;
145    let mut moment_of_inertia_center = 0.0;
146
147    for body in microstate.bodies() {
148        if !should_zero_body(body) {
149            continue;
150        }
151
152        let position = body.item.properties.position();
153        let momentum = body.item.properties.momentum();
154        let mass = body.item.properties.mass();
155
156        let r = *position - center_of_mass;
157
158        angular_momentum_center += r.wedge(momentum);
159
160        moment_of_inertia_center += r.norm_squared() * mass;
161    }
162
163    if moment_of_inertia_center > 0.0 {
164        let angular_velocity_center = angular_momentum_center / moment_of_inertia_center;
165
166        for body_index in 0..microstate.bodies().len() {
167            let body = &microstate.bodies()[body_index];
168            if !should_zero_body(body) {
169                continue;
170            }
171
172            let mut body_properties = body.item.properties.clone();
173
174            let position = body_properties.position();
175            let mass = body_properties.mass();
176
177            let r = *position - center_of_mass;
178
179            *body_properties.momentum_mut() -= r.perpendicular() * angular_velocity_center * mass;
180
181            microstate
182                .update_body_properties(body_index, body_properties)
183                .expect("Bodies and sites should remain in simulation boundary.");
184        }
185    }
186}
187
188impl<R, S, X, C> ZeroCenterAngularMomentum<DynamicOrientedPoint<Cartesian<3>, R>, S>
189    for Microstate<DynamicOrientedPoint<Cartesian<3>, R>, S, X, C>
190where
191    R: RotationalMotionTypes,
192    DynamicOrientedPoint<Cartesian<3>, R>: Transform<S> + Clone,
193    S: Position<Position = Cartesian<3>> + Default,
194    X: PointUpdate<Cartesian<3>, SiteKey>,
195    C: Wrap<DynamicOrientedPoint<Cartesian<3>, R>> + Wrap<S> + GenerateGhosts<S>,
196{
197    #[inline]
198    fn zero_center_angular_momentum_with_filter<
199        F: Fn(&Tagged<Body<DynamicOrientedPoint<Cartesian<3>, R>, S>>) -> bool,
200    >(
201        &mut self,
202        should_zero_body: F,
203    ) {
204        zero_angular_momentum_3d(self, should_zero_body);
205    }
206}
207
208impl<S, X, C> ZeroCenterAngularMomentum<DynamicPoint<Cartesian<3>>, S>
209    for Microstate<DynamicPoint<Cartesian<3>>, S, X, C>
210where
211    DynamicPoint<Cartesian<3>>: Transform<S>,
212    S: Position<Position = Cartesian<3>> + Default,
213    X: PointUpdate<Cartesian<3>, SiteKey>,
214    C: Wrap<DynamicPoint<Cartesian<3>>> + Wrap<S> + GenerateGhosts<S>,
215{
216    #[inline]
217    fn zero_center_angular_momentum_with_filter<
218        F: Fn(&Tagged<Body<DynamicPoint<Cartesian<3>>, S>>) -> bool,
219    >(
220        &mut self,
221        should_zero_body: F,
222    ) {
223        zero_angular_momentum_3d(self, should_zero_body);
224    }
225}
226
227impl<R, S, X, C> ZeroCenterAngularMomentum<DynamicOrientedPoint<Cartesian<2>, R>, S>
228    for Microstate<DynamicOrientedPoint<Cartesian<2>, R>, S, X, C>
229where
230    R: RotationalMotionTypes,
231    DynamicOrientedPoint<Cartesian<2>, R>: Transform<S> + Clone,
232    S: Position<Position = Cartesian<2>> + Default,
233    X: PointUpdate<Cartesian<2>, SiteKey>,
234    C: Wrap<DynamicOrientedPoint<Cartesian<2>, R>> + Wrap<S> + GenerateGhosts<S>,
235{
236    #[inline]
237    fn zero_center_angular_momentum_with_filter<
238        F: Fn(&Tagged<Body<DynamicOrientedPoint<Cartesian<2>, R>, S>>) -> bool,
239    >(
240        &mut self,
241        should_zero_body: F,
242    ) {
243        zero_angular_momentum_2d(self, should_zero_body);
244    }
245}
246
247impl<S, X, C> ZeroCenterAngularMomentum<DynamicPoint<Cartesian<2>>, S>
248    for Microstate<DynamicPoint<Cartesian<2>>, S, X, C>
249where
250    DynamicPoint<Cartesian<2>>: Transform<S>,
251    S: Position<Position = Cartesian<2>> + Default,
252    X: PointUpdate<Cartesian<2>, SiteKey>,
253    C: Wrap<DynamicPoint<Cartesian<2>>> + Wrap<S> + GenerateGhosts<S>,
254{
255    #[inline]
256    fn zero_center_angular_momentum_with_filter<
257        F: Fn(&Tagged<Body<DynamicPoint<Cartesian<2>>, S>>) -> bool,
258    >(
259        &mut self,
260        should_zero_body: F,
261    ) {
262        zero_angular_momentum_2d(self, should_zero_body);
263    }
264}
hoomd_md/modify/zero_center_angular_momentum.rs

hoomd_md/modify/
zero_center_angular_momentum.rs
