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
somebody...)
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
--
========================================
Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
========================================
--
gmx-users mailing list [email protected]
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to [email protected].
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
