freeH partition function calculation

[Dempsey]

Dear Users,

I have been trying to reproduce ln(qtrans), ln(qrot), and ln(qvib) as given by the freeH script. I can reproduce the vibrational partition function well, but I cannot reproduce the translational/rotational values using the equations printed in the output. I have also tried using equations given in "Molecular Thermodynamics" by McQuarrie and Simon. I'm not sure what I am missing, but I cannot get the equations for qtrans and qrot to work.  In fact, for qrot freeH says to use A*B*C whereas McQuarrie and Simon has a very similar equation but it uses the rotational temperatures instead, that equation also doesn't reproduce ln(qrot).

FreeH says it is using the following:
qtrans=(m*k*T/(2*pi))^1.5*v h/(2*pi)=1a.u.(omitted)  m=mol. mass in a.u. v=vol. per mol. of ideal gas
qrot=((2*pi*kT)^3*A*B*C)^0.5/(sigma*pi) A,B,C=moments of inertia in a.u. 

For qrot, I calculate the same values of A, B, and C from the inertia tensor but fail to reproduce qrot from those. Which is strange since it is simple maths from there, maybe I am missing some unit conversion?

Could you share more details about freeH, there is nothing in the manual about it.

I have an example below for a water molecule showing that I am using the equations in freeH but getting a different result. For clarity, my mass is in a.u. and I took volume to be 22.4 Lmol-1 converted to borh^3  (1.5125e+29 bohr^3).

Code: [Select]

   
freeH output:

          ----------------------------------------
          eigenvalues & -vectors of inertia tensor
          ----------------------------------------

         2.17654    6.32312    4.14658


    T        p       ln(qtrans) ln(qrot) ln(qvib)
  (K)      (MPa)                                 

 298.15   0.1000000      14.94     3.77     0.00


Code: [Select]

My equations and results:
q_trans = ((mass * k * T) / (2 * pi))**1.5 * volume
q_rot = np.sqrt((2*pi*k*T)**3*A*B*C)/(sigma*pi)

BV = 0.9914       
v = np.array(v[6:])*BV
top = c*h*v
bot = k*T
vexp = 1 - np.exp(-top/bot)

q_vib = np.prod(1/(1-np.exp(-(v*h*c)/(k*T))))


Rotational Moments
A =  2.17654
B =  4.14658
C =  6.32312
Partition Functions
ln(qtrans) =  -1.64
ln(qrot) =  -67.47
ln(qvib) =  0.0

Best,
Dempsey

[Dempsey]

Dear Users,

Unsurprisingly it was a mistake on my part with units, I am able to reproduce all partition functions now. There is a minor difference in qrot and qtrans possibly due to my choice of atomic masses and constants being to lower precision than those used in Turbomole but I cannot say for sure.

I would like to know how the zero-point energy is calculated. Textbooks typically quote the equation E₀ = ½ ℏ sqrt(^k/_µ) or E₀ = E_elec +  Σ(modes)½ ℏ sqrt(^k/_µ) or some equivalent. How is it done in freeH? I'm not sure if there is information about k and µ in my output or whether those equations only hold for simple diatomics, should I be able to calculate the zero-point correction from $vibspectrum, or do I need more information?

Could somebody point me towards a good textbook on the topic?

Best,
Dempsey

[uwe]

Hi,

sqrt(k/mu) is the wave number (or angular frequency as Wikipedia names it, see: https://en.wikipedia.org/wiki/Zero-point_energy) which is printed by freeh, but easier to get from the keyword $vibrational spectrum in the control file. Note that the unit is 1/cm as shown in the header of this keyword.

The conversion to kJ/mol can be done in different ways (all hopefully resulting in the same zero point energy value), but the perhaps easiest or clearest way is this:

• Conversion 1: from 1/cm to atomic units (Hartree) for each individual wave number. Conversion factor is 4.55634E-6 or 1/219474.63.
  
  
• Sum up ZPE in atomic units: sum up all wave numbers in atomic units and multiply the sum with 0.5. As we use atomic units here, there is no additional factor for h_bar as it is 1 by definition, see https://en.wikipedia.org/wiki/Hartree_atomic_units
  The number you get is the zero point vibrational energy in atomic units, which is Hartree per molecule.
  
  
• Conversion 2: To convert Hartree/molecule to kJ/mol you need two conversion factors: Hartree to kJ and 1 particle to 1 mol (Avogardo). The combination of both is needed so often, that many converters do that automatically and one has to be very careful if conversions factors are from Hartree to kJ or from Hartree to kJ/mol.
  A nice web site to convert energies is from Prof. Martin (Weizmann Institute) https://www.weizmann.ac.il/oc/martin/tools/hartree.html where you can add a 1 to any of the input fields to get the conversion factors.
  
  If you enter the ZPE value from step 2 on the Weizmann Institute site, you will get kJ/mol directly.

 

[Dempsey]

Dear Uwe,

Thank you for the help. I have managed to calculate everything now, it was a good learning experience. I do have one further question, in order to calculate the internal thermal energy ("energy" in freeH output) I found that the zero point energy is not included. According to the output freeH calculates "energy" like this:

Code: [Select]

energy=ZPE+3RT +sum(i) e(i)*(1+exp(-e(i)/kT))/2*(1-exp(-e(i)/kT)=3/2RT for atoms etc.  
But I could not get this to work, I recognize that this is the ZPE + translational (3/2RT) + rotational (3/2RT) + vibrational contributions to the internal energy. The vibrational term is not in the units per mol and including avogadros constant did not fix the result. Instead, I calculated "energy" according to how it is written in "Thermochemistry in Gaussian" https://gaussian.com/wp-content/uploads/dl/thermo.pdf where the vibrational contributions are as follows:

Code: [Select]

Ev = R * sum(k) theta(k) (0.5 + 1 / e(theta(k)/T) - 1) where theta(k) is the characteristic vibrational temperature theta(k) = hv(k)/k_b
The equations listed there reproduce the values calculated by freeH but "energy" does not include the zero point energy. I'm not sure if this is an error in the freeH output or a missunderstanding. Is the zero point included in "energy" and hence enthalpy? Or is it not included?