diff --git a/src/core/src/consts.rs b/src/core/src/consts.rs
index a1de8a7a0..b4fa8cfdb 100644
--- a/src/core/src/consts.rs
+++ b/src/core/src/consts.rs
@@ -8,6 +8,6 @@ pub const K_BOLTZMANN: f64 = 8.31446284161522e-7;
 /// Bohr radius
 pub const BOHR_RADIUS: f64 = 0.52917720859;
 /// Avogadro number
-pub const NA: f64 = 6.02214179e23;
-/// 4 * pi * epsilon_0
-pub const ELCC: f64 = 7.197589831304046;
+pub const AVOGADRO_NUMBER: f64 = 6.02214179e23;
+/// 4 * π * epsilon_0
+pub const FOUR_PI_EPSILON_0: f64 = 7.197589831304046;
diff --git a/src/core/src/energy/global/ewald.rs b/src/core/src/energy/global/ewald.rs
index f0ffb8551..f83736d48 100644
--- a/src/core/src/energy/global/ewald.rs
+++ b/src/core/src/energy/global/ewald.rs
@@ -12,7 +12,7 @@ use ndarray::Zip;
 use math::*;
 use sys::{Configuration, UnitCell, CellShape};
 use types::{Matrix3, Vector3D, Array3, Complex, Zero, One};
-use consts::ELCC;
+use consts::FOUR_PI_EPSILON_0;
 use energy::{PairRestriction, RestrictionInfo};
 use parallel::ThreadLocalStore;
 use parallel::prelude::*;
@@ -360,7 +360,7 @@ impl Ewald {
         for charge in configuration.particles().charge {
             q2 += charge * charge;
         }
-        q2 /= ELCC;
+        q2 /= FOUR_PI_EPSILON_0;
 
         let natoms = configuration.size() as f64;
         let lengths = configuration.cell.lengths();
@@ -429,11 +429,11 @@ impl Ewald {
         }
 
         if !info.excluded {
-            qiqj / ELCC * erfc(self.alpha * r) / r
+            qiqj / FOUR_PI_EPSILON_0 * erfc(self.alpha * r) / r
         } else {
             // use a correction for excluded interaction, removing the energy
             // from kspace
-            - qiqj / ELCC * erf(self.alpha * r) / r
+            - qiqj / FOUR_PI_EPSILON_0 * erf(self.alpha * r) / r
         }
     }
 
