Get hydrogen-bond sites#

Recognizing attributed donor–hydrogen pairs and acceptor atoms.

Prepare chemical candidates independently of coordinates. Complete declared connectivity and element symbols are required. Candidate membership depends on the chosen method; it does not establish a hydrogen bond or favorable energy.

Added in version 1.0.0.

Basic usage#

Let’s show how this function works with a declared urea/water chemical graph. RDKit is used here to construct the graph with indexed hydrogen atoms; the public tool also accepts native Topology, MolSys and chemical H5MSM forms with the required information.

import molsysmt as msm
from rdkit import Chem
molsys = Chem.AddHs(Chem.MolFromSmiles('NC(=O)N.O.S.F'))
sites = msm.physchem.get_hbond_sites(molsys)
print(sites['donor_hydrogen_pairs'])
print(sites['acceptor_atom_indices'], sites['chemical_state_index'])
print(sites['method_reference'])
[[ 0  7]
 [ 0  8]
 [ 3  9]
 [ 3 10]
 [ 4 11]
 [ 4 12]]
[0 2 3 4] 0
{'software': 'MDTraj', 'commit': '80f7cf2ddb43dd0d32905e490464f2d8765329df', 'pi_pi': 'https://github.com/mdtraj/mdtraj/blob/80f7cf2ddb43dd0d32905e490464f2d8765329df/mdtraj/geometry/pi_stacking.py', 'hbonds': 'https://github.com/mdtraj/mdtraj/blob/80f7cf2ddb43dd0d32905e490464f2d8765329df/mdtraj/geometry/hbond.py', 'paper': 'https://doi.org/10.1016/j.bpj.2015.08.015', 'hydrogen_bond_definitions': {'baker_hubbard': 'https://doi.org/10.1016/0079-6107(84)90007-5', 'wernet_nilsson': 'https://doi.org/10.1126/science.1096205'}}

Choosing chemical definitions#

Method

Donor element bonded to indexed H

Acceptor rule

elemental_nitrogen_oxygen

N/O

Every N/O

elemental_fluorine_oxygen_nitrogen

F/O/N

Every F/O/N

smarts_donor_acceptor

ProLIF 2.2.2 donor SMARTS

ProLIF 2.2.2 acceptor SMARTS

The elemental rules include amide N as an acceptor; ProLIF excludes it. These are reproduced definitions, not interchangeable assignments of universal chemistry. No partial-charge guesses or residue/water filters are added. The following loop demonstrates the retained software-name aliases.

profiles = {method: msm.physchem.get_hbond_sites(molsys, method=method)
            for method in ('mdtraj', 'cpptraj', 'prolif')}
for method, result in profiles.items():
    print(method, result['donor_hydrogen_pairs'].tolist(), result['acceptor_atom_indices'].tolist())
mdtraj [[0, 7], [0, 8], [3, 9], [3, 10], [4, 11], [4, 12]] [0, 2, 3, 4]
cpptraj [[0, 7], [0, 8], [3, 9], [3, 10], [4, 11], [4, 12], [6, 15]] [0, 2, 3, 4, 6]
prolif [[0, 7], [0, 8], [3, 9], [3, 10], [4, 11], [4, 12], [5, 13], [5, 14]] [2, 4]

Selecting indexed participants#

Both donor and hydrogen must be selected for a pair to survive. Indices refer to the input graph; they are not atom IDs. Recognition uses the full graph before filtering.

donor, hydrogen = sites['donor_hydrogen_pairs'][0]
selected = msm.physchem.get_hbond_sites(molsys, selection=[int(donor), int(hydrogen)])
print(selected['donor_hydrogen_pairs'])
print(selected['selected_atom_indices'])
[[0 7]]
[0 7]

Distinguishing hydrogen annotations#

Bracket-declared or implicit hydrogen counts support chemical interpretation but do not provide atom indices or coordinates. This tool never constructs hydrogen atoms.

implicit = Chem.MolFromSmiles('O')
print(msm.physchem.get_hbond_sites(implicit)['donor_hydrogen_pairs'].shape)
assert msm.physchem.get_hbond_sites(implicit)['donor_hydrogen_pairs'].shape == (0, 2)
(0, 2)

Selecting stored chemical states#

Use the reference state or an explicit state index. State assignment by structure must resolve to one known state across the request. Native conversion preserves site indices. ProLIF also needs formal charges, bond orders and explicit atom/bond aromaticity.

native = msm.convert(molsys, to_form='molsysmt.Topology')
other = msm.physchem.get_hbond_sites(native, chemical_state=0, method='smarts_donor_acceptor')
print(other['chemical_state_index'], other['software'])
assert other['donor_hydrogen_pairs'].tolist() == profiles['prolif']['donor_hydrogen_pairs'].tolist()
0 {'molsysmt': '0.21.0+606.ga03eb4bf6', 'rdkit': '2025.09.5'}

Making limits explicit#

Unknown or partial connectivity fails unless you explicitly set assume_complete_connectivity=True. The assumption is recorded and does not repair the source. max_matches limits each ProLIF SMARTS query before selection; exceeding it raises rather than silently truncating the site set. Non-ProLIF profiles do not require RDKit at calculation time. Returned integer indices carry no physical units.

Apply geometry with Get attributed hydrogen bonds.

Preserving the reference#

The returned dictionary includes the canonical method, its versioned definition, producer versions and attribution bibliography. A reference implementation is distinct from an executed dependency. With Ackredit installed, completed recognition also contributes to the enclosing workflow; without it these metadata are identical. See Methods and attribution.

print(sites['method'], sites['method_definition'])
print([(item['title'], item['roles']) for item in sites['attribution']['items']])
elemental_nitrogen_oxygen molsysmt.hbond_sites.elemental_nitrogen_oxygen.default@1
[('MDTraj: A Modern Open Library for the Analysis of Molecular Dynamics Trajectories', ['reference_implementation']), ('MolSysMT', ['executed_software'])]