take aromatic rings into account for hydrogen bond calculation

src/mol-model-props/computed/interactions/hydrogen-bonds.ts

@@ -68,8 +68,11 @@ function addUnitHydrogenDonors(structure: Structure, unit: Unit.Atomic, builder:
     const { totalH } = getUnitValenceModel(structure, unit)
     const { elements } = unit
     const { x, y, z } = unit.model.atomicConformation
+    const { elementAromaticRingIndices } = unit.rings
 
     for (let i = 0 as StructureElement.UnitIndex, il = elements.length; i < il; ++i) {
+        if (elementAromaticRingIndices.has(i)) continue;
+
         const element = typeSymbol(unit, i)
         if ((
             // include both nitrogen atoms in histidine due to
@@ -99,7 +102,7 @@ function addUnitWeakHydrogenDonors(structure: Structure, unit: Unit.Atomic, buil
             (
                 bondToElementCount(structure, unit, i, Elements.N) > 0 ||
                 bondToElementCount(structure, unit, i, Elements.O) > 0 ||
-                inAromaticRingWithElectronNegativeElement(structure, unit, i)
+                inAromaticRingWithElectronNegativeElement(unit, i)
             )
         ) {
             builder.add(FeatureType.WeakHydrogenDonor, FeatureGroup.None, x[elements[i]], y[elements[i]], z[elements[i]], i)
@@ -107,27 +110,21 @@ function addUnitWeakHydrogenDonors(structure: Structure, unit: Unit.Atomic, buil
     }
 }
 
-function inAromaticRingWithElectronNegativeElement(structure: Structure, unit: Unit.Atomic, index: StructureElement.UnitIndex) {
-    return false // TODO
-    // if (!a.isAromatic()) return false
-
-    // const ringData = a.residueType.getRings()
-    // if (!ringData) return false
-
-    // let hasElement = false
-    // const rings = ringData.rings
-    // rings.forEach(ring => {
-    //     if (hasElement) return  // already found one
-    //     if (ring.some(idx => (a.index - a.residueAtomOffset) === idx)) {  // in ring
-    //         hasElement = ring.some(idx => {
-    //             const atomTypeId = a.residueType.atomTypeIdList[ idx ]
-    //             const number = a.atomMap.get(atomTypeId).number
-    //             return number === Elements.N || number === Elements.O
-    //         })
-    //     }
-    // })
-
-    // return hasElement
+function inAromaticRingWithElectronNegativeElement(unit: Unit.Atomic, index: StructureElement.UnitIndex) {
+    const { elementAromaticRingIndices, all } = unit.rings
+    const ringIndices = elementAromaticRingIndices.get(index)
+    if (ringIndices === undefined) return false
+
+    for (let i = 0, il = ringIndices.length; i < il; ++i) {
+        const ring = all[ringIndices[i]]
+        for (let j = 0, jl = ring.length; j < jl; ++j) {
+            const element = typeSymbol(unit, ring[j])
+            if (element === Elements.N || element === Elements.O) {
+                return true
+            }
+        }
+    }
+    return false
 }
 
 /**
@@ -137,12 +134,15 @@ function addUnitHydrogenAcceptors(structure: Structure, unit: Unit.Atomic, build
     const { charge, implicitH, idealGeometry } = getUnitValenceModel(structure, unit)
     const { elements } = unit
     const { x, y, z } = unit.model.atomicConformation
+    const { elementAromaticRingIndices } = unit.rings
 
     const add = (i: StructureElement.UnitIndex) => {
         builder.add(FeatureType.HydrogenAcceptor, FeatureGroup.None, x[elements[i]], y[elements[i]], z[elements[i]], i)
     }
 
     for (let i = 0 as StructureElement.UnitIndex, il = elements.length; i < il; ++i) {
+        if (elementAromaticRingIndices.has(i)) continue;
+
         const element = typeSymbol(unit, i)
         if (element === Elements.O) {
             // Basically assume all oxygen atoms are acceptors!

src/mol-model/structure/structure/unit/rings.ts

@@ -25,6 +25,7 @@ class UnitRings {
     private _byFingerprint?: ReadonlyMap<UnitRing.Fingerprint, ReadonlyArray<UnitRings.Index>>;
     private _index?: {
         readonly elementRingIndices: ReadonlyMap<StructureElement.UnitIndex, UnitRings.Index[]>,
+        readonly elementAromaticRingIndices: ReadonlyMap<StructureElement.UnitIndex, UnitRings.Index[]>,
         readonly ringComponentIndex: ReadonlyArray<UnitRings.ComponentIndex>,
         readonly ringComponents: ReadonlyArray<ReadonlyArray<UnitRings.Index>>
     };
@@ -32,7 +33,7 @@ class UnitRings {
 
     private get index() {
         if (this._index) return this._index;
-        this._index = createIndex(this.all);
+        this._index = createIndex(this.all, this.aromaticRings);
         return this._index;
     }
 
@@ -42,11 +43,15 @@ class UnitRings {
         return this._byFingerprint;
     }
 
-    /** Maps atom index inside a Unit to the smallest ring index (an atom can be part of more than one ring) */
+    /** Maps atom index inside a Unit to ring indices (an atom can be part of more than one ring) */
     get elementRingIndices() {
         return this.index.elementRingIndices;
     }
 
+    get elementAromaticRingIndices() {
+        return this.index.elementAromaticRingIndices;
+    }
+
     /** Maps UnitRings.Index to index to ringComponents */
     get ringComponentIndex() {
         return this.index.ringComponentIndex;

src/mol-model/structure/structure/unit/rings/compute.ts

@@ -253,8 +253,9 @@ function buildFinderprint(elements: string[], offset: number) {
 type RingIndex = import('../rings').UnitRings.Index
 type RingComponentIndex = import('../rings').UnitRings.ComponentIndex
 
-export function createIndex(rings: ArrayLike<SortedArray<StructureElement.UnitIndex>>) {
+export function createIndex(rings: ArrayLike<SortedArray<StructureElement.UnitIndex>>, aromaticRings: ReadonlyArray<RingIndex>) {
     const elementRingIndices: Map<StructureElement.UnitIndex, RingIndex[]> = new Map();
+    const elementAromaticRingIndices: Map<StructureElement.UnitIndex, RingIndex[]> = new Map();
 
     // for each ring atom, assign all rings that it is present in
     for (let rI = 0 as RingIndex, _rI = rings.length; rI < _rI; rI++) {
@@ -266,6 +267,17 @@ export function createIndex(rings: ArrayLike<SortedArray<StructureElement.UnitIn
         }
     }
 
+    // for each ring atom, assign all aromatic rings that it is present in
+    for (let aI = 0, _aI = aromaticRings.length; aI < _aI; aI++) {
+        const rI = aromaticRings[aI]
+        const r = rings[rI];
+        for (let i = 0, _i = r.length; i < _i; i++) {
+            const e = r[i];
+            if (elementAromaticRingIndices.has(e)) elementAromaticRingIndices.get(e)!.push(rI);
+            else elementAromaticRingIndices.set(e, [rI]);
+        }
+    }
+
     // create a graph where vertices are rings, edge if two rings share at least one atom
     const graph = new IntAdjacencyGraph.UniqueEdgeBuilder(rings.length);
     for (let rI = 0 as RingIndex, _rI = rings.length; rI < _rI; rI++) {
@@ -296,5 +308,5 @@ export function createIndex(rings: ArrayLike<SortedArray<StructureElement.UnitIn
         ringComponents[ringComponentIndex[rI]].push(rI);
     }
 
-    return { elementRingIndices, ringComponentIndex, ringComponents };
+    return { elementRingIndices, elementAromaticRingIndices, ringComponentIndex, ringComponents };
 }