i was trying to search hbonds intermolecular in packing crystals. I use Molecule.HBondCriterion() to define hbonds.

        but the results are different with that shown in Mercury.

        for example:

        firstly, I define the HBondCriterion()

        def HBondCriterion():
            criterion = Molecule.HBondCriterion()
            criterion.distance_range = (-5.0, 0.0)
            criterion.intermolecular = 'intermolecular'
            criterion.angle_tolerance = 120.0
            criterion.donor_types['sp C'] = True
            criterion.donor_types['sp2 C'] = True
            criterion.donor_types['sp3 C'] = True
            criterion.donor_types['aromatic C'] = True
            criterion.donor_types['carboncationic C'] = True
            criterion.acceptor_types['unclassified F'] = True
            return criterion

criterion = HBondCriterion()

i searched hbond in ABABOW

csd_reader = io.EntryReader('CSD')

cry = csd_reader.crystal('ABABOW')


the result is nothing; while I used cry.hbonds(hbond_criterion=criterion), the result is : 


I used Mercury to show this packing crystal, and the setting about hbond is the same with pythonAPI. It displayed some Hbonds (see attachment).

there are some thing different between pythonAPI and Mercury? why the result is nothing after packing?




the return value of the packing() function on a Crystal is an instance of a Molecule, not a Crystal, and the 'intermolecular' attribute of the HBondCriterion does not apply: all the HBonds are considered to be intramolecular, even if they are in disconnected components of the molecule.  Additionally, the HBonds will not be found since the path_length_range attribute if the HBondCriterion is left at its default setting of (4, 999), requiring the HBond to be within a single connected component of the molecule.

You can discover the HBonds by setting the path_length_range to (0, 999), allowing HBonds to be between disconnected components of the molecule, and you can determine which HBonds are between disconnected components by testing the shortest_path() method of the atoms of the HBond, that is to say

criterion.path_length_range = (0, 999)
criterion.intermolecular = 'ANY'

packed = cry.packing()
hbonds = packed.hbonds(hbond_criterion=criterion)
inter_component_hbonds = [hb for hb in hbonds if packed.shortest_path(hb.atoms[0], hb.atoms[-1]) == 0]

I hope this is helpful,  If anything is unclear please get back to us.

Best wishes



You must be signed in to post in this forum.