Re: [gmx-users] TPI and chemical potential

[João M . Damas]

I think the Widom particle insertion method can also apply to NPT as long
as the volume fluctuations are taken into account (see equation 2.69b in
Allen and Tildesley's Computer Simulation of Liquids ). From what I recall
from the code, I think they are.

Do not forget the result is the excess chemical potential and not the
chemical potential. I don't know if this is relevant for you or not.

TPI is not going to work for a big molecule. I would say impossible to
converge, but someone correct me if I am wrong. TPI is very limited in that
sense. You can try to explore more recent methods.

João

On Wed, Aug 17, 2016 at 10:18 AM, Gmx QA  wrote:

> Thanks again
>
> In my hands at least, -pbc mol -ur compact seems to be necessary, as
> without it I get a lot of  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 :
>
> > 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  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  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 , while having the right magnitude, is still off by about
> a
> > > factor of 2.
> > > After 2.5 ns I get  = -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 :
> > >
> > > > Well, 5 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  as you increase sampling. But to be inf after 5
> > > > 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 
> > wrote:
> > > >
> > > > > Hi Joao
> > > > >
> > > > > Thank you for your reply. My mdp-file is pasted below:
> > > > >
> > > > > integrator  = tpi
> > > > > emtol = 1000.0
> > > > > emstep  = 0.01
> > > > > nsteps= 5
> > > > >
> > > > > 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 

Re: [gmx-users] TPI and chemical potential

[Gmx QA]

Thanks again

In my hands at least, -pbc mol -ur compact seems to be necessary, as
without it I get a lot of  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 :

> 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  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  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 , while having the right magnitude, is still off by about a
> > factor of 2.
> > After 2.5 ns I get  = -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 :
> >
> > > Well, 5 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  as you increase sampling. But to be inf after 5
> > > 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 
> wrote:
> > >
> > > > Hi Joao
> > > >
> > > > Thank you for your reply. My mdp-file is pasted below:
> > > >
> > > > integrator  = tpi
> > > > emtol = 1000.0
> > > > emstep  = 0.01
> > > > nsteps= 5
> > > >
> > > > 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  (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 :
> > > >
> > > > > 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 

Re: [gmx-users] TPI and chemical potential

[João M . Damas]

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  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  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 , while having the right magnitude, is still off by about a
> factor of 2.
> After 2.5 ns I get  = -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 :
>
> > Well, 5 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  as you increase sampling. But to be inf after 5
> > 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  wrote:
> >
> > > Hi Joao
> > >
> > > Thank you for your reply. My mdp-file is pasted below:
> > >
> > > integrator  = tpi
> > > emtol = 1000.0
> > > emstep  = 0.01
> > > nsteps= 5
> > >
> > > 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  (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 :
> > >
> > > > 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 
> > 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 :
> > > > >
> > > > > > Dear list
> > > > > >
> > > > > > I am trying to teach myself how the test particle method for
> excess
> > > > > > chemical potential calculations work. To this end, I 

Re: [gmx-users] TPI and chemical potential

[Gmx QA]

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 , while having the right magnitude, is still off by about a
factor of 2.
After 2.5 ns I get  = -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 :

> Well, 5 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  as you increase sampling. But to be inf after 5
> 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  wrote:
>
> > Hi Joao
> >
> > Thank you for your reply. My mdp-file is pasted below:
> >
> > integrator  = tpi
> > emtol = 1000.0
> > emstep  = 0.01
> > nsteps= 5
> >
> > 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  (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 :
> >
> > > 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 
> 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 :
> > > >
> > > > > 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:
> > > > > 
> > > > >   884SOL OW 2650   2.466   0.788   2.747 -0.5599 -0.0387
> 0.0570
> > > > >   884SOLHW1 2651   2.529   0.835   2.802  0.0024 -0.9901
> 0.2261
> > > > >   884SOLHW2 2652   2.506   0.703   2.732 -1.4472 -0.3598
> -0.5022
> > > > >   885SOL OW 2653   0.000   0.000   0.000  0.  0.
> 0.
> > > > >   885SOLHW1 2654   0.000   0.000   0.000  0.  0.
> 0.
> > > > >   885SOLHW2 2655   0.000   0.000   0.000  0.  0.
> 0.
> > > > >3.01188   3.01188   3.01188
> > > > >

Re: [gmx-users] TPI and chemical potential

[João M . Damas]

Well, 5 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  as you increase sampling. But to be inf after 5
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  wrote:

> Hi Joao
>
> Thank you for your reply. My mdp-file is pasted below:
>
> integrator  = tpi
> emtol = 1000.0
> emstep  = 0.01
> nsteps= 5
>
> 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  (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 :
>
> > 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  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 :
> > >
> > > > 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:
> > > > 
> > > >   884SOL OW 2650   2.466   0.788   2.747 -0.5599 -0.0387  0.0570
> > > >   884SOLHW1 2651   2.529   0.835   2.802  0.0024 -0.9901  0.2261
> > > >   884SOLHW2 2652   2.506   0.703   2.732 -1.4472 -0.3598 -0.5022
> > > >   885SOL OW 2653   0.000   0.000   0.000  0.  0.  0.
> > > >   885SOLHW1 2654   0.000   0.000   0.000  0.  0.  0.
> > > >   885SOLHW2 2655   0.000   0.000   0.000  0.  0.  0.
> > > >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   7.240e+00
> > > > mu  7.671e+00   7.242e+00
> > > > mu  6.987e+00   7.241e+00
> > > > mu  7.010e+00   7.239e+00
> > > > mu  7.691e+00   7.241e+00
> > > > mu  7.439e+00   7.242e+00
> > > > mu  6.757e+00   7.240e+00
> > > > mu  8.114e+00   7.243e+00
> > > > mu  7.173e+00   7.243e+00
> > > > mu  8.768e+00   7.249e+00
> > > > Reading frame 190 time  950.000   mu  7.799e+00   7.252e+00
> > > > mu  9.284e+00   7.259e+00
> > > > mu  8.103e+00   7.262e+00
> > > > mu  7.835e+00   7.265e+00
> > > > mu  7.321e+00   7.265e+00
> > > > mu  7.576e+00   7.267e+00
> > > > mu  9.348e+00   7.274e+00
> > > > mu  7.021e+00   7.272e+00
> > > > mu  7.162e+00   7.272e+00
> > > > mu  

Re: [gmx-users] TPI and chemical potential

[Gmx QA]

Hi Joao

Thank you for your reply. My mdp-file is pasted below:

integrator  = tpi
emtol = 1000.0
emstep  = 0.01
nsteps= 5

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

> 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  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 :
> >
> > > 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:
> > > 
> > >   884SOL OW 2650   2.466   0.788   2.747 -0.5599 -0.0387  0.0570
> > >   884SOLHW1 2651   2.529   0.835   2.802  0.0024 -0.9901  0.2261
> > >   884SOLHW2 2652   2.506   0.703   2.732 -1.4472 -0.3598 -0.5022
> > >   885SOL OW 2653   0.000   0.000   0.000  0.  0.  0.
> > >   885SOLHW1 2654   0.000   0.000   0.000  0.  0.  0.
> > >   885SOLHW2 2655   0.000   0.000   0.000  0.  0.  0.
> > >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   7.240e+00
> > > mu  7.671e+00   7.242e+00
> > > mu  6.987e+00   7.241e+00
> > > mu  7.010e+00   7.239e+00
> > > mu  7.691e+00   7.241e+00
> > > mu  7.439e+00   7.242e+00
> > > mu  6.757e+00   7.240e+00
> > > mu  8.114e+00   7.243e+00
> > > mu  7.173e+00   7.243e+00
> > > mu  8.768e+00   7.249e+00
> > > Reading frame 190 time  950.000   mu  7.799e+00   7.252e+00
> > > mu  9.284e+00   7.259e+00
> > > mu  8.103e+00   7.262e+00
> > > mu  7.835e+00   7.265e+00
> > > mu  7.321e+00   7.265e+00
> > > mu  7.576e+00   7.267e+00
> > > mu  9.348e+00   7.274e+00
> > > mu  7.021e+00   7.272e+00
> > > mu  7.162e+00   7.272e+00
> > > mu  7.695e+00   7.274e+00
> > > Reading frame 200 time 1000.000   mu  7.104e+00   7.273e+00
> > >
> > > Now, is  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
> > ...
> > >  > cd=21=0ahUKEwj-7-2gqMPOAhVFWiwKHdp4DYY4FBAWCBww
> > 

Re: [gmx-users] TPI and chemical potential

[João M . Damas]

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  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 :
>
> > 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:
> > 
> >   884SOL OW 2650   2.466   0.788   2.747 -0.5599 -0.0387  0.0570
> >   884SOLHW1 2651   2.529   0.835   2.802  0.0024 -0.9901  0.2261
> >   884SOLHW2 2652   2.506   0.703   2.732 -1.4472 -0.3598 -0.5022
> >   885SOL OW 2653   0.000   0.000   0.000  0.  0.  0.
> >   885SOLHW1 2654   0.000   0.000   0.000  0.  0.  0.
> >   885SOLHW2 2655   0.000   0.000   0.000  0.  0.  0.
> >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   7.240e+00
> > mu  7.671e+00   7.242e+00
> > mu  6.987e+00   7.241e+00
> > mu  7.010e+00   7.239e+00
> > mu  7.691e+00   7.241e+00
> > mu  7.439e+00   7.242e+00
> > mu  6.757e+00   7.240e+00
> > mu  8.114e+00   7.243e+00
> > mu  7.173e+00   7.243e+00
> > mu  8.768e+00   7.249e+00
> > Reading frame 190 time  950.000   mu  7.799e+00   7.252e+00
> > mu  9.284e+00   7.259e+00
> > mu  8.103e+00   7.262e+00
> > mu  7.835e+00   7.265e+00
> > mu  7.321e+00   7.265e+00
> > mu  7.576e+00   7.267e+00
> > mu  9.348e+00   7.274e+00
> > mu  7.021e+00   7.272e+00
> > mu  7.162e+00   7.272e+00
> > mu  7.695e+00   7.274e+00
> > Reading frame 200 time 1000.000   mu  7.104e+00   7.273e+00
> >
> > Now, is  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
> ...
> >  cd=21=0ahUKEwj-7-2gqMPOAhVFWiwKHdp4DYY4FBAWCBww
> AA=https%3A%2F%2Fopus4.kobv.de%2Fopus4-zib%2Ffiles%
> 2F5097%2FZIB-Report_14-25.pdf=AFQjCNECRk9dif8TGxKJeeztHV6T_
> FmVuw=-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-requ...@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-requ...@gromacs.org.

Re: [gmx-users] TPI and chemical potential

[Gmx QA]

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 :

> 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:
> 
>   884SOL OW 2650   2.466   0.788   2.747 -0.5599 -0.0387  0.0570
>   884SOLHW1 2651   2.529   0.835   2.802  0.0024 -0.9901  0.2261
>   884SOLHW2 2652   2.506   0.703   2.732 -1.4472 -0.3598 -0.5022
>   885SOL OW 2653   0.000   0.000   0.000  0.  0.  0.
>   885SOLHW1 2654   0.000   0.000   0.000  0.  0.  0.
>   885SOLHW2 2655   0.000   0.000   0.000  0.  0.  0.
>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   7.240e+00
> mu  7.671e+00   7.242e+00
> mu  6.987e+00   7.241e+00
> mu  7.010e+00   7.239e+00
> mu  7.691e+00   7.241e+00
> mu  7.439e+00   7.242e+00
> mu  6.757e+00   7.240e+00
> mu  8.114e+00   7.243e+00
> mu  7.173e+00   7.243e+00
> mu  8.768e+00   7.249e+00
> Reading frame 190 time  950.000   mu  7.799e+00   7.252e+00
> mu  9.284e+00   7.259e+00
> mu  8.103e+00   7.262e+00
> mu  7.835e+00   7.265e+00
> mu  7.321e+00   7.265e+00
> mu  7.576e+00   7.267e+00
> mu  9.348e+00   7.274e+00
> mu  7.021e+00   7.272e+00
> mu  7.162e+00   7.272e+00
> mu  7.695e+00   7.274e+00
> Reading frame 200 time 1000.000   mu  7.104e+00   7.273e+00
>
> Now, is  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 ...
> 
>
-- 
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.

[gmx-users] TPI and chemical potential

[Gmx QA]

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:

  884SOL OW 2650   2.466   0.788   2.747 -0.5599 -0.0387  0.0570
  884SOLHW1 2651   2.529   0.835   2.802  0.0024 -0.9901  0.2261
  884SOLHW2 2652   2.506   0.703   2.732 -1.4472 -0.3598 -0.5022
  885SOL OW 2653   0.000   0.000   0.000  0.  0.  0.
  885SOLHW1 2654   0.000   0.000   0.000  0.  0.  0.
  885SOLHW2 2655   0.000   0.000   0.000  0.  0.  0.
   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   7.240e+00
mu  7.671e+00   7.242e+00
mu  6.987e+00   7.241e+00
mu  7.010e+00   7.239e+00
mu  7.691e+00   7.241e+00
mu  7.439e+00   7.242e+00
mu  6.757e+00   7.240e+00
mu  8.114e+00   7.243e+00
mu  7.173e+00   7.243e+00
mu  8.768e+00   7.249e+00
Reading frame 190 time  950.000   mu  7.799e+00   7.252e+00
mu  9.284e+00   7.259e+00
mu  8.103e+00   7.262e+00
mu  7.835e+00   7.265e+00
mu  7.321e+00   7.265e+00
mu  7.576e+00   7.267e+00
mu  9.348e+00   7.274e+00
mu  7.021e+00   7.272e+00
mu  7.162e+00   7.272e+00
mu  7.695e+00   7.274e+00
Reading frame 200 time 1000.000   mu  7.104e+00   7.273e+00

Now, is  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 ...

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