Thermodynamic quantities from vibrational frequencies using "freeh"

[hannesb]

Hello Everyone,

I have a question regarding the determination of thermodynamic quantities using the internal program freeh. I was trying to use frequency dependent scaling factors and wanted to recalculate thermodynamic properties by using the frequencies from my aoforce calculations.
I used the equations given in the output of the freeh calculation.
However, I was only able to recalculate:

qvib          = product(i) 1/(1-exp(-e(i)/kT)     and
chem.pot. = ZPE-RT*ln(qtrans*qrot*qvib)   and
the ZPE itself.

Unfortunately I could not get the same results for the internal energy E:

energy=ZPE + 3RT + sum(i) e(i)*(1+exp(-e(i)/kT))/2*(1-exp(-e(i)/kT)

Frequencies:

mode     symmetry     wave number     intensity     selection rules
                          cm**(-1)         km/mol   
     1                        0.00       0.0000000        -       -
     2                        0.00       0.0000000        -       -
     3                        0.00       0.0000000        -       -
     4                        0.00       0.0000000        -       -
     5                        0.00       0.0000000        -       -
     6                        0.00       0.0000000        -       -
     7        a              76.29       0.4158000       YES     YES
     8        a             430.57       4.6595700       YES     YES
     9        a             545.01      29.8828900       YES     YES
    10        a             590.45      35.8843700       YES     YES
    11        a             661.64      **********       YES     YES
    12        a             856.86       4.4155700       YES     YES
    13        a            1000.43      81.9322100       YES     YES
    14        a            1070.80       7.6460300       YES     YES
    15        a            1203.96      **********       YES     YES
    16        a            1337.60      32.6408900       YES     YES
    17        a            1411.57      50.2645300       YES     YES
    18        a            1474.99      15.7249400       YES     YES
    19        a            1478.70      10.8687900       YES     YES
    20        a            1820.67      **********       YES     YES
    21        a            3051.63       1.5803000       YES     YES
    22        a            3112.18       4.2956400       YES     YES
    23        a            3161.39       4.9184500       YES     YES
    24        a            3736.08      55.9035900       YES

My calculations:

ei = waveNumber(mode_i)*c*h;                                                                     % [1/cm --> J]
E  = ZPE + 3*R*T + sum(ei.*(1+exp(-ei./(k*T)))/2.*(1-exp(-ei./(k*T))));        % calculation of E

When I tried to compare the equation above with those from Turbomole theory I somehow could not see the connection of both equations.
(http://hincklab.uthscsa.edu/html/soft_packs/msi_docs/insight980/turbo/2_Theory.html       -    equation 44)

Can anyone help? I tried several ways as well via the helmholtz energy E = -kT (d(A/T)/dT). Maybe these are all due to wrong calculations but I checked the code several times.
But maybe one of you can find a general problem in my approach.

Thanks for your help!

Hannes

[diedenho]

Hi Hannes,

I think there might be a problem with a bracket

old
energy=ZPE + 3RT + sum(i) e(i)*(1+exp(-e(i)/kT))/2*(1-exp(-e(i)/kT)
new
energy=ZPE + 3RT + sum(i) e(i)*(1+exp(-e(i)/kT))/ [ 2*(1-exp(-e(i)/kT) ]

the upper formula is for non-linear polyatomic molecules
3RT = 3/2 RT (trans) + 3/2 RT (rot)
for linear molecules this would be
3/2 RT (trans) + RT (rot)

dos this solve your problem?
Michael

[hannesb]

Hey Michael,

first of all thank you for your reply!
I changed the equation but still there are big differences between the results from freeh and the discussed equation:

      freeh   eq. above           abs Diff
T [K]      E [kJ/mol]   E [kJ/mol]   [-]   
0      160,23        160,23   0,00       <--- starting from the same ZPE
50      161,59        160,60   0,99
100      163,20        161,85   1,35
200      167,07        164,34   2,73
300      172,35        166,84   5,51
400      179,16        169,33   9,83
500       187,38        171,82   15,56
600       196,79        174,32   22,47
700       207,21        176,81   30,40
800       218,45        179,31   39,14
900       230,41        181,80   48,61
1000   242,97        184,30   58,67
1100   256,06        186,79   69,27
1200   269,61        189,28   80,33
1300   283,55        191,78   91,77
1400   297,84        194,27   103,57
1500   312,44        196,77   115,67

Somehow I have the feeling that they calculate the internal energy in a different way. Maybe via the partition function which i could not yet reproduce..
Thank you anyway! If you have another idea, please let me know

Cheers,

Hannes

[hannesb]

Hello again,

I solved it in a different way and calculate the free energy via the partition function (chem.pot. = ZPE-RT*ln(qtrans*qrot*qvib)) and then using the gibbs-helmholtz relation:
-h/(T^2) = d(G/T)/dT
I don't think that TURBOMOLE uses -- energy=ZPE + 3RT + sum(i) e(i)*(1+exp(-e(i)/kT))/ [ 2*(1-exp(-e(i)/kT) ] --  to calculate the internal energy as hereby the rotational and translational contributions are simplified to 3RT..

Best regards,

Hannes

[Arnim]

Hi Hannes,

freeh actually uses the formular as printed in the output with the backets clearified Michael.

Cheers,

Arnim

[uwe]

Hi,

which results did you reproduce using the partition function: The freeh output of the free energy?

Uwe

[hannesb]

Hi Uwe,

First of all I could reproduce the partition function. With the calculated partition function I got the same results for the chemical potential.
Then I used the Gibbs-Helmholtz relation -h/(T^2) = d(G/T)/dT to get the enthalpy (I differentiated numerically G/T with respect to T).
The results were were identical with the results from Turbomole. Then I calculated cp (which was my main goal) as the derivative of h with respect to T.
However, I run into numerical issues but got good results using a very fine grid.

When Armin said that Turbomole uses the formula printed in the output I checked again and couldn't reproduce the results from Turbomole. However, maybe I use wrong dimensions in this case.

energy=ZPE + 3RT + sum(i) e(i)*(1+exp(-e(i)/kT))/ [ 2*(1-exp(-e(i)/kT) ]

For e(i) I used [J] coming from the wavenumbers:   e(i) = vibSpec*c*h  (vibSpec is a vector including all frequencies in wavenumbers cm-1, c = speed of light in cm/s, h is planck's constant)

Cheers,

Hannes

[uwe]

Hi,

I think the ZPE is already included in the (1+a)/(1-a) part, where a is exp(-E/kt). About the conversion factors: I am pretty sure that freeh in Turbomole does it correctly since it was used and modified by several generations of PhD students and Post-Docs... especially if you were able to reproduce the numbers coming from Gibbs-Helmholtz.

Regards,

Uwe

[hannesb]

Hi Uwe,

please don't get me wrong, I don't doubt that the results from TURBOMOLE are correct. I was just wondering why the I couldn't get the same results using the aformentioned equation for the internal energy. Maybe you have an idea where I am missing something.

Anyway thank you for your help!

Hannes