[gmx-users] Different Cv and Cp

Thu Jun 11 19:26:31 CEST 2015

Dear Andre,

Many thanks for the message, was very useful. I followed your comments
carefully. Below;

>
>
> (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.
>

I tried this, but it did not help. The average volume is so close to the
last configuration volume.

>
> (2) are you sure that 10 ns is enough to properly sample heat capacities?
> Regarding this topic, I could ask you several questions:
>
> Actually not only 10 ns, but I have also tried a run up to 80 ns. I check
that all energy parameters are well equilibriated

> - 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?
>
> Yes, I have done.

> - 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).

That part I did not know, thanks for the comment. So I just checked them,
and I see the correlations also get equlibrated.

>
>
- 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).
>
> This also I just checked, but it does not improve that much.

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.
>
> Yes, I read the paper you gave and they can be different in some
molecules. But still you are right and I am not convinced about my
simulation yet. So I am trying to solve the problem. What I realized during
the checking is, when I change thermostat from V-rescale to Nose-hoover, I
get different results. Moreover, tau_t is also crucial in my simulation. I
had before 0.1, but when I change it to 0.2 I get a bigger C_v.

> 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
> --
> 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.
>