[gmx-users] TPI and chemical potential
João M. Damas
jmdamas at itqb.unl.pt
Tue Aug 16 16:39:16 CEST 2016
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
More information about the gromacs.org_gmx-users
mailing list