[gmx-users] several questions about g_hbond -contact
Justin A. Lemkul
jalemkul at vt.edu
Tue Feb 21 21:39:48 CET 2012
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 g_hbond
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