tentatively think this is working

cctk/parse_gaussian.py

@@ -967,7 +967,7 @@ def parse_modes(freq_block):
 
     modes = list()
     for f, m, k, d in zip(freqs, masses, force_ks, displacements):
-        k *= 143.836 # mdyne Å**-1 to kcal/mol Å**-2
-        modes.append(cctk.VibrationalMode(f, m, k, d))
+        k *= 143.9326 # mdyne Å**-1 to kcal/mol Å**-2
+        modes.append(cctk.VibrationalMode(frequency=f, reduced_mass=m, force_constant=k, displacements=d))
 
     return modes

cctk/quasiclassical.py

@@ -51,7 +51,7 @@ def apply_vibration(molecule, mode, min_freq=50, temperature=298):
 
     # here we could compute atom velocities if we wanted to! initializer lines 440-480
 
-    print(f"Mode {mode.frequency:.2f}\t QC Level {level}\t Shift {rel_shift:.2%} of a potential {max_shift:.2f} Å\tPE = {potential_energy:.2f} kcal/mol")
+    print(f"Mode {mode.frequency:.2f} ({mode.energy():.2f} kcal/mol)\t QC Level {level}\t Shift {rel_shift:.2%} of a potential {max_shift:.2f} Å\tPE = {potential_energy:.2f} kcal/mol\tk = {mode.force_constant:.2f} kcal/mol Å^-2")
 
     return potential_energy, kinetic_energy
 

cctk/vibrational_mode.py

@@ -141,9 +141,8 @@ def quantum_distribution_value(self, x, level=0):
         H = get_hermite_polynomial(n)
 
         # following https://github.com/ekwan/Jprogdyn/blob/master/src/main/java/edu/harvard/chemistry/ekwan/Jprogdyn/HarmonicOscillatorDistribution.java, line 109
-
-        # 4 * pi * 3 * 10**10 / (1000 * 10**23 * 6.022 * 10**23 * 6.626 * 10^-34) = 9.448 * 10 ** -6, take it or leave it
-        omega_term = 9.448 * 10 ** -6 * self.reduced_mass * freq
+        # 4 * pi * 3 * 10**8 / (1000 * 10**20 * 6.022 * 10**23 * 6.626 * 10^-34) = 0.000094411, take it or leave it
+        omega_term = 9.4411e-5 * self.reduced_mass * freq
         val = math.sqrt(omega_term) * math.exp(-1 * omega_term * math.pi * x ** 2 ) * (H(math.sqrt(omega_term * math.pi) * x) ** 2) / (2 ** n * math.factorial(n))
         return val