I have encountered an interesting glitch, which is not necessarily a bug, but potentially could be improved on the MDAnalysis side. In AMBER, waters might have 3 bonds per molecule in the topology (including H-H bond). This can be easily tested with the following code:
import MDAnalysis as mda
from MDAnalysis.tests.datafiles import PRM_NCBOX, TRJ_NCBOX
u = mda.Universe(PRM_NCBOX, TRJ_NCBOX)
waters = u.select_atoms("water")
n_waters = waters.residues.n_residues
n_water_bonds = len(waters.bonds)
print(f"Number of water residues: {n_waters}")
print(f"Total number of bonds associated with water molecules: {n_water_bonds}")
print(f"Average number of bonds per water residue: {n_water_bonds / n_waters:.1f}")
# Number of water residues: 464
# Total number of bonds associated with water molecules: 1392
# Average number of bonds per water residue: 3.0The thing is that while this bond is in the topology, it is not a real chemical bond, and this will create some issues downstream, when, for example, trying to convert a water molecule to rdkit:
from rdkit import Chem
u.select_atoms("water").residues[0].atoms.convert_to("RDKIT")The output will be:
[16:41:30] Explicit valence for atom # 1 H, 2, is greater than permitted
---------------------------------------------------------------------------
AtomValenceException Traceback (most recent call last)
[/tmp/ipykernel_676/1174561067.py](https://localhost:8080/#) in <cell line: 0>()
1 from rdkit import Chem
----> 2 u.select_atoms("water").residues[0].atoms.convert_to("RDKIT")
4 frames
[/usr/local/lib/python3.12/dist-packages/MDAnalysis/converters/RDKitInferring.py](https://localhost:8080/#) in _standardize_patterns(self, mol, max_iter)
401 # fragment mol (reactions must have single reactant and product)
402 fragments = list(
--> 403 Chem.GetMolFrags(mol, asMols=True, sanitizeFrags=self.sanitize)
404 )
405 for reactant in fragments:
AtomValenceException: Explicit valence for atom # 1 H, 2, is greater than permitted
Which is exactly what you expect for a water molecule with 3 bonds instead of 2.
I suggest adding an optional post-parsing for amber loaders, which will ensure that no H-H bonds in waters are there in the generated topology. If the plan is approved, I can start working on corresponding PR
I have encountered an interesting glitch, which is not necessarily a bug, but potentially could be improved on the
MDAnalysisside. In AMBER, waters might have 3 bonds per molecule in the topology (includingH-Hbond). This can be easily tested with the following code:The thing is that while this bond is in the topology, it is not a real chemical bond, and this will create some issues downstream, when, for example, trying to convert a water molecule to
rdkit:The output will be:
Which is exactly what you expect for a water molecule with 3 bonds instead of 2.
I suggest adding an optional post-parsing for amber loaders, which will ensure that no H-H bonds in waters are there in the generated topology. If the plan is approved, I can start working on corresponding PR