diff --git a/src/mol-model-props/computed/interactions/halogen-bonds.ts b/src/mol-model-props/computed/interactions/halogen-bonds.ts new file mode 100644 index 0000000000000000000000000000000000000000..5487be38e05b7249de7096230c180c8b369cabe0 --- /dev/null +++ b/src/mol-model-props/computed/interactions/halogen-bonds.ts @@ -0,0 +1,212 @@ +/** + * 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 diff --git a/src/mol-model-props/computed/interactions/hydrogen-bonds.ts b/src/mol-model-props/computed/interactions/hydrogen-bonds.ts index dda69658d4253615aeb6cc8c8edb92972b598b28..67177f5eec5e0c3f3017f6035f2dd1e4496963a7 100644 --- a/src/mol-model-props/computed/interactions/hydrogen-bonds.ts +++ b/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'; diff --git a/src/mol-model-props/computed/interactions/interactions.ts b/src/mol-model-props/computed/interactions/interactions.ts index fe0a53513f37b222ce33668a2a9f374a03b745e3..8c69fe5c6107688f741c4c28f74ea5c2cf7ee4f2 100644 --- a/src/mol-model-props/computed/interactions/interactions.ts +++ b/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) } } }