TURBOMOLE Users Forum
Forum General => Miscellaneous => Topic started by: hannesb on February 18, 2014, 05:43:21 PM
-
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
-
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
-
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
-
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
-
Hi Hannes,
freeh actually uses the formular as printed in the output with the backets clearified Michael.
Cheers,
Arnim
-
Hi,
which results did you reproduce using the partition function: The freeh output of the free energy?
Uwe
-
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
-
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
-
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