[gmx-users] Transition difficulties: version 3.3.3 to 4.0.3 regarding pull

Discussion:

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

Steve Fiedler

2009-01-22 18:22:24 UTC

Dear all,

I have encountered an odd behavior with use of the "pull_geometry =
distance" option of the pull code, upon transitioning from Gromacs
version 3.3.3 to version 4.0.3. It appears to be related to the center
of mass distances of the two pull groups, which has an effect of
abruptly displacing the coordinates of the less massive group. A
diagnostic is a discrepancy between the distances between the pull
groups from the preprocessor output in version 4.0.3, and the distance
between the groups as calculated using the difference of the groups'
centers of mass from the g_traj utility. For example, using the
coordinates of a system previously equilibrated with the constraint
force approach from version 3.3.3, the grompp output from version 4.0.3 is:

Pull group natoms pbc atom distance at start reference at t=0
0 2672 1336
1 60 6818 2.673 0.400

Using g_traj (4.0.3 version), the difference of the distance between the
center of masses of the two groups is: 0.39911 nm versus the 2.673 value
from above.

This issue does not exist in previous versions of Gromacs including
version 3.3.3. In version 4.0.3, this behavior occurs for both
pull=umbrella and pull = constraint, on 32 and 64 bit architecture
systems, and in both single and double precision calculations. A test
of a two atom system determined that the pull_start option was not
appropriate. The pull options used in the mdp file are listed below, as
well as the contents of the ppa file which has worked previously.

Suggestions would be appreciated,

Thank you,

Steve Fiedler

Gromacs 4.0.0 mdp pull options
pull = constraint
pull_geometry = distance
pull_dim = N N Y
pull_start = no
pull_init1 = .4
pull_group0 = Ref
pull_group1 = Buk

Gromacs 3.3.3 ppa file:
runtype = constraint
group_1 = BUK
reference_group = Ref
reftype = com
constraint_direction = 0 0 1
constraint_distance1 = 0.400

Berk Hess

2009-01-22 21:55:16 UTC

Permalink

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

Date: Thu, 22 Jan 2009 13:22:24 -0500
From: fiedler at umich.edu
To: gmx-users at gromacs.org
Subject: [gmx-users] Transition difficulties: version 3.3.3 to 4.0.3 regarding pull_geometry=distance
Dear all,
I have encountered an odd behavior with use of the "pull_geometry =
distance" option of the pull code, upon transitioning from Gromacs
version 3.3.3 to version 4.0.3. It appears to be related to the center
of mass distances of the two pull groups, which has an effect of
abruptly displacing the coordinates of the less massive group. A
diagnostic is a discrepancy between the distances between the pull
groups from the preprocessor output in version 4.0.3, and the distance
between the groups as calculated using the difference of the groups'
centers of mass from the g_traj utility. For example, using the
coordinates of a system previously equilibrated with the constraint
Pull group natoms pbc atom distance at start reference at t=0
0 2672 1336
1 60 6818 2.673 0.400
Using g_traj (4.0.3 version), the difference of the distance between the
center of masses of the two groups is: 0.39911 nm versus the 2.673 value
from above.
This issue does not exist in previous versions of Gromacs including
version 3.3.3. In version 4.0.3, this behavior occurs for both
pull=umbrella and pull = constraint, on 32 and 64 bit architecture
systems, and in both single and double precision calculations. A test
of a two atom system determined that the pull_start option was not
appropriate. The pull options used in the mdp file are listed below, as
well as the contents of the ppa file which has worked previously.
Suggestions would be appreciated,
Thank you,
Steve Fiedler
Gromacs 4.0.0 mdp pull options
pull = constraint
pull_geometry = distance
pull_dim = N N Y
pull_start = no
pull_init1 = .4
pull_group0 = Ref
pull_group1 = Buk
runtype = constraint
group_1 = BUK
reference_group = Ref
reftype = com
constraint_direction = 0 0 1
constraint_distance1 = 0.400
_______________________________________________
gmx-users mailing list gmx-users at gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php

_________________________________________________________________
Express yourself instantly with MSN Messenger! Download today it's FREE!
http://messenger.msn.click-url.com/go/onm00200471ave/direct/01/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20090122/cd70d952/attachment.html>

Chris Neale

2009-01-22 22:12:29 UTC

Permalink

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

Steve Fiedler

2009-01-23 17:14:43 UTC

Permalink

Hi,

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

-----------------
conf.gro
Three atoms
3
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
-----------------
index.ndx
[ System ]
1 2 3
[ Ref ]
1 2
[ Pulled ]
3
-----------------
grompp.mdp
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.top
; 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 ]
AAAA 1

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

