[gmx-users] Transition difficulties: version 3.3.3 to 4.0.3 regarding pull_geometry=distance

chris.neale at utoronto.ca chris.neale at utoronto.ca
Fri Jan 23 20:03:40 CET 2009

Hi Steve,

what I intended to suggest was actually something different (and much easier).

The idea is not that you need some special system to be able to  
utilize the pull code, but that the pull code is correct whereas the  
g_dist and g_traj programs are not as good at treating pbc in the way  
that one desires.

I suggest the following.

1. Take your original system and run the pull code for a very short  
simulation. Use the last line of the output to calculate the relevant  

2. Now use trjconv -b -e to get the last frame of the .xtc that  
resulted from that short MD run as a .gro file, call it final.gro. I  
suspect that your groups are not entirely in the same simulation box  
in final.gro.

3. Now make a new .ndx file from that .gro and give it a single  
residue that is near your binding pocket, call it R_1

4. Now apply trjconv -center -pbc mol -ur compact while selecting R_1  
for centering, call the new .gro file final_center.gro

5. Visualize final_center.gro and ensure that all of your relevant  
atoms are in the same image in the way that puts the minimum distance  
between them along a path that is entirely contained within the unit  
cell. If not, go back to step 3 and try making a group R_2, etc.   
until this process works. NOTE: you might think that giving trjconv  
-center the relevant groups that you use for pulling will be a good  
idea here, but it is not. The problem there is that the atoms may be  
"centered" by placing half on the left boundary and half on the right  
boundary. I find using one logically selected residue or atom is the  
best method here.

6. Assuming that you got what you wanted in step 5, now run g_traj and  
g_dist on final_center.gro. In my case, I found that g_traj and g_dist  
give the same answer as the pull code output when I am using  
final_center.gro, but not always when I am using final.gro.

*** I always laugh when these problems arise because, in an important  
sense, the protein *did* jump out of the simulation box... at least as  
far as g_traj and g_dist are concerned. This, we must hope, is  
correctly treated in the pull code even though it is incorrectly (or  
at least unintuitively) treated by g_traj and g_dist.


-- original message --


Thank you Berk and Chris for the suggestions.

To address the possibility that this issue is related to periodic
boundaries, I used two approaches:
1.  The pull group of interest (permeant) was centered in the x-y plane
of the box using Chris' approach.  I then used the genconf utility to
replicate my lipid box to a 9x9 grid in the x-y plane and removed all
but the center box.  This generated the coordinates for a bilayer system
with all lipid molecules inside a box and intact.  The discrepancy
between the grompp (version 4.0.3) output and distances as calculated by
g_traj (version 4.0.3) persist, 2.667 vs. 0.3996 nm.
2.  I constructed a three atom system containing 2 reference atoms of
type A, and a "pull" atom of type B.  Proper output from grompp was
observed for all coordinates of both the reference and pulled atoms,
include coordinates for atoms moved outside the box in the x-y plane.
The coordinate, topology, and run control parameter file are given below.

If there are additional suggestions, I would be greatly appreciative.

Thank you,

    Steve Fiedler

Three atoms
     1AAA      A    1   1.500   1.500   1.000
     2AAA      A    2   0.500   1.500   1.000
     3BBB      B    3  -1.500   1.500   1.700
    3.00000   3.00000   3.00000
[ System ]
    1    2    3
[ Ref ]
    1    2
[ Pulled ]
title                    = ThreeAtoms
integrator               = md
dt                       = 0.001
nsteps                   = 1
ns_type                  = grid
pbc                      = xyz
coulombtype              = shift
rlist                    = 1.4
rcoulomb                 = 1.4
rvdw                     = 1.4
tcoupl                   = no
pcoupl                   = no
constraint_algorithm     = shake
shake_tol                = 1e-4
gen-vel                  = no
gen-temp                 = 0

nstxout                  = 1
nstvout                  = 0
nstfout                  = 0

pull = umbrella
pull_geometry = distance
pull_dim = N N Y
pull_start = no
pull_init1 = 0.7
pull_group0 = Ref
pull_group1 = Pulled
pull_k1 = 10000
; topology for two partially charged atoms

[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1     3     yes      0.125  0.5

[ atomtypes ]
;name     mass     charge ptype  sig           eps
    A   1000.0000    0.000  A     0.50000      9.90000
    B      9.0000    0.000  A     0.30000      9.00000

[ nonbond_params ]
   ; i    j    func    sig          eps

[ moleculetype ]

[ atoms ]
;   nr   type  resnr residue  atom   cgnr  charge       mass
     1     A     1      AAA    A      1      0.000   1000.0000

[ moleculetype ]

[ atoms ]
;   nr   type  resnr residue  atom   cgnr  charge     mass
     1     B     1      BBB    B      1      0.000      9.00

[ system ]
; name
Three atoms

[ molecules ]
; name   number
AAAA      2
BBBB      1

Chris Neale wrote:
> I just checked similar simulations of mine and Berk's suggestion  
> accounts for similar discrepancies that I notice on a quick  
> evaluation where g_traj and g_dist fail to give me the same distance  
> as I obtain from the pull pos.xvg file. As Berk suggests, once I  
> first trjconv -center -pbc mol -ur compact (giving an appropriate  
> residue for centering that puts all relevant pulled atoms in the  
> same box) then g_traj and g_dist both give me the exact same answer  
> as I calculate based on pull pos.xvg. Chris -- original message --  
> Hi, There could be a problem with periodic boundary conditions. Do  
> you have multiple molecules in a pull group, or broken molecules? In  
> that case the COM position of 3.3.3 and g_traj are both incorrect.  
> The pull code in 4.0 grompp and mdrun are (as far as I know) always  
> correct. Berk

More information about the gromacs.org_gmx-users mailing list