wip, charged interactions

src/mol-math/linear-algebra/3d/vec3.ts

@@ -505,6 +505,17 @@ namespace Vec3 {
         return add(out, scale(out, copy(out, vector), scalar), origin);
     }
 
+    export function projectOnVector(out: Vec3, p: Vec3, vector: Vec3 ) {
+        const scalar = dot(vector, p) / squaredMagnitude(vector);
+        return scale(out, vector, scalar);
+    }
+
+    const tmpProject = Vec3()
+    export function projectOnPlane(out: Vec3, p: Vec3, normal: Vec3) {
+        projectOnVector(tmpProject, p, normal);
+        return sub(out, p, tmpProject);
+    }
+
     /** Get a vector that is similar to `b` but orthogonal to `a` */
     export function orthogonalize(out: Vec3, a: Vec3, b: Vec3) {
         return normalize(out, cross(out, cross(out, a, b), a));

src/mol-model-props/computed/chemistry/util.ts

@@ -8,6 +8,7 @@ import { Structure, Unit } from '../../../mol-model/structure';
 import { StructureElement } from '../../../mol-model/structure/structure';
 import { Elements } from '../../../mol-model/structure/model/properties/atomic/types';
 import { BondType } from '../../../mol-model/structure/model/types';
+import { SortedArray } from '../../../mol-data/int';
 
 export function typeSymbol(unit: Unit.Atomic, index: StructureElement.UnitIndex) {
     return unit.model.atomicHierarchy.atoms.type_symbol.value(unit.elements[index])
@@ -26,7 +27,7 @@ export function altLoc(unit: Unit.Atomic, index: StructureElement.UnitIndex) {
 }
 
 export function compId(unit: Unit.Atomic, index: StructureElement.UnitIndex) {
-    return unit.model.atomicHierarchy.residues.label_comp_id.value(unit.elements[index])
+    return unit.model.atomicHierarchy.residues.label_comp_id.value(unit.getResidueIndex(index))
 }
 
 //
@@ -94,4 +95,16 @@ export function eachIntraBondedAtom(unit: Unit.Atomic, index: StructureElement.U
 export function eachBondedAtom(structure: Structure, unit: Unit.Atomic, index: StructureElement.UnitIndex, cb: (unit: Unit.Atomic, index: StructureElement.UnitIndex) => void): void {
     eachInterBondedAtom(structure, unit, index, cb)
     eachIntraBondedAtom(unit, index, cb)
+}
+
+//
+
+export function eachResidueAtom(unit: Unit.Atomic, index: StructureElement.UnitIndex, cb: (index: StructureElement.UnitIndex) => void): void {
+    const { offsets } = unit.model.atomicHierarchy.residueAtomSegments
+    const rI = unit.getResidueIndex(index)
+    for (let i = offsets[rI], il = offsets[rI + 1]; i < il; ++i) {
+        // TODO optimize, avoid search with .indexOf
+        const idx = SortedArray.indexOf(unit.elements, i)
+        if (idx !== -1) cb(idx as StructureElement.UnitIndex)
+    }
 }

src/mol-model-props/computed/interactions/charged.ts

@@ -0,0 +1,361 @@
+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ * @author Fred Ludlow <Fred.Ludlow@astx.com>
+ *
+ * based in part on NGL (https://github.com/arose/ngl)
+ */
+
+import { ParamDefinition as PD } from '../../../mol-util/param-definition';
+import { Structure, Unit, StructureElement } from '../../../mol-model/structure';
+import { FeaturesBuilder, Features } from './features';
+import { ProteinBackboneAtoms, PolymerNames, BaseNames, ElementSymbol } from '../../../mol-model/structure/model/types';
+import { typeSymbol, atomId, altLoc, eachBondedAtom } from '../chemistry/util';
+import { Elements } from '../../../mol-model/structure/model/properties/atomic/types';
+import { ValenceModelProvider } from '../valence-model';
+import { degToRad } from '../../../mol-math/misc';
+import { FeatureType, FeatureGroup, InteractionType } from './common';
+import { LinkProvider } from './links';
+import { Segmentation, SortedArray } from '../../../mol-data/int';
+import { isGuanidine, isAcetamidine, isPhosphate, isSulfonicAcid, isSulfate, isCarboxylate } from '../chemistry/functional-group';
+import { PrincipalAxes } from '../../../mol-math/linear-algebra/matrix/principal-axes';
+import { getPositions } from '../../../mol-model/structure/util';
+import { Vec3 } from '../../../mol-math/linear-algebra';
+
+export const ChargedParams = {
+    piStackingDistanceMax: PD.Numeric(5.5, { min: 1, max: 8, step: 0.1 }),
+    piStackingOffsetMax: PD.Numeric(2.0, { min: 0, max: 4, step: 0.1 }),
+    piStackingAngleDevMax: PD.Numeric(30, { min: 0, max: 180, step: 1 }),
+    cationPiDistanceMax: PD.Numeric(6.0, { min: 1, max: 8, step: 0.1 }),
+    cationPiOffsetMax: PD.Numeric(2.0, { min: 0, max: 4, step: 0.1 }),
+    ionicDistanceMax: PD.Numeric(5.0, { min: 0, max: 8, step: 0.1 }),
+}
+export type ChargedParams = typeof ChargedParams
+export type ChargedProps = PD.Values<ChargedParams>
+
+//
+
+const PositvelyCharged = ['ARG', 'HIS', 'LYS']
+const NegativelyCharged = ['GLU', 'ASP']
+
+function getUnitValenceModel(structure: Structure, unit: Unit.Atomic) {
+    const valenceModel = ValenceModelProvider.getValue(structure).value
+    if (!valenceModel) throw Error('expected valence model to be available')
+    const unitValenceModel = valenceModel.get(unit.id)
+    if (!unitValenceModel) throw Error('expected valence model for unit to be available')
+    return unitValenceModel
+}
+
+export function addUnitPositiveCharges(structure: Structure, unit: Unit.Atomic, builder: FeaturesBuilder) {
+    const { charge } = getUnitValenceModel(structure, unit)
+    const { elements } = unit
+    const { x, y, z } = unit.model.atomicConformation
+
+    const addedElements = new Set<StructureElement.UnitIndex>()
+
+    const { label_comp_id } = unit.model.atomicHierarchy.residues
+    const residueIt = Segmentation.transientSegments(unit.model.atomicHierarchy.residueAtomSegments, elements)
+
+    while (residueIt.hasNext) {
+        const { index: residueIndex, start, end } = residueIt.move();
+        const compId = label_comp_id.value(residueIndex)
+
+        if (PositvelyCharged.includes(compId)) {
+            builder.startState()
+            for (let j = start as StructureElement.UnitIndex; j < end; ++j) {
+                if (typeSymbol(unit, j) === Elements.N && !ProteinBackboneAtoms.has(atomId(unit, j))) {
+                    builder.pushMember(x[elements[j]], y[elements[j]], z[elements[j]], j)
+                }
+            }
+            builder.finishState(FeatureType.PositiveCharge, FeatureGroup.None)
+        } else if (!PolymerNames.has(compId)) {
+            for (let j = start as StructureElement.UnitIndex; j < end; ++j) {
+                builder.startState()
+                if (typeSymbol(unit, j) === Elements.N && !ProteinBackboneAtoms.has(atomId(unit, j))) {
+                    builder.pushMember(x[elements[j]], y[elements[j]], z[elements[j]], j)
+                }
+                builder.finishState(FeatureType.PositiveCharge, FeatureGroup.None)
+
+                let group = FeatureGroup.None
+                if (isGuanidine(structure, unit, j)) {
+                    group = FeatureGroup.Guanidine
+                } else if (isAcetamidine(structure, unit, j)) {
+                    group = FeatureGroup.Acetamidine
+                }
+                if (group) {
+                    builder.startState()
+                    eachBondedAtom(structure, unit, j, (_, k) => {
+                        if (typeSymbol(unit, k) === Elements.N) {
+                            addedElements.add(k)
+                            builder.pushMember(x[elements[k]], y[elements[k]], z[elements[k]], k)
+                        }
+                    })
+                    builder.finishState(FeatureType.PositiveCharge, group)
+                }
+            }
+
+            for (let j = start as StructureElement.UnitIndex; j < end; ++j) {
+                if (charge[j] > 0 && !addedElements.has(j)) {
+                    builder.add(FeatureType.PositiveCharge, FeatureGroup.None, x[elements[j]], y[elements[j]], z[elements[j]], j)
+                }
+            }
+        }
+    }
+}
+
+export function addUnitNegativeCharges(structure: Structure, unit: Unit.Atomic, builder: FeaturesBuilder) {
+    const { charge } = getUnitValenceModel(structure, unit)
+    const { elements } = unit
+    const { x, y, z } = unit.model.atomicConformation
+
+    const addedElements = new Set<StructureElement.UnitIndex>()
+
+    const { label_comp_id } = unit.model.atomicHierarchy.residues
+    const residueIt = Segmentation.transientSegments(unit.model.atomicHierarchy.residueAtomSegments, elements)
+
+    while (residueIt.hasNext) {
+        const { index: residueIndex, start, end } = residueIt.move();
+        const compId = label_comp_id.value(residueIndex)
+
+        if (NegativelyCharged.includes(compId)) {
+            builder.startState()
+            for (let j = start as StructureElement.UnitIndex; j < end; ++j) {
+                if (typeSymbol(unit, j) === Elements.O && !ProteinBackboneAtoms.has(atomId(unit, j))) {
+                    builder.pushMember(x[elements[j]], y[elements[j]], z[elements[j]], j)
+                }
+            }
+            builder.finishState(FeatureType.NegativeCharge, FeatureGroup.None)
+        } else if (BaseNames.has(compId)) {
+            for (let j = start as StructureElement.UnitIndex; j < end; ++j) {
+                if (isPhosphate(structure, unit, j)) {
+                    builder.startState()
+                    eachBondedAtom(structure, unit, j, (_, k) => {
+                        if (typeSymbol(unit, k) === Elements.O) {
+                            builder.pushMember(x[elements[k]], y[elements[k]], z[elements[k]], k)
+                        }
+                    })
+                    builder.finishState(FeatureType.NegativeCharge, FeatureGroup.Phosphate)
+                }
+            }
+        } else if (!PolymerNames.has(compId)) {
+            for (let j = start as StructureElement.UnitIndex; j < end; ++j) {
+                builder.startState()
+                if (typeSymbol(unit, j) === Elements.N && !ProteinBackboneAtoms.has(atomId(unit, j))) {
+                    builder.pushMember(x[elements[j]], y[elements[j]], z[elements[j]], j)
+                }
+                builder.finishState(FeatureType.PositiveCharge, FeatureGroup.None)
+
+                let group = FeatureGroup.None
+                if (isSulfonicAcid(structure, unit, j)) {
+                    group = FeatureGroup.SulfonicAcid
+                } else if (isPhosphate(structure, unit, j)) {
+                    group = FeatureGroup.Phosphate
+                } else if (isSulfate(structure, unit, j)) {
+                    group = FeatureGroup.Sulfate
+                } else if (isCarboxylate(structure, unit, j)) {
+                    group = FeatureGroup.Carboxylate
+                }
+                if (group) {
+                    builder.startState()
+                    eachBondedAtom(structure, unit, j, (_, k) => {
+                        if (typeSymbol(unit, k) === Elements.O) {
+                            addedElements.add(k)
+                            builder.pushMember(x[elements[k]], y[elements[k]], z[elements[k]], k)
+                        }
+                    })
+                    builder.finishState(FeatureType.PositiveCharge, group)
+                }
+            }
+
+            for (let j = start as StructureElement.UnitIndex; j < end; ++j) {
+                if (charge[j] < 0 && !addedElements.has(j)) {
+                    builder.add(FeatureType.NegativeCharge, FeatureGroup.None, x[elements[j]], y[elements[j]], z[elements[j]], j)
+                }
+            }
+        }
+    }
+}
+
+const AromaticRingElements = [
+    Elements.B, Elements.C, Elements.N, Elements.O,
+    Elements.SI, Elements.P, Elements.S,
+    Elements.GE, Elements.AS,
+    Elements.SN, Elements.SB,
+    Elements.BI
+] as ElementSymbol[]
+const AromaticRingPlanarityThreshold = 0.05
+
+function isRingAromatic(unit: Unit.Atomic, ring: SortedArray<StructureElement.UnitIndex>) {
+    // TODO also check `chem_comp_bond.pdbx_aromatic_flag`
+    let hasAromaticRingElement = false
+    for (let i = 0, il = ring.length; i < il; ++i) {
+        if (AromaticRingElements.includes(typeSymbol(unit, ring[i]))) {
+            hasAromaticRingElement = true
+            break
+        }
+    }
+    if (!hasAromaticRingElement) return
+
+    const ma = PrincipalAxes.calculateMomentsAxes(getPositions(unit, ring))
+    return Vec3.magnitude(ma.dirC) < AromaticRingPlanarityThreshold
+}
+
+export function addUnitAromaticRings(structure: Structure, unit: Unit.Atomic, builder: FeaturesBuilder) {
+    const { elements } = unit
+    const { x, y, z } = unit.model.atomicConformation
+
+    for (const ring of unit.rings.all) {
+        if (!isRingAromatic(unit, ring)) continue
+
+        builder.startState()
+        for (let i = 0, il = ring.length; i < il; ++i) {
+            const j = ring[i]
+            builder.pushMember(x[elements[j]], y[elements[j]], z[elements[j]], j)
+        }
+        builder.finishState(FeatureType.AromaticRing, FeatureGroup.None)
+    }
+}
+
+function isIonicInteraction(ti: FeatureType, tj: FeatureType) {
+    return (
+        (ti === FeatureType.NegativeCharge && tj === FeatureType.PositiveCharge) ||
+        (ti === FeatureType.PositiveCharge && tj === FeatureType.NegativeCharge)
+    )
+}
+
+function isPiStacking(ti: FeatureType, tj: FeatureType) {
+    return ti === FeatureType.AromaticRing && tj === FeatureType.AromaticRing
+}
+
+function isCationPi(ti: FeatureType, tj: FeatureType) {
+    return (
+        (ti === FeatureType.AromaticRing && tj === FeatureType.PositiveCharge) ||
+        (ti === FeatureType.PositiveCharge && tj === FeatureType.AromaticRing)
+    )
+}
+
+function areFeaturesWithinDistanceSq(infoA: Features.Info, infoB: Features.Info, distanceSq: number): boolean {
+    // TODO
+    // const sn = atomSet1.length
+    // const sm = atomSet2.length
+    // for (let si = 0; si < sn; ++si) {
+    //   ap1.index = atomSet1[ si ]
+    //   for (let sj = 0; sj < sm; ++sj) {
+    //     ap2.index = atomSet2[ sj ]
+    //     if (ap1.distanceTo(ap2) <= maxDist) {
+    //       return true
+    //     }
+    //   }
+    // }
+    return false
+}
+
+
+const tmpVecA = Vec3()
+const tmpVecB = Vec3()
+const tmpVecC = Vec3()
+const tmpVecD = Vec3()
+
+function getNormal(out: Vec3, info: Features.Info) {
+    const { unit, feature, offsets} = info
+    const { elements } = unit
+    const { x, y, z } = unit.model.atomicConformation
+
+    const i = offsets[feature]
+    Vec3.set(tmpVecA, x[elements[i]], y[elements[i]], z[elements[i]])
+    Vec3.set(tmpVecB, x[elements[i + 1]], y[elements[i + 1]], z[elements[i + 1]])
+    Vec3.set(tmpVecC, x[elements[i + 2]], y[elements[i + 2]], z[elements[i + 2]])
+
+    return Vec3.triangleNormal(out, tmpVecA, tmpVecB, tmpVecC)
+}
+
+const getOffset = function (infoA: Features.Info, infoB: Features.Info, normal: Vec3) {
+    Vec3.set(tmpVecA, infoA.x[infoA.feature], infoA.y[infoA.feature], infoA.z[infoA.feature])
+    Vec3.set(tmpVecB, infoB.x[infoB.feature], infoB.y[infoB.feature], infoB.z[infoB.feature])
+
+    Vec3.sub(tmpVecC, tmpVecA, tmpVecB)
+
+    Vec3.projectOnPlane(tmpVecD, tmpVecC, normal)
+    Vec3.add(tmpVecD, tmpVecD, tmpVecB)
+    return Vec3.distance(tmpVecD, tmpVecB)
+}
+
+function getOptions(props: ChargedProps) {
+    return {
+        piStackingDistanceMaxSq: props.piStackingDistanceMax * props.piStackingDistanceMax,
+        piStackingOffsetMax: props.piStackingOffsetMax,
+        piStackingAngleDevMax: degToRad(props.piStackingAngleDevMax),
+        cationPiDistanceMaxSq: props.cationPiDistanceMax * props.cationPiDistanceMax,
+        cationPiOffsetMax: props.cationPiOffsetMax,
+        ionicDistanceMaxSq: props.ionicDistanceMax * props.ionicDistanceMax,
+    }
+}
+type Options = ReturnType<typeof getOptions>
+
+const deg180InRad = degToRad(180)
+const deg90InRad = degToRad(90)
+
+const tmpNormalA = Vec3()
+const tmpNormalB = Vec3()
+
+function testCharged(structure: Structure, infoA: Features.Info, infoB: Features.Info, distanceSq: number, opts: Options): InteractionType | undefined {
+    const typeA = infoA.types[infoA.feature]
+    const typeB = infoB.types[infoB.feature]
+
+    const indexA = infoA.members[infoA.offsets[infoA.feature]]
+    const indexB = infoB.members[infoB.offsets[infoB.feature]]
+
+    if (indexA === indexB) return // to self
+
+    const altA = altLoc(infoA.unit, indexA)
+    const altB = altLoc(infoB.unit, indexB)
+
+    if (altA && altB && altA !== altB) return // incompatible alternate location id
+    if (infoA.unit.residueIndex[infoA.unit.elements[indexA]] === infoB.unit.residueIndex[infoB.unit.elements[indexB]]) return // same residue
+
+    if (isIonicInteraction(typeA, typeB)) {
+        if (areFeaturesWithinDistanceSq(infoA, infoB, opts.ionicDistanceMaxSq)) {
+            return InteractionType.IonicInteraction
+        }
+    } else if (isPiStacking(typeA, typeB)) {
+        if (distanceSq <= opts.piStackingDistanceMaxSq) {
+            getNormal(tmpNormalA, infoA)
+            getNormal(tmpNormalB, infoB)
+
+            const angle = Vec3.angle(tmpNormalA, tmpNormalB)
+            const offset = Math.min(getOffset(infoA, infoB, tmpNormalB), getOffset(infoB, infoA, tmpNormalA))
+            if (offset <= opts.piStackingOffsetMax) {
+                if (angle <= opts.piStackingAngleDevMax || angle >= deg180InRad - opts.piStackingAngleDevMax) {
+                    return InteractionType.PiStacking  // parallel
+                } else if (angle <= opts.piStackingAngleDevMax + deg90InRad && angle >= deg90InRad - opts.piStackingAngleDevMax) {
+                    return InteractionType.PiStacking  // t-shaped
+                }
+            }
+        }
+    } else if (isCationPi(typeA, typeB)) {
+        if (distanceSq <= opts.cationPiDistanceMaxSq) {
+            const [infoR, infoC] = typeA === FeatureType.AromaticRing ? [infoA, infoB] : [infoB, infoA]
+
+            getNormal(tmpNormalA, infoR)
+            const offset = getOffset(infoC, infoR, tmpNormalA)
+            if (offset <= opts.cationPiOffsetMax) {
+                return InteractionType.CationPi
+            }
+        }
+    }
+
+}
+
+export const ChargedProvider: LinkProvider<ChargedParams> = {
+    name: 'charged',
+    params: ChargedParams,
+    createTester: (props: ChargedProps) => {
+        const maxDistance = Math.max(props.ionicDistanceMax + 2, props.piStackingDistanceMax, props.cationPiDistanceMax)
+        const opts = getOptions(props)
+        return {
+            maxDistanceSq: maxDistance * maxDistance,
+            getType: (structure, infoA, infoB, distanceSq) => testCharged(structure, infoA, infoB, distanceSq, opts)
+        }
+    }
+}

src/mol-model-props/computed/interactions/features.ts

@@ -103,21 +103,28 @@ namespace Features {
         unit: Unit.Atomic,
         types: ArrayLike<FeatureType>,
         feature: number,
+        x: ArrayLike<number>
+        y: ArrayLike<number>
+        z: ArrayLike<number>
         members: ArrayLike<StructureElement.UnitIndex>,
         offsets: ArrayLike<number>,
         idealGeometry: Int8Array
     }
-    export function Info(structure: Structure, unit: Unit.Atomic, features: Features) {
+    export function Info(structure: Structure, unit: Unit.Atomic, features: Features): Info {
         const valenceModel = ValenceModelProvider.getValue(structure).value
         if (!valenceModel || !valenceModel.has(unit.id)) throw new Error('valence model required')
 
         return {
             unit,
             types: features.types,
+            feature: -1,
+            x: features.x,
+            y: features.y,
+            z: features.z,
             members: features.members,
             offsets: features.offsets,
             idealGeometry: valenceModel.get(unit.id)!.idealGeometry
-        } as Info
+        }
     }
 
     export interface Provider {
@@ -129,10 +136,10 @@ namespace Features {
 export { FeaturesBuilder }
 
 interface FeaturesBuilder {
-    clearState: () => void
+    startState: () => void
     pushMember: (x: number, y: number, z: number, member: StructureElement.UnitIndex) => void
-    addState: (type: FeatureType, group: FeatureGroup) => void
-    addOne: (type: FeatureType, group: FeatureGroup, x: number, y: number, z: number, member: StructureElement.UnitIndex) => void
+    finishState: (type: FeatureType, group: FeatureGroup) => void
+    add: (type: FeatureType, group: FeatureGroup, x: number, y: number, z: number, member: StructureElement.UnitIndex) => void
     getFeatures: (elementsCount: number) => Features
 }
 
@@ -151,14 +158,21 @@ namespace FeaturesBuilder {
         const state: State = { x: 0, y: 0, z: 0, offset: 0, count: 0 }
 
         return {
-            clearState: () => {
-                state.x = 0, state.y = 0, state.z = 0, state.offset = members.elementCount, state.count = 0
+            startState: () => {
+                state.x = 0
+                state.y = 0
+                state.z = 0
+                state.offset = members.elementCount
+                state.count = 0
             },
             pushMember: (x: number, y: number, z: number, member: StructureElement.UnitIndex) => {
                 ChunkedArray.add(members, member)
-                state.x += x, state.y += y, state.z += z
+                state.x += x
+                state.y += y
+                state.z += z
+                state.count += 1
             },
-            addState: (type: FeatureType, group: FeatureGroup) => {
+            finishState: (type: FeatureType, group: FeatureGroup) => {
                 const { count } = state
                 if (count === 0) return
                 ChunkedArray.add(types, type)
@@ -168,7 +182,7 @@ namespace FeaturesBuilder {
                 ChunkedArray.add(zCenters, state.z / count)
                 ChunkedArray.add(offsets, state.offset)
             },
-            addOne: (type: FeatureType, group: FeatureGroup, x: number, y: number, z: number, member: StructureElement.UnitIndex) => {
+            add: (type: FeatureType, group: FeatureGroup, x: number, y: number, z: number, member: StructureElement.UnitIndex) => {
                 ChunkedArray.add(types, type)
                 ChunkedArray.add(groups, group)
                 ChunkedArray.add(xCenters, x)

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

@@ -37,7 +37,7 @@ export function addUnitHalogenDonors (structure: Structure, unit: Unit.Atomic, b
     for (let i = 0 as StructureElement.UnitIndex, il = elements.length; i < il; ++i) {
         const element = typeSymbol(unit, i)
         if (halBondElements.includes(element)) {
-            builder.addOne(FeatureType.HalogenDonor, FeatureGroup.None, x[elements[i]], y[elements[i]], z[elements[i]], i)
+            builder.add(FeatureType.HalogenDonor, FeatureGroup.None, x[elements[i]], y[elements[i]], z[elements[i]], i)
         }
     }
 }
@@ -62,17 +62,17 @@ export function addUnitHalogenAcceptors (structure: Structure, unit: Unit.Atomic
                 }
             })
             if (flag) {
-                builder.addOne(FeatureType.HalogenAcceptor, FeatureGroup.None, x[elements[i]], y[elements[i]], z[elements[i]], i)
+                builder.add(FeatureType.HalogenAcceptor, FeatureGroup.None, x[elements[i]], y[elements[i]], z[elements[i]], i)
             }
         }
     }
 }
 
 function isHalogenBond (ti: FeatureType, tj: FeatureType) {
-  return (
-    (ti === FeatureType.HalogenAcceptor && tj === FeatureType.HalogenDonor) ||
-    (ti === FeatureType.HalogenDonor && tj === FeatureType.HalogenAcceptor)
-  )
+    return (
+        (ti === FeatureType.HalogenAcceptor && tj === FeatureType.HalogenDonor) ||
+        (ti === FeatureType.HalogenDonor && tj === FeatureType.HalogenAcceptor)
+    )
 }
 
 // http://www.pnas.org/content/101/48/16789.full

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

@@ -66,7 +66,7 @@ export function addUnitHydrogenDonors(structure: Structure, unit: Unit.Atomic, b
                 (element === Elements.N || element === Elements.O || element === Elements.S)
             )
         ) {
-            builder.addOne(FeatureType.HydrogenDonor, FeatureGroup.None, x[elements[i]], y[elements[i]], z[elements[i]], i)
+            builder.add(FeatureType.HydrogenDonor, FeatureGroup.None, x[elements[i]], y[elements[i]], z[elements[i]], i)
         }
     }
 }
@@ -89,7 +89,7 @@ export function addUnitWeakHydrogenDonors(structure: Structure, unit: Unit.Atomi
                 inAromaticRingWithElectronNegativeElement(structure, unit, i)
             )
         ) {
-            builder.addOne(FeatureType.WeakHydrogenDonor, FeatureGroup.None, x[elements[i]], y[elements[i]], z[elements[i]], i)
+            builder.add(FeatureType.WeakHydrogenDonor, FeatureGroup.None, x[elements[i]], y[elements[i]], z[elements[i]], i)
         }
     }
 }
@@ -126,7 +126,7 @@ export function addUnitHydrogenAcceptors(structure: Structure, unit: Unit.Atomic
     const { x, y, z } = unit.model.atomicConformation
 
     function add(i: StructureElement.UnitIndex) {
-        builder.addOne(FeatureType.HydrogenAcceptor, FeatureGroup.None, x[elements[i]], y[elements[i]], z[elements[i]], i)
+        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) {

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

@@ -16,6 +16,7 @@ import { Vec3 } from '../../../mol-math/linear-algebra';
 import { addUnitLinks, LinkTester, addStructureLinks } from './links';
 import { addUnitHalogenDonors, addUnitHalogenAcceptors, HalogenBondsProvider } from './halogen-bonds';
 import { addUnitHydrogenDonors, addUnitWeakHydrogenDonors, addUnitHydrogenAcceptors, HydrogenBondsProvider } from './hydrogen-bonds';
+import { addUnitPositiveCharges, addUnitNegativeCharges, addUnitAromaticRings, ChargedProvider } from './charged';
 
 export { Interactions }
 
@@ -101,6 +102,7 @@ namespace Interactions {
 export const InteractionsParams = {
     hydrogenBonds: PD.Group(HydrogenBondsProvider.params),
     halogenBonds: PD.Group(HalogenBondsProvider.params),
+    charged: PD.Group(ChargedProvider.params),
 }
 export type InteractionsParams = typeof InteractionsParams
 export type InteractionsProps = PD.Values<InteractionsParams>
@@ -111,7 +113,8 @@ export async function computeInteractions(runtime: RuntimeContext, structure: St
 
     const linkTesters: LinkTester[] = [
         HydrogenBondsProvider.createTester(p.hydrogenBonds),
-        HalogenBondsProvider.createTester(p.halogenBonds)
+        HalogenBondsProvider.createTester(p.halogenBonds),
+        ChargedProvider.createTester(p.charged),
     ]
 
     const unitsFeatures = IntMap.Mutable<Features>()
@@ -139,6 +142,10 @@ const FeatureProviders: Features.Provider[] = [
 
     { name: 'halogen-donors', add: addUnitHalogenDonors },
     { name: 'halogen-acceptors', add: addUnitHalogenAcceptors },
+
+    { name: 'positive-charges', add: addUnitPositiveCharges },
+    { name: 'negative-charges', add: addUnitNegativeCharges },
+    { name: 'aromatic-rings', add: addUnitAromaticRings },
 ]
 
 function findIntraUnitLinksAndFeatures(structure: Structure, unit: Unit, linkTesters: ReadonlyArray<LinkTester>) {

src/mol-repr/structure/visual/aromatic-ring-mesh.ts

@@ -0,0 +1 @@
+// TODO