Print molecular dipole moment

Author: mlund

src/main.rs

@@ -182,6 +182,14 @@ fn do_scan(cmd: &Commands) -> Result<()> {
         ref_a.net_charge(),
         ref_b.net_charge(),
     );
+
+    const ELECTRON_ANGSTROM_TO_DEBYE: f64 = 4.80320425;
+    info! {
+        "Molecular dipole moments: [{:.2} D, {:.2} D]",
+        ref_a.dipole_moment().norm() * ELECTRON_ANGSTROM_TO_DEBYE,
+        ref_b.dipole_moment().norm() * ELECTRON_ANGSTROM_TO_DEBYE,
+    };
+
     info!(
         "Molecular masses (g/mol): [{:.2}e, {:.2}e]",
         ref_a.total_mass(),

src/structure.rs

@@ -152,6 +152,21 @@ impl Structure {
         self.charges.iter().sum()
     }
 
+    /// Molecular dipole moment of the structure with respect to the mass center
+    ///
+    /// The dipole moment is calculated as
+    /// 𝒎 = ∑ 𝒓ᵢ𝑞ᵢ where 𝒓ᵢ = 𝒑ᵢ - 𝑪 and 𝑪 is the mass center.
+    /// Note that if the molecule is not neutral, the dipole moment depends on the
+    /// choice of origin, here the mass center.
+    pub fn dipole_moment(&self) -> Vector3 {
+        let center = self.mass_center();
+        self.pos
+            .iter()
+            .zip(&self.charges)
+            .map(|(pos, charge)| (pos - center).scale(*charge))
+            .fold(Vector3::zeros(), |sum, i| sum + i)
+    }
+
     /// Total mass of the structure
     pub fn total_mass(&self) -> f64 {
         self.masses.iter().sum()