[gmx-users] Different Cv and Cp

Wed Jun 10 16:40:02 CEST 2015

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 <fpoosaneh at 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 <mrshirts at 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 <fpoosaneh at 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 <mrshirts at 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 <
> 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.
> > > > >
> > > > --
> > > > 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.
>

-- 
_____________

Prof. Dr. André Farias de Moura
Department of Chemistry
Federal University of São Carlos
São Carlos - Brazil
phone: +55-16-3351-8090