Dear Faezeh, Michael is absolutely right, both methods should yield exactly the same Cp and Cv values, but that is only true in the limit of infinite sampling. If sampling is not infinite but is long enough, Cp and Cv values should agree within the statistical uncertainty of each method. So it seems to me that you have improve your simulations until you get this agreement.
regardless of which method you are using to compute Cp and Cv and based on your last messages, I think that you should try some more "best practice testing" before deciding whether Cv and Cp should be the same or not (that's your decision, not ours): (1) I would suggest you not to take the last frame of the NPT simulation as the "correct volume". You should instead take an average volume from the equilibrated NPT simulation and then edit the gro file to reach that average volume. (2) are you sure that 10 ns is enough to properly sample heat capacities? Regarding this topic, I could ask you several questions: - Have you checked that all relevant structural and energetic patterns are fully relaxed within this time scale? For instance, energy components, density, radial distribution functions, relative orientation between molecules and so on. - Have you discarded the frames before structural and energetic relaxations were achieved? - Considering that lutidine is slightly bulky and has a well defined dipole moment, it may take a while for a full rotational relaxation. Have you computed the rotational correlation function to be sure that it relaxed long before 10 ns? If it did not relax, then 10 ns is not enough. - I like to compute these fluctuation-based properties as cumulative averages along the simulation (after structural and energetic relaxations, of course), it is easier to see if the average is approaching some plateau value (if not, longer simulation time is required). - I have never tested myself (has anyone tested?), but fluctuation-based calculations are obviously sensitive to the quality of the fluctuation, and double precision integration yields less noisier raw data than single precision integration does, so I would guess that double precision mdrun might provide a better-quality sampling of the fluctuation and would then improve the quality of your Cp and Cv estimates (just guessing). As I already stated (and provided you published experimental data in support of my statement), Cp and Cv need not be the same for molecular liquids. But it is up to you to prove that your simulations are reliable (my questions above are just a fraction of the checks you need to perform in order to convince yourself and your audience that you have obtained reliable data and not just simulation artifacts). In principle, you should include these checks in the supporting information of your manuscripts to convince the editors and referees (and the general audience after publication) that you have been thorough in your investigation. best Andre On Wed, Jun 10, 2015 at 9:54 AM, Faezeh Pousaneh <fpoosa...@gmail.com> wrote: > Dear Michael, > > Could you please comment on my last question, thanks a lot. > I have noticed that when I run both NVT and NPT simulations from a same > .gro file (obtained from energy minimization) I obtain same Cv and Cp for > two ensembles. However, so far I was running NVT after NPT (meaning I used > .gro file obtained from NPT as initial configuration for NVT, since I > wanted to have correct volume). The second way gives that big difference in > Cv and Cp. > I run a long simulations in all cases, but why both ways produce such > difference? > > > > Best regards > > > On Tue, Jun 9, 2015 at 5:45 PM, Michael Shirts <mrshi...@gmail.com> wrote: > > > The variance formula is derived from the derivative formula + the > > assumption the distribution in Boltzmann, so they must agree if the > > distribution is Boltzmann. > > > > On Tue, Jun 9, 2015 at 11:43 AM, Faezeh Pousaneh <fpoosa...@gmail.com> > > wrote: > > > > > Thank you so much for the reply. > > > > > > Yes, I use contact applied pressure and I am careful about units. I > > checked > > > and average enthalpy and U are close, meaning that PV is negligible. > But > > > the point is variance of enthalpy in NPT differs from variance of > energy > > in > > > NVT and that causes the difference. > > > > > > You had given me an article showing that for example for Benzene the > > > difference in Cv and Cp is 25%, and here Ii get similar ( my molecule > is > > > also carbon ring). But still I can answer why two both ways does not > give > > > same Cv-Cp. > > > I will follow your suggestion. > > > > > > > > > > > > Best regards > > > > > > > > > On Tue, Jun 9, 2015 at 5:29 PM, Michael Shirts <mrshi...@gmail.com> > > 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). > > > > > > > > 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. > > > > > > > > > -- > > > > 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. > -- _____________ Prof. Dr. André Farias de Moura Department of Chemistry Federal University of São Carlos São Carlos - Brazil phone: +55-16-3351-8090 -- 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.
