dynamics/
com_zero.rs

1//! Relates to zeroing center-of-mass drift adn possibly rotation.
2
3use lin_alg::f32::Vec3;
4
5use crate::MdState;
6
7impl MdState {
8    /// Remove center-of-mass drift for "atoms" only (exclude water).
9    pub fn zero_linear_momentum_atoms(&mut self) {
10        let mut m_sum = 0.0;
11        let mut p_sum = Vec3::new_zero(); // Σ m v
12
13        for a in &self.atoms {
14            let m = a.mass;
15            m_sum += m;
16            p_sum += a.vel * m;
17        }
18        if m_sum <= 0.0 {
19            return;
20        }
21
22        let v_cm = p_sum / m_sum; // COM velocity
23        // Subtract uniformly so Σ m v' = 0
24        for a in &mut self.atoms {
25            a.vel -= v_cm;
26        }
27    }
28
29    // todo: Assess if you want this for multi-molecule systems.
30    /// Remove rigid-body rotation for "atoms" only (exclude water).
31    /// Computes ω from I ω = L about the atoms' COM, then sets v' = v - ω × (r - r_cm).
32    pub fn zero_angular_momentum_atoms(&mut self) {
33        // COM position for atoms (wrapped is fine; all r use the same frame)
34        let mut m_sum = 0.0;
35        let mut m_r_sum = Vec3::new_zero();
36        for a in &self.atoms {
37            let m = a.mass;
38            m_sum += m;
39            m_r_sum += a.posit * m;
40        }
41        if m_sum <= 0.0 {
42            return;
43        }
44        let r_cm = m_r_sum / m_sum;
45
46        // Build inertia tensor I and angular momentum L about r_cm
47        let mut i_xx = 0.0f32;
48        let mut i_xy = 0.0f32;
49        let mut i_xz = 0.0f32;
50        let mut i_yy = 0.0f32;
51        let mut i_yz = 0.0f32;
52        let mut i_zz = 0.0f32;
53        let mut L = Vec3::new_zero();
54
55        for a in &self.atoms {
56            let m = a.mass;
57            // r relative to COM (use minimum-image if you prefer: self.cell.min_image(a.posit - r_cm))
58            let r = a.posit - r_cm;
59            let vxr = r.cross(a.vel);
60            L += vxr * m;
61
62            let rx = r.x;
63            let ry = r.y;
64            let rz = r.z;
65            let r2 = rx * rx + ry * ry + rz * rz;
66            i_xx += m * (r2 - rx * rx);
67            i_yy += m * (r2 - ry * ry);
68            i_zz += m * (r2 - rz * rz);
69            i_xy -= m * (rx * ry);
70            i_xz -= m * (rx * rz);
71            i_yz -= m * (ry * rz);
72        }
73
74        // Inertia tensor (symmetric)
75        let I = [[i_xx, i_xy, i_xz], [i_xy, i_yy, i_yz], [i_xz, i_yz, i_zz]];
76
77        // Solve I * ω = L  (3x3). Add tiny Tikhonov for degeneracy.
78        let eps = 1.0e-6f32;
79        let Ireg = [
80            [I[0][0] + eps, I[0][1], I[0][2]],
81            [I[1][0], I[1][1] + eps, I[1][2]],
82            [I[2][0], I[2][1], I[2][2] + eps],
83        ];
84
85        // Inverse of 3x3 (cofactor / det)
86        let det = Ireg[0][0] * (Ireg[1][1] * Ireg[2][2] - Ireg[1][2] * Ireg[2][1])
87            - Ireg[0][1] * (Ireg[1][0] * Ireg[2][2] - Ireg[1][2] * Ireg[2][0])
88            + Ireg[0][2] * (Ireg[1][0] * Ireg[2][1] - Ireg[1][1] * Ireg[2][0]);
89
90        if det.abs() < 1e-12 {
91            return;
92        }
93
94        let inv_det = 1.0 / det;
95        let inv = [
96            [
97                (Ireg[1][1] * Ireg[2][2] - Ireg[1][2] * Ireg[2][1]) * inv_det,
98                (Ireg[0][2] * Ireg[2][1] - Ireg[0][1] * Ireg[2][2]) * inv_det,
99                (Ireg[0][1] * Ireg[1][2] - Ireg[0][2] * Ireg[1][1]) * inv_det,
100            ],
101            [
102                (Ireg[1][2] * Ireg[2][0] - Ireg[1][0] * Ireg[2][2]) * inv_det,
103                (Ireg[0][0] * Ireg[2][2] - Ireg[0][2] * Ireg[2][0]) * inv_det,
104                (Ireg[0][2] * Ireg[1][0] - Ireg[0][0] * Ireg[1][2]) * inv_det,
105            ],
106            [
107                (Ireg[1][0] * Ireg[2][1] - Ireg[1][1] * Ireg[2][0]) * inv_det,
108                (Ireg[0][1] * Ireg[2][0] - Ireg[0][0] * Ireg[2][1]) * inv_det,
109                (Ireg[0][0] * Ireg[1][1] - Ireg[0][1] * Ireg[1][0]) * inv_det,
110            ],
111        ];
112
113        let omega = Vec3 {
114            x: inv[0][0] * L.x + inv[0][1] * L.y + inv[0][2] * L.z,
115            y: inv[1][0] * L.x + inv[1][1] * L.y + inv[1][2] * L.z,
116            z: inv[2][0] * L.x + inv[2][1] * L.y + inv[2][2] * L.z,
117        };
118
119        // If ω is tiny, nothing to do
120        if !omega.x.is_finite() || !omega.y.is_finite() || !omega.z.is_finite() {
121            return;
122        }
123        if omega.magnitude_squared() < 1e-16 {
124            return;
125        }
126
127        // v' = v - ω × (r - r_cm)
128        for a in &mut self.atoms {
129            let r = a.posit - r_cm;
130            a.vel -= omega.cross(r);
131        }
132
133        // Clean up any translation introduced by roundoff
134        self.zero_linear_momentum_atoms();
135    }
136}
dynamics/com_zero.rs

dynamics/
com_zero.rs
