[gmx-users] Different Cv and Cp

Faezeh Pousaneh fpoosaneh at gmail.com
Thu Jun 18 17:24:57 CEST 2015 

Just for those who may face with similar problem as I had. I found where my
mistake was.

The problem: I have been obtaining Cv and Cp differently for Lutidine. The
obtained heat capacities in Groamcs were consistent with

kB T^2 c_P = Var (Enthalpy)
kB T^2 c_V = Var (Energy)

but not with

Cv = (d<U>/dT)_V

Cp = (d<H>/dT)_P.

Answer: The answer was for (d<U>/dT)_V I had simulated (NVT) at two
different temperatures and was obtaining energy derivatives, while the
volumes were different at two temperatures. Now I kept the volume the same
at two temperature and I get the same results as approach 1.

In fact there is a difference in Cv and Cp in my molecule, as it is seen in
the link for Benzene:

https://books.google.se/books?id=eH_1dIZr-zMC&pg=PA171&dq=liquid+benzene+Cp-Cv&hl=en&sa=X&ved=0CB8Q6AEwAGoVChMIzoPJ0diWxgIVRI4sCh0wCwl0#v=onepage&q=liquid%20benzene%20Cp-Cv&f=false

Best regards


On Thu, Jun 11, 2015 at 8:10 PM, David van der Spoel <spoel at xray.bmc.uu.se>
wrote:

