[gmx-developers] details about Heat capacity caculation
David van der Spoel
spoel at xray.bmc.uu.se
Wed Nov 2 09:20:55 CET 2011
On 2011-11-01 22:09, Haiqing Zhao wrote:
>
> Hi David,
>
> Yes, after normalize the Enthalpy, I get the exactly same result with Gromacs. Thank you for patience.
> I do add the flag (-nmol 465) when using g_energy_d. But why the total energy can get normalized output, not Enthalpy?
This is a technicality, and it should be fixed.
Another question is that in my case, with the organic molecule's rate
increasing, the Cv value becomes bigger, than water's? Actually the real
water's only 74.15 when @25ºc.
> Maybe that's why QM is necessary.(?)
(Experimental) Cv for organic molecules in the paper we submitted range
from 60-300.
>
> Best,
> Haiqing
>
> On Nov 1, 2011, at 3:02 PM, David van der Spoel wrote:
>
>> On 2011-11-01 21:00, Haiqing Zhao wrote:
>>>
>>> Hi,
>>>
>>> Energy Average Err.Est. RMSD Tot-Drift
>>> -------------------------------------------------------------------------------
>>> Potential -44.6187 0.007 0.293318 0.0281861 (kJ/mol)
>>> Kinetic En. 7.59226 0.0032 0.19886 0.0201858 (kJ/mol)
>>> Total Energy -37.0264 0.0091 0.359776 0.0483719 (kJ/mol)
>>> pV 0.438734 9.5e-05 0.00136276 -0.000679764 (kJ/mol)
>>> Enthalpy -17216.8 4.2 167.296 22.4922 (kJ/mol)
>>>
>>> Yes, I did not doubt Cp=Var(Enthalpy)/k T^2.
>>> Actually my question is only about Enthalpy. As we see in the example, the total energy= -37.0264KJ/mol, pv = 0.4387 Kj/mol, How can enthalpy= -17216.8 KJ/mol?
>>>
>>> sorry for the time.
>> You did not say how many molecules. Enthalpy is not normalized in the g_energy output.
>>
>>>
>>>
>>> On Nov 1, 2011, at 2:35 PM, David van der Spoel wrote:
>>>
>>>> On 2011-11-01 20:15, Haiqing Zhao wrote:
>>>>>
>>>>> Hi, David,
>>>>>
>>>>> thanks for reply.
>>>>>
>>>>> I checked it. Version 4.5.5 really makes a difference from 4.5.3 which I used.
>>>>> In 4.5.5, Cv works exactly well by nmol*(RMSD(E))^2/RT^2. (Even though I dont know here why need to multiply the # of molecular.)
>>>>>
>>>> This is supposed to work for liquids, hence you need to normalize.
>>>> How many molecules do you have, 465?
>>>> Heat capacity of 80 resembles water. I'm pretty sure it works correctly, since we just submitted a paper about it and I have looked into this in detail. We used:
>>>>
>>>> kB T^2 c_P = Var (Enthalpy)
>>>>
>>>> which comes straight from Allen and Tildesley.
>>>>
>>>>> But still I'm sure the Cp is not from simple Enthalpy Fluctuation, as you see from the example below--167296 J/mol, and there is several orders mistake.
>>>>>
>>>>> Please help me check it.
>>>>> Thanks.
>>>>>
>>>>> e.g.
>>>>>
>>>>> Enthalpy -17216.8 4.2 167.296 22.4922 (kJ/mol)
>>>>>
>>>>> You may want to use the -driftcorr flag in order to correct
>>>>> for spurious drift in the graphs. Note that this is not
>>>>> a substitute for proper equilibration and sampling!
>>>>>
>>>>> Temperature dependent fluctuation properties at T = 300.076.
>>>>>
>>>>> Heat capacities obtained from fluctuations do *not* include
>>>>> quantum corrections. If you want to get a more accurate estimate
>>>>> please use the g_dos program.
>>>>>
>>>>> WARNING: Please verify that your simulations are converged and perform
>>>>> a block-averaging error analysis (not implemented in g_energy yet)
>>>>> Volume = 1.8548e-05 m^3/mol
>>>>> Enthalpy = -37.0255 kJ/mol
>>>>> Coefficient of Thermal Expansion Alpha_P = 0.0129106 (1/K)
>>>>> Isothermal Compressibility Kappa = 3.00555e-10 (J/m^3)
>>>>> Adiabatic bulk modulus = 3.32718e+09 (m^3/J)
>>>>> Heat capacity at constant pressure Cp = 80.3937 J/mol K
>>>>>
>>>>>
>>>>>
>>>>> On Nov 1, 2011, at 2:05 PM, David van der Spoel wrote:
>>>>>
>>>>>> On 2011-11-01 16:54, Haiqing Zhao wrote:
>>>>>>>
>>>>>>> Dear Gmxers,
>>>>>>>
>>>>>>> I have one question about how Gromacs calculate the heat capacity. I dont think it is only Enthalpy fluctuation÷(kT^2) in NPT ensemble. And also I checked it in NVT, it is not simply Total energy(/internal energy) fluctuation÷(kT^2).
>>>>>>>
>>>>>> Which gmx version? It's more or less correct in 4.5.4 but you really need qm correctiions. Check g_dos for that
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>> One data for NPT is: ( -nmol 465 -nconstr 3 )
>>>>>>>
>>>>>>> Energy Average Err.Est. RMSD Tot-Drift
>>>>>>> -------------------------------------------------------------------------------
>>>>>>> Potential -44.6187 0.007 0.293318 0.0281861 (kJ/mol)
>>>>>>> Kinetic En. 7.59226 0.0032 0.19886 0.0201858 (kJ/mol)
>>>>>>> Total Energy -37.0264 0.0091 0.359776 0.0483719 (kJ/mol)
>>>>>>> Temperature 300.076 0.13 7.85971 0.797822 (K)
>>>>>>> Pressure 0.663728 0.16 447.392 -0.388627 (bar)
>>>>>>> Volume 14.3219 0.0093 0.133543 -0.0667073 (nm^3)
>>>>>>> pV 0.438734 9.5e-05 0.00136276 -0.000679764 (kJ/mol)
>>>>>>> Enthalpy -17216.8 4.2 167.296 22.4922 (kJ/mol)
>>>>>>>
>>>>>>> Temperature dependent fluctuation properties at T = 300.076. #constr/mol = 3
>>>>>>> Isothermal Compressibility: 3.00555e-05 /bar
>>>>>>> Adiabatic bulk modulus: 33271.8 bar
>>>>>>> Heat capacity at constant pressure Cp: 67.9219 J/mol K
>>>>>>> Thermal expansion coefficient alphaP: 9.66233e-05 1/K
>>>>>>>
>>>>>>> And data for NVT, ( -nmol 7631 -nconstr 3 )
>>>>>>>
>>>>>>> Energy Average Err.Est. RMSD Tot-Drift
>>>>>>> -------------------------------------------------------------------------------
>>>>>>> Total Energy -33.2264 0.0003 0.0900307 -0.000219237 (kJ/mol)
>>>>>>> Temperature 300.008 0.0063 1.95324 -0.0148441 (K)
>>>>>>>
>>>>>>> Temperature dependent fluctuation properties at T = 300.008. #constr/mol = 3
>>>>>>> Heat capacity at constant volume Cv: 70.1813 J/mol K
>>>>>>>
>>>>>>> (BTW, I found that (1)if you choose Enthalpy as output, g_energy will always give Cp, no matter what ensemble it is. (2)And only you choose tot energy and temp as output, no other terms, g_energy will give you Cv. If not, you cannot get Cv by g_energy.)
>>>>>>>
>>>>>>> I appreciate if anybody could help me to check it.
>>>>>>> Thanks!!
>>>>>>>
>>>>>>> ---------------------------------------------
>>>>>>> Haiqing Zhao
>>>>>>> PH.D. candidate
>>>>>>> in Computational Biophysics
>>>>>>> Michigan Technological University
>>>>>>> http://www.phy.mtu.edu/~haizhao
>>>>>>>
>>>>>>
>>>>>>
>>>>>> --
>>>>>> David van der Spoel, Ph.D., Professor of Biology
>>>>>> Dept. of Cell& Molec. Biol., Uppsala University.
>>>>>> Box 596, 75124 Uppsala, Sweden. Phone: +46184714205.
>>>>>> spoel at xray.bmc.uu.se http://folding.bmc.uu.se
>>>>>> --
>>>>>> gmx-developers mailing list
>>>>>> gmx-developers at gromacs.org
>>>>>> http://lists.gromacs.org/mailman/listinfo/gmx-developers
>>>>>> Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-developers-request at gromacs.org.
>>>>>
>>>>
>>>>
>>>> --
>>>> David van der Spoel, Ph.D., Professor of Biology
>>>> Dept. of Cell& Molec. Biol., Uppsala University.
>>>> Box 596, 75124 Uppsala, Sweden. Phone: +46184714205.
>>>> spoel at xray.bmc.uu.se http://folding.bmc.uu.se
>>>> --
>>>> gmx-developers mailing list
>>>> gmx-developers at gromacs.org
>>>> http://lists.gromacs.org/mailman/listinfo/gmx-developers
>>>> Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-developers-request at gromacs.org.
>>>
>>
>>
>> --
>> David van der Spoel, Ph.D., Professor of Biology
>> Dept. of Cell& Molec. Biol., Uppsala University.
>> Box 596, 75124 Uppsala, Sweden. Phone: +46184714205.
>> spoel at xray.bmc.uu.se http://folding.bmc.uu.se
>> --
>> gmx-developers mailing list
>> gmx-developers at gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-developers
>> Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-developers-request at gromacs.org.
>
--
David van der Spoel, Ph.D., Professor of Biology
Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone: +46184714205.
spoel at xray.bmc.uu.se http://folding.bmc.uu.se
More information about the gromacs.org_gmx-developers
mailing list