Gromacs stores matrix data in the xpm file format. This implementation of a Python reader is based on Tsjerk Wassenaar’s post to gmx-users numerical matrix from xpm file (Mon Oct 4 13:05:26 CEST 2010). This version returns a NumPy array and can guess an appropriate dtype for the array.
Class to make a Gromacs XPM matrix available as a NumPy numpy.ndarray.
The data is available in the attribute XPM.array.
Note
By default, the rows (2nd dimension) in the XPM.array are re-ordered so that row 0 (i.e. array[:,0] corresponds to the first residue/hydrogen bond/etc. The original xpm matrix is obtained for reverse = False. The XPM reader always reorders the XPM.yvalues (obtained from the xpm file) to match the order of the rows.
Initialize xpm structure.
Arguments : |
|
---|
Values of on the x-axis, extracted from the xpm file.
Values of on the y-axis, extracted from the xpm file. These are in the same order as the rows in the xpm matrix. If reverse = False then this is typically a descending list of numbers (highest to lowest residue number, index number, etc). For reverse = True it is resorted accordingly.
compiled regular expression to parse the colors in the xpm file:
static char *gromacs_xpm[] = {
"14327 9 2 1",
" c #FFFFFF " /* "None" */,
"o c #FF0000 " /* "Present" */,
Matches are named “symbol”, “color” (hex string), and “value”. “value” is typically autoconverted to appropriate values with gromacs.fileformats.convert.Autoconverter. The symbol is matched as a printable ASCII character in the range 0x20 (space) to 0x7E (~).
XPM matrix as a numpy.ndarray.
The attribute itself cannot be assigned a different array but the contents of the array can be modified.
Parse colour specification
Read and parse mdp file filename.
Return string s with C-style comments /* ... */ removed.
Return string s with quotes " removed.
Run gromacs.g_hbond() to produce the existence map (and the log file for the atoms involved in the bonds; the ndx file is also useful):
gromacs.g_hbond(s=TPR, f=XTC, g="hbond.log", hbm="hb.xpm", hbn="hb.ndx")
Load the XPM:
hb = XPM("hb.xpm", reverse=True)
Calculate the fraction of time that each H-bond existed:
hb_fraction = hb.array.mean(axis=0)
Get the descriptions of the bonds:
desc = [line.strip() for line in open("hbond.log") if not line.startswith('#')]
Note
It is important that reverse=True is set so that the rows in the xpm matrix are brought in the same order as the H-bond labels.
Show the results:
print "\n".join(["%-40s %4.1f%%" % p for p in zip(desc, 100*hb_fraction)])
See also
gromacs.analysis.plugins.hbonds