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 <sp...@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 <fpoosa...@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 <mrshi...@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 <fpoosa...@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 <fpoosa...@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 < >>>>>> >>>>> fpoosa...@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 <mrshi...@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 < >>>>>>>> >>>>>>> fpoosa...@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-requ...@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-requ...@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-requ...@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-requ...@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-requ...@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. > sp...@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-requ...@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-requ...@gromacs.org.