Merge branch 'master' of https://github.com/molstar/molstar into anim-rock

CHANGELOG.md

@@ -9,6 +9,10 @@ Note that since we don't clearly distinguish between a public and private interf
 - Add PDBj as a ``pdb-provider`` option
 - Move Viewer APP to a separate file to allow use without importing light theme & index.html
 - Add symmetry support for mol2 files (only spacegroup setting 1)
+- Fix mol2 files element symbol assignment
+- Improve bond assignment from ``IndexPairBonds``
+    - Add ``key`` field for mapping to source data
+    - Fix assignment of bonds with unphysical length
 - Fix label/stats of single atom selection in multi-chain units
 - [Breaking] Add rock animation to trackball controls
     - Add ``animate`` to ``TrackballControlsParams``, remove ``spin`` and ``spinSpeed``

src/mol-model-formats/structure/mol2.ts

@@ -31,8 +31,17 @@ async function getModels(mol2: Mol2File, ctx: RuntimeContext) {
         const A = Column.ofConst('A', atoms.count, Column.Schema.str);
 
         const type_symbol = new Array<string>(atoms.count);
+        let hasAtomType = false;
         for (let i = 0; i < atoms.count; ++i) {
-            type_symbol[i] = guessElementSymbolString(atoms.atom_name.value(i));
+            if (atoms.atom_type.value(i).includes('.')) {
+                hasAtomType = true;
+                break;
+            }
+        }
+        for (let i = 0; i < atoms.count; ++i) {
+            type_symbol[i] = hasAtomType
+                ? atoms.atom_type.value(i).split('.')[0].toUpperCase()
+                : guessElementSymbolString(atoms.atom_name.value(i));
         }
 
         const atom_site = Table.ofPartialColumns(BasicSchema.atom_site, {
@@ -77,6 +86,7 @@ async function getModels(mol2: Mol2File, ctx: RuntimeContext) {
         if (_models.frameCount > 0) {
             const indexA = Column.ofIntArray(Column.mapToArray(bonds.origin_atom_id, x => x - 1, Int32Array));
             const indexB = Column.ofIntArray(Column.mapToArray(bonds.target_atom_id, x => x - 1, Int32Array));
+            const key = bonds.bond_id;
             const order = Column.ofIntArray(Column.mapToArray(bonds.bond_type, x => {
                 switch (x) {
                     case 'ar': // aromatic
@@ -103,7 +113,7 @@ async function getModels(mol2: Mol2File, ctx: RuntimeContext) {
                         return BondType.Flag.Covalent;
                 }
             }, Int8Array));
-            const pairBonds = IndexPairBonds.fromData({ pairs: { indexA, indexB, order, flag }, count: atoms.count });
+            const pairBonds = IndexPairBonds.fromData({ pairs: { key, indexA, indexB, order, flag }, count: atoms.count });
 
             const first = _models.representative;
             IndexPairBonds.Provider.set(first, pairBonds);

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

 /**
- * Copyright (c) 2019-2021 Mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2019-2022 Mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
  */
@@ -10,9 +10,9 @@ 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';
-import { DefaultBondMaxRadius } from '../../../../mol-model/structure/structure/unit/bonds/common';
 
 export type IndexPairsProps = {
+    readonly key: ArrayLike<number>
     readonly order: ArrayLike<number>
     readonly distance: ArrayLike<number>
     readonly flag: ArrayLike<BondType.Flag>
@@ -22,17 +22,19 @@ export type IndexPairBonds = { bonds: IndexPairs, maxDistance: number }
 
 function getGraph(indexA: ArrayLike<ElementIndex>, indexB: ArrayLike<ElementIndex>, props: Partial<IndexPairsProps>, count: number): IndexPairs {
     const builder = new IntAdjacencyGraph.EdgeBuilder(count, indexA, indexB);
+    const key = new Int32Array(builder.slotCount);
     const order = new Int8Array(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(key, props.key ? props.key[i] : -1);
         builder.assignProperty(order, props.order ? props.order[i] : 1);
         builder.assignProperty(distance, props.distance ? props.distance[i] : -1);
         builder.assignProperty(flag, props.flag ? props.flag[i] : BondType.Flag.Covalent);
     }
 
-    return builder.createGraph({ order, distance, flag });
+    return builder.createGraph({ key, order, distance, flag });
 }
 
 export namespace IndexPairBonds {
@@ -45,15 +47,33 @@ export namespace IndexPairBonds {
     export type Data = {
         pairs: {
             indexA: Column<number>,
-            indexB: Column<number>
+            indexB: Column<number>,
+            key?: Column<number>,
             order?: Column<number>,
+            /**
+             * Useful for bonds in periodic cells. That is, only bonds within the given
+             * distance are added. This allows for bond between periodic image but
+             * avoids unwanted bonds with wrong distances. If negative, test using the
+             * `maxDistance` option from `Props`.
+             */
             distance?: Column<number>,
             flag?: Column<BondType.Flag>,
         },
         count: number
     }
 
-    export const DefaultProps = { maxDistance: DefaultBondMaxRadius };
+    export const DefaultProps = {
+        /**
+         * If negative, test using element-based threshold, otherwise distance in Angstrom.
+         *
+         * This option exists to handle bonds in periodic cells. For systems that are
+         * made from beads (as opposed to atomic elements), set to a specific distance.
+         *
+         * Note that `Data` has a `distance` field which allows specifying a distance
+         * for each bond individually which takes precedence over this option.
+         */
+        maxDistance: -1
+    };
     export type Props = typeof DefaultProps
 
     export function fromData(data: Data, props: Partial<Props> = {}): IndexPairBonds {
@@ -61,11 +81,12 @@ export namespace IndexPairBonds {
         const { pairs, count } = data;
         const indexA = pairs.indexA.toArray() as ArrayLike<ElementIndex>;
         const indexB = pairs.indexB.toArray() as ArrayLike<ElementIndex>;
+        const key = pairs.key && pairs.key.toArray();
         const order = pairs.order && pairs.order.toArray();
         const distance = pairs.distance && pairs.distance.toArray();
         const flag = pairs.flag && pairs.flag.toArray();
         return {
-            bonds: getGraph(indexA, indexB, { order, distance, flag }, count),
+            bonds: getGraph(indexA, indexB, { key, order, distance, flag }, count),
             maxDistance: p.maxDistance
         };
     }

src/mol-model/structure/structure/unit/bonds/common.ts

 /**
- * Copyright (c) 2017-2021 Mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2017-2022 Mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author David Sehnal <david.sehnal@gmail.com>
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
@@ -110,6 +110,15 @@ export function getElementThreshold(i: number) {
     return r;
 }
 
+export function getPairingThreshold(elementIndexA: number, elementIndexB: number, thresholdA: number, thresholdB: number) {
+    const thresholdAB = getElementPairThreshold(elementIndexA, elementIndexB);
+    return thresholdAB > 0
+        ? thresholdAB
+        : elementIndexB < 0
+            ? thresholdA
+            : (thresholdA + thresholdB) / 1.95; // not sure if avg or min but max is too big
+}
+
 const H_ID = __ElementIndex['H']!;
 export function isHydrogen(i: number) {
     return i === H_ID;

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

 /**
- * Copyright (c) 2017-2021 Mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2017-2022 Mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author David Sehnal <david.sehnal@gmail.com>
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
@@ -8,7 +8,7 @@
 import { BondType, MoleculeType } from '../../../model/types';
 import { Structure } from '../../structure';
 import { Unit } from '../../unit';
-import { getElementIdx, getElementPairThreshold, getElementThreshold, isHydrogen, BondComputationProps, MetalsSet, DefaultBondComputationProps } from './common';
+import { getElementIdx, getElementThreshold, isHydrogen, BondComputationProps, MetalsSet, DefaultBondComputationProps, getPairingThreshold } from './common';
 import { InterUnitBonds, InterUnitEdgeProps } from './data';
 import { SortedArray } from '../../../../../mol-data/int';
 import { Vec3, Mat4 } from '../../../../../mol-math/linear-algebra';
@@ -82,11 +82,31 @@ function findPairBonds(unitA: Unit.Atomic, unitB: Unit.Atomic, props: BondComput
 
                 const _bI = SortedArray.indexOf(unitB.elements, bI) as StructureElement.UnitIndex;
                 if (_bI < 0) continue;
-                if (type_symbolA.value(aI) === 'H' && type_symbolB.value(bI) === 'H') continue;
+
+                const aeI = getElementIdx(type_symbolA.value(aI));
+                const beI = getElementIdx(type_symbolA.value(bI));
 
                 const d = distance[i];
                 const dist = getDistance(unitA, aI, unitB, bI);
-                if ((d !== -1 && equalEps(dist, d, 0.5)) || dist < maxDistance) {
+
+                let add = false;
+                if (d >= 0) {
+                    add = equalEps(dist, d, 0.3);
+                } else if (maxDistance >= 0) {
+                    add = dist < maxDistance;
+                } else {
+                    const pairingThreshold = getPairingThreshold(
+                        aeI, beI, getElementThreshold(aeI), getElementThreshold(beI)
+                    );
+                    add = dist < pairingThreshold;
+
+                    if (isHydrogen(aeI) && isHydrogen(beI)) {
+                        // TODO handle molecular hydrogen
+                        add = false;
+                    }
+                }
+
+                if (add) {
                     builder.add(_aI, _bI, { order: order[i], flag: flag[i] });
                 }
             }
@@ -155,13 +175,7 @@ function findPairBonds(unitA: Unit.Atomic, unitB: Unit.Atomic, props: BondComput
             const dist = Math.sqrt(squaredDistances[ni]);
             if (dist === 0) continue;
 
-            const thresholdAB = getElementPairThreshold(aeI, beI);
-            const pairingThreshold = thresholdAB > 0
-                ? thresholdAB
-                : beI < 0
-                    ? thresholdA
-                    : (thresholdA + getElementThreshold(beI)) / 1.95; // not sure if avg or min but max is too big
-
+            const pairingThreshold = getPairingThreshold(aeI, beI, thresholdA, getElementThreshold(beI));
             if (dist <= pairingThreshold) {
                 const atomIdB = label_atom_idB.value(bI);
                 const compIdB = label_comp_idB.value(residueIndexB[bI]);

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

 /**
- * Copyright (c) 2017-2021 Mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2017-2022 Mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author David Sehnal <david.sehnal@gmail.com>
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
@@ -9,7 +9,7 @@ import { BondType } from '../../../model/types';
 import { IntraUnitBonds } from './data';
 import { Unit } from '../../unit';
 import { IntAdjacencyGraph } from '../../../../../mol-math/graph';
-import { BondComputationProps, getElementIdx, MetalsSet, getElementThreshold, isHydrogen, getElementPairThreshold, DefaultBondComputationProps } from './common';
+import { BondComputationProps, getElementIdx, MetalsSet, getElementThreshold, isHydrogen, DefaultBondComputationProps, getPairingThreshold } from './common';
 import { SortedArray } from '../../../../../mol-data/int';
 import { getIntraBondOrderFromTable } from '../../../model/properties/atomic/bonds';
 import { StructureElement } from '../../element';
@@ -62,7 +62,8 @@ function findIndexPairBonds(unit: Unit.Atomic) {
 
     for (let _aI = 0 as StructureElement.UnitIndex; _aI < atomCount; _aI++) {
         const aI = atoms[_aI];
-        const isHa = type_symbol.value(aI) === 'H';
+        const aeI = getElementIdx(type_symbol.value(aI));
+        const isHa = isHydrogen(aeI);
 
         const srcA = sourceIndex.value(aI);
 
@@ -72,11 +73,30 @@ function findIndexPairBonds(unit: Unit.Atomic) {
 
             const _bI = SortedArray.indexOf(unit.elements, bI) as StructureElement.UnitIndex;
             if (_bI < 0) continue;
-            if (isHa && type_symbol.value(bI) === 'H') continue;
+
+            const beI = getElementIdx(type_symbol.value(bI));
 
             const d = distance[i];
             const dist = getDistance(unit, aI, bI);
-            if ((d !== -1 && equalEps(dist, d, 0.5)) || dist < maxDistance) {
+
+            let add = false;
+            if (d >= 0) {
+                add = equalEps(dist, d, 0.3);
+            } else if (maxDistance >= 0) {
+                add = dist < maxDistance;
+            } else {
+                const pairingThreshold = getPairingThreshold(
+                    aeI, beI, getElementThreshold(aeI), getElementThreshold(beI)
+                );
+                add = dist < pairingThreshold;
+
+                if (isHa && isHydrogen(beI)) {
+                    // TODO handle molecular hydrogen
+                    add = false;
+                }
+            }
+
+            if (add) {
                 atomA[atomA.length] = _aI;
                 atomB[atomB.length] = _bI;
                 orders[orders.length] = order[i];
@@ -214,13 +234,7 @@ function findBonds(unit: Unit.Atomic, props: BondComputationProps): IntraUnitBon
             const dist = Math.sqrt(squaredDistances[ni]);
             if (dist === 0) continue;
 
-            const thresholdAB = getElementPairThreshold(aeI, beI);
-            const pairingThreshold = thresholdAB > 0
-                ? thresholdAB
-                : beI < 0
-                    ? thresholdA
-                    : (thresholdA + getElementThreshold(beI)) / 1.95; // not sure if avg or min but max is too big
-
+            const pairingThreshold = getPairingThreshold(aeI, beI, thresholdA, getElementThreshold(beI));
             if (dist <= pairingThreshold) {
                 atomA[atomA.length] = _aI;
                 atomB[atomB.length] = _bI;