index-pairs bonds fixes & improvements, better cif-core bond handling

src/mol-model-formats/structure/cif-core.ts

@@ -6,7 +6,7 @@
 
 import { Column, Table } from '../../mol-data/db';
 import { Model, Symmetry } from '../../mol-model/structure/model';
-import { MoleculeType } from '../../mol-model/structure/model/types';
+import { MoleculeType, BondType } from '../../mol-model/structure/model/types';
 import { RuntimeContext, Task } from '../../mol-task';
 import { createModels } from './basic/parser';
 import { BasicSchema, createBasic } from './basic/schema';
@@ -22,6 +22,7 @@ import { IndexPairBonds } from './property/bonds/index-pair';
 import { AtomSiteAnisotrop } from './property/anisotropic';
 import { guessElementSymbolString } from './util';
 import { Trajectory } from '../../mol-model/structure';
+import { cantorPairing } from '../../mol-data/util';
 
 function getSpacegroupNameOrNumber(space_group: CifCore_Database['space_group']) {
     const groupNumber = space_group.IT_number.value(0);
@@ -36,7 +37,7 @@ function getSymmetry(db: CifCore_Database): Symmetry {
     const nameOrNumber = getSpacegroupNameOrNumber(space_group);
     const spaceCell = SpacegroupCell.create(nameOrNumber,
         Vec3.create(cell.length_a.value(0), cell.length_b.value(0), cell.length_c.value(0)),
-        Vec3.scale(Vec3.zero(), Vec3.create(cell.angle_alpha.value(0), cell.angle_beta.value(0), cell.angle_gamma.value(0)), Math.PI / 180));
+        Vec3.scale(Vec3(), Vec3.create(cell.angle_alpha.value(0), cell.angle_beta.value(0), cell.angle_gamma.value(0)), Math.PI / 180));
 
     return {
         spacegroup: Spacegroup.create(spaceCell),
@@ -166,28 +167,59 @@ async function getModels(db: CifCore_Database, format: CifCoreFormat, ctx: Runti
                 labelIndexMap[label.value(i)] = i;
             }
 
+            const bond_type = format.data.frame.categories.ccdc_geom_bond_type?.getField('');
+
             const indexA: number[] = [];
             const indexB: number[] = [];
             const order: number[] = [];
-            const symmetryA: string[] = [];
-            const symmetryB: string[] = [];
+            const dist: number[] = [];
+            const flag: number[] = [];
+
+            const included = new Set<number>();
+            let j = 0;
 
-            const { atom_site_label_1, atom_site_label_2, valence, site_symmetry_1, site_symmetry_2 } = db.geom_bond;
+            const { atom_site_label_1, atom_site_label_2, valence, distance } = db.geom_bond;
             for (let i = 0; i < bondCount; ++i) {
-                indexA[i] = labelIndexMap[atom_site_label_1.value(i)];
-                indexB[i] = labelIndexMap[atom_site_label_2.value(i)];
-                // TODO derive order from bond length if undefined
-                order[i] = valence.isDefined ? valence.value(i) : 1;
-                symmetryA[i] = site_symmetry_1.value(i) || '1_555';
-                symmetryB[i] = site_symmetry_2.value(i) || '1_555';
+                const iA = labelIndexMap[atom_site_label_1.value(i)];
+                const iB = labelIndexMap[atom_site_label_2.value(i)];
+                const id = iA < iB ? cantorPairing(iA, iB) : cantorPairing(iB, iA);
+                if (included.has(id)) continue;
+                included.add(id);
+
+                indexA[j] = iA;
+                indexB[j] = iB;
+                dist[j] = distance.value(i) || -1;
+
+                if (bond_type) {
+                    const t = bond_type.str(i);
+                    if (t === 'D') {
+                        order[j] = 2;
+                        flag[j] = BondType.Flag.Covalent;
+                    } else if (t === 'A') {
+                        order[j] = 1;
+                        flag[j] = BondType.Flag.Covalent | BondType.Flag.Aromatic;
+                    } else if (t === 'S') {
+                        order[j] = 1;
+                        flag[j] = BondType.Flag.Covalent;
+                    } else {
+                        order[j] = 1;
+                        flag[j] = BondType.Flag.Covalent;
+                    }
+                } else {
+                    flag[j] = BondType.Flag.Covalent;
+                    // TODO derive order from bond length if undefined
+                    order[j] = valence.isDefined ? valence.value(i) : 1;
+                }
+
+                j += 1;
             }
 
             IndexPairBonds.Provider.set(first, IndexPairBonds.fromData({ pairs: {
                 indexA: Column.ofIntArray(indexA),
                 indexB: Column.ofIntArray(indexB),
                 order: Column.ofIntArray(order),
-                symmetryA: Column.ofStringArray(symmetryA),
-                symmetryB: Column.ofStringArray(symmetryB)
+                distance: Column.ofFloatArray(dist),
+                flag: Column.ofIntArray(flag)
             }, count: indexA.length }));
         }
     }

