[gmx-users] R: Re: several questions about g_hbond -contact
anna.marabotti at isa.cnr.it
Wed Feb 22 13:04:35 CET 2012
thakn you for your answer, although it is quite discouraging...
g_mindist is only partially suitable for my scope: it gives me the number of
contacts, and also the residues that contact the ligand at least once in the
simulation (if I use -or in addition to -od), but with g_hbond I obtained a
matrix of contacts and I could do another analysis using your plot_hbmap.pl
script. g_mindist does not produce the matrix, therefore I don't know the %
of time that a contact is present between protein and ligand (which was
rather useful for my scope).
Any other comment about this subject? (especially from developers) As I told
before, I'm not able to go into the code and see what -r2 does, but people
who wrote the code should be the ones who know this answer, I presume...
Date: Tue, 21 Feb 2012 15:39:48 -0500
From: "Justin A. Lemkul" <jalemkul at vt.edu>
Subject: Re: [gmx-users] several questions about g_hbond -contact
To: Discussion list for GROMACS users <gmx-users at gromacs.org>
Message-ID: <4F440114.9010607 at vt.edu>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Anna Marabotti wrote:
> Dear users,
> I'm using g_hbond with option -contact in order to find the contacts
> between my protein and a series of ligands, in a radius of 0.5 nm. I
> made several calculations but I'm quite uncertain about the results,
> therefore I'm asking you some questions.
> The command I used is:
> g_hbond -f my_prot.xtc -s my_prot.tpr -n index.ndx -g logfile.log -num
> contact_num.xvg -hbn contact.ndx -hbm contact.xpm -contact (for options
> -r and -r2 see below).
> Here my questions:
> 1) I have to indicate two groups when prompted. If I select first "1"
> (protein) and then "15" (ligand), the program does not find contacts (it
> calculates, but says "No contacts found"). On the contrary, if I first
> select "15" and then "1", the program finds contacts. Why? I'm expecting
> that if an atom of the protein contacts an atom of the ligand, it is a
> reversible thing, especially because this is not a H-bond in which the
> different order D-A could affect the results.
> To proceed with my test, I always selected first "15" and then "1".
> 2) If I fix only -r 0.5, the program starts calculating. It finds a
> number of "donor" and "acceptor" atoms and starts doing a grid search on
> a 18x18x13 grid, with rcut=0.5. It calculates for all the frames of the
> trajectory, then at the end it finds a number of "hbonds" and 0
> different atom-pairs within H-bond distance. At the end, I obtain a
> range checking error: "Variable y has value 0. It should have been
> within [ 0 .. 0 ] and the files .xvg, .log and .ndx are created, but the
> file .xpm is not.
> If I understand well, the -r2 flag is then absolutely necessary to tell
> the program that it does not have to treat this as an H-bond search, but
> as a contact search, right? The -contact flag is not sufficient, alone?
> 3) If I fix -r 0.5 -r2 0.6, the program starts calculating. It finds the
> same number as above of "donor" and "acceptor" atoms and starts doing a
> grid search on a 15x15x11 grid with rcut=0.6. It calculates for all the
> frames of the trajectory, then at the end it finds a number of
> "contacts" and a number of different atom-pairs within second cutoff
> distance. It produces all the files (.xvg, .ndx, .log and .xpm).
> In this case, the program has indeed calculated the "contacts", instead
> of the "hbonds". However, I don't understand if it calculates the
> contacts within a cutoff distance of 0.5 or within a cutoff distance of
> 0.6 nm. If I change this second value to, say, 1 nm, the number of
> contacts is still the same, but the number of different atom-pairs
> increases to 15951, and the grid has different dimensions (9x9x6)
> In this case, my question is: how the -r2 flag affects the number of
> contacts to find? (I would like to ask: what is the function of flag
> -r2? but looking at the gmx-users archive, it seems to me that nobody
> knows it...and I'm not able to look into the code, as suggested by
All of the above are valid points. These features have little or no
documentation (there is an open redmine issue about this). I find the
code hopelessly confusing, myself.
> My final question is, obviously: what is the correct command to provide
> to obtain what I want, i.e. the number of contacts between the protein
> and the ligand within a cutoff radius of 0.5 nm?
Don't use g_hbond at all.
g_mindist -d 0.5 -on
Justin A. Lemkul
ICTAS Doctoral Scholar
Department of Biochemistry
jalemkul[at]vt.edu | (540) 231-9080
More information about the gromacs.org_gmx-users