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