> On 09/06/15 17:29, Michael Shirts wrote:
>
>> If the simulation are generating configurations with the Boltzmann
>> probability distribution, the results should the same up to error.
>>
>> Cv and Cp should not be exactly the same, though for liquids at room
>> temperature, they are pretty close (look up the precise numbers for the
>> fluid you are interested in).
>>
> This is not correct. cP and cV can easily differ 20-30%. Not many cV have
> been measured but you can compute the difference from other fluctuation
> formulae, see e.g. Caleman et al. J. Chem. Theory Comput. 2012, 8, 61–74,
> http://dx.doi.org/10.1021/ct200731v.
> Note that to get accurate numbers you need quantum corrections as well. My
> group are working on implementing a method for that in gmx dos.
>
>
>
>> Are you calculating enthalpy as U + PV, where P is the constant APPLIED
>> pressure, not the instantaneous pressure, and PV is in same units?  Since
>> the PV term should be quite low for liquids, the two heat capacities
>> should
>> be relatively close (within noise -- fluctuation based calculations I
>> think
>> are noisier).
>>
>> Else you need to check if the Boltzmann distributions are being correctly
>> generated: See the code and paper describing it linked here:
>> https://github.com/shirtsgroup/checkensemble
>>
>>
>>
>>
>> On Tue, Jun 9, 2015 at 11:23 AM, Faezeh Pousaneh <fpoosaneh at gmail.com>
>> wrote:
>>
>>  Dear Michael,
>>>
>>> Can I ask a question concerning your previous email,
>>> I followed
>>>
>>> Cv = (d<U>/dT)_V
>>>
>>> Cp = (d<H>/dT)_P
>>>
>>> for my lutidine molecule, and I get same values for Cv and Cp. But when I
>>> test with
>>>
>>> kB T^2 c_P = Var (Enthalpy)
>>> kB T^2 c_V = Var (Energy)
>>>
>>> I get 40 J/mol.K difference in Cv and Cp.
>>>
>>> Mean that fluctuation play big role. Which way of checking I can rely?
>>>
>>>
>>> Best regards
>>>
>>>
>>> On Tue, May 26, 2015 at 3:38 PM, Michael Shirts <mrshirts at gmail.com>
>>> wrote:
>>>
>>>  By definition (more fundamental that fluctuation formulas)
>>>>
>>>> Cv = (d<U>/dT)_V
>>>>
>>>> Cp = (d<H>/dT)_P
>>>>
>>>> Run two simulations at different T and estimate the derivatives.
>>>>
>>>> On Tue, May 26, 2015 at 5:12 AM, Faezeh Pousaneh <fpoosaneh at gmail.com>
>>>> wrote:
>>>>
>>>>  Dear Michael,
>>>>>
>>>>> I still would like to know what was your method you mentioned on last
>>>>> paragraph, just for learning:
>>>>>
>>>>> ''Also, to be sure, you should double check by calculating both heat
>>>>> capacities by finite difference formulas as well with two simulations
>>>>>
>>>> at
>>>
>>>> T+dt/2 and T-dt/2 -- if the fluctuation and finite difference resutls
>>>>>
>>>> don't
>>>>
>>>>> agree within propagated error, then something is off.''
>>>>>
>>>>> ?
>>>>> thanks
>>>>>
>>>>>
>>>>> Best regards
>>>>>
>>>>>
>>>>> On Mon, May 25, 2015 at 5:10 PM, Faezeh Pousaneh <fpoosaneh at gmail.com>
>>>>> wrote:
>>>>>
>>>>>  Dear Andre,
>>>>>>
>>>>>> thank you for the link, you are probably right, It seems that my
>>>>>>
>>>>> molecule
>>>>
>>>>> has the difference Cp-Cv in the same range as benzene (since it has
>>>>>>
>>>>> also
>>>>
>>>>> ring structure).
>>>>>>
>>>>>>
>>>>>> Best regards
>>>>>>
>>>>>>
>>>>>> On Mon, May 25, 2015 at 4:44 PM, Faezeh Pousaneh <
>>>>>>
>>>>> fpoosaneh at gmail.com>
>>>
>>>> wrote:
>>>>>>
>>>>>>  Dear Michael,
>>>>>>>
>>>>>>> I use  Parrinello-Rahman for barostat and v-rescale for thermostat.
>>>>>>>
>>>>>>> Sorry, could you explain more the second paragraph please? I did not
>>>>>>>
>>>>>> get
>>>>
>>>>> the method. What I checked so far is checking if gromacs correctly
>>>>>>>
>>>>>> gives
>>>>
>>>>> Cv,Cp= Var(Energy or Enthalpy)/kBT^2 , and I find that it gives.
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> Best regards
>>>>>>>
>>>>>>>
>>>>>>> On Mon, May 25, 2015 at 4:11 PM, Michael Shirts <mrshirts at gmail.com
>>>>>>>
>>>>>>
>>>>  wrote:
>>>>>>>
>>>>>>>  Are you running with the Berendsen thermostat or barostat?  The
>>>>>>>>
>>>>>>> gromacs
>>>>
>>>>> g_energy functions for heat capacity use the fluctuation formula,
>>>>>>>>
>>>>>>> and
>>>
>>>> the
>>>>>
>>>>>> fluctuations with both of these algorithms are wrong (as should be
>>>>>>>> printed
>>>>>>>> in the log file warning message). Make sure you use
>>>>>>>>
>>>>>>> ensemble-preserving
>>>>
>>>>> thermostats if you want fluctuation properties.
>>>>>>>>
>>>>>>>> Also, to be sure, you should double check by calculating both heat
>>>>>>>> capacities by finite difference formulas as well with two
>>>>>>>>
>>>>>>> simulations
>>>
>>>> at
>>>>>
>>>>>> T+dt/2 and T-dt/2 -- if the fluctuation and finite difference
>>>>>>>>
>>>>>>> resutls
>>>
>>>> don't
>>>>>>>> agree within propagated error, then something is off.
>>>>>>>>
>>>>>>>>
>>>>>>>> On Mon, May 25, 2015 at 5:59 AM, Faezeh Pousaneh <
>>>>>>>>
>>>>>>> fpoosaneh at gmail.com>
>>>>
>>>>> wrote:
>>>>>>>>
>>>>>>>>  Hi,
>>>>>>>>>
>>>>>>>>> I do not know why I obtain two difference cp and cv from NVT and
>>>>>>>>>
>>>>>>>> NPT
>>>>
>>>>> simulations.
>>>>>>>>> What I do is, I take 1000 lutidne molecules, and I do firstly an
>>>>>>>>>
>>>>>>>> energy
>>>>>
>>>>>> minimization with steep integrator, then NPT simulation at T=300
>>>>>>>>>
>>>>>>>> and
>>>>
>>>>> P=1
>>>>>>>>
>>>>>>>>> atm for 10ns, (I obtain Cp= 230), then I run NVT for 10 ns with
>>>>>>>>>
>>>>>>>> same
>>>>
>>>>> mdp
>>>>>>>>
>>>>>>>>> file except no pressure coupling, and with initial .gro file
>>>>>>>>>
>>>>>>>> obtained
>>>>
>>>>> from
>>>>>>>>
>>>>>>>>> NPT run, (I obtain Cv=180).
>>>>>>>>> Does some one know where is my mistake? (In both runs, I obtain
>>>>>>>>>
>>>>>>>> Cv
>>>
>>>> and
>>>>>
>>>>>> Cp
>>>>>>>>
>>>>>>>>> from g_energy and in different time intervals and after
>>>>>>>>>
>>>>>>>> equilibrited
>>>>
>>>>> time)
>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> Best regards
>>>>>>>>> --
>>>>>>>>> Gromacs Users mailing list
>>>>>>>>>
>>>>>>>>> * Please search the archive at
>>>>>>>>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List
>>>>>>>>>
>>>>>>>> before
>>>
>>>> posting!
>>>>>>>>>
>>>>>>>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>>>>>>>
>>>>>>>>> * For (un)subscribe requests visit
>>>>>>>>>
>>>>>>>>>
>>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users
>>>
>>>> or
>>>>>
>>>>>> send a mail to gmx-users-request at gromacs.org.
>>>>>>>>>
>>>>>>>>>  --
>>>>>>>> Gromacs Users mailing list
>>>>>>>>
>>>>>>>> * Please search the archive at
>>>>>>>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>>>>>>>> posting!
>>>>>>>>
>>>>>>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>>>>>>
>>>>>>>> * For (un)subscribe requests visit
>>>>>>>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users
>>>>>>>>
>>>>>>> or
>>>>
>>>>> send a mail to gmx-users-request at gromacs.org.
>>>>>>>>
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>  --
>>>>> Gromacs Users mailing list
>>>>>
>>>>> * Please search the archive at
>>>>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>>>>> posting!
>>>>>
>>>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>>>
>>>>> * For (un)subscribe requests visit
>>>>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>>>>> send a mail to gmx-users-request at gromacs.org.
>>>>>
>>>>>  --
>>>> Gromacs Users mailing list
>>>>
>>>> * Please search the archive at
>>>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>>>> posting!
>>>>
>>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>>
>>>> * For (un)subscribe requests visit
>>>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>>>> send a mail to gmx-users-request at gromacs.org.
>>>>
>>>>  --
>>> Gromacs Users mailing list
>>>
>>> * Please search the archive at
>>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>>> posting!
>>>
>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>
>>> * For (un)subscribe requests visit
>>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>>> send a mail to gmx-users-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
>
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>

More information about the gromacs.org_gmx-users mailing list