[gmx-users] Re: boundaries in PMF
Thomas Schlesier
schlesi at uni-mainz.de
Fri Jun 1 15:33:51 CEST 2012
One comment:
If the channel is horizontal orientated, and the COM of MOL lies in the
middle, you have two directions to go out of the channel: two the left
side (quasi negative distance) and to the right side (quasi positive
distance).
What happens with 'pull_geometry=distance' is, that only the distance
between MOL and Na matters. So you get the same (positive) distance if
you go to the left or go to the right. The resulting PMF will be for
only half of the channel (since the are no negative distances), but it
will be evaluated from going into both directions (since the distances
for left and right are the same).
To distingiush if you go left or right, you would need to use
'pull_geometry=direction/position'.
If you are 100% sure, that the system is symmetric and that the COM of
MOL is in the middle of the channel, then your results are right. But if
the COM of MOL is a little more to the right or left (relative to the
middle of the channel) you would get problems. Since then the left and
right direction are not equal.
One thing you could check is.
Make a analysis with g_wham on the windows which are on the left side of
the center and one for the right side of the center. If everything is
converged and the system is symmetric + COM of MOL is in the middle of
the channel, both PMF should look the same. If they differ one or more
of the above things are not correct.
Am 31.05.2012 22:03, schrieb gmx-users-request at gromacs.org:
> OK,
> however, just one point about your last comment:
>
>> > I suspect this is why g_wham is finding a range of values only equal to half of
>> > your expected reaction coordinate; it is considering only the positive
>> > displacement along the reaction coordinate.
> It seems like all the channel is explored, not only one half. If g_wham was only considering the positive displacement I should see a profile for just one half of the channel, shouldn?t I? and I can see a profile typical for the entire channel.
>
> Rebeca.
>
>> > Date: Thu, 31 May 2012 15:42:06 -0400
>> > From:jalemkul at vt.edu
>> > To:gmx-users at gromacs.org
>> > Subject: Re: [gmx-users] Re: boundaries in PMF
>> >
>> >
>> >
>> > On 5/31/12 3:37 PM, Rebeca Garc?a Fandi?o wrote:
>>> > > Hi,
>>> > > the center of mass of my channel is at the middle of the ion channel, and it is
>>> > > a symmetric system, so I suppose these results would be OK. Anyway, I will check
>>> > > the options you propose.
>> >
>> > If you are sampling regions above and below the channel/membrane, then the
>> > "distance" geometry is not appropriate; you need either "direction" or (perhaps
>> > the most flexible option) "position." There are a number of useful discussions
>> > on such topics in the list archive and in my umbrella sampling tutorial for you
>> > to consider. You can likely create a series of .tpr files from new .mdp files
>> > with correct options to run the analysis.
>> >
>> > I suspect this is why g_wham is finding a range of values only equal to half of
>> > your expected reaction coordinate; it is considering only the positive
>> > displacement along the reaction coordinate.
>> >
>> > -Justin
>> >
>>> > > Thanks a lot for all your comments!!
>>> > > Best wishes,
>>> > > Rebeca.
>>> > >
>>> > > > Date: Thu, 31 May 2012 20:08:26 +0200
>>> > > > From:schlesi at uni-mainz.de
>>> > > > To:gmx-users at gromacs.org
>>> > > > Subject: [gmx-users] Re: boundaries in PMF
>>> > > >
>>> > > > Where is the center of mass of reference group (MOL) located?
>>> > > >
>>> > > > It seems that the COM is near the middle of the ion channel. Since you
>>> > > > use 'pull_geometry=distance', g_wham will look only for the distance
>>> > > > between 'MOL' and 'Na' and that leads to problem.
>>> > > > If the com of 'MOL' sits in the center of the channel (which is around
>>> > > > 5nm long), one has 2.5nm in both directions. Since g_wham uses only the
>>> > > > distance you get the PMF for half of the channel, but with the data of
>>> > > > both parts.
>>> > > > If the channel would be symmetric and the com of 'MOL' would lie exactly
>>> > > > in the middle of the channel, this could be ok. But if even one of both
>>> > > > assumptions is wrong, the results would be errorous.
>>> > > >
>>> > > > A better approach would be to use 'pull_geometry=direction', since the
>>> > > > you define a vector along which the windows lie and do not have the
>>> > > > problem that the distance is in both directions (along positive and
>>> > > > negative vector) the same.
>>> > > > Only problem could be that g_wham supports 'pull_geometry=direction'
>>> > > > only from gromacs 4.5.x (don't know this, since instead of umbrella
>>> > > > smapling i use thermodynamik integration, where one uses constraints
>>> > > > (instead of restraints) and integrates the constraint-force; but the
>>> > > > conceptual problem with 'distance/direction' is the same).
>>> > > >
>>> > > > Another approach (with 'pull_geometry=distance') would be to use a
>>> > > > reference group which is just outside of the channel (like going some
>>> > > > steps away from the channel, along the vector which goes through the
>>> > > > channel). If there is only water, it would be bad, because then the
>>> > > > reference group would be move away.
>>> > > > But then on could use the entry and exit of the channel as a reference
>>> > > > point for two simulations. In the case that the entry is the reference
>>> > > > group, the PMF would be ill defined near the entry, but from the
>>> > > > simulation with the exit as reference you would get the right PMF for
>>> > > > the entry region and vice versa.
>>> > > >
>>> > > > Greetings
>>> > > > Thomas
>>> > > >
>>> > > >
>>> > > > Am 31.05.2012 19:39, schriebgmx-users-request at gromacs.org:
>>> > > > > Thanks a lot for your quick answer. The mdp file I have used is copied
>>> > > > > below. What is strange is that when I look at the *gro files for the
>>> > > > > different windows (100 windows in total), i. e: window 1: 8704Na Na56458
>>> > > > > 5.134 5.085 5.824 window 50: 8704Na Na56458 5.053 5.081 3.459 window
>>> > > > > 100: 8704Na Na56458 4.990 5.042 0.951 you can see that the z-coordinate
>>> > > > > goes from 0.951 to 5.824 nm I should have a total distance in the x-axis
>>> > > > > of about 5 nm, and however, the boundaries calculated by g_wham are
>>> > > > > "Determined boundaries to 0.000035 to 2.603290 " Can you see anything in
>>> > > > > the mdp file which is causing this behaviour...? Thanks again for your
>>> > > > > help! MDP FILE USED: title = Umbrella pulling simulation define =
>>> > > > > -DPOSRES_B 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 Tcoupl =
>>> > > > > V-rescale tau_t = 0.1 0.1 0.1 tc-grps = MOL dop WAT_Cl_Na ref_t = 300
>>> > > > > 300 300 ; Pressure coupling Pcoupl = Parrinello-Rahman Pcoupltype =
>>> > > > > Semiisotropic ; Time constant (ps), compressibility (1/bar) and
>>> > > > > reference P (bar) tau-p = 1.0 1.0 compressibility = 4.6E-5 4.6E-5 ref-p
>>> > > > > = 1.0 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 = MOL pull_group1 = Na 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
>>> > > > >> > Date: Thu, 31 May 2012 13:25:26 -0400
>>> > > > >> > From:jalemkul at vt.edu
>>> > > > >> > To:gmx-users at gromacs.org
>>> > > > >> > Subject: Re: [gmx-users] boundaries in PMF
>>> > > > >> >
>>> > > > >> >
>>> > > > >> >
>>> > > > >> > On 5/31/12 1:20 PM, Rebeca Garc?a Fandi?o wrote:
>>> > > > >>> > > Hi,
>>> > > > >>> > > I am trying to calculate a PMF for an ion along a channel. Everything
>>> > > went OK,
>>> > > > >>> > > but when I used g_wham I got a profile with strange dimensions for
>>> > > the x-axis.
>>> > > > >>> > > What are the boundaries g_wham is using for calculating the units of
>>> > > x-axis?
>>> > > > >>> > >
>>> > > > >> >
>>> > > > >> > The values are based on the restraint distances along the reaction
>>> > > coordinate.
>>> > > > >> >
>>> > > > >>> > > I have used:
>>> > > > >>> > > g_wham -it tpr-files.dat -if pullf-files.dat -o file_output.xvg -hist
>>> > > > >>> > > file_histo_output.xvg -unit kCal -cycl yes
>>> > > > >>> > >
>>> > > > >>> > > (version 4.0.7)
>>> > > > >>> > >
>>> > > > >>> > > and the units I got in the x-axis are determined by the boundaries:
>>> > > > >>> > >
>>> > > > >>> > > "Determined boundaries to 0.000035 to 2.603290"
>>> > > > >>> > >
>>> > > > >>> > > Which are these units? nm?
>>> > > > >>> > >
>>> > > > >> >
>>> > > > >> > Yes.
>>> > > > >> >
>>> > > > >>> > > The z coordinate for my ion explores at least 5 nm!!
>>> > > > >>> > >
>>> > > > >>> > > I am a bit confused about that.
>>> > > > >>> > >
>>> > > > >> >
>>> > > > >> > The exact output depends on how you set up the umbrella sampling (in the
>>> > > .mdp
>>> > > > >> > file). The range of values corresponds to whatever the distances are
>>> > > that are
>>> > > > >> > sampled in the various windows. Perhaps there is a sign issue here? Do you
>>> > > > >> > have some restraints that are at negative displacement and others at
>>> > > positive?
>>> > > > >> > Did you set up the .mdp files properly to account for this behavior?
>>> > > > >> >
>>> > > > >> > -Justin
More information about the gromacs.org_gmx-users
mailing list