You'll need to create an index of triplets (Hydrogen-Donor-Acceptor) for all the catalytic residues and the ligand. Then you can use gmx hbond which which will use the angle and distance to compute the number of hydrogen bonds.

gmx hbond: http://manual.gromacs.org/programs/gmx-hbond.html

MDTraj also has several different ways of computing hydrogen bonds. For that you'd use the MDTraj Python module to load the trajectory. You'd also have to create an index of the catalytic residues and ligand, but I find indexing in MDTraj more intuitive.

Baker-Hubbard: http://mdtraj.org/latest/api/generated/mdtraj.baker_hubbard.html#mdtraj.baker_hubbard

Kabsch-Sander: http://mdtraj.org/latest/api/generated/mdtraj.kabsch_sander.html#mdtraj.kabsch_sander

Wernet-Nilsson: http://mdtraj.org/latest/api/generated/mdtraj.wernet_nilsson.html#mdtraj.wernet_nilsson

Thank you so much.. I will try to use MDtraj python module.