wip, interactions, halogen bonds

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

+/**
+ * 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 { calcAngles } from '../chemistry/geometry';
+import { FeaturesBuilder, Features } from './features';
+import { ElementSymbol } from '../../../mol-model/structure/model/types';
+import { typeSymbol, altLoc, eachBondedAtom } from '../chemistry/util';
+import { Elements } from '../../../mol-model/structure/model/properties/atomic/types';
+import { degToRad } from '../../../mol-math/misc';
+import { FeatureType, FeatureGroup, InteractionType } from './common';
+import { IntraLinksBuilder, InterLinksBuilder } from './builder';
+import { Mat4, Vec3 } from '../../../mol-math/linear-algebra';
+
+export const HalogenBondsParams = {
+    distanceMax: PD.Numeric(4.0, { min: 1, max: 5, step: 0.1 }),
+    angleMax: PD.Numeric(30, { min: 0, max: 60, step: 1 }),
+}
+export type HalogenBondsParams = typeof HalogenBondsParams
+export type HalogenBondsProps = PD.Values<HalogenBondsParams>
+
+const halBondElements = [Elements.CL, Elements.BR, Elements.I, Elements.AT] as ElementSymbol[]
+
+/**
+ * Halogen bond donors (X-C, with X one of Cl, Br, I or At) not F!
+ */
+export function addUnitHalogenDonors (structure: Structure, unit: Unit.Atomic, builder: FeaturesBuilder) {
+    const { elements } = unit
+    const { x, y, z } = unit.model.atomicConformation
+
+    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)
+        }
+    }
+}
+
+const X = [Elements.N, Elements.O, Elements.S] as ElementSymbol[]
+const Y = [Elements.C, Elements.N, Elements.P, Elements.S] as ElementSymbol[]
+
+/**
+ * Halogen bond acceptors (Y-{O|N|S}, with Y=C,P,N,S)
+ */
+export function addUnitHalogenAcceptors (structure: Structure, unit: Unit.Atomic, builder: FeaturesBuilder) {
+    const { elements } = unit
+    const { x, y, z } = unit.model.atomicConformation
+
+    for (let i = 0 as StructureElement.UnitIndex, il = elements.length; i < il; ++i) {
+        const element = typeSymbol(unit, i)
+        if (X.includes(element)) {
+            let flag = false
+            eachBondedAtom(structure, unit, i, (unitB, indexB) => {
+                if (Y.includes(typeSymbol(unitB, indexB))) {
+                    flag = true
+                }
+            })
+            if (flag) {
+                builder.addOne(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)
+  )
+}
+
+// http://www.pnas.org/content/101/48/16789.full
+const OptimalHalogenAngle = degToRad(180)  // adjusted from 165 to account for spherical statistics
+const OptimalAcceptorAngle = degToRad(120)
+
+interface Info {
+    unit: Unit.Atomic,
+    types: ArrayLike<FeatureType>,
+    feature: number,
+    members: ArrayLike<StructureElement.UnitIndex>,
+    offsets: ArrayLike<number>,
+}
+function Info(structure: Structure, unit: Unit.Atomic, features: Features) {
+    return {
+        unit,
+        types: features.types,
+        members: features.members,
+        offsets: features.offsets,
+    } as Info
+}
+
+function getOptions(props: HalogenBondsProps) {
+    return {
+        distanceMax: props.distanceMax,
+        angleMax: degToRad(props.angleMax),
+    }
+}
+type Options = ReturnType<typeof getOptions>
+
+function testHalogenBond(structure: Structure, infoA: Info, infoB: Info, opts: Options): InteractionType | undefined {
+    const typeA = infoA.types[infoA.feature]
+    const typeB = infoB.types[infoB.feature]
+
+    if (!isHalogenBond(typeA, typeB)) return
+
+    const [don, acc] = typeA === FeatureType.HalogenDonor ? [infoA, infoB] : [infoB, infoA]
+
+    const donIndex = don.members[don.offsets[don.feature]]
+    const accIndex = acc.members[acc.offsets[acc.feature]]
+
+    if (accIndex === donIndex) return // DA to self
+
+    const altD = altLoc(don.unit, donIndex)
+    const altA = altLoc(acc.unit, accIndex)
+
+    if (altD && altA && altD !== altA) return // incompatible alternate location id
+    if (don.unit.residueIndex[don.unit.elements[donIndex]] === acc.unit.residueIndex[acc.unit.elements[accIndex]]) return // same residue
+
+    const halogenAngles = calcAngles(structure, don.unit, donIndex, acc.unit, accIndex)
+    // Singly bonded halogen only (not bromide ion for example)
+    if (halogenAngles.length !== 1) return
+    if (OptimalHalogenAngle - halogenAngles[0] > opts.angleMax) return
+
+    const acceptorAngles = calcAngles(structure, acc.unit, accIndex, don.unit, donIndex)
+    // Angle must be defined. Excludes water as acceptor. Debatable
+    if (acceptorAngles.length === 0) return
+    if (acceptorAngles.some(acceptorAngle => OptimalAcceptorAngle - acceptorAngle > opts.angleMax)) return
+
+    return InteractionType.HalogenBond
+}
+
+/**
+ * All intra-unit pairs of halogen donor and acceptor atoms
+ */
+export function addUnitHalogenBonds(structure: Structure, unit: Unit.Atomic, features: Features, builder: IntraLinksBuilder, props: HalogenBondsProps) {
+    const opts = getOptions(props)
+    const { x, y, z, count, lookup3d } = features
+
+    const infoA = Info(structure, unit, features)
+    const infoB = { ...infoA }
+
+    for (let i = 0; i < count; ++i) {
+        const { count, indices } = lookup3d.find(x[i], y[i], z[i], opts.distanceMax)
+        if (count === 0) continue
+
+        infoA.feature = i
+
+        for (let r = 0; r < count; ++r) {
+            const j = indices[r]
+            if (j <= i) continue
+
+            infoB.feature = j
+            const bondType = testHalogenBond(structure, infoA, infoB, opts)
+            if (bondType) builder.add(i, j, bondType)
+        }
+    }
+}
+
+const _imageTransform = Mat4()
+
+//
+
+/**
+ * All inter-unit pairs of halogen donor and acceptor atoms
+ */
+export function addStructureHalogenBonds(structure: Structure, unitA: Unit.Atomic, featuresA: Features, unitB: Unit.Atomic, featuresB: Features, builder: InterLinksBuilder, props: HalogenBondsProps) {
+    const opts = getOptions(props)
+
+    const { count, x: xA, y: yA, z: zA } = featuresA;
+    const { lookup3d } = featuresB;
+
+    // the lookup queries need to happen in the "unitB space".
+    // that means imageA = inverseOperB(operA(i))
+    const imageTransform = Mat4.mul(_imageTransform, unitB.conformation.operator.inverse, unitA.conformation.operator.matrix)
+    const isNotIdentity = !Mat4.isIdentity(imageTransform)
+    const imageA = Vec3()
+
+    const { center: bCenter, radius: bRadius } = lookup3d.boundary.sphere;
+    const testDistanceSq = (bRadius + opts.distanceMax) * (bRadius + opts.distanceMax);
+
+    const infoA = Info(structure, unitA, featuresA)
+    const infoB = Info(structure, unitB, featuresB)
+
+    builder.startUnitPair(unitA, unitB)
+
+    for (let i = 0; i < count; ++i) {
+        Vec3.set(imageA, xA[i], yA[i], zA[i])
+        if (isNotIdentity) Vec3.transformMat4(imageA, imageA, imageTransform)
+        if (Vec3.squaredDistance(imageA, bCenter) > testDistanceSq) continue
+
+        const { indices, count } = lookup3d.find(imageA[0], imageA[1], imageA[2], opts.distanceMax)
+        if (count === 0) continue
+
+        infoA.feature = i
+
+        for (let r = 0; r < count; ++r) {
+            const j = indices[r]
+            infoB.feature = j
+            const bondType = testHalogenBond(structure, infoA, infoB, opts)
+            if (bondType) builder.add(i, j, bondType)
+        }
+    }
+
+    builder.finishUnitPair()
+}
\ No newline at end of file

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

@@ -2,6 +2,9 @@
  * 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';

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

@@ -14,6 +14,7 @@ import { InteractionsIntraLinks, InteractionsInterLinks } from './common';
 import { IntraLinksBuilder, InterLinksBuilder } from './builder';
 import { IntMap } from '../../../mol-data/int';
 import { Vec3 } from '../../../mol-math/linear-algebra';
+import { addUnitHalogenBonds, HalogenBondsParams, addUnitHalogenDonors, addUnitHalogenAcceptors, addStructureHalogenBonds } from './halogen-bonds';
 
 export { Interactions }
 
@@ -98,6 +99,7 @@ namespace Interactions {
 
 export const InteractionsParams = {
     hydrogenBonds: PD.Group(HydrogenBondsParams),
+    halogenBonds: PD.Group(HalogenBondsParams),
 }
 export type InteractionsParams = typeof InteractionsParams
 export type InteractionsProps = PD.Values<InteractionsParams>
@@ -131,12 +133,16 @@ function findIntraUnitLinksAndFeatures(structure: Structure, unit: Unit, props:
         addUnitHydrogenDonors(structure, unit, featuresBuilder)
         addUnitWeakHydrogenDonors(structure, unit, featuresBuilder)
         addUnitHydrogenAcceptors(structure, unit, featuresBuilder)
+
+        addUnitHalogenDonors(structure, unit, featuresBuilder)
+        addUnitHalogenAcceptors(structure, unit, featuresBuilder)
     }
     const features = featuresBuilder.getFeatures(count)
 
     const linksBuilder = IntraLinksBuilder.create(features, count)
     if (Unit.isAtomic(unit)) {
         addUnitHydrogenBonds(structure, unit, features, linksBuilder, props.hydrogenBonds)
+        addUnitHalogenBonds(structure, unit, features, linksBuilder, props.halogenBonds)
     }
 
     return { features, links: linksBuilder.getLinks() }
@@ -167,8 +173,10 @@ function findInterUnitLinks(structure: Structure, unitsFeatures: IntMap<Features
 
             if (unitB.elements.length >= unitA.elements.length) {
                 addStructureHydrogenBonds(structure, unitA, featuresA, unitB, featuresB, builder, props.hydrogenBonds)
+                addStructureHalogenBonds(structure, unitA, featuresA, unitB, featuresB, builder, props.halogenBonds)
             } else {
                 addStructureHydrogenBonds(structure, unitB, featuresB, unitA, featuresA, builder, props.hydrogenBonds)
+                addStructureHalogenBonds(structure, unitB, featuresB, unitA, featuresA, builder, props.halogenBonds)
             }
         }
     }