src/mol-model-formats/structure/property/bonds/index-pair.ts

@@ -8,27 +8,29 @@ import { CustomPropertyDescriptor } from '../../../../mol-model/custom-property'
 import { IntAdjacencyGraph } from '../../../../mol-math/graph';
 import { Column } from '../../../../mol-data/db';
 import { FormatPropertyProvider } from '../../common/property';
+import { BondType } from '../../../../mol-model/structure/model/types';
+import { ElementIndex } from '../../../../mol-model/structure';
 
 export type IndexPairBondsProps = {
     readonly order: ArrayLike<number>
-    readonly symmetryA: ArrayLike<string>
-    readonly symmetryB: ArrayLike<string>
+    readonly distance: ArrayLike<number>
+    readonly flag: ArrayLike<BondType.Flag>
 }
-export type IndexPairBonds = IntAdjacencyGraph<number, IndexPairBondsProps>
+export type IndexPairBonds = IntAdjacencyGraph<ElementIndex, IndexPairBondsProps>
 
-function getGraph(indexA: ArrayLike<number>, indexB: ArrayLike<number>, props: Partial<IndexPairBondsProps>, count: number): IndexPairBonds {
+function getGraph(indexA: ArrayLike<ElementIndex>, indexB: ArrayLike<ElementIndex>, props: Partial<IndexPairBondsProps>, count: number): IndexPairBonds {
     const builder = new IntAdjacencyGraph.EdgeBuilder(count, indexA, indexB);
     const order = new Int8Array(builder.slotCount);
-    const symmetryA = new Array(builder.slotCount);
-    const symmetryB = new Array(builder.slotCount);
+    const distance = new Array(builder.slotCount);
+    const flag = new Array(builder.slotCount);
     for (let i = 0, _i = builder.edgeCount; i < _i; i++) {
         builder.addNextEdge();
         builder.assignProperty(order, props.order ? props.order[i] : 1);
-        builder.assignProperty(symmetryA, props.symmetryA ? props.symmetryA[i] : '');
-        builder.assignProperty(symmetryB, props.symmetryB ? props.symmetryB[i] : '');
+        builder.assignProperty(distance, props.distance ? props.distance[i] : -1);
+        builder.assignProperty(flag, props.flag ? props.flag[i] : BondType.Flag.Covalent);
     }
 
-    return builder.createGraph({ order, symmetryA, symmetryB });
+    return builder.createGraph({ order, distance, flag });
 }
 
 export namespace IndexPairBonds {
@@ -43,19 +45,19 @@ export namespace IndexPairBonds {
             indexA: Column<number>,
             indexB: Column<number>
             order?: Column<number>,
-            symmetryA?: Column<string>,
-            symmetryB?: Column<string>,
+            distance?: Column<number>,
+            flag?: Column<BondType.Flag>,
         },
         count: number
     }
 
     export function fromData(data: Data) {
         const { pairs, count } = data;
-        const indexA = pairs.indexA.toArray();
-        const indexB = pairs.indexB.toArray();
+        const indexA = pairs.indexA.toArray() as ArrayLike<ElementIndex>;
+        const indexB = pairs.indexB.toArray() as ArrayLike<ElementIndex>;
         const order = pairs.order && pairs.order.toArray();
-        const symmetryA = pairs.symmetryA && pairs.symmetryA.toArray();
-        const symmetryB = pairs.symmetryB && pairs.symmetryB.toArray();
-        return getGraph(indexA, indexB, { order, symmetryA, symmetryB }, count);
+        const distance = pairs.distance && pairs.distance.toArray();
+        const flag = pairs.flag && pairs.flag.toArray();
+        return getGraph(indexA, indexB, { order, distance, flag }, count);
     }
 }

