[gmx-users] large error bars in PMF

chris.neale at utoronto.ca chris.neale at utoronto.ca
Fri Jul 22 14:00:19 CEST 2011


I don't see why the box-type makes any difference whatsoever. It is  
possible that if you use a rhombic dodecahedron, you may reduce the  
system size, thus simulate more ns/day, thus converge faster, but that  
should be the only effect. I would be interested to hear more from  
Justin about how the box-shape is expected to affect peptide  
rotation... perhaps I misunderstand this point.

I have a few other comments.

1. If you allow the peptide to rotate freely, then you do indeed need  
to converge all of their different rotational interactions. An  
alternative is to apply orientational restraints during the pulling  
(assuming that you know the bound state) and then to correct to an  
unrestrained state at large separations. You can see, for instance,   
D. L. Mobley, J. D. Chodera, K. A. Dill. "On the use of orientational  
restraints and symmetry number corrections in alchemical free energy  
calculations", ...

2. You are using "pull_dim = N N Y" which, if I haven't entirely  
forgotten how the pull-code works, means that the distance along Z is  
restrained but the distance along X and Y is free to change. What you  
end up with by sampling in this way is pretty strange and will require  
a really weird volumetric correction in the absence of infinite  
sampling time. You must decide to either: (i) pull to a spherical  
distance:

pull_dim = Y Y Y
pull_geometry = distance

or (ii) to pull along a defined vector

pull_dim = Y Y Y
pull_geometry = direction

What you have done:

pull_dim = N N Y
pull_geometry = distance

is only really useful when the system is isotropic along the XY plane  
(at least in the time averaged sense), such as for a lipid bilayer, or  
when the freedom in X and Y is very low, such as in a channel.

3. Finally, just because you sampled *more* doesn't mean that your  
values are converged. Look into block averaging and test to see if  
your binding free energy is drifting over time.

Good luck,
Chris.

-- original message --

Hi again,
I have one doub about the suggestion of using a dodecahedral box for  
my umbrella sampling to remove the problems I am having with the  
peptides rotating. I cannot see why a dodecaheral box is going to  
avoid this. Would it be better a truncated octahedron?
Thanks a lot for your help.
Best wishes,
Rebeca.

[Hide Quoted Text]
Date: Thu, 21 Jul 2011 15:16:52 -0400
From: jalemkul at vt.edu
To: gmx-users at gromacs.org
Subject: Re: FW: [gmx-users] large error bars in PMF



Rebeca García Fandiño wrote:


I am trying to achieve the binding energy of the dimer composed by the
two small cyclic peptides, to compare it with experimental. What
advantages would I have using 3D PMF instead only 1D for this calculation?
Intuitively, two molecules diffuse through solution until they find  
one another,
which to me sounds a lot like a 3D path.  Further, using a  
dodecahedral box for
your umbrella sampling removes the problems you're having with the peptides
rotating.  It sounds like you're trying to pull in one direction along a
rectangular box, but the peptides are not playing nice.  I feel like this
discussion has come up at least once or twice before, though...

-Justin
Thanks a lot!
Rebeca.

   > Date: Thu, 21 Jul 2011 14:14:44 -0400
   > From: jalemkul at vt.edu
   > To: gmx-users at gromacs.org
   > Subject: Re: [gmx-users] large error bars in PMF
   >
   >
   >
   > Rebeca García Fandiño wrote:
   > > Hi,
   > > thanks a lot for your quick answer.
   > > What I am trying to pull are two small peptides one from another (r_1
   > > and r_2).
   > > I did not understand very well your last suggestion: "...if you want
   > > reasonable error bars you will not lots of well-converged data".
   >
   > Oops, that should have read "you will *need* lots of well-converged
data."
   >
   > > Do you mean I will need also more windows besides extending the
simulations?
   >
   > I doubt you need more windows. Likely you just need more time in each.
   >
   > > I think the problem could be also that the peptides I am using
rotate in
   > > the box and they do not remain flat one respect to the other. They
   > > gyrate freely and some parts of their structure interact along the
   > > pulling...
   >
   > Interactions are part of the dissociation process and are not
problematic per
   > se. But if you're trying to obtain only a one-dimensional PMF then your
   > rotation could be a problem. Is there some reason you need a
one-dimensional
   > PMF and not a three-dimensional PMF? What are you trying to achieve?
   >
   > -Justin
   >
   > > Thanks a lot again for your help.
   > > Best wishes,
   > > Rebeca.
   > >
   > >
   > >
   > >
------------------------------------------------------------------------
   > > From: regafan at hotmail.com
   > > To: gmx-users at gromacs.org
   > > Date: Thu, 21 Jul 2011 16:36:59 +0000
   > > Subject: [gmx-users] large error bars in PMF
   > >
   > >
   > > Hi,
   > > I am trying to calculate the binding energy of two molecules using the
   > > PMF (Umbrella Sampling method) and Gromacs 4.0.
   > > Some weeks ago I have written to the list because changing the
