[gmx-users] g_hbond results inconsistent between v4.0.7 and v4.5.1
chris.neale at utoronto.ca
chris.neale at utoronto.ca
Fri Oct 15 16:55:46 CEST 2010
Has anybody checked g_hbond results for consistency between v4.0.7 and v4.5.1?
My usage is:
g_hbond -f my.xtc -r 0.68 -contact -life -s my.tpr -ac
and I find that the output files hbnum.xvg, hbac.xvg, hblife.xvg, and
the ac-reported lifetime of interactions are all different. To show
the simplest metric, here are the differences for hbnum.xvg:
paste hbnum.4.0.7.xvg hbnum.4.5.1.xvg|grep -v '[#|@]'|head
0 135 0 0 97 0
50 143 0 50 101 0
100 135 0 100 95 0
150 143 0 150 100 0
200 142 0 200 100 0
250 133 0 250 96 0
300 143 0 300 100 0
350 135 0 350 91 0
400 144 0 400 97 0
450 140 0 450 95 0
If I redo this, selecting only a single acceptor atom and 64 donor
atoms, with a large cutoff distance of 3 nm, v4.0.7 gives me 44
contacts in one frame while v4.5.1 gives me only 1 contact. g_dist
confirms that there are 44 contacts within 4 nm (using both versions
of g_dist).
Then I looked at the more standard actual h-bond usage for a single
frame and gave SOL SOL as the selection. The numbers of h-bonds vary
widely: 3915 vs. 7836. These numbers are nearly separated by 2x, but
not exactly 2x. I can find no record of any intentional changes.
gpc-f109n001-$ echo -e "SOL\nSOL\n" |
/project/pomes/cneale/GPC/exe/intel/gromacs-4.0.5/exec/bin/g_hbond -f
../../8ns/RESULTS/3.2/0/0/0-205ns.xtc -n index.ndx -s
make_nrexl10_tpr/nrexl10.gromacs4.0.5.tpr -e 0
:-) G R O M A C S (-:
Gromacs Runs One Microsecond At Cannonball Speeds
:-) VERSION 4.0.5 (-:
Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
Copyright (c) 1991-2000, University of Groningen, The Netherlands.
Copyright (c) 2001-2008, The GROMACS development team,
check out http://www.gromacs.org for more information.
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.
:-) /project/pomes/cneale/GPC/exe/intel/gromacs-4.0.5/exec/bin/g_hbond (-:
Option Filename Type Description
------------------------------------------------------------
-f ../../8ns/RESULTS/3.2/0/0/0-205ns.xtc Input Trajectory: xtc trr
trj gro g96 pdb cpt
-s make_nrexl10_tpr/nrexl10.gromacs4.0.5.tpr Input Run input file:
tpr tpb tpa
-n index.ndx Input, Opt! Index file
-num hbnum.xvg Output xvgr/xmgr file
-g hbond.log Output, Opt. Log file
-ac hbac.xvg Output, Opt. xvgr/xmgr file
-dist hbdist.xvg Output, Opt. xvgr/xmgr file
-ang hbang.xvg Output, Opt. xvgr/xmgr file
-hx hbhelix.xvg Output, Opt. xvgr/xmgr file
-hbn hbond.ndx Output, Opt. Index file
-hbm hbmap.xpm Output, Opt. X PixMap compatible matrix file
-don donor.xvg Output, Opt. xvgr/xmgr file
-dan danum.xvg Output, Opt. xvgr/xmgr file
-life hblife.xvg Output, Opt. xvgr/xmgr file
-nhbdist nhbdist.xvg Output, Opt. xvgr/xmgr file
Option Type Value Description
------------------------------------------------------
-[no]h bool no Print help info and quit
-nice int 19 Set the nicelevel
-b time 0 First frame (ps) to read from trajectory
-e time 0 Last frame (ps) to read from trajectory
-dt time 0 Only use frame when t MOD dt = first time (ps)
-[no]xvgr bool yes Add specific codes (legends etc.) in the output
xvg files for the xmgrace program
-[no]ins bool no Analyze solvent insertion
-a real 30 Cutoff angle (degrees, Acceptor - Donor -
Hydrogen)
-r real 0.35 Cutoff radius (nm, X - Acceptor, see next option)
-[no]da bool yes Use distance Donor-Acceptor (if TRUE) or
Hydrogen-Acceptor (FALSE)
-r2 real 0 Second cutoff radius. Mainly useful with -contact
and -ac
-abin real 1 Binwidth angle distribution (degrees)
-rbin real 0.005 Binwidth distance distribution (nm)
-[no]nitacc bool yes Regard nitrogen atoms as acceptors
-[no]contact bool no Do not look for hydrogen bonds, but merely for
contacts within the cut-off distance
-shell real -1 when > 0, only calculate hydrogen bonds within #
nm shell around one particle
-fitstart real 1 Time (ps) from which to start fitting the
correlation functions in order to obtain the
forward and backward rate constants for HB
breaking and formation
-temp real 298.15 Temperature (K) for computing the Gibbs energy
corresponding to HB breaking and reforming
-smooth real -1 If >= 0, the tail of the ACF will be smoothed by
fitting it to an exponential function: y = A
exp(-x/tau)
-dump int 0 Dump the first N hydrogen bond ACFs in a single
xvg file for debugging
-max_hb real 0 Theoretical maximum number of hydrogen bonds used
for normalizing HB autocorrelation function. Can
be useful in case the program estimates it wrongly
-[no]merge bool yes H-bonds between the same donor and acceptor, but
with different hydrogen are treated as a single
H-bond. Mainly important for the ACF.
-acflen int -1 Length of the ACF, default is half the number of
frames
-[no]normalize bool yes Normalize ACF
-P enum 0 Order of Legendre polynomial for ACF (0 indicates
none): 0, 1, 2 or 3
-fitfn enum none Fit function: none, exp, aexp, exp_exp, vac,
exp5, exp7 or exp9
-ncskip int 0 Skip N points in the output file of correlation
functions
-beginfit real 0 Time where to begin the exponential fit of the
correlation function
-endfit real -1 Time where to end the exponential fit of the
correlation function, -1 is till the end
No option -sel
Reading file make_nrexl10_tpr/nrexl10.gromacs4.0.5.tpr, VERSION 4.0.5
(single precision)
Specify 2 groups to analyze:
Group 0 ( System) has 17135 elements
Group 1 ( LEUJ) has 14 elements
Group 2 ( DOPC) has 3456 elements
Group 3 ( SOL) has 13665 elements
Group 4 ( P) has 64 elements
Group 5 ( N) has 64 elements
Group 6 ( SOL_&_OW) has 4555 elements
Group 7 ( SOL_&_!OW) has 9110 elements
Select a group: Selected 3: 'SOL'
Select a group: Selected 3: 'SOL'
Calculating hydrogen bonds in SOL (13665 atoms)
Found 4555 donors and 4555 acceptors
Reading frame 0 time 0.000
Will do grid-seach on 10x10x27 grid, rcut=0.35
Last frame 0 time 0.000
Back Off! I just backed up hbnum.xvg to ./#hbnum.xvg.10#
Average number of hbonds per timeframe 7836.000 out of 1.0374e+07 possible
gcq#337: "Act like Prometheus would" (Gogol Bordello)
gpc-f109n001-$ echo -e "SOL\nSOL\n" |
/project/pomes/cneale/GPC/exe/intel/gromacs-4.5.1/exec/bin/g_hbond -f
../../8ns/RESULTS/3.2/0/0/0-205ns.xtc -n index.ndx -s
make_nrexl10_tpr/nrexl10.gromacs4.5.1.tpr -e 0
:-) G R O M A C S (-:
Gravel Rubs Often Many Awfully Cauterized Sores
:-) VERSION 4.5.1 (-:
Written by Emile Apol, Rossen Apostolov, Herman J.C. Berendsen,
Aldert van Buuren, Pär Bjelkmar, Rudi van Drunen, Anton Feenstra,
Gerrit Groenhof, Peter Kasson, Per Larsson, Peiter Meulenhoff,
Teemu Murtola, Szilard Pall, Sander Pronk, Roland Schultz,
Michael Shirts, Alfons Sijbers, Peter Tieleman,
Berk Hess, David van der Spoel, and Erik Lindahl.
Copyright (c) 1991-2000, University of Groningen, The Netherlands.
Copyright (c) 2001-2010, The GROMACS development team at
Uppsala University & The Royal Institute of Technology, Sweden.
check out http://www.gromacs.org for more information.
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.
:-) /project/pomes/cneale/GPC/exe/intel/gromacs-4.5.1/exec/bin/g_hbond (-:
Option Filename Type Description
------------------------------------------------------------
-f ../../8ns/RESULTS/3.2/0/0/0-205ns.xtc Input Trajectory: xtc trr
trj gro g96 pdb cpt
-s make_nrexl10_tpr/nrexl10.gromacs4.5.1.tpr Input Run input file:
tpr tpb tpa
-n index.ndx Input, Opt! Index file
-num hbnum.xvg Output xvgr/xmgr file
-g hbond.log Output, Opt. Log file
-ac hbac.xvg Output, Opt. xvgr/xmgr file
-dist hbdist.xvg Output, Opt. xvgr/xmgr file
-ang hbang.xvg Output, Opt. xvgr/xmgr file
-hx hbhelix.xvg Output, Opt. xvgr/xmgr file
-hbn hbond.ndx Output, Opt. Index file
-hbm hbmap.xpm Output, Opt. X PixMap compatible matrix file
-don donor.xvg Output, Opt. xvgr/xmgr file
-dan danum.xvg Output, Opt. xvgr/xmgr file
-life hblife.xvg Output, Opt. xvgr/xmgr file
-nhbdist nhbdist.xvg Output, Opt. xvgr/xmgr file
Option Type Value Description
------------------------------------------------------
-[no]h bool no Print help info and quit
-[no]version bool no Print version info and quit
-nice int 19 Set the nicelevel
-b time 0 First frame (ps) to read from trajectory
-e time 0 Last frame (ps) to read from trajectory
-dt time 0 Only use frame when t MOD dt = first time (ps)
-xvg enum xmgrace xvg plot formatting: xmgrace, xmgr or none
-a real 30 Cutoff angle (degrees, Acceptor - Donor -
Hydrogen)
-r real 0.35 Cutoff radius (nm, X - Acceptor, see next option)
-[no]da bool yes Use distance Donor-Acceptor (if TRUE) or
Hydrogen-Acceptor (FALSE)
-r2 real 0 Second cutoff radius. Mainly useful with -contact
and -ac
-abin real 1 Binwidth angle distribution (degrees)
-rbin real 0.005 Binwidth distance distribution (nm)
-[no]nitacc bool yes Regard nitrogen atoms as acceptors
-[no]contact bool no Do not look for hydrogen bonds, but merely for
contacts within the cut-off distance
-shell real -1 when > 0, only calculate hydrogen bonds within #
nm shell around one particle
-fitstart real 1 Time (ps) from which to start fitting the
correlation functions in order to obtain the
forward and backward rate constants for HB
breaking and formation. With -gemfit we suggest
-fitstart 0
-fitstart real 1 Time (ps) to which to stop fitting the
correlation functions in order to obtain the
forward and backward rate constants for HB
breaking and formation (only with -gemfit)
-temp real 298.15 Temperature (K) for computing the Gibbs energy
corresponding to HB breaking and reforming
-smooth real -1 If >= 0, the tail of the ACF will be smoothed by
fitting it to an exponential function: y = A
exp(-x/tau)
-dump int 0 Dump the first N hydrogen bond ACFs in a single
xvg file for debugging
-max_hb real 0 Theoretical maximum number of hydrogen bonds used
for normalizing HB autocorrelation function. Can
be useful in case the program estimates it wrongly
-[no]merge bool yes H-bonds between the same donor and acceptor, but
with different hydrogen are treated as a single
H-bond. Mainly important for the ACF.
-geminate enum none Use reversible geminate recombination for the
kinetics/thermodynamics calclations. See
Markovitch et al., J. Chem. Phys 129, 084505
(2008) for details.: none, dd, ad, aa or a4
-diff real -1 Dffusion coefficient to use in the rev. gem.
recomb. kinetic model. If non-positive, then it
will be fitted to the ACF along with ka and kd.
-acflen int -1 Length of the ACF, default is half the number of
frames
-[no]normalize bool yes Normalize ACF
-P enum 0 Order of Legendre polynomial for ACF (0 indicates
none): 0, 1, 2 or 3
-fitfn enum none Fit function: none, exp, aexp, exp_exp, vac,
exp5, exp7 or exp9
-ncskip int 0 Skip N points in the output file of correlation
functions
-beginfit real 0 Time where to begin the exponential fit of the
correlation function
-endfit real -1 Time where to end the exponential fit of the
correlation function, -1 is until the end
Reading file make_nrexl10_tpr/nrexl10.gromacs4.5.1.tpr, VERSION 4.5.1
(single precision)
Specify 2 groups to analyze:
Group 0 ( System) has 17135 elements
Group 1 ( LEUJ) has 14 elements
Group 2 ( DOPC) has 3456 elements
Group 3 ( SOL) has 13665 elements
Group 4 ( P) has 64 elements
Group 5 ( N) has 64 elements
Group 6 ( SOL_&_OW) has 4555 elements
Group 7 ( SOL_&_!OW) has 9110 elements
Select a group: Selected 3: 'SOL'
Select a group: Selected 3: 'SOL'
Calculating hydrogen bonds in SOL (13665 atoms)
Found 4555 donors and 4555 acceptors
Reading frame 0 time 0.000
Will do grid-seach on 10x10x27 grid, rcut=0.35
Last frame 0 time 0.000
Back Off! I just backed up hbnum.xvg to ./#hbnum.xvg.11#
Average number of hbonds per timeframe 3915.000 out of 1.0374e+07 possible
gcq#317: "It's Unacceptable That Choclate Makes You Fat" (MI 3)
##########################
Finally, I checked to see that the default options had not changed,
and they have not:
gpc-f101n084-$ diff a b
5c5
< -s make_nrexl10_tpr/nrexl10.gromacs4.0.5.tpr Input Run input file:
---
> -s make_nrexl10_tpr/nrexl10.gromacs4.5.1.tpr Input Run input file:
23a24
> -[no]version bool no Print version info and quit
28,30c29
< -[no]xvgr bool yes Add specific codes (legends etc.) in the output
< xvg files for the xmgrace program
< -[no]ins bool no Analyze solvent insertion
---
> -xvg enum xmgrace xvg plot formatting: xmgrace, xmgr or none
48c47,52
< breaking and formation
---
> breaking and formation. With -gemfit we suggest
> -fitstart 0
> -fitstart real 1 Time (ps) to which to stop fitting the
> correlation functions in order to obtain the
> forward and backward rate constants for HB
> breaking and formation (only with -gemfit)
61a66,72
> -geminate enum none Use reversible geminate recombination for the
> kinetics/thermodynamics calclations. See
> Markovitch et al., J. Chem. Phys 129, 084505
> (2008) for details.: none, dd, ad, aa or a4
> -diff real -1 Dffusion coefficient to use in the rev. gem.
> recomb. kinetic model. If non-positive, then it
> will be fitted to the ACF along with ka and kd.
74c85
< correlation function, -1 is till the end
---
> correlation function, -1 is until the end
So perhaps v4.0.7 double counts the hydrogen bonds and 4.5.1 does not
(by maxing out at 1), but then v4.5.1 misses true double donor
situations (at least double donor according to the geometric
definitions)?
Thank you,
Chris.
More information about the gromacs.org_gmx-users
mailing list