I am attempting to count the number of hydrogen bonds around molecules in several hundred structures from the CSD. I would like to count these on a "per molecule" basis -- i.e. ignoring crystal symmetry equivalence, simple count the number of hydrogen bonds made by each chemically-independent species in a crystal or co-crystal.
Crystal.hbonds function seems to behave fairly sensibly, returning several hydrogen bonds for a structure like e.g. CATKIT, although it is not counting in the fashion that I would like (i.e. it at least partly accounts for crystal symmetry and tries to return only completely unique H-bonds). Given that I want to consider the immediate environment around each molecule regardless of the crystal symmetry, this seems like it would be the perfect use case for the
molecular_shell function in the
However, when I build the molecular shell for CATKIT and use the
Molecule.hbonds function on the resulting object (which I understand is itself a collection of
Molecule() objects for each "component"), I get no hydrogen bonds (i.e. the function returns an empty tuple). I have attempted to create both molecular and packing shells and use
Molecule.hbonds on the resulting objects to no avail (both methods return empty tuples).
Similarly, even calling
Molecule.hbonds on the original
Crystal.Molecule object (i.e. the "base" molecule of this entry in the CSD) returns an empty tuple, despite there clearly being at least one intermolecular hydrogen bond present in the asymmetric unit. This seems counter-intuitive to me -- surely even if a
Molecule object is itself composed of several further
Molecule.hbonds() function should return all contacts that match the hydrogen bond criteria in the
Molecule class instance it is called on?
Any help you could provide with getting the
Molecule.hbonds function to correctly count hydrogen bonds would be greatly appreciated. Alternatively, if you could suggest another means of counting all hydrogen bonds around distinct, independent molecules using the API, that would be greatly appreciated. (To give you an idea of what I'm looking for, in CATKIT, the bi-napthalene-diol species I would count as participating in 3 hydrogen bonds, while the other benzamidazole species participates in 2.)
Thanks in advance.