Path A - Module 38: Advanced H-Bond Algorithms#
In the previous module, we used a general criterion for Hydrogen Bonds. However, in high-end Biophysics, the definition of an H-bond depends on the specific physical model you want to follow.
MolSysMT implements classic algorithms from the literature (like Buch or Luzar-Chandler) and gives you total control over the rules of engagement.
import molsysmt as msm
from molsysmt import systems
# Load our solvated complex
fibril = msm.convert('pdb:2BEG', to_form='molsysmt.MolSys')
msm.build.add_missing_hydrogens(fibril)
peptide = msm.build.build_peptide('KLVFF')
msm.build.add_missing_hydrogens(peptide)
molsys = msm.merge([fibril, peptide])
1. The Luzar-Chandler Algorithm#
This is one of the most widely used criteria in molecular dynamics. It uses a distance cutoff and an angle cutoff to define the H-bond.
lc_hbonds = msm.hbonds.get_luzard_chandler_hbonds(molsys, selection='molecule_type=="peptide"',
selection_2='molecule_type=="protein"')
print(f"Luzar-Chandler H-bonds: {len(lc_hbonds)}")
2. Customizing Donors and Acceptors#
By default, MolSysMT knows which atoms can form H-bonds. But you can refine these lists using Inclusion and Exclusion rules.
# Get all potential acceptors but EXCLUDE those from the Phenylalanine sidechains
my_acceptors = msm.hbonds.get_acceptor_atoms(molsys, selection='molecule_type=="peptide"',
exclusion_rules=['group_name=="PHE"'])
print(f"Found {len(my_acceptors)} custom acceptor atoms in the peptide.")
3. The Buch Algorithm#
The Buch criterion is often used for specific hydration studies and has slightly different geometric requirements.
buch_hbonds = msm.hbonds.get_buch_hbonds(molsys, selection='molecule_type=="peptide"',
selection_2='molecule_type=="protein"')
print(f"Buch H-bonds: {len(buch_hbonds)}")
🏆 Path A Challenge: The Interaction Physicist#
Compare the number of H-bonds between the peptide and the fibril using Luzar-Chandler vs Buch.
Identify the atom indices of the 3 strongest H-bonds (those with the shortest distance).
Use
msm.hbonds.get_donor_atoms()to find if there are any Nitrogen atoms in your peptide that are NOT currently acting as donors.
You have reached the peak of interaction analysis! In Module 39, we will finish this phase by looking at the Physicochemical Properties of the complex: mass, charge, and surface area.