[gmx-users] TPI and chemical potential
Gmx QA
gmxquestions at gmail.com
Tue Aug 16 14:12:44 CEST 2016
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.
>
More information about the gromacs.org_gmx-users
mailing list