Calculating the hydrogen bond occupancy is a frequent task in the analysis of molecular dynamics simulations. Therefore, it has been implemented in packages like VMD hbonds or MDAnalysis.analysis.hbonds. However, these implementations do not fully treat periodic boundary conditions, as hydrogen bonds crossing the simulation box boundary are not incorporated into the analysis.
The attached python script takes a PSF file and a DCD file (or any other combination of topology and trajectory that is supported by MDAnalysis and contains cell informations) and calculates the correct hydrogen bond occupancy. This works for arbitrary simulation cell geometries (not only orthorhombic ones).
python hbondpbc.py PSFFILE DCDFILE 'donor selector' 'intermediate selector' 'acceptor selector' cutoff angle
The list of supported selectors is given in the MDAnalysis documentation. Intermediate atoms are hydrogen atoms. The selection can be used to speed up the analysis. The script will print the current frame number as progress information followed by the atom names that form a hydrogen bond and the occupancy of this bond over the course of the trajectory in percent. The cutoff is the maximum distance of the donor and acceptor that is accepted for a hydrogen bond while angle is the maximum deviation from straight 180 degrees. Cutoff in angstroms, angle in degrees. Recommended values: 3 Angstrom and 20 degrees (vmd default).