[gmx-users] PMF using umbrella sampling and Gromacs 5.0
Christopher Neale
chris.neale at alum.utoronto.ca
Tue Jul 28 15:27:07 CEST 2015
Dear Eudes:
then I return to being stumped. It sure looks like it might be an issue with a constant somewhere. I'll note that v5.0.5 g_wham uses what I think is an incorrect constant to convert to kcal (http://redmine.gromacs.org/issues/1787), but you're using kj so that shouldn't be it. The gas constant defined by g_wham on line 1086 in v5.0.5 is also a little bit different from the one used in core gromacs (for e.g. replica exchange to get the boltzmann constant), though again I don't see how this could do anything but affect the overall scale by a tiny amount.
However, I do wonder about "8.314e-3" in lines like this where it possibly could affect the value:
denom += invg*window[j].N[k]*exp(-U/(8.314e-3*opt->Temperature) + window[j].z[k]);
(line 966 in calc_profile() in v5.0.5)
Are you using g_wham? Can you try a different wham software and see if you get the same thing?
Also, are you setting -temp properly in g_wham ?
Beyond that, I've run out of ideas. If you do ever sort this out, I'd appreciate it if you could send me an email off list to let me know (and post it here of course, I just don't always check the list).
Thank you,
Chris.
________________________________________
From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se <gromacs.org_gmx-users-bounces at maillist.sys.kth.se> on behalf of Eudes Fileti <fileti at gmail.com>
Sent: 28 July 2015 07:59
To: gromacs.org_gmx-users at maillist.sys.kth.se
Subject: Re: [gmx-users] PMF using umbrella sampling and Gromacs 5.0
Dear Chris, thank you for help. I performed the test you suggested, with no
pressure coupling. Note that the behavior of the profile remained the same
even with no pressure coupling.
https://goo.gl/photos/2d8tsYFwp2eHCuCWA
The simulation parameters were exactly the same (except the sampling in the
test was 5 ns window and by coupling pressure, which was turned off).
Bests
eef
_______________________________________
Eudes Eterno Fileti
Instituto de Ciência e Tecnologia da UNIFESP
Rua Talim, 330, São José dos Campos - SP
Página: sites.google.com/site/fileti/
> ------------------------------
>
> Message: 2
> Date: Sun, 26 Jul 2015 02:45:50 +0000
> From: Christopher Neale <chris.neale at alum.utoronto.ca>
> To: "gmx-users at gromacs.org" <gmx-users at gromacs.org>
> Subject: Re: [gmx-users] PMF using umbrella sampling and Gromacs 5.0
> (Eudes Fileti)
> Message-ID:
> <
> BLUPR03MB184037E6EF2EC6D63BEBE21C58F0 at BLUPR03MB184.namprd03.prod.outlook.com
> >
>
> Content-Type: text/plain; charset="iso-8859-1"
>
> Dear Eudes:
>
> Glad that you solved one of the two issues. As for the bumps in the PMF, I
> have a new theory: the bumps come from pressure coupling. When the sampled
> distance, d, between the two molecules fluctuate a little closer than the
> center of restraint, d0, that adds a repulsive force that contributes to
> the virial and the box gets a little larger. Conversely, slight
> fluctuations of d that are larger than d0 will add a small bias to box
> contraction. This should be more noticeable when the restrained distance
> involves larger masses.
>
> It is at the moment unclear to me whether this might exert an effect
> indirectly due to overall system density or more directly as coordinate
> scaling impacts the instantaneous value of d. If it is the latter, then
> semi-isotropic pressure coupling, may also enhance the effect since the
> virial will be computed independently along the order parameter (I presume)
> and hence there is less noise from other dimensions.
>
> Can you please try again without pressure coupling (single precision
> should be fine for this test). Hopefully this is not the source of the
> bumps because, if it is affecting the PMF noticeably and the underlying
> free energy surface has a large gradient, then d will always be on one side
> of d0 and the effect will not be merely bumps but also some type of bias in
> the PMF. Whether this bias is accurate or artifactual falls outside of my
> mathematical abilities. The thing is, the force is a real force between
> real atoms so it seems like it really should be included in the virial (as
> it certainly is... I checked). I can tell you one thing for sure: the
> effect on box volume is real and noticeable. That is, if you look at the
> average system volume when d<<d0, it differs in a statistically significant
> manner from the average system volume when d>>d0 (something that I also
> checked).
>
> Thank you for looking into this further,
> Chris.
>
> ________________________________________
> From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se <
> gromacs.org_gmx-users-bounces at maillist.sys.kth.se> on behalf of Eudes
> Fileti <fileti at gmail.com>
> Sent: 24 July 2015 16:47
> To: gromacs.org_gmx-users at maillist.sys.kth.se
> Subject: Re: [gmx-users] PMF using umbrella sampling and Gromacs 5.0
> (Eudes Fileti)
>
> Hello Chris, I write to report the results of the tests you suggested.
> To recap, I have two problems to solve. 1) The bad sampling around z = 0
> and 2) the bumps along to the profile.
>
> I solved the first discarding all the my initial configurations and
> performing a new pulling (SMD). Only this time I used a higher force
> constant (5000 kJ/mol nm2). Thus I got configurations close to z = 0, where
> they were not generated.
>
> For the second problem, you suggested recalculate the PMF using double
> precision. The results of this test showed that it does not solve the
> problem, on the contrary the bumps were even more pronounced, as indicated
> by the plot in this link. https://goo.gl/photos/wGNhdyNG9pdqGfcB6
>
> All tests were performed with prototypes simulations, with 40 windows of 2
> ns spaced by 0.1 nm. A larger sample is obviously necessary for reliable
> profiles, but this was enough to show the trend that I wanted to watch.
>
> As I have mentioned before, I've done several tests aiming to eliminate
> these bumps: use of higher sampling, up to 20ns per window; reducing the
> spacing between the windows (from 0.1 to 0.05 nm); changing the spring
> constant from 1000 up to 5000 kJ/mol nm2, use of two different versions of
> the Gromacs (4.6 e 5.0) and also I tested it with two different sets of
> initial settings.
>
> None of this attempts solved the bumps problem.
>
> If you (or someone else) have any other tips please let me know.
>
> Thank you
> eef
>
> _______________________________________
> Eudes Eterno Fileti
> Instituto de Ci?ncia e Tecnologia da UNIFESP
> Rua Talim, 330, S?o Jos? dos Campos - SP
> P?gina: sites.google.com/site/fileti/
>
> >
> > ------------------------------
> >
> > Message: 4
> > Date: Wed, 22 Jul 2015 14:56:50 +0000
> > From: Christopher Neale <chris.neale at alum.utoronto.ca>
> > To: "gmx-users at gromacs.org" <gmx-users at gromacs.org>
> > Subject: Re: [gmx-users] PMF using umbrella sampling and Gromacs 5.0
> > Message-ID:
> > <
> >
> BLUPR03MB18470A33A4E3FF149C95963C5830 at BLUPR03MB184.namprd03.prod.outlook.com
> > >
> >
> > Content-Type: text/plain; charset="iso-8859-1"
> >
> > Dear Eudes:
> >
> > There are two issues: The first issue is the fact that you've got a
> > sampling problem near the bilayer center. The second issue is the
> periodic
> > bumps that you see in your PMFs. I'll take the second part first.
> >
> > The source of the periodic bumps in PMFs from umbrella sampling is, to
> me,
> > a million dollar question. I've seen them in my own work. I've seen then
> in
> > the literature (as you noted for Fig. 4 in
> > http://www.mdpi.com/1422-0067/13/11/14451/htm ). I've seen them when
> > people use g_wham and Alan Grossfield's version of WHAM. What I don't yet
> > know is if people also see this when running US simulations with AMBER,
> > CHARMM, or NAMD. Frankly, either answer scares me a little. If you see
> wild
> > oscillations of the PMF, at large distances then one possible source is
> > that you are near a distance that is half your box length along the order
> > parameter. Note that with constant pressure simulations you will get
> > oscillations in the length of the box and also when your membrane changes
> > shape it could flip the vector of closest approach by 180 deg if you are
> > close to the half-box limit. However, that should not be relevant to the
> > bumps (out of interest, how large is your box along z?). Although I do
> > n't know what is going on with these PMF bumps, I also note that
> > membranes as a whole always tend to migrate toward positive z and lipids
> > tend to flow to positive x in gromacs simulations, so I wonder if it is a
> > rounding issue. Could you try again with double precision and see if you
> > get the bumps?
> >
> > As for the problem with sampling near the bilayer center... my first
> guess
> > is that you've got some of your replicas on the wrong side of the
> bilayer's
> > center. Did you intend to go across the center and sample also in the
> lower
> > leaflet? I've only ever used gromacs 3 and 4 to do pulling simulations,
> > never gmx 5. You can see my concerns about gmx5 for this type of issue
> near
> > the bilayer center here:
> >
> https://www.mail-archive.com/gromacs.org_gmx-users%40maillist.sys.kth.se/msg11324.html
> > One simply has to test it. However, if gmx5 does indeed work properly
> with
> > this new US philosophy, then I presume that you simply need to move your
> > starting coordinates such that your solute is in the positive leaflet for
> > all replicas. You can test this by visualization or g_dist, which reports
> > the sign of the displacement (look in the 5th column I think for the z
> > value). Note that you are using pull_start=yes, so which side of the
> > bilayer your solute is initially on should make a difference here.
> >
> > Chris.
> >
> > -- original message:
> >
> > Hi Chris, a few days ago I posted a question on GMX list but
> unfortunately
> > I not received an answer yet. So I write to you for help at your
> > convenience.
> > http://permalink.gmane.org/gmane.science.biology.gromacs.user/78439
> >
> > Besides the problem exposed the link above I'm trying to understand the
> > importance of the relative distance between the COM. Due to fluctuation
> in
> > the position of the center of mass position of the membrane, the minimum
> > distance is never zero, as shown in the graphs of my pulling below. What
> is
> > the effect of this difference in the calculation of the PMF? Could this
> be
> > the source of my problem?
> > https://goo.gl/photos/cwwLzog6wXAgVJG88
> >
> > Bests
> > eef
> >
> > PS. I found this paper below where authors have published this unusual
> > behavior (Figure 4). They did not explain the reason for the oscillations
> > at large distances and I also do not know why this occurs.
> > http://www.mdpi.com/1422-0067/13/11/14451/htm
> >
> >
> >
> > ________________________________________
> > From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se <
> > gromacs.org_gmx-users-bounces at maillist.sys.kth.se> on behalf of Eudes
> > Fileti <fileti at gmail.com>
> > Sent: 18 July 2015 21:06
> > To: gmx-users at gromacs.org
> > Subject: [gmx-users] PMF using umbrella sampling and Gromacs 5.0
> >
> > Hello guys, I'm trying to calculate the PMF using umbrella sampling for a
> > small molecule to penetrate a lipid membrane. The 1D reaction coordinate
> is
> > along the z axis, which corresponds to the bilayer normal. The umbrella
> > potential acts on the center of mass of the molecule. The initial
> > configurations for each window, separated by a distance of 0.1nm, in a
> > range of 4 nm (distance between the centers of mass of the membrane and
> the
> > molecule), were obtained by a SMD.
> >
> > The umbrella potential was applied according to the parameters below (for
> > Gromacs version 5.0.3). Each window was sampled by 5 ns.
> >
> > ; COM PULLING
> > pull = umbrella
> > pull_geometry = direction
> > pull_dim = N N Y
> > pull_start = yes
> > pull-print-reference = no
> > pull_nstxout = 1000
> > pull_nstfout = 1000
> > pull_ngroups = 2
> > pull-ncoords = 1
> > ; Group name, weight (default all 1), vector, init, rate (nm/ps),
> > kJ/(mol*nm^2)
> > pull-group1-name = DPPC ; ref
> > pull-group2-name = LIG ; pulled
> > pull-coord1-groups = 1 2
> > pull-coord1-origin = 0 0 0
> > pull-coord1-vec = 0 0 1
> > pull-coord1-init = 0
> > pull-coord1-rate = 0.0
> > pull-coord1-k = 3000
> >
> > Histograms and the PMF obtained are shown in figure the link below.
> > https://goo.gl/photos/RkW9gbWrgeKEfV3R7
> >
> > I repeated the same test using other parameters (larger and small k
> values)
> > and options (cylinder and position, for this later I have used gmx 4.6)
> and
> > so far I could not get a satisfactory profile. In all tests, a problem
> that
> > I have observed is that the region near z = 0 is not sampled (there is a
> > gap between 0.0 and 0.5 nm). In addition the profile presents not smooth,
> > but somehow oscillating mainly for large z.
> >
> > As the gmx distance give a value different of zero for the smallest
> > distance of separation between the centers of mass (around 0.1-0.2 nm) I
> > believe that this weird behavior is related to the reference distance
> that
> > I have used.
> >
> > Could someone give me a light? Most of the tips I read in GMX list are
> > related to previous versions of the Gromacs and in a way the other tips
> > were already included in my tests.
> >
> > Thank you
> > eef
> >
> > _______________________________________
> > Eudes Eterno Fileti
> > Instituto de Ci?ncia e Tecnologia da UNIFESP
> > Rua Talim, 330, S?o Jos? dos Campos - SP
> > P?gina: sites.google.com/site/fileti/
> > --
> > 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.
>
> ------------------------------
>
> Message: 3
> Date: Sun, 26 Jul 2015 10:15:51 +0100
> From: Dawid das <addiw7 at googlemail.com>
> To: gmx-users at gromacs.org
> Subject: Re: [gmx-users] Default cut-off scheme.
> Message-ID:
> <
> CAKSLqn60ybuXhuDjbgNDFDsm1caG4OqUFgeBcz4hg7xigJZijQ at mail.gmail.com>
> Content-Type: text/plain; charset=UTF-8
>
> Actually, I am working with Gromacs 4.6.7, but thank you for a tip ;).
> I forgot about the actual manual.
>
> All the best,
>
> Dawid Grabarek
>
> 2015-07-25 23:44 GMT+01:00 Mark Abraham <mark.j.abraham at gmail.com>:
>
> > On Sun, Jul 26, 2015 at 12:32 AM Tamisra Pal <tamisra1985 at gmail.com>
> > wrote:
> >
> > > Hi ,
> > >
> > > The cutoff-scheme by default is Verlet , unless you specify it as
> group.
> > > Yes , in this case , the ns_type is grid
> > >
> >
> > ... assuming this is version 5. But Dawid will get an accurate answer
> much
> > faster by consulting the manual for the version he is using ;-)
> >
> > Mark
> >
> >
> > > On Sun, Jul 26, 2015 at 1:33 AM, Dawid das <addiw7 at googlemail.com>
> > wrote:
> > >
> > > > Dear Gromacs Experts,
> > > >
> > > > I was issued with following minimize.mdp file:
> > > >
> > > > ; Minimization of bovine pancreatic trypsin inhibitor
> > > > title = minimization of BPTI
> > > > ; this is a comment
> > > > integrator = l-bfgs ; (other: cg, l-bfgs)
> > > > emtol = 500 ; convergence criterium- maximum value
> > of
> > > > force [kJ/(mol*nm)]
> > > > emstep = 0.01 ; initial step size [nm]
> > > > nsteps = 2000 ; number of steps
> > > > nstlist = 1 ; frequency update of neighbour list
> > > > ns_type = grid ; algorithm of updating neighbour
> list.
> > > It
> > > > seems to be cell list method.
> > > > nstenergy = 50 ;
> > > > rlist = 1.0
> > > > coulombtype = PME
> > > > rcoulomb = 1.0
> > > > rvdw = 1.0
> > > > pbc = xyz
> > > >
> > > > my question is what is the default cutoff-scheme here? group or
> > Verlet? I
> > > > assume there has to be default option because otherwhise nstlist and
> > > > ns_type will not make sense, won't they?
> > > > The other question is whether I understand ns_type=grid correctly. Is
> > it
> > > > cell list method as I commented?
> > > >
> > > > Best wishes,
> > > >
> > > > Dawid Grabarek
> > > > --
> > > > 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.
> > >
> > --
> > 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.
>
> End of gromacs.org_gmx-users Digest, Vol 135, Issue 136
> *******************************************************
>
--
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