[ moleculetype ]
BBBB 1

[ 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

Post by Chris Neale
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

Date: Thu, 22 Jan 2009 13:22:24 -0500
From: fiedler at umich.edu
To: gmx-users at gromacs.org
Subject: [gmx-users] Transition difficulties: version 3.3.3 to

4.0.3 regarding pull_geometry=distance

Post by Steve Fiedler
Dear all,
I have encountered an odd behavior with use of the "pull_geometry

= > distance" option of the pull code, upon transitioning from
Gromacs > version 3.3.3 to version 4.0.3. It appears to be related
to the center > of mass distances of the two pull groups, which has
an effect of > abruptly displacing the coordinates of the less
massive group. A > diagnostic is a discrepancy between the distances
between the pull > groups from the preprocessor output in version
4.0.3, and the distance > between the groups as calculated using the
difference of the groups' > centers of mass from the g_traj utility.
For example, using the > coordinates of a system previously
equilibrated with the constraint > force approach from version 3.3.3,

Post by Steve Fiedler
Pull group natoms pbc atom distance at start reference at t=0

0 2672 1336 > 1 60 6818

2.673 0.400 > > Using g_traj (4.0.3 version), the
difference of the distance between the > center of masses of the two
groups is: 0.39911 nm versus the 2.673 value > from above. > > This
issue does not exist in previous versions of Gromacs including >
version 3.3.3. In version 4.0.3, this behavior occurs for both >
pull=umbrella and pull = constraint, on 32 and 64 bit architecture >
systems, and in both single and double precision calculations. A
test > of a two atom system determined that the pull_start option was
not > appropriate. The pull options used in the mdp file are listed
below, as > well as the contents of the ppa file which has worked
previously. > > Suggestions would be appreciated,

_______________________________________________
gmx-users mailing list gmx-users at gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php

--
Steve Fiedler, Ph.D.
Research Fellow
Department of Mechanical Engineering
The University of Michigan
2024 G.G. Brown
2350 Hayward St.
Ann Arbor, MI 48109-2125

Berk Hess

2009-01-24 17:58:15 UTC

Permalink

Hi,

I just found out that I introduced a bug in 4.0.3, which could cause
the pull code to crash or give wrong results when running single processor.
In parallel it is correct.
I committed fixed for 4.0.4.

If you run in parallel things should be correct,
or change in src/kernel/md.c:
if (ir->pull) {
dd_make_local_pull_groups(NULL,ir->pull,mdatoms);
to
if (ir->pull && PAR(cr)) {
dd_make_local_pull_groups(NULL,ir->pull,mdatoms);

Berk

Date: Fri, 23 Jan 2009 12:14:43 -0500
From: fiedler at umich.edu
To: gmx-users at gromacs.org
Subject: Re: [gmx-users] Transition difficulties: version 3.3.3 to 4.0.3 regarding pull_geometry=distance
Hi,
Thank you Berk and Chris for the suggestions.
To address the possibility that this issue is related to periodic
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
-----------------
conf.gro
Three atoms
3
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
-----------------
index.ndx
[ System ]
1 2 3
[ Ref ]
1 2
[ Pulled ]
3
-----------------
grompp.mdp
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.top
; 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 ]
AAAA 1
[ atoms ]
; nr type resnr residue atom cgnr charge mass
1 A 1 AAA A 1 0.000 1000.0000
[ moleculetype ]
BBBB 1
[ 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

Date: Thu, 22 Jan 2009 13:22:24 -0500
From: fiedler at umich.edu
To: gmx-users at gromacs.org
Subject: [gmx-users] Transition difficulties: version 3.3.3 to

4.0.3 regarding pull_geometry=distance

Post by Steve Fiedler
Dear all,
I have encountered an odd behavior with use of the "pull_geometry

Post by Steve Fiedler
Pull group natoms pbc atom distance at start reference at t=0

0 2672 1336 > 1 60 6818

--
Steve Fiedler, Ph.D.
Research Fellow
Department of Mechanical Engineering
The University of Michigan
2024 G.G. Brown
2350 Hayward St.
Ann Arbor, MI 48109-2125
_______________________________________________
gmx-users mailing list gmx-users at gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php

Steve Fiedler

2009-01-26 04:45:29 UTC

Permalink

Hi,

I rebuilt the code using both sets of Berk's corrections. On an
isolated bilayer system with no periodic boundaries in the x-y plane
prepared in the manner previously explained, the recompiled grompp
output yielded no improvement in the discrepancies of the position of a
buckyball permeant. By visual inspection, it appears that the 4.0.x
versions still attempt to distort the position of the permeant
significant upward from the bilayer center, i.e., the center of mass of
the reference group. Interestingly, the 3.3.x and 4.0.x version
identically treated a test bilayer system constructed with a single
dummy particle permeant.

I appreciate the help of Berk, and Chris' detailed suggestions with this
interesting problem.

Sincerely,

Steve Fiedler

Post by Berk Hess
Hi,
I just found out that I introduced a bug in 4.0.3, which could cause
the pull code to crash or give wrong results when running single processor.
In parallel it is correct.
I committed fixed for 4.0.4.
If you run in parallel things should be correct,
if (ir->pull) {
dd_make_local_pull_groups(NULL,ir->pull,mdatoms);
to
if (ir->pull && PAR(cr)) {
dd_make_local_pull_groups(NULL,ir->pull,mdatoms);
Berk

Date: Fri, 23 Jan 2009 12:14:43 -0500
From: fiedler at umich.edu
To: gmx-users at gromacs.org
Subject: Re: [gmx-users] Transition difficulties: version 3.3.3 to

4.0.3 regarding pull_geometry=distance

Hi,
Thank you Berk and Chris for the suggestions.
To address the possibility that this issue is related to periodic
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
-----------------
conf.gro
Three atoms
3
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
-----------------
index.ndx
[ System ]
1 2 3
[ Ref ]
1 2
[ Pulled ]
3
-----------------
grompp.mdp
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.top
; 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 ]
AAAA 1
[ atoms ]
; nr type resnr residue atom cgnr charge mass
1 A 1 AAA A 1 0.000 1000.0000
[ moleculetype ]
BBBB 1
[ 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

Post by Chris Neale
I just checked similar simulations of mine and Berk's suggestion
accounts for similar discrepancies that I notice on a quick

evaluation

Post by Chris Neale
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

Post by Chris Neale
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

Date: Thu, 22 Jan 2009 13:22:24 -0500
From: fiedler at umich.edu
To: gmx-users at gromacs.org
Subject: [gmx-users] Transition difficulties: version 3.3.3 to

4.0.3 regarding pull_geometry=distance

Post by Steve Fiedler
Dear all,
I have encountered an odd behavior with use of the

"pull_geometry

Post by Chris Neale

3.3.3,

Post by Chris Neale

Post by Steve Fiedler
Pull group natoms pbc atom distance at start reference at t=0

0 2672 1336 > 1 60 6818

was

Post by Chris Neale

not > appropriate. The pull options used in the mdp file are listed
below, as > well as the contents of the ppa file which has worked
previously. > > Suggestions would be appreciated,

_______________________________________________
gmx-users mailing list gmx-users at gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php

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/mailing_lists/users.php

------------------------------------------------------------------------
Express yourself instantly with MSN Messenger! MSN Messenger
<http://clk.atdmt.com/AVE/go/onm00200471ave/direct/01/>
------------------------------------------------------------------------
_______________________________________________
gmx-users mailing list gmx-users at gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php

--
Steve Fiedler, Ph.D.
Research Fellow
Department of Mechanical Engineering
The University of Michigan
2024 G.G. Brown
2350 Hayward St.
Ann Arbor, MI 48109-2125

Berk Hess

2009-01-26 08:59:09 UTC

Permalink

Hi,

Strange, I think I have really fixed all problems now.
All test systems work fine for me.

Can you file a bugzilla entry?
Please include all the files needed to run grompp
as well as the correct pull distance and the pull distance you obtained.

Berk

Date: Sun, 25 Jan 2009 23:45:29 -0500
From: fiedler at umich.edu
To: gmx-users at gromacs.org
Subject: Re: [gmx-users] Transition difficulties: version 3.3.3 to 4.0.3 regarding pull_geometry=distance
Hi,
I rebuilt the code using both sets of Berk's corrections. On an
isolated bilayer system with no periodic boundaries in the x-y plane
prepared in the manner previously explained, the recompiled grompp
output yielded no improvement in the discrepancies of the position of a
buckyball permeant. By visual inspection, it appears that the 4.0.x
versions still attempt to distort the position of the permeant
significant upward from the bilayer center, i.e., the center of mass of
the reference group. Interestingly, the 3.3.x and 4.0.x version
identically treated a test bilayer system constructed with a single
dummy particle permeant.
I appreciate the help of Berk, and Chris' detailed suggestions with this
interesting problem.
Sincerely,
Steve Fiedler

Date: Fri, 23 Jan 2009 12:14:43 -0500
From: fiedler at umich.edu
To: gmx-users at gromacs.org
Subject: Re: [gmx-users] Transition difficulties: version 3.3.3 to

4.0.3 regarding pull_geometry=distance

calculated by

below.

Post by Chris Neale
I just checked similar simulations of mine and Berk's suggestion
accounts for similar discrepancies that I notice on a quick

evaluation

could be

Post by Berk Hess

Date: Thu, 22 Jan 2009 13:22:24 -0500
From: fiedler at umich.edu
To: gmx-users at gromacs.org
Subject: [gmx-users] Transition difficulties: version 3.3.3 to

4.0.3 regarding pull_geometry=distance

Post by Steve Fiedler
Dear all,
I have encountered an odd behavior with use of the

"pull_geometry

Post by Chris Neale

Post by Berk Hess
= > distance" option of the pull code, upon transitioning from
Gromacs > version 3.3.3 to version 4.0.3. It appears to be related
to the center > of mass distances of the two pull groups, which has
an effect of > abruptly displacing the coordinates of the less
massive group. A > diagnostic is a discrepancy between the distances
between the pull > groups from the preprocessor output in version
4.0.3, and the distance > between the groups as calculated using the
difference of the groups' > centers of mass from the g_traj utility.
For example, using the > coordinates of a system previously
equilibrated with the constraint > force approach from version

3.3.3,

Post by Chris Neale

Post by Berk Hess

Post by Steve Fiedler
Pull group natoms pbc atom distance at start reference at t=0

0 2672 1336 > 1 60 6818

was

Post by Chris Neale

Post by Berk Hess
not > appropriate. The pull options used in the mdp file are listed
below, as > well as the contents of the ppa file which has worked
previously. > > Suggestions would be appreciated,

_______________________________________________
gmx-users mailing list gmx-users at gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php

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/mailing_lists/users.php

chris.neale

2009-01-23 19:03:40 UTC

Permalink

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
displacement

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.

Chris.

-- original message --

Hi,

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

-----------------
conf.gro
Three atoms
3
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
-----------------
index.ndx
[ System ]
1 2 3
[ Ref ]
1 2
[ Pulled ]
3
-----------------
grompp.mdp
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.top
; 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 ]
AAAA 1

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

[ moleculetype ]
BBBB 1

[ 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

Post by Chris Neale
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

Steve Fiedler

2009-01-26 22:46:26 UTC

Permalink

Hi Chris,

Your suggestions were helpful. Due to the nature of this problem, there
were limitations in the application of the sequence of steps that you
recommended. Specifically, the difference in the equilibrated structure
from version 3.3.x to the calculated reference center of mass from
version 4.0.x causes an unreasonable "jump" with constraint force
calculation. Resultant bad contacts prohibit generation of a pull
output file.

More generally though, I believe I understand your broader point: the
sensitive nature of the periodic box with respect to the pull code and
possible discrepancies in the pbc treatment by various utilities, may
cause confusion with the diagnostic process. The test bilayer system I
created however, centered distantly from the x and y box edges, may
effectively remove this set of concerns. For example, this
configuration is relatively unchanged upon "recentering", however it is
still susceptible to the problem that is the subject of this discussion.

If I can reproduce this phenomenon with a publicly available topology
file, I can provide the necessary input files for inspection.

Thank you again,

Steve Fiedler

Post by chris.neale
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
displacement
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.
Chris.
-- original message --
Hi,
Thank you Berk and Chris for the suggestions.
To address the possibility that this issue is related to periodic
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
-----------------
conf.gro
Three atoms
3
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
-----------------
index.ndx
[ System ]
1 2 3
[ Ref ]
1 2
[ Pulled ]
3
-----------------
grompp.mdp
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.top
; 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 ]
AAAA 1
[ atoms ]
; nr type resnr residue atom cgnr charge mass
1 A 1 AAA A 1 0.000 1000.0000
[ moleculetype ]
BBBB 1
[ 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

Post by Chris Neale
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

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

--
Steve Fiedler, Ph.D.
Research Fellow
Department of Mechanical Engineering
The University of Michigan
2024 G.G. Brown
2350 Hayward St.
Ann Arbor, MI 48109-2125

Berk Hess

2009-01-27 08:37:26 UTC

Permalink

Hi,

Just to be sure, the pull code is producing the correct results now?
The only problems should be in several analysis tools.

The analysis tools have been written only with analysis of molecules
or parts of molecules in mind. When you want to analyze groups
that cover two molecules or more there will be pbc problems.
These are not simple to fix. If the distances between all the atoms
in the group are less than half a box length there is a unique COM.
The procedure should be similar to what trjconv -pbc cluster does.

I can try to implement this for g_traj and g_dist.
Are there other tools that cause problems?

Berk

Date: Mon, 26 Jan 2009 17:46:26 -0500
From: fiedler at umich.edu
To: gmx-users at gromacs.org
Subject: Re: [gmx-users] Transition difficulties: version 3.3.3 to 4.0.3 regarding pull_geometry=distance
Hi Chris,
Your suggestions were helpful. Due to the nature of this problem, there
were limitations in the application of the sequence of steps that you
recommended. Specifically, the difference in the equilibrated structure
from version 3.3.x to the calculated reference center of mass from
version 4.0.x causes an unreasonable "jump" with constraint force
calculation. Resultant bad contacts prohibit generation of a pull
output file.
More generally though, I believe I understand your broader point: the
sensitive nature of the periodic box with respect to the pull code and
possible discrepancies in the pbc treatment by various utilities, may
cause confusion with the diagnostic process. The test bilayer system I
created however, centered distantly from the x and y box edges, may
effectively remove this set of concerns. For example, this
configuration is relatively unchanged upon "recentering", however it is
still susceptible to the problem that is the subject of this discussion.
If I can reproduce this phenomenon with a publicly available topology
file, I can provide the necessary input files for inspection.
Thank you again,
Steve Fiedler

Post by Chris Neale
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

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

_________________________________________________________________
What can you do with the new Windows Live? Find out
http://www.microsoft.com/windows/windowslive/default.aspx
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20090127/48a1cf26/attachment.html>

Ran Friedman

2009-01-27 08:58:32 UTC

Permalink

Dear Berk,

I think that g_clustsize has a similar problem. IIRC I fixed it in a
similar way on my own copy.
There are probably other tools that will suffer from this as well, but
for g_clustsize it's important because it may deal with such groups.

Same goes for trjconv, e.g., for a protein and ligand.

Ran.

Post by Berk Hess
Hi,
Just to be sure, the pull code is producing the correct results now?
The only problems should be in several analysis tools.
The analysis tools have been written only with analysis of molecules
or parts of molecules in mind. When you want to analyze groups
that cover two molecules or more there will be pbc problems.
These are not simple to fix. If the distances between all the atoms
in the group are less than half a box length there is a unique COM.
The procedure should be similar to what trjconv -pbc cluster does.
I can try to implement this for g_traj and g_dist.
Are there other tools that cause problems?
Berk

Date: Mon, 26 Jan 2009 17:46:26 -0500
From: fiedler at umich.edu
To: gmx-users at gromacs.org
Subject: Re: [gmx-users] Transition difficulties: version 3.3.3 to

4.0.3 regarding pull_geometry=distance

Hi Chris,
Your suggestions were helpful. Due to the nature of this problem, there
were limitations in the application of the sequence of steps that you
recommended. Specifically, the difference in the equilibrated structure
from version 3.3.x to the calculated reference center of mass from
version 4.0.x causes an unreasonable "jump" with constraint force
calculation. Resultant bad contacts prohibit generation of a pull
output file.
More generally though, I believe I understand your broader point: the
sensitive nature of the periodic box with respect to the pull code and
possible discrepancies in the pbc treatment by various utilities, may
cause confusion with the diagnostic process. The test bilayer system I
created however, centered distantly from the x and y box edges, may
effectively remove this set of concerns. For example, this
configuration is relatively unchanged upon "recentering", however it is
still susceptible to the problem that is the subject of this discussion.
If I can reproduce this phenomenon with a publicly available topology
file, I can provide the necessary input files for inspection.
Thank you again,
Steve Fiedler

and

Post by chris.neale
g_dist on final_center.gro. In my case, I found that g_traj and

g_dist

Post by chris.neale
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

Post by chris.neale
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.
Chris.
-- original message --
Hi,
Thank you Berk and Chris for the suggestions.
To address the possibility that this issue is related to periodic
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

Post by chris.neale
with all lipid molecules inside a box and intact. The discrepancy
between the grompp (version 4.0.3) output and distances as

calculated by

Post by chris.neale
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.

Post by chris.neale
If there are additional suggestions, I would be greatly appreciative.
Thank you,
Steve Fiedler
-----------------
conf.gro
Three atoms
3
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
-----------------
index.ndx
[ System ]
1 2 3
[ Ref ]
1 2
[ Pulled ]
3
-----------------
grompp.mdp
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.top
; 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 ]
AAAA 1
[ atoms ]
; nr type resnr residue atom cgnr charge mass
1 A 1 AAA A 1 0.000 1000.0000
[ moleculetype ]
BBBB 1
[ 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

same