[gmx-developers] measuring parallel vs. anti-parallel hydrogen bonding between polymers

[Andrew Jewett]

> I wrote a simple program in python to help interpret hydrogen bond
> information from gromacs trajectories.  It generates a "connectivity matrix".
> There are many ways this matrix could be used, and I will send a couple
> of programs as examples.

I have attached a program (xpm2correlation) which attempts
to determine whether two polymers are bonded together in
a parallel or anti-parallel fashion.  (It could also work on
two different regions of the same polymer.)
If the two strands are parallel, this program returns a positive
number close to +1.0.  If they are anti-parallel, it returns a
negative number close to -1.0).  The formula used to compute
this number is included at the end of this message.  This can be useful
if the two polymers are not straight, but are long and flexible.  Again,
this code is just an example.  There may be better ways of calculating
this.

This code has been tested using the CVS version of g_hbond (version
3.3.1).  You will need to have python installed (preferably in /usr/bin)
to run the script.  Windows users can change the name of the file to
add a .py extension and run it in through python interpreter.

Syntax (example):

xpm2correlation hbond.ndx residues1.ndx residues2.ndx [DA,AD] < hbmap.xpm

   As an input this program requires a pair of files ("residues1.ndx" and
"residues2.ndx" in this example), containing a list of which
atoms belonging to each of the residues from either polymer.
(There should be one group per residue, and one file for either polymer.)
Input:

 hbmap.xpm An .xpm file created by "g_hbond -hbm" that you want to analyze
              (See documentation for g_hbond.)

Arguments:

 hbond.ndx    The ndx file created by "g_hbond -hbn ..."
 residues1.ndx  a list of the atoms in each residue (1 group per
residue) for polymer1
 residues2.ndx  a list of the atoms in each residue (1 group per
residue) for polymer2
 DA or AD     The 4th argument DA / AD is optional.

             By using the "DA" and "AD" arguments, the user can limit
             the atoms that the program considers from groups1.ndx
             and groups2.ndx.  The DA command line argument causes
             xpm2matrix to ignore all atoms from the first list of
             groups which are not potential donor atoms, and all atoms
             from the second list which are not potential acceptors.
             (The reverse is true for the AD command.)

     ---- formula ----

Let "K" and "L", denote the length of the two polymers, respectively.
And let's assume that "L" is the number of residues in the shorter of
the two.  Then the number calculated and returned by xpm2correlation is:

       =  A / B

where A = \sum_{over bonded pairs i,j}  (i - <i>)(j - <j>)

          B = \sum_{j=1...L} (j-L/2)*(j-L/2)

where i and j are integers ranging from 1...K and 1...L
indicate which two residues are hydrogen
bonded together, <i> and <j> denote the average position of a residue
participating in a hydrogen bond (from either molecule,
respectively).  "A" measures the correlation between the
locations of bonded residues from either polymer, and "B"
measures the same thing for two (maximally bonded)
parallel polymers of length L.

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

A side note:  There is some controversy over the order of h-bonds in
the "hbond.ndx"
file generated by g_hbond -hbn.  If the gromacs developers eventually decide
to reverse the order of the bonds in the "hbond.ndx" file generated by g_hbond,
then to make this program work all you have to do is uncomment lines 171 and
172 of this xpm2correlation code (the code included in this attachment.)
-------------- next part --------------
A non-text attachment was scrubbed...
Name: xpm2correlation
Type: application/octet-stream
Size: 13457 bytes
Desc: not available
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20070907/4bd0010b/attachment.obj>