@@ -450,14 +450,14 @@ impl Ewald {
         }
 
         if !info.excluded {
-            qiqj / (ELCC * r * r) * (
+            qiqj / (FOUR_PI_EPSILON_0 * r * r) * (
                 self.alpha * FRAC_2_SQRT_PI * exp(-self.alpha * self.alpha * r * r)
                 + erfc(self.alpha * r) / r
             )
         } else {
             // use a correction for excluded interaction, removing the force
             // from kspace
-            qiqj / (ELCC * r * r) * (
+            qiqj / (FOUR_PI_EPSILON_0 * r * r) * (
                 self.alpha * FRAC_2_SQRT_PI * exp(-self.alpha * self.alpha * r * r)
                 - erf(self.alpha * r) / r
             )
@@ -655,7 +655,7 @@ impl Ewald {
                               .iter()
                               .map(|q| q * q)
                               .sum::<f64>();
-        return -self.alpha / sqrt(PI) * q2 / ELCC;
+        return -self.alpha / sqrt(PI) * q2 / FOUR_PI_EPSILON_0;
     }
 }
 
@@ -713,7 +713,7 @@ impl Ewald {
             .and(&self.rho)
             .par_map(|(factor, rho)| factor * rho.norm2())
             .sum();
-        return energy / ELCC;
+        return energy / FOUR_PI_EPSILON_0;
     }
 
     /// k-space contribution to the forces
@@ -744,7 +744,7 @@ impl Ewald {
 
         let charges = configuration.particles().charge;
         for (force, &charge, field) in izip!(&mut *forces, charges, &self.efield) {
-            *force += charge * field / ELCC;
+            *force += charge * field / FOUR_PI_EPSILON_0;
         }
     }
 
@@ -757,7 +757,7 @@ impl Ewald {
             .par_map(|(factor, rho)| rho.norm2() * factor)
             .sum::<Matrix3>();
 
-        return virial / ELCC;
+        return virial / FOUR_PI_EPSILON_0;
     }
 
     /// k-space contribution to the molecular virial
@@ -845,7 +845,7 @@ impl Ewald {
         for (factor, &rho) in izip!(&self.factors.energy, &self.rho) {
             old_energy += factor * rho.norm2();
         }
-        old_energy /= ELCC;
+        old_energy /= FOUR_PI_EPSILON_0;
 
         let delta_rho = self.delta_rho_move_rigid_molecules(
             configuration, molecule_id, new_positions
@@ -855,7 +855,7 @@ impl Ewald {
         for (factor, &rho, &delta) in izip!(&self.factors.energy, &self.rho, &delta_rho) {
             new_energy += factor * (rho + delta).norm2();
         }
-        new_energy /= ELCC;
+        new_energy /= FOUR_PI_EPSILON_0;
 
         self.updater = Some(Box::new(move |ewald: &mut Ewald| {
             for (rho, &delta) in izip!(&mut ewald.rho, &delta_rho) {
diff --git a/src/core/src/energy/global/wolf.rs b/src/core/src/energy/global/wolf.rs
index 133e9e748..c19b09665 100644
--- a/src/core/src/energy/global/wolf.rs
+++ b/src/core/src/energy/global/wolf.rs
@@ -2,7 +2,7 @@
 // Copyright (C) Lumol's contributors — BSD license
 use std::f64::consts::PI;
 
-use consts::ELCC;
+use consts::FOUR_PI_EPSILON_0;
 use energy::{PairRestriction, RestrictionInfo};
 use math::*;
 use parallel::ThreadLocalStore;
@@ -87,13 +87,13 @@ impl Wolf {
         if rij > self.cutoff || info.excluded {
             return 0.0;
         }
-        info.scaling * qiqj * (erfc(self.alpha * rij) / rij - self.energy_cst) / ELCC
+        info.scaling * qiqj * (erfc(self.alpha * rij) / rij - self.energy_cst) / FOUR_PI_EPSILON_0
     }
 
     /// Compute the energy for self interaction of a particle with charge `qi`
     #[inline]
     fn energy_self(&self, qi: f64) -> f64 {
-        qi * qi * (self.energy_cst / 2.0 + self.alpha / sqrt(PI)) / ELCC
+        qi * qi * (self.energy_cst / 2.0 + self.alpha / sqrt(PI)) / FOUR_PI_EPSILON_0
     }
 
     /// Compute the force for self the pair of particles with charge `qi` and
@@ -107,7 +107,7 @@ impl Wolf {
         }
         let factor = erfc(self.alpha * d) / (d * d)
             + 2.0 * self.alpha / sqrt(PI) * exp(-self.alpha * self.alpha * d * d) / d;
-        return info.scaling * qiqj * (factor - self.force_cst) * rij.normalized() / ELCC;
+        return info.scaling * qiqj * (factor - self.force_cst) * rij.normalized() / FOUR_PI_EPSILON_0;
     }
 }
 
@@ -118,8 +118,8 @@ impl GlobalCache for Wolf {
         molecule_id: usize,
         new_positions: &[Vector3D],
     ) -> f64 {
-        let mut old_energynergy = 0.0;
-        let mut new_energynergy = 0.0;
+        let mut old_energy = 0.0;
+        let mut new_energy = 0.0;
 
         let charges = configuration.particles().charge;
         let positions = configuration.particles().position;
@@ -146,13 +146,13 @@ impl GlobalCache for Wolf {
                     let path = configuration.bond_path(part_i, part_j);
                     let info = self.restriction.information(path);
 
-                    old_energynergy += self.energy_pair(info, qi * qj, old_r);
-                    new_energynergy += self.energy_pair(info, qi * qj, new_r);
+                    old_energy += self.energy_pair(info, qi * qj, old_r);
+                    new_energy += self.energy_pair(info, qi * qj, new_r);
                 }
             }
         }
 
-        return new_energynergy - old_energynergy;
+        return new_energy - old_energy;
     }
 
     fn update(&self) {
diff --git a/src/core/src/units.rs b/src/core/src/units.rs
index bc836dafa..7f6756746 100644
--- a/src/core/src/units.rs
+++ b/src/core/src/units.rs
@@ -21,7 +21,7 @@ use std::num;
 
 // Using a separated module because lazy_static does not support pub(crate)
 mod detail {
-    use consts::{BOHR_RADIUS, NA};
+    use consts::{BOHR_RADIUS, AVOGADRO_NUMBER};
     use std::collections::BTreeMap;
     use std::f64::consts::PI;
 
@@ -55,7 +55,7 @@ mod detail {
             // Temperature units.
             assert!(map.insert("K", 1.0).is_none());
             // Quantity of matter units.
-            assert!(map.insert("mol", NA).is_none());
+            assert!(map.insert("mol", AVOGADRO_NUMBER).is_none());
 
             // Angle units.
             assert!(map.insert("rad", 1.0).is_none());