Add degeneracy

Author: mlund

scripts/cppm.sh

@@ -9,7 +9,7 @@ RUST_LOG="Debug" cargo run --release -- \
     -2 cppm-p18.xyz \
     --rmin 40.0 --rmax 150 --dr 1.0 \
     --top topology.yaml \
-    --resolution 0.5 \
+    --resolution 0.8 \
     --cutoff 1000 \
     --molarity 0.005 \
     #--savetable table.dat.gz \

src/icoscan.rs

@@ -171,9 +171,10 @@ pub fn do_icoscan(
     let calc_partition_func = |r: f64, omega: f64| {
         table.get_icospheres(r, omega).unwrap().flat_iter().fold(
             Sample::default(),
-            |sum, (_, _, data_b)| {
+            |sum, (vertex_i, vertex_j, data_b)| {
+                let degeneracy = vertex_i.norm() * vertex_j.norm();
                 let energy = data_b.get().unwrap();
-                sum + Sample::new(*energy, *temperature)
+                sum + Sample::new(*energy, *temperature, degeneracy)
             },
         )
     };

src/sample.rs

@@ -23,30 +23,33 @@ pub struct Sample {
     n: u64,
     /// Thermal energy, RT in kJ/mol
     thermal_energy: f64,
-    /// Boltzmann weighted energy, U * exp(-U/kT)
+    /// Boltzmann weighted energy, U * g x exp(-U/kT)
     mean_energy: f64,
-    /// Boltzmann weighted squared energy, U^2 * exp(-U/kT)
+    /// Boltzmann weighted squared energy, U^2 * g x exp(-U/kT)
     mean_squared_energy: f64,
-    /// Boltzmann factored energy, exp(-U/kT)
+    /// Boltzmann factored energy, g x exp(-U/kT)
     exp_energy: f64,
-    /// Boltzmann factored energy minus one, exp(-U/kT) - 1
+    /// Boltzmann factored energy minus one, g x exp(-U/kT) - 1
     exp_energy_m1: f64,
+    /// Degeneracy of the state, g
+    degeneracy: f64,
 }
 
 impl Sample {
     /// New from energy in kJ/mol and temperature in K
-    pub fn new(energy: f64, temperature: f64) -> Self {
+    pub fn new(energy: f64, temperature: f64, degeneracy: f64) -> Self {
         const KJ_PER_J: f64 = 1e-3;
         let thermal_energy = MOLAR_GAS_CONSTANT * temperature * KJ_PER_J; // kJ/mol
-        let exp_energy = (-energy / thermal_energy).exp();
-        let exp_energy_m1 = (-energy / thermal_energy).exp_m1();
+        let exp_energy = degeneracy * (-energy / thermal_energy).exp();
+        let exp_energy_m1 = degeneracy * (-energy / thermal_energy).exp_m1();
         Self {
             n: 1,
             thermal_energy,
             mean_energy: energy * exp_energy,
             mean_squared_energy: energy.powi(2) * exp_energy,
             exp_energy,
             exp_energy_m1,
+            degeneracy,
         }
     }
     /// Thermal energy (kJ/mol)
@@ -67,11 +70,11 @@ impl Sample {
     }
     /// Free energy (kJ / mol)
     pub fn free_energy(&self) -> f64 {
-        (self.exp_energy / self.n as f64).ln().neg() * self.thermal_energy
+        (self.exp_energy / self.degeneracy as f64).ln().neg() * self.thermal_energy
     }
     /// Mean <exp(-U/kT)-1>
     pub fn mean_exp_energy_m1(&self) -> f64 {
-        self.exp_energy_m1 / self.n as f64
+        self.exp_energy_m1 / self.degeneracy as f64
     }
     /// Number of samples
     pub fn n(&self) -> u64 {
@@ -83,12 +86,13 @@ impl std::fmt::Display for Sample {
     fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
         write!(
             f,
-            "n: {}, ⟨U⟩: {:.3} kJ/mol, C/𝑘𝐵: {:.3}, -𝑘𝑇⟨exp(-𝛽U)⟩: {:.3} kJ/mol, ⟨exp(-𝛽U)-1⟩: {:.3}",
+            "n: {}, ⟨U⟩: {:.3} kJ/mol, C/𝑘𝐵: {:.3}, -𝑘𝑇⟨exp(-𝛽U)⟩: {:.3} kJ/mol, ⟨exp(-𝛽U)-1⟩: {:.3}, g = {:.3}",
             self.n,
             self.mean_energy(),
             self.heat_capacity(),
             self.free_energy(),
-            self.mean_exp_energy_m1()
+            self.mean_exp_energy_m1(),
+            self.degeneracy,
         )
     }
 }
@@ -110,6 +114,7 @@ impl Add for Sample {
             mean_squared_energy: self.mean_squared_energy + other.mean_squared_energy,
             exp_energy: self.exp_energy + other.exp_energy,
             exp_energy_m1: self.exp_energy_m1 + other.exp_energy_m1,
+            degeneracy: self.degeneracy + other.degeneracy,
         }
     }
 }
@@ -122,5 +127,6 @@ impl AddAssign for Sample {
         self.exp_energy += other.exp_energy;
         self.exp_energy_m1 += other.exp_energy_m1;
         self.thermal_energy = f64::max(self.thermal_energy, other.thermal_energy);
+        self.degeneracy += other.degeneracy;
     }
 }