src/mol-model/structure/structure/unit/bonds/inter-compute.ts

@@ -18,6 +18,7 @@ import { getInterBondOrderFromTable } from '../../../model/properties/atomic/bon
 import { IndexPairBonds } from '../../../../../mol-model-formats/structure/property/bonds/index-pair';
 import { InterUnitGraph } from '../../../../../mol-math/graph/inter-unit-graph';
 import { StructConn } from '../../../../../mol-model-formats/structure/property/bonds/struct_conn';
+import { equalEps } from '../../../../../mol-math/linear-algebra/3d/common';
 
 const MAX_RADIUS = 4;
 
@@ -29,7 +30,8 @@ function getDistance(unitA: Unit.Atomic, indexA: ElementIndex, unitB: Unit.Atomi
     return Vec3.distance(tmpDistVecA, tmpDistVecB);
 }
 
-const _imageTransform = Mat4.zero();
+const _imageTransform = Mat4();
+const _imageA = Vec3();
 
 function findPairBonds(unitA: Unit.Atomic, unitB: Unit.Atomic, props: BondComputationProps, builder: InterUnitGraph.Builder<Unit.Atomic, StructureElement.UnitIndex, InterUnitEdgeProps>) {
     const { elements: atomsA, residueIndex: residueIndexA } = unitA;
@@ -50,36 +52,34 @@ function findPairBonds(unitA: Unit.Atomic, unitB: Unit.Atomic, props: BondComput
     const indexPairs = unitA.model === unitB.model && IndexPairBonds.Provider.get(unitA.model);
 
     // the lookup queries need to happen in the "unitB space".
-    // that means imageA = inverseOperB(operA(aI))
+    // that means _imageA = inverseOperB(operA(aI))
     const imageTransform = Mat4.mul(_imageTransform, unitB.conformation.operator.inverse, unitA.conformation.operator.matrix);
     const isNotIdentity = !Mat4.isIdentity(imageTransform);
-    const imageA = Vec3.zero();
 
     const { center: bCenter, radius: bRadius } = lookup3d.boundary.sphere;
     const testDistanceSq = (bRadius + MAX_RADIUS) * (bRadius + MAX_RADIUS);
 
     builder.startUnitPair(unitA, unitB);
-    const symmUnitA = unitA.conformation.operator.name;
-    const symmUnitB = unitB.conformation.operator.name;
 
     for (let _aI = 0 as StructureElement.UnitIndex; _aI < atomCount; _aI++) {
         const aI = atomsA[_aI];
-        Vec3.set(imageA, xA[aI], yA[aI], zA[aI]);
-        if (isNotIdentity) Vec3.transformMat4(imageA, imageA, imageTransform);
-        if (Vec3.squaredDistance(imageA, bCenter) > testDistanceSq) continue;
+        Vec3.set(_imageA, xA[aI], yA[aI], zA[aI]);
+        if (isNotIdentity) Vec3.transformMat4(_imageA, _imageA, imageTransform);
+        if (Vec3.squaredDistance(_imageA, bCenter) > testDistanceSq) continue;
 
         if (!props.forceCompute && indexPairs) {
-            const { order, symmetryA, symmetryB } = indexPairs.edgeProps;
+            const { order, distance, flag } = indexPairs.edgeProps;
             for (let i = indexPairs.offset[aI], il = indexPairs.offset[aI + 1]; i < il; ++i) {
                 const bI = indexPairs.b[i];
                 if (aI >= bI) continue;
 
-                const _bI = SortedArray.indexOf(unitA.elements, bI) as StructureElement.UnitIndex;
+                const _bI = SortedArray.indexOf(unitB.elements, bI) as StructureElement.UnitIndex;
                 if (_bI < 0) continue;
-                if (symmetryA[i] === symmetryB[i]) continue;
-                if (type_symbolA.value(aI) === 'H' && type_symbolB.value(indexPairs.b[i]) === 'H') continue;
-                if (symmUnitA === symmetryA[i] && symmUnitB === symmetryB[i]) {
-                    builder.add(_aI, _bI, { order: order[i], flag: BondType.Flag.Covalent });
+                if (type_symbolA.value(aI) === 'H' && type_symbolB.value(bI) === 'H') continue;
+
+                const d = distance[i];
+                if (d === -1 || equalEps(getDistance(unitA, aI, unitB, bI), d, 0.5)) {
+                    builder.add(_aI, _bI, { order: order[i], flag: flag[i] });
                 }
             }
             continue; // assume `indexPairs` supplies all bonds
@@ -107,7 +107,7 @@ function findPairBonds(unitA: Unit.Atomic, unitB: Unit.Atomic, props: BondComput
 
         const occA = occupancyA.value(aI);
 
-        const { indices, count, squaredDistances } = lookup3d.find(imageA[0], imageA[1], imageA[2], MAX_RADIUS);
+        const { indices, count, squaredDistances } = lookup3d.find(_imageA[0], _imageA[1], _imageA[2], MAX_RADIUS);
         if (count === 0) continue;
 
         const aeI = getElementIdx(type_symbolA.value(aI));

src/mol-model/structure/structure/unit/bonds/intra-compute.ts

@@ -16,6 +16,9 @@ import StructureElement from '../../element';
 import { IndexPairBonds } from '../../../../../mol-model-formats/structure/property/bonds/index-pair';
 import { ComponentBond } from '../../../../../mol-model-formats/structure/property/bonds/chem_comp';
 import { StructConn } from '../../../../../mol-model-formats/structure/property/bonds/struct_conn';
+import { Vec3 } from '../../../../../mol-math/linear-algebra';
+import { ElementIndex } from '../../../model/indexing';
+import { equalEps } from '../../../../../mol-math/linear-algebra/3d/common';
 
 function getGraph(atomA: StructureElement.UnitIndex[], atomB: StructureElement.UnitIndex[], _order: number[], _flags: number[], atomCount: number): IntraUnitBonds {
     const builder = new IntAdjacencyGraph.EdgeBuilder(atomCount, atomA, atomB);
@@ -30,6 +33,14 @@ function getGraph(atomA: StructureElement.UnitIndex[], atomB: StructureElement.U
     return builder.createGraph({ flags, order });
 }
 
+const tmpDistVecA = Vec3();
+const tmpDistVecB = Vec3();
+function getDistance(unit: Unit.Atomic, indexA: ElementIndex, indexB: ElementIndex) {
+    unit.conformation.position(indexA, tmpDistVecA);
+    unit.conformation.position(indexB, tmpDistVecB);
+    return Vec3.distance(tmpDistVecA, tmpDistVecB);
+}
+
 const __structConnAdded = new Set<StructureElement.UnitIndex>();
 
 function _computeBonds(unit: Unit.Atomic, props: BondComputationProps): IntraUnitBonds {
@@ -69,12 +80,15 @@ function _computeBonds(unit: Unit.Atomic, props: BondComputationProps): IntraUni
 
                 const _bI = SortedArray.indexOf(unit.elements, bI) as StructureElement.UnitIndex;
                 if (_bI < 0) continue;
-                if (edgeProps.symmetryA[i] !== edgeProps.symmetryB[i]) continue;
-                if (type_symbol.value(aI) === 'H' && type_symbol.value(indexPairs.b[i]) === 'H') continue;
-                atomA[atomA.length] = _aI;
-                atomB[atomB.length] = _bI;
-                order[order.length] = edgeProps.order[i];
-                flags[flags.length] = BondType.Flag.Covalent;
+                if (type_symbol.value(aI) === 'H' && type_symbol.value(bI) === 'H') continue;
+
+                const d = edgeProps.distance[i];
+                if (d === -1 || equalEps(getDistance(unit, aI, bI), d, 0.5)) {
+                    atomA[atomA.length] = _aI;
+                    atomB[atomB.length] = _bI;
+                    order[order.length] = edgeProps.order[i];
+                    flags[flags.length] = edgeProps.flag[i];
+                }
             }
             continue; // assume `indexPairs` supplies all bonds
         }