[gmx-users] TPI and chemical potential
Gmx QA
gmxquestions at gmail.com
Wed Aug 17 10:19:02 CEST 2016
Thanks again
In my hands at least, -pbc mol -ur compact seems to be necessary, as
without it I get a lot of <infs> instead of actual numbers.
I tried to find a reference for the chemical potential of TIP3P, but
haven't succeeded yet. Related to gromacs, I find one thread in the
archives that claims that for TIP3p as well as SPC/E the agreement with
experiment is reasonable, but with problems for TIP5P.
https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2014-May/089540.html
Also, does it matter if the rerun is performed on a NVT or NPT trajectory?
This I could of course easily test, but it seems that formally at least the
Widow formula is applicable to NVT only?
Anyway, I think I now know enough to proceed. How reliable are TPI results
if I want to measure the chemical potential of a drug molecule? How hard
are they to converge?
Best
/PK
2016-08-16 16:38 GMT+02:00 João M. Damas <jmdamas at itqb.unl.pt>:
> I don't recall any necessary post-processing of the trajectory being
> necessary, as it should deal well with PBC conditions.
>
> I really don't know what is happening, but it seems like the fixes are
> going in the right direction.
>
> Try to plot the <mu> along time, to see if it's converging or not.
> Depending on the trend, you may want to extend the simulation time. Other
> reasons for incorrect excess chemical potential may be that the TIP3p model
> may not reproduce well that quantity.
>
> And yeah, it's kJ/mol.
>
> João
>
> On Tue, Aug 16, 2016 at 2:12 PM, Gmx QA <gmxquestions at gmail.com> wrote:
>
> > Thanks Joao
> >
> > The water box is equilibrated at 298K (I checked density and temperature
> > with gmx energy, and they show no drift from reasonable averages).
> >
> > However, I am a bit confused as I just tried the approach on the same
> > trajectory, but extended to 2.5 ns and after one pass through trjconv
> with
> > -pbc mol -ur compact options.
> >
> > This seems to work in the sense that I no longer get any infs, but the
> > value of <mu>, while having the right magnitude, is still off by about a
> > factor of 2.
> > After 2.5 ns I get <mu> = -45 (I assume the unit is kJ/mol), whereas if i
> > read correctly the chemical potential should be -23.5 or so. Can it be
> that
> > I now need to increase the run time of my water simulation even more?
> >
> > Thanks
> > /PK
> >
> >
> >
> > 2016-08-16 13:47 GMT+02:00 João M. Damas <jmdamas at itqb.unl.pt>:
> >
> > > Well, 50000 insertions per frame, doesn't look bad. 1 ns should not be
> an
> > > issue.
> > >
> > > Well, actually that value normally starts at inf and should go to the
> > > equilibrium <mu> as you increase sampling. But to be inf after 50000
> > > inserts, it does sound weird.
> > >
> > > How did you build your box of water? Is it equilibrated at that
> > temperature
> > > (298 K)?
> > > Please make make that the water molecule to be inserted has the correct
> > > geometrical coordinates and that it is centered at 0,0,0.
> > >
> > > João
> > >
> > > On Tue, Aug 16, 2016 at 12:45 PM, Gmx QA <gmxquestions at gmail.com>
> wrote:
> > >
> > > > Hi Joao
> > > >
> > > > Thank you for your reply. My mdp-file is pasted below:
> > > >
> > > > integrator = tpi
> > > > emtol = 1000.0
> > > > emstep = 0.01
> > > > nsteps = 50000
> > > >
> > > > nstlist = 1
> > > > cutoff-scheme = group
> > > > ns_type = grid
> > > > coulombtype = pme
> > > > rcoulomb = 1.0
> > > > rvdw = 1.0
> > > > pbc = xyz
> > > > nstxtcout = 5000
> > > >
> > > > ; Temperature coupling is on
> > > > tcoupl = V-rescale
> > > > tc-grps = system
> > > > tau_t = 0.1
> > > > ref_t = 298
> > > >
> > > > ; Pressure coupling is on
> > > > pcoupl = no
> > > > pcoupltype = isotropic
> > > > tau_p = 2.0
> > > > ref_p = 1.0
> > > > compressibility = 4.5e-5
> > > >
> > > > ld_seed = -1
> > > >
> > > > I am using the amber99sb forcefield for tip3p, so I believe the 1.0
> nm
> > > > cutoffs to be correct?
> > > >
> > > > I have tried a couple different values for the number of insertions,
> > and
> > > > for some frames different values sometimes give non-inf results, but
> > for
> > > > the most part it does not seem to make a difference.
> > > >
> > > > 1 ns is short, I agree, but since the value of <mu> (which I take to
> be
> > > the
> > > > excess chemical potential I am after) immediately goes to inf when
> the
> > > > first mu=inf appears, extending does not make sense to me? Is the
> > > chemical
> > > > potential reported elsewhere as well?
> > > >
> > > > Thanks
> > > > /PK
> > > >
> > > > 2016-08-16 12:35 GMT+02:00 João M. Damas <jmdamas at itqb.unl.pt>:
> > > >
> > > > > Yes, all water atoms can't be with coordinates 0,0,0. I don't know
> > how
> > > > did
> > > > > that even work... Weird.
> > > > >
> > > > > 1 ns is very short time scale. Also, how many insertions are you
> > trying
> > > > per
> > > > > frame? That could also help improve the sampling.
> > > > >
> > > > > The inf values are normal as it just means that the water is too
> > > > structured
> > > > > and it's hard to insert a new water from the sampling point of
> view.
> > > > >
> > > > > Can you post the .mdp you used for the tpi?
> > > > >
> > > > > João
> > > > >
> > > > > On Tue, Aug 16, 2016 at 9:52 AM, Gmx QA <gmxquestions at gmail.com>
> > > wrote:
> > > > >
> > > > > > Please - anyone?
> > > > > >
> > > > > > I gathered from the mail-list that maybe the entire molecule
> (water
> > > in
> > > > > this
> > > > > > case) should be centered on origo (rather than having all atoms
> at
> > > > > (0,0,0),
> > > > > > but this gives me only a lot of inf for the chemical potential.
> > > > > >
> > > > > > How can I calculate the excess chemical potential of water using
> > TPI
> > > in
> > > > > > gromacs?
> > > > > >
> > > > > >
> > > > > >
> > > > > > 2016-08-15 13:31 GMT+02:00 Gmx QA <gmxquestions at gmail.com>:
> > > > > >
> > > > > > > Dear list
> > > > > > >
> > > > > > > I am trying to teach myself how the test particle method for
> > excess
> > > > > > > chemical potential calculations work. To this end, I created a
> > > small
> > > > > > system
> > > > > > > of about 900 water molecules (TIP3p), and simulated for 1 ns.
> > > > > > >
> > > > > > > I then set up files for inserting an extra water molecule to
> > > > calculate
> > > > > > mu,
> > > > > > > the chemical potential.
> > > > > > >
> > > > > > > The end of my tpi_start.gro looks like this:
> > > > > > > <snip>
> > > > > > > 884SOL OW 2650 2.466 0.788 2.747 -0.5599 -0.0387
> > > 0.0570
> > > > > > > 884SOL HW1 2651 2.529 0.835 2.802 0.0024 -0.9901
> > > 0.2261
> > > > > > > 884SOL HW2 2652 2.506 0.703 2.732 -1.4472 -0.3598
> > > -0.5022
> > > > > > > 885SOL OW 2653 0.000 0.000 0.000 0.0000 0.0000
> > > 0.0000
> > > > > > > 885SOL HW1 2654 0.000 0.000 0.000 0.0000 0.0000
> > > 0.0000
> > > > > > > 885SOL HW2 2655 0.000 0.000 0.000 0.0000 0.0000
> > > 0.0000
> > > > > > > 3.01188 3.01188 3.01188
> > > > > > >
> > > > > > > So I think this last water molecule is the one to be inserted.
> > The
> > > > > > > topology was also updated accordingly.
> > > > > > >
> > > > > > > I then made a rerun like this:
> > > > > > >
> > > > > > > $ gmx mdrun -v -deffnm tpi -rerun md.xtc
> > > > > > >
> > > > > > > And the output is a series of lines like this:
> > > > > > >
> > > > > > > Reading frame 180 time 900.000 mu 8.924e+00 <mu>
> > 7.240e+00
> > > > > > > mu 7.671e+00 <mu> 7.242e+00
> > > > > > > mu 6.987e+00 <mu> 7.241e+00
> > > > > > > mu 7.010e+00 <mu> 7.239e+00
> > > > > > > mu 7.691e+00 <mu> 7.241e+00
> > > > > > > mu 7.439e+00 <mu> 7.242e+00
> > > > > > > mu 6.757e+00 <mu> 7.240e+00
> > > > > > > mu 8.114e+00 <mu> 7.243e+00
> > > > > > > mu 7.173e+00 <mu> 7.243e+00
> > > > > > > mu 8.768e+00 <mu> 7.249e+00
> > > > > > > Reading frame 190 time 950.000 mu 7.799e+00 <mu>
> > 7.252e+00
> > > > > > > mu 9.284e+00 <mu> 7.259e+00
> > > > > > > mu 8.103e+00 <mu> 7.262e+00
> > > > > > > mu 7.835e+00 <mu> 7.265e+00
> > > > > > > mu 7.321e+00 <mu> 7.265e+00
> > > > > > > mu 7.576e+00 <mu> 7.267e+00
> > > > > > > mu 9.348e+00 <mu> 7.274e+00
> > > > > > > mu 7.021e+00 <mu> 7.272e+00
> > > > > > > mu 7.162e+00 <mu> 7.272e+00
> > > > > > > mu 7.695e+00 <mu> 7.274e+00
> > > > > > > Reading frame 200 time 1000.000 mu 7.104e+00 <mu>
> > 7.273e+00
> > > > > > >
> > > > > > > Now, is <mu> here the average computed chemical potential? What
> > is
> > > > the
> > > > > > TPI
> > > > > > > energy distribution being reported in the tpi.xvg-file?
> > > > > > >
> > > > > > > I tried to google for the chemical potential of water and
> > according
> > > > to
> > > > > > > e.g. [1] (page 13) it should be around -23.5 kJ/mol. I would
> > really
> > > > > > > appreciate if someone could comment on this, as I need to
> > > understand
> > > > it
> > > > > > > before moving on to more relevant (for me) systems….
> > > > > > >
> > > > > > > Thanks
> > > > > > > /PK
> > > > > > >
> > > > > > >
> > > > > > >
> > > > > > >
> > > > > > >
> > > > > > > [1] Chemical potential of liquids and mixtures via Adaptive
> > > > Resolution
> > > > > > ...
> > > > > > > <https://www.google.se/url?sa=t&rct=j&q=&esrc=s&source=web&
> > > > > > cd=21&ved=0ahUKEwj-7-2gqMPOAhVFWiwKHdp4DYY4FBAWCBww
> > > > > > AA&url=https%3A%2F%2Fopus4.kobv.de%2Fopus4-zib%2Ffiles%
> > > > > > 2F5097%2FZIB-Report_14-25.pdf&usg=AFQjCNECRk9dif8TGxKJeeztHV6T_
> > > > > > FmVuw&sig2=-vsjUB9HRdlcnt4vyKWcbw>
> > > > > > >
> > > > > > --
> > > > > > 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.
> > > > > >
> > > > >
> > > > >
> > > > >
> > > > > --
> > > > > João M. Damas
> > > > > PhD Student
> > > > > Protein Modelling Group
> > > > > ITQB-UNL, Oeiras, Portugal
> > > > > Tel:+351-214469613
> > > > > --
> > > > > 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.
> > > >
> > >
> > >
> > >
> > > --
> > > João M. Damas
> > > PhD Student
> > > Protein Modelling Group
> > > ITQB-UNL, Oeiras, Portugal
> > > Tel:+351-214469613
> > > --
> > > 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.
> >
>
>
>
> --
> João M. Damas
> PhD Student
> Protein Modelling Group
> ITQB-UNL, Oeiras, Portugal
> Tel:+351-214469613
> --
> 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.
>
More information about the gromacs.org_gmx-users
mailing list