the easiest way is to use g_contacts. Unfortunately it only works with GROMACS 4.6.7. Another way is to use gmx hbond -contact -r2 0.35 to get the total number of contacts within each frame and then have indexes only for the residues of interest.
The plot you have shared is generated by "Simulation Interaction Diagram" module in Schrodinger Maestro. You can obtain the Free Maestro from Schrodinger or the Desmond Academic release from DEShaw Research Group. Import the trajectory (may need to convert to desmond trajectroy; VMD probably can do that) and run the analysis. It's GUI based, a few clicks will do the trick.