number of
   > > windows used in the Umbrella Sampling calculations different results
   > > were obtained, and I was suggested to extend my simulations since the
   > > error bars associated to each windows were too high.
   > > I have now extended my simulations from 1 ns to 8 ns, however, I
cannot
   > > see much different from the shorter calculations. I send you the
   > > comparison of the two PMF including the error bars (attached).
   > > Now I am using 50 windows, but the shorter simulations were done using
   > > 100 windows, so I don't think increasing the number of windows
could help.
   > > My system has about 29200 atoms (where 29000 are chloroform atoms).
The
   > > *mdp file I am using is copied below.
   > > Would you have any suggestion to improve the results and decrease the
   > > error bars in the calculations?
   > >
   > > ----------------------------MDP file---------------------------
   > > title = Umbrella pulling simulation
   > > define =
   > > define =
   > > ; Run parameters
   > > integrator = md
   > > dt = 0.002
   > > tinit = 0
   > > nsteps = 500000 ; 1 ns
   > > nstcomm = 10
   > > ; Output parameters
   > > nstxout = 5000 ; every 10 ps
   > > nstvout = 5000
   > > nstfout = 5000
   > > nstxtcout = 5000 ; every 10 ps
   > > nstenergy = 5000
   > > ; Bond parameters
   > > constraint_algorithm = lincs
   > > constraints = all-bonds
   > > continuation = yes
   > > ; Single-range cutoff scheme
   > > nstlist = 5
   > > ns_type = grid
   > > rlist = 1.4
   > > rcoulomb = 1.4
   > > rvdw = 1.4
   > > ; PME electrostatics parameters
   > > coulombtype = PME
   > > fourierspacing = 0.12
   > > fourier_nx = 0
   > > fourier_ny = 0
   > > fourier_nz = 0
   > > pme_order = 4
   > > ewald_rtol = 1e-5
   > > optimize_fft = yes
   > > ; Berendsen temperature coupling is on in two groups
   > > Tcoupl = Nose-Hoover
   > > tc_grps = ACH CL3
   > > tau_t = 0.5 0.5
   > > ref_t = 300 300
   > > ; Pressure coupling is on
   > > Pcoupl = Parrinello-Rahman
   > > pcoupltype = isotropic
   > > tau_p = 1.0
   > > compressibility = 4.5e-5
   > > ref_p = 1.0
   > > ; Generate velocities is off
   > > gen_vel = no
   > > ; Periodic boundary conditions are on in all directions
   > > pbc = xyz
   > > ; Long-range dispersion correction
   > > DispCorr = EnerPres
   > > ; Pull code
   > > pull = umbrella
   > > pull_geometry = distance
   > > pull_dim = N N Y
   > > pull_start = yes
   > > pull_ngroups = 1
   > > pull_group0 = r_1
   > > pull_group1 = r_2
   > > pull_init1 = 0
   > > pull_rate1 = 0.0
   > > pull_k1 = 1000 ; kJ mol^-1 nm^-2
   > > pull_nstxout = 1000 ; every 2 ps
   > > pull_nstfout = 1000 ; every 2 ps
   > >
   > > -----------------------------------------------
   > >
   > >
   > > Thanks a lot in advance.
   > >
   > > Best wishes,
   > >
   > > Dr. Rebeca Garcia
   > > Santiago de Compostela University
   > > Spain
   > >
   > >
   > >
   > > -- gmx-users mailing list gmx-users at gromacs.org
   > > http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the
   > > archive at http://www.gromacs.org/Support/Mailing_Lists/Search before
   > > posting! Please don't post (un)subscribe requests to the list. Use the
   > > www interface or send it to gmx-users-request at gromacs.org. Can't post?
   > > Read http://www.gromacs.org/Support/Mailing_Lists
   > >
   >
   > --
   > ========================================
   >
   > Justin A. Lemkul
   > Ph.D. Candidate
   > ICTAS Doctoral Scholar
   > MILES-IGERT Trainee
   > Department of Biochemistry
   > Virginia Tech
   > Blacksburg, VA
   > jalemkul[at]vt.edu | (540) 231-9080
   > http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
   >
   > ========================================
   > --
   > gmx-users mailing list gmx-users at gromacs.org
   > http://lists.gromacs.org/mailman/listinfo/gmx-users
   > Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
   > Please don't post (un)subscribe requests to the list. Use the
   > www interface or send it to gmx-users-request at gromacs.org.
   > Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================
-- 
gmx-users mailing list    gmx-users at gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at  
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-request at gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists





More information about the gromacs.org_gmx-users mailing list