[gmx-users] Re: boundaries in PMF

[Thomas Schlesier]

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, schrieb gmx-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