[gmx-users] Can g_wham support using different temperature for different windows?
Jianguo Li
ljggmx at yahoo.com.sg
Tue Feb 22 09:12:19 CET 2011
Sorry I forgot to attach my mdp files.
Here is the mdp file for pulling simulaition:
-------------------------------------------------------------------------
title = Martini
cpp = /usr/bin/cpp
define = -DPOSRES_LIP
integrator = md
; start time and timestep in ps
tinit = 0.0
dt = 0.02
nsteps = 5000000
; number of steps for center of mass motion removal =
nstcomm = 10
comm-grps =
nstxout = 5000
nstvout = 500000
nstfout = 0
; Output frequency for energies to log file and energy file =
nstlog = 1000
nstenergy = 1000
; Output frequency and precision for xtc file =
nstxtcout = 1000
xtc_precision = 100
nstlist = 10
; ns algorithm (simple or grid) =
ns_type = grid
; Periodic boundary conditions: xyz or none =
pbc = xyz
; nblist cut-off =
rlist = 1.4
coulombtype = PME
rcoulomb = 1.4
fourierspacing = 0.12
pme_order = 4
ewald_rtol = 1e-5
optimize_fft = yes
; Dielectric constant (DC) for cut-off or DC of reaction field =
epsilon_r = 15
; Method for doing Van der Waals =
vdw_type = Shift
; cut-off lengths =
rvdw_switch = 0.9
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure =
DispCorr = No
; Temperature coupling =
tcoupl = V-rescale
; Groups to couple separately =
tc-grps = Protein_lipid Sol_Ion
; Time constant (ps) and reference temperature (K) =
tau_t = 1.5 1.5
ref_t = 310 310
; Pressure coupling =
Pcoupl = Parrinello-Rahman
Pcoupltype = semiisotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar) =
tau_p = 10.0 10.0
compressibility = 3e-5 3e-5
ref_p = 1.0 1.0
constraints = none
; Type of constraint algorithm =
constraint_algorithm = Lincs
; Do not constrain the start configuration =
unconstrained_start = no
; Highest order in the expansion of the constraint coupling matrix =
lincs_order = 4
; Lincs will write a warning to the stderr if in one step a bond =
; rotates over more degrees than =
lincs_warnangle = 90
; pull staff
; pull staff
pull = umbrella
pull_geometry = position
pull_vec1 = 0 0 -1
pull_dim = N N Y
pull_start = no ; define initial COM distance > 0
pull_ngroups = 1
pull_group0 = lipid1
pull_group1 = Protein
pull_init1 = 0.0 0.0 4.50
pull_rate1 = 0.001 ; 0.01 nm per ps = 10 nm per ns
pull_k1 = 1000 ; kJ mol^-1 nm^-2
================================================================
Here is the pull part of the mpd file for the windowed umbrella sampling
simulation, other part of the mdp file are same as that of pulling simulation.
--------------------------------------------------------------------------------
pull = umbrella
pull_geometry = position
pull_vec1 = 0 0 -1
pull_dim = N N Y
pull_start = no ; define initial COM distance > 0
pull_ngroups = 1
pull_group0 = lipid1
pull_group1 = Protein
pull_init1 = 0.0 0.0 1.2
pull_k1 = 1000 ; kJ mol^-1 nm^-2
Cheers,
Jianguo
________________________________
From: Jianguo Li <ljggmx at yahoo.com.sg>
To: jalemkul at vt.edu; Discussion list for GROMACS users <gmx-users at gromacs.org>
Sent: Tuesday, 22 February 2011 14:27:34
Subject: Re: [gmx-users] Can g_wham support using different temperature for
different windows?
Thanks Justin and Chris and sorry for confusing interpretation.
Let me make it more clear. My peptide is flexible Martini beads, and highly
positively charged. My membrane is a mixture of negatively charged lipids (25%)
and zitterionic lipids(75%). So there is strong electrostatic attraction
between peptide and membrane. To get the PMF, I did the following:
(1) I did pulling simulation along (0 0 -1) direction to pull my peptide across
the membrane. Then I got different configurations corresponding to different
windows along the reaction coordinates, which is the z-distance
between peptide and membrane. This figure
(http://www.flickr.com/photos/lijg/5467080971/) shows some of the configurations
at certain reaction coordinates.
(2) In each window, I used the corresponding configuration that generated by the
pulling simulation as initial input and run umbrella sampling. The size of each
window is 0.15 nm, but close to the bilayer cneter (e.g., -0.6<d<0.6), I have
increased number of windows so that the width of the window is to be 0.05 or 0.1
nm, I also tried to use different force constant in these windows.
From the figure (http://www.flickr.com/photos/lijg/5467080971/) , we can
classify the peptide conformation to be either extended (interacting with two
bilayers) or compact (interacting with only one bilayer). Ideally, the peptide
conformation should be similar for d=x and d=-x. The problem is that the
configuration of peptide is not symmetric with respect to the bilayer center.
For example, the peptide configuration is compact at d=0.6 and d=0.9, but the
peptide is extended at d=-0.6 and d=-0.9. This leads Hysteresis. If I use g_wham
to generate PMF, then the PMF is not symmetric with respect to the bilayer
center. Using more number of windows and different force constant did not remove
the problem.
In my opinion, at least in some windows, the peptide should sample both compact
and extended structure. But what I found is that the windowed umbrella
simulation depends on the initial peptide conformation. If the initial peptide
conformation is compact, then after 100 ns, it is still compact; if the initial
peptide in that window is extended, the final configuration is also extended. I
also tried to run longer equilibrium time (e.g., 200 ns), but the problem still
exists.
My question is how to increase sampling of the peptide conformation? I just
think of two choices:
(1) use high temperature (e.g., 500K) at those bad windows. As I mentioned, I am
wondering if g_wham can unbias the effect of using different temperatures in
different windows.
(2) use REMD in those bad windows. These need a lot of computational resources.
Is there any other method to deal with the insufficient sampling?
Any suggestions are welcome, thanks for your time reading this email!
Cheers,
Jianguo
________________________________
From: Justin A. Lemkul <jalemkul at vt.edu>
To: Gromacs Users' List <gmx-users at gromacs.org>
Sent: Tuesday, 22 February 2011 11:13:05
Subject: Re: [gmx-users] Can g_wham support using different temperature for
different windows?
Jianguo Li wrote:
> Thanks for your comments, Justin.
>
> Using timestep of 20 fs, in each window the simulation runs for 100 ns CG time.
>The pulling rate is 0.001 nm/ps. Is it too fast?
>
Let me clarify things, since I'm not convinced I understand your procedure.
You generate a series of configurations with 0.001 nm/ps pulling, but then how
many windows do you generate for independent simulations? What are your .mdp
parameters during those windows? The pull rate should be 0 during the actual
umbrella sampling, to restrain the peptide within the window. What force
constant(s) do you use?
> My system is a little different. My peptide is highly positively charged. The
>NMR experiments show that the conformation of the peptide in water is very
>dynamic, so I make it flexible without fixing any secondary structure in Martini
>model.
As was discussed in the last few days, do not interpret changes in structure too
directly when using MARTINI. It is not designed to faithfully mimic secondary
structure changes.
> In the membrane, 25% of the lipids are negatively charged, so there are very
>strong electrostatic attraction between peptides and membrane.
>
> During the peptide approaching the membrane from the top, peptide can take
>different configurations at different reaction coordinates. When pulling the
>peptide into the membrane, the peptide takes relatively compact structure and
>interacts with only the top leaflet until the distance becomes smaller than 0.45
>nm, after that the peptide becomes extended structure and interacts with both
>leaflets. This extended structure remains until the distance becomes -1.05 nm.
>Further pulling leads to compact structure and interacts only with the lower
>leaflet. So the comformation of the peptide is not symmetric between the center
>of the bilayer, which leads to Hysteresis. It seems that there is a huge
>
I guess I'm confused here, too. The peptide is compact when interacting with
the top leaflet, extended further in the membrane, then compact again when
interacting with the lower leaflet. What's strange about that?
> energy barrier for the peptide to translocate across the membrane because if
>the initial conformation in a certain window is extended (interacting with both
>leaflets), then it remains extended. Similarly, it the initial conformation in a
>certain window is compact (interacting with only one leaflet), it will remain
>compact.
>
I don't see how that is necessarily unexpected or problematic. Peptides will
change conformation depending on their environment. If you want a static
structure to cross the membrane (which may or may not represent reality) you'll
have to introduce some kind of intramolecular restraints.
-Justin
> Any Suggestions of dealing with the highly charged system?
>
> Cheers,
> Jianguo
>
>
> ------------------------------------------------------------------------
> *From:* Justin A. Lemkul <jalemkul at vt.edu>
> *To:* Gromacs Users' List <gmx-users at gromacs.org>
> *Sent:* Tuesday, 22 February 2011 09:58:36
> *Subject:* Re: [gmx-users] Can g_wham support using different temperature for
>different windows?
>
>
>
> Jianguo Li wrote:
> > Thanks Justin.
> > I tried your suggestions by either increase more windows and change the
>force constant, but it seems the samplings are still bad in some windows. When I
>did pulling in (0 0 1) direction and a reverse pulling in (0 0 -1) direction, I
>got different configurations at certain reaction coordinates. And the windowed
>umbrella sampling seems depends strongly on the initial configurations in that
>window. Therefore I got different PMFs using pulling in (0 0 1) direction and
>reverse pulling in （0 0 -1) direction.
> >
>
> How long are each of the simulations in each window? Sufficient sampling
>should eliminate any configurational bias and/or hysteresis. Also, if the
>pulling that sets up the initial configurations is done slowly enough, you won't
>see these problems. Sounds to me like you're pulling too fast or hard, such
>that the system is not stable.
>
> > In my simulation, I exert constraints on phosphate atoms in z direction, so
>there is no lipid flip-flop and the membrane will be stable at high
>temperatures. Then I am thinking of increasing temperature in those bad windows
>to enhance sampling...
> >
>
> I don't know if I can make a convincing argument here, but intuitively, these
>windows would be sampling in a different ensemble, so the free energy landscape
>in these windows would be discontinuous with any adjacent windows that are done
>at different temperatures, and perhaps the forces required to restrain your
>peptide at a given COM distance will still result in a discontinuous PMF. I
>would also suspect that g_wham can't handle this situation; it has a -temp flag,
>but it only takes one value. So if you construct your PMF curve using WHAM, but
>supply incorrect or inconsistent information, I certainly wouldn't believe the
>result.
>
> I guess the main point is, there are tons of published demonstrations of
>peptides and other molecules crossing a membrane with SMD and umbrella sampling,
>so it should be possible to generate stable configurations without any funny
>tricks.
>
> -Justin
>
> > best regards,
> > Jianguo
> >
> >
> >
> > ------------------------------------------------------------------------
> > *From:* Justin A. Lemkul <jalemkul at vt.edu <mailto:jalemkul at vt.edu>>
> > *To:* Discussion list for GROMACS users <gmx-users at gromacs.org
><mailto:gmx-users at gromacs.org>>
> > *Sent:* Tuesday, 22 February 2011 09:35:37
> > *Subject:* Re: [gmx-users] Can g_wham support using different temperature
>for different windows?
> >
> >
> >
> > Jianguo Li wrote:
> > > Dear all,
> > >
> > > I want to get the PMF of my peptide across the membrane bilayer. First I
>pulled my peptide across the membrane and then did windowed umbrella sampling
>along the reaction coordinates which is the z-distance between peptide and
>membrane. However, I found that sampling is not sufficient in some windows(e.g.,
>around the center of the membrane). To enhance the sampling, I am thinking to
>run the simulation in those windows at higher temperature (e.g., 500K), but this
>will introduce a bias. My question is: can g_wham remove the bias due to using
>different temperatures in different windows?
> > >
> > > If g_wham cannot deal with the bias due to using different T, I may need
>to do REMD in those windows. But that will be very expensive computationally.
>Anybody have an idea of enhancing sampling in those windows?
> > >
> > > Btw, I am using Martini CG model.
> > >
> > > Any suggestions will be highly appreciated, thank you!
> > >
> >
> > A more straightforward approach is to (1) add more sampling windows or (2)
>increase the force constant in regions where there's poor sampling, or perhaps
>both.
> >
> > -Justin
> >
> > > Cheers,
> > > Jianguo
> > >
> >
> > -- ========================================
> >
> > 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
><mailto:gmx-users at gromacs.org> <mailto:gmx-users at gromacs.org
><mailto: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
><mailto:gmx-users-request at gromacs.org> <mailto:gmx-users-request at gromacs.org
><mailto: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
><mailto: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 <mailto: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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20110222/d4a2cd88/attachment.html>
More information about the gromacs.org_gmx-users
mailing list