[gmx-users] R: Re: several questions about g_hbond -contact

[Anna Marabotti]

Dear Justin,
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...

Thanks
Anna

----------------------------------------------------------------------

Message: 1
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 
> 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

========================================