improved accessible-surface-area computation

src/mol-model-props/computed/accessible-surface-area.ts

@@ -7,7 +7,7 @@
 import { ParamDefinition as PD } from '../../mol-util/param-definition'
 import { ShrakeRupleyComputationParams, AccessibleSurfaceArea } from './accessible-surface-area/shrake-rupley';
 import { Structure, CustomPropertyDescriptor } from '../../mol-model/structure';
-import { Task } from '../../mol-task';
+import { Task, RuntimeContext } from '../../mol-task';
 import { idFactory } from '../../mol-util/id-factory';
 
 
@@ -27,9 +27,13 @@ export namespace ComputedAccessibleSurfaceArea {
     }
 
     export function createAttachTask(params: Partial<AccessibleSurfaceAreaComputationProps> = {}) {
-        return (structure: Structure) => Task.create('Compute Accessible Surface Area', async ctx => {
-            if (get(structure)) return true;
-            return await attachFromCifOrCompute(structure, params)
+        return (structure: Structure) => attachTask(structure, params)
+    }
+
+    export function attachTask(structure: Structure, params: Partial<AccessibleSurfaceAreaComputationProps> = {}) {
+        return Task.create('Compute Accessible Surface Area', async ctx => {
+            if (get(structure)) return;
+            return await attachFromCifOrCompute(ctx, structure, params)
         });
     }
 
@@ -39,14 +43,13 @@ export namespace ComputedAccessibleSurfaceArea {
         // TODO `cifExport` and `symbol`
     });
 
-    export async function attachFromCifOrCompute(structure: Structure, params: Partial<AccessibleSurfaceAreaComputationProps> = {}) {
-        if (structure.customPropertyDescriptors.has(Descriptor)) return true;
+    export async function attachFromCifOrCompute(ctx: RuntimeContext, structure: Structure, params: Partial<AccessibleSurfaceAreaComputationProps> = {}) {
+        if (structure.customPropertyDescriptors.has(Descriptor)) return;
 
-        const compAccessibleSurfaceArea = await computeAccessibleSurfaceArea(structure, params)
+        const compAccessibleSurfaceArea = await computeAccessibleSurfaceArea(ctx, structure, params)
 
         structure.customPropertyDescriptors.add(Descriptor);
         set(structure, compAccessibleSurfaceArea);
-        return true;
     }
 }
 
@@ -56,9 +59,9 @@ export const AccessibleSurfaceAreaComputationParams = {
 export type AccessibleSurfaceAreaComputationParams = typeof AccessibleSurfaceAreaComputationParams
 export type AccessibleSurfaceAreaComputationProps = PD.Values<AccessibleSurfaceAreaComputationParams>
 
-async function computeAccessibleSurfaceArea(structure: Structure, params: Partial<AccessibleSurfaceAreaComputationProps>): Promise<ComputedAccessibleSurfaceArea.Property> {
+async function computeAccessibleSurfaceArea(ctx: RuntimeContext, structure: Structure, params: Partial<AccessibleSurfaceAreaComputationProps>): Promise<ComputedAccessibleSurfaceArea.Property> {
     const p = { ...PD.getDefaultValues(AccessibleSurfaceAreaComputationParams), params };
 
-    const accessibleSurfaceArea = await AccessibleSurfaceArea.compute(structure, p);
-    return { id: nextAccessibleSurfaceAreaId(), asa: accessibleSurfaceArea }
+    const asa = await AccessibleSurfaceArea.compute(structure, p).runInContext(ctx);
+    return { id: nextAccessibleSurfaceAreaId(), asa }
 }

src/mol-model-props/computed/accessible-surface-area/shrake-rupley.ts

@@ -2,68 +2,69 @@
  * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author Sebastian Bittrich <sebastian.bittrich@rcsb.org>
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
  */
 
 import { Task, RuntimeContext } from '../../../mol-task';
-import { BitFlags } from '../../../mol-util';
+// import { BitFlags } from '../../../mol-util';
 import { ParamDefinition as PD } from '../../../mol-util/param-definition'
 import { Vec3 } from '../../../mol-math/linear-algebra';
 import { Structure } from '../../../mol-model/structure';
 import { assignRadiusForHeavyAtoms } from './shrake-rupley/radii';
-import { ShrakeRupleyContext, VdWLookup } from './shrake-rupley/common';
-import { computePerResidue } from './shrake-rupley/per-residue';
-import { normalizeAccessibleSurfaceArea } from './shrake-rupley/normalize';
+import { ShrakeRupleyContext, VdWLookup, MaxAsa, DefaultMaxAsa } from './shrake-rupley/common';
+import { computeArea } from './shrake-rupley/area';
 
 export const ShrakeRupleyComputationParams = {
-    numberOfSpherePoints: PD.Numeric(92, {}, { description: 'number of sphere points to sample per atom: 92 (original paper), 960 (BioJava), 3000 (EPPIC) - see Shrake A, Rupley JA: Environment and exposure to solvent of protein atoms. Lysozyme and insulin. J Mol Biol 1973.' }),
-    probeSize: PD.Numeric(1.4, {}, { description: 'corresponds to the size of a water molecule: 1.4 (original paper), 1.5 (occassionally used)' }),
-    buriedRasaThreshold: PD.Numeric(0.16, { min: 0.0, max: 1.0 }, { description: 'below this cutoff of relative accessible surface area a residue will be considered buried - see: Rost B, Sander C: Conservation and prediction of solvent accessibility in protein families. Proteins 1994.' }),
-    nonPolymer: PD.Boolean(false, { description: 'Include non-polymer atoms in computation.' })
+    numberOfSpherePoints: PD.Numeric(92, { min: 12, max: 120, step: 1 }, { description: 'number of sphere points to sample per atom: 92 (original paper), 960 (BioJava), 3000 (EPPIC) - see Shrake A, Rupley JA: Environment and exposure to solvent of protein atoms. Lysozyme and insulin. J Mol Biol 1973.' }),
+    probeSize: PD.Numeric(1.4, { min: 0.1, max: 4, step: 0.01 }, { description: 'corresponds to the size of a water molecule: 1.4 (original paper), 1.5 (occassionally used)' }),
+    // buriedRasaThreshold: PD.Numeric(0.16, { min: 0.0, max: 1.0 }, { description: 'below this cutoff of relative accessible surface area a residue will be considered buried - see: Rost B, Sander C: Conservation and prediction of solvent accessibility in protein families. Proteins 1994.' }),
+    nonPolymer: PD.Boolean(false, { description: 'Include non-polymer atoms as occluders.' })
 }
 export type ShrakeRupleyComputationParams = typeof ShrakeRupleyComputationParams
+export type ShrakeRupleyComputationProps = PD.Values<ShrakeRupleyComputationParams>
+
+// TODO
+// - add back buried and relative asa
 
 namespace AccessibleSurfaceArea {
     /**
      * Adapts the BioJava implementation by Jose Duarte. That implementation is based on the publication by Shrake, A., and
      * J. A. Rupley. "Environment and Exposure to Solvent of Protein Atoms. Lysozyme and Insulin." JMB (1973).
      */
-    export function compute(structure: Structure,
-        params: Partial<PD.Values<ShrakeRupleyComputationParams>> = {}) {
-        params = { ...PD.getDefaultValues(ShrakeRupleyComputationParams), ...params };
-        return Task.create('Compute Accessible Surface Area', async rtctx => {
-            return await _compute(rtctx, structure, params);
-        }).run();
+    export function compute(structure: Structure, props: Partial<ShrakeRupleyComputationProps> = {}) {
+        const p = { ...PD.getDefaultValues(ShrakeRupleyComputationParams), ...props };
+        return Task.create('Compute Accessible Surface Area', async runtime => {
+            return await calculate(runtime, structure, p);
+        });
     }
 
-    async function _compute(rtctx: RuntimeContext, structure: Structure, params: Partial<PD.Values<ShrakeRupleyComputationParams>> = {}): Promise<AccessibleSurfaceArea> {
-        const ctx = initialize(rtctx, structure, params);
+    async function calculate(runtime: RuntimeContext, structure: Structure, props: ShrakeRupleyComputationProps): Promise<AccessibleSurfaceArea> {
+        const ctx = initialize(structure, props);
 
         assignRadiusForHeavyAtoms(ctx);
-        computePerResidue(ctx);
-        normalizeAccessibleSurfaceArea(ctx);
+        await computeArea(runtime, ctx);
 
+        const { accessibleSurfaceArea, serialResidueIndex } = ctx
         return {
-            accessibleSurfaceArea: ctx.accessibleSurfaceArea,
-            relativeAccessibleSurfaceArea: ctx.relativeAccessibleSurfaceArea,
-            buried: (index: number) => ctx.relativeAccessibleSurfaceArea[index] < 0.16
+            serialResidueIndex,
+            accessibleSurfaceArea
         };
     }
 
-    function initialize(rtctx: RuntimeContext, structure: Structure, params: Partial<PD.Values<ShrakeRupleyComputationParams>>): ShrakeRupleyContext {
-        const { elementCount, polymerResidueCount } = structure;
+    function initialize(structure: Structure, props: ShrakeRupleyComputationProps): ShrakeRupleyContext {
+        const { elementCount, atomicResidueCount } = structure;
+        const { probeSize, nonPolymer, numberOfSpherePoints } = props
 
         return {
-            rtctx: rtctx,
-            structure: structure,
-            probeSize: params.probeSize!,
-            nonPolymer: params.nonPolymer!,
-            spherePoints: generateSpherePoints(params.numberOfSpherePoints!),
-            cons: 4.0 * Math.PI / params.numberOfSpherePoints!,
-            maxLookupRadius: 2 * params.probeSize! + 2 * VdWLookup[2], // 2x probe size + 2x largest VdW
-            atomRadius: new Int8Array(elementCount),
-            accessibleSurfaceArea: new Float32Array(polymerResidueCount),
-            relativeAccessibleSurfaceArea: new Float32Array(polymerResidueCount),
-            updateChunk: 25000
+            structure,
+            probeSize,
+            nonPolymer,
+            spherePoints: generateSpherePoints(numberOfSpherePoints!),
+            scalingConstant: 4.0 * Math.PI / numberOfSpherePoints!,
+            maxLookupRadius: 2 * props.probeSize + 2 * VdWLookup[2], // 2x probe size + 2x largest VdW
+            atomRadiusType: new Int8Array(elementCount),
+            serialResidueIndex: new Int32Array(elementCount),
+            accessibleSurfaceArea: new Float32Array(atomicResidueCount)
         }
     }
 
@@ -76,26 +77,31 @@ namespace AccessibleSurfaceArea {
             const y = k * offset - 1.0 + (offset / 2.0);
             const r = Math.sqrt(1.0 - y * y);
             const phi = k * inc;
-            points[points.length] = [Math.cos(phi) * r, y, Math.sin(phi) * r] as Vec3;
+            points[points.length] = Vec3.create(Math.cos(phi) * r, y, Math.sin(phi) * r);
         }
         return points;
     }
 
-    export namespace SolventAccessibility {
-        export const is: (t: number, f: Flag) => boolean = BitFlags.has
-        export const create: (f: Flag) => number = BitFlags.create
-        export const enum Flag {
-            _ = 0x0,
-            BURIED = 0x1,
-            ACCESSIBLE = 0x2
-        }
+    // export namespace SolventAccessibility {
+    //     export const is: (t: number, f: Flag) => boolean = BitFlags.has
+    //     export const create: (f: Flag) => number = BitFlags.create
+    //     export const enum Flag {
+    //         _ = 0x0,
+    //         BURIED = 0x1,
+    //         ACCESSIBLE = 0x2
+    //     }
+    // }
+
+    /** Get relative area for a given component id */
+    export function normalize(compId: string, asa: number) {
+        const maxAsa = MaxAsa[compId] || DefaultMaxAsa;
+        return asa / maxAsa
     }
 }
 
 interface AccessibleSurfaceArea {
+    readonly serialResidueIndex: ArrayLike<number>
     readonly accessibleSurfaceArea: ArrayLike<number>
-    readonly relativeAccessibleSurfaceArea: ArrayLike<number>
-    buried(index: number): boolean
 }
 
 export { AccessibleSurfaceArea }

src/mol-model-props/computed/accessible-surface-area/shrake-rupley/area.ts

@@ -0,0 +1,103 @@
+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Sebastian Bittrich <sebastian.bittrich@rcsb.org>
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import { ShrakeRupleyContext, VdWLookup } from './common';
+import { Vec3 } from '../../../../mol-math/linear-algebra';
+import { RuntimeContext } from '../../../../mol-task';
+import { StructureProperties, StructureElement, Structure } from '../../../../mol-model/structure/structure';
+
+// TODO
+// - iterate over units and elements
+// - avoid using serial-element index whenever possible
+// - calculate atomRadiusType only for invariant units
+// - factor serialResidueIndex out
+
+const updateChunk = 5000
+export async function computeArea(runtime: RuntimeContext, ctx: ShrakeRupleyContext) {
+    const { atomRadiusType: atomRadius } = ctx;
+    for (let i = 0; i < atomRadius.length; i += updateChunk) {
+        if (runtime.shouldUpdate) {
+            await runtime.update({ message: 'Computing per residue surface accessibility...', current: i, max: atomRadius.length });
+        }
+
+        computeRange(ctx, i, Math.min(i + updateChunk, atomRadius.length));
+    }
+}
+
+const aPos = Vec3();
+const bPos = Vec3();
+const testPoint = Vec3();
+const aLoc = StructureElement.Location.create()
+const bLoc = StructureElement.Location.create()
+
+function setLocation(l: StructureElement.Location, structure: Structure, serialIndex: number) {
+    l.unit = structure.units[structure.serialMapping.unitIndices[serialIndex]]
+    l.element = structure.serialMapping.elementIndices[serialIndex]
+    return l
+}
+
+function computeRange(ctx: ShrakeRupleyContext, begin: number, end: number) {
+    const { structure, atomRadiusType, serialResidueIndex, accessibleSurfaceArea, spherePoints, scalingConstant, maxLookupRadius, probeSize } = ctx;
+    const { x, y, z } = StructureProperties.atom
+    const { lookup3d, serialMapping, unitIndexMap } = structure;
+    const { cumulativeUnitElementCount } = serialMapping
+
+    for (let aI = begin; aI < end; ++aI) {
+        const radius1 = VdWLookup[atomRadiusType[aI]];
+        if (radius1 === VdWLookup[0]) continue;
+
+        setLocation(aLoc, structure, aI)
+        const aX = x(aLoc)
+        const aY = y(aLoc)
+        const aZ = z(aLoc)
+
+        // pre-filter by lookup3d (provides >10x speed-up compared to naive evaluation)
+        const { count, units, indices, squaredDistances } = lookup3d.find(aX, aY, aZ, maxLookupRadius);
+
+        // collect neighbors for each atom
+        const cutoff1 = probeSize + probeSize + radius1;
+        const neighbors = []; // TODO reuse
+        for (let iI = 0; iI < count; ++iI) {
+            const bUnit = units[iI]
+            StructureElement.Location.set(bLoc, bUnit, bUnit.elements[indices[iI]])
+            const bI = cumulativeUnitElementCount[unitIndexMap.get(bUnit.id)] + indices[iI]
+
+            const radius2 = VdWLookup[atomRadiusType[bI]];
+            if (StructureElement.Location.areEqual(aLoc, bLoc) || radius2 === VdWLookup[0]) continue;
+
+            const cutoff2 = (cutoff1 + radius2) * (cutoff1 + radius2);
+            if (squaredDistances[iI] < cutoff2) {
+                neighbors[neighbors.length] = bI;
+            }
+        }
+
+        // for all neighbors: test all sphere points
+        Vec3.set(aPos, aX, aY, aZ)
+        const scale = probeSize + radius1;
+        let accessiblePointCount = 0;
+        for (let sI = 0; sI < spherePoints.length; ++sI) {
+            Vec3.scaleAndAdd(testPoint, aPos, spherePoints[sI], scale)
+            let accessible = true;
+
+            for (let _nI = 0; _nI < neighbors.length; ++_nI) {
+                const nI = neighbors[_nI];
+                setLocation(bLoc, structure, nI)
+                Vec3.set(bPos, x(bLoc), y(bLoc), z(bLoc))
+                const radius3 = VdWLookup[atomRadiusType[nI]];
+                const cutoff3 = (radius3 + probeSize) * (radius3 + probeSize);
+                if (Vec3.squaredDistance(testPoint, bPos) < cutoff3) {
+                    accessible = false;
+                    break;
+                }
+            }
+
+            if (accessible) ++accessiblePointCount;
+        }
+
+        accessibleSurfaceArea[serialResidueIndex[aI]] += scalingConstant * accessiblePointCount * scale * scale;
+    }
+}

src/mol-model-props/computed/accessible-surface-area/shrake-rupley/common.ts

@@ -4,25 +4,22 @@
  * @author Sebastian Bittrich <sebastian.bittrich@rcsb.org>
  */
 
-import { RuntimeContext } from '../../../../mol-task';
 import { Structure } from '../../../../mol-model/structure';
 import { Vec3 } from '../../../../mol-math/linear-algebra';
 
 export interface ShrakeRupleyContext {
-    rtctx: RuntimeContext,
     structure: Structure,
     spherePoints: Vec3[],
     probeSize: number,
     nonPolymer: boolean,
-    cons: number,
+    scalingConstant: number,
     maxLookupRadius: number,
-    atomRadius: Int8Array,
-    accessibleSurfaceArea: Float32Array,
-    relativeAccessibleSurfaceArea: Float32Array,
-    updateChunk: number
+    atomRadiusType: Int8Array,
+    serialResidueIndex: Int32Array,
+    accessibleSurfaceArea: Float32Array
 }
 
-// Chothia's amino acid and nucleotide atom vdw radii
+/** Chothia's amino acid and nucleotide atom vdw radii */
 export const VdWLookup = [
     -1.0, // 0: missing
     1.76, // 1: trigonal C
@@ -36,9 +33,8 @@ export const VdWLookup = [
     1.40  // 9: P (nucleic)
 ] // can still be appended on-the-fly for rare elements like selenium
 
-
 /** Maximum accessible surface area observed for amino acids. Taken from: http://dx.doi.org/10.1371/journal.pone.0080635 */
-export const MaxAsa = {
+export const MaxAsa: { [k: string]: number } = {
     'ALA': 121.0,
     'ARG': 265.0,
     'ASN': 187.0,

src/mol-model-props/computed/accessible-surface-area/shrake-rupley/normalize.ts

@@ -1,32 +0,0 @@
-/**
- * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
- *
- * @author Sebastian Bittrich <sebastian.bittrich@rcsb.org>
- */
-
-import { ShrakeRupleyContext } from './common';
-import { isPolymer, MaxAsa, DefaultMaxAsa } from '../../../../mol-model/structure/model/types';
-
-export async function normalizeAccessibleSurfaceArea(ctx: ShrakeRupleyContext) {
-    const updateChunk = ctx.updateChunk / 10;
-    const residueCount = ctx.relativeAccessibleSurfaceArea.length;
-    for (let i = 0; i < residueCount; i += updateChunk) {
-        computeRange(ctx, i, Math.min(i + updateChunk, residueCount));
-    }
-}
-
-function computeRange(ctx: ShrakeRupleyContext, begin: number, end: number) {
-    const { accessibleSurfaceArea, relativeAccessibleSurfaceArea, structure } = ctx;
-    const { residues, derived } = structure.model.atomicHierarchy;
-
-    for (let i = begin; i < end; ++i) {
-        // skip entities not part of a polymer chain
-        if (!ctx.nonPolymer) {
-            if (!isPolymer(derived.residue.moleculeType[i])) continue;
-        }
-
-        const maxAsa = (MaxAsa as any)[residues.label_comp_id.value(i)];
-        const rasa = accessibleSurfaceArea[i] / (maxAsa === undefined ? DefaultMaxAsa : maxAsa);
-        relativeAccessibleSurfaceArea[i] = rasa;
-    }
-}

src/mol-model-props/computed/accessible-surface-area/shrake-rupley/per-residue.ts

@@ -1,86 +0,0 @@
-/**
- * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
- *
- * @author Sebastian Bittrich <sebastian.bittrich@rcsb.org>
- */
-
-import { ShrakeRupleyContext, VdWLookup } from './common';
-import { Vec3 } from '../../../../mol-math/linear-algebra';
-
-export async function computePerResidue(ctx: ShrakeRupleyContext) {
-    const { rtctx, updateChunk, atomRadius } = ctx;
-    for (let i = 0; i < atomRadius.length; i += updateChunk) {
-        computeRange(ctx, i, Math.min(i + updateChunk, atomRadius.length));
-
-        if (rtctx.shouldUpdate) {
-            rtctx.update({ message: 'Computing per residue surface accessibility...', current: i, max: atomRadius.length });
-        }
-    }
-    // async yields 4x speed-up
-}
-
-const aPos = Vec3.zero();
-const bPos = Vec3.zero();
-let testPoint = Vec3.zero();
-function computeRange(ctx: ShrakeRupleyContext, begin: number, end: number) {
-    const { structure, atomRadius, accessibleSurfaceArea, spherePoints, cons, maxLookupRadius, probeSize } = ctx;
-        const { model } = structure;
-        const { x, y, z } = model.atomicConformation;
-        const { residueAtomSegments } = model.atomicHierarchy;
-        const { lookup3d } = structure;
-
-    const position = (i: number, v: Vec3) => Vec3.set(v, x[i], y[i], z[i]);
-
-    // console.log(`computing ASA chunk ${ begin } to ${ end }`);
-    for (let aI = begin; aI < end; aI++) {
-        const radius1 = VdWLookup[atomRadius[aI]];
-        if (radius1 === VdWLookup[0]) continue;
-
-        // pre-filter by lookup3d
-        // lookup provides >10x speed-up compared to naive evaluation
-        const { count, units, indices, squaredDistances } = lookup3d.find(x[aI], y[aI], z[aI], maxLookupRadius);
-        // we could keep track of all found neighbors of each atom, however slow and memory intensive
-
-        // collect neighbors for each atom
-        const cutoff1 = probeSize + probeSize + radius1;
-        const neighbors = [];
-        for (let iI = 0; iI < count; ++iI) {
-            const bI = units[iI].elements[indices[iI]];
-
-            const radius2 = VdWLookup[atomRadius[bI]];
-            if (aI === bI || radius2 === VdWLookup[0]) continue;
-
-            const cutoff2 = (cutoff1 + radius2) * (cutoff1 + radius2);
-            if (squaredDistances[iI] < cutoff2) {
-                neighbors[neighbors.length] = bI;
-            }
-        }
-
-        // for all neighbors: test all sphere points
-        position(aI, aPos);
-        const scalar = probeSize + radius1;
-        let accessiblePointCount = 0;
-        for (let sI = 0; sI < spherePoints.length; ++sI) {
-            const spherePoint = spherePoints[sI];
-            testPoint = [spherePoint[0] * scalar + aPos[0],
-                spherePoint[1] * scalar + aPos[1],
-                spherePoint[2] * scalar + aPos[2]] as Vec3;
-            let accessible = true;
-
-            for (let _nI = 0; _nI < neighbors.length; ++_nI) {
-                const nI = neighbors[_nI];
-                position(nI, bPos);
-                const radius3 = VdWLookup[atomRadius[nI]];
-                const cutoff3 = (radius3 + probeSize) * (radius3 + probeSize);
-                if (Vec3.squaredDistance(testPoint, bPos) < cutoff3) {
-                    accessible = false;
-                    break;
-                }
-            }
-
-            if (accessible) ++accessiblePointCount;
-        }
-
-        accessibleSurfaceArea[residueAtomSegments.index[aI]] += cons * accessiblePointCount * scalar * scalar;
-    }
-}

src/mol-model-props/computed/accessible-surface-area/shrake-rupley/radii.ts

@@ -2,61 +2,75 @@
  * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author Sebastian Bittrich <sebastian.bittrich@rcsb.org>
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
  */
 
 import { ShrakeRupleyContext, VdWLookup } from './common';
 import { getElementIdx, isHydrogen } from '../../../../mol-model/structure/structure/unit/links/common';
 import { isPolymer, isNucleic, MoleculeType, ElementSymbol } from '../../../../mol-model/structure/model/types';
 import { VdwRadius } from '../../../../mol-model/structure/model/properties/atomic';
+import { StructureElement, StructureProperties } from '../../../../mol-model/structure/structure';
+import { getElementMoleculeType } from '../../../../mol-model/structure/util';
 
-export async function assignRadiusForHeavyAtoms(ctx: ShrakeRupleyContext) {
-    const { updateChunk, atomRadius } = ctx;
-    for (let i = 0; i < atomRadius.length; i += updateChunk) {
-        computeRange(ctx, i, Math.min(i + updateChunk, atomRadius.length));
-    }
-}
+const l = StructureElement.Location.create()
 
-function computeRange(ctx: ShrakeRupleyContext, begin: number, end: number) {
-    const { structure } = ctx;
-    const { model } = structure;
-    const { atoms: atomInfo, derived, residues, residueAtomSegments } = model.atomicHierarchy;
-    const { label_comp_id } = residues;
-    const { moleculeType } = derived.residue;
-    const { type_symbol, label_atom_id } = atomInfo;
-
-    for (let aI = begin; aI < end; ++aI) {
-        const rI = residueAtomSegments.index[aI];
-        const element = type_symbol.value(aI);
-        const elementIdx = getElementIdx(element);
-        // skip hydrogen atoms
-        if (isHydrogen(elementIdx)) {
-            ctx.atomRadius[aI] = VdWLookup[0];
-            continue;
-        }
+export function assignRadiusForHeavyAtoms(ctx: ShrakeRupleyContext) {
+    const { label_comp_id, key } = StructureProperties.residue
+    const { type_symbol, label_atom_id } = StructureProperties.atom
+    const { structure, atomRadiusType, serialResidueIndex } = ctx;
+
+    let prevResidueIdx = 0
+    let residueIdx = 0
+    let serialResidueIdx = -1
+
+    for (let i = 0, m = 0, il = structure.units.length; i < il; ++i) {
+        const unit = structure.units[i]
+        const { elements } = unit
+        l.unit = unit
+
+        prevResidueIdx = -1
+
+        for (let j = 0, jl = elements.length; j < jl; ++j) {
+            const eI = elements[j]
+            const mj = m + j
+
+            l.element = eI
+            residueIdx = key(l)
+
+            if (prevResidueIdx !== residueIdx) ++serialResidueIdx
+            prevResidueIdx = residueIdx
 
-        const residueType = moleculeType[rI];
-        // skip non-polymer groups
-        if (!ctx.nonPolymer) {
-            if (!isPolymer(residueType)) {
-                ctx.atomRadius[aI] = VdWLookup[0];
+            const element = type_symbol(l);
+            const elementIdx = getElementIdx(element);
+
+            // skip hydrogen atoms
+            if (isHydrogen(elementIdx)) {
+                atomRadiusType[mj] = VdWLookup[0];
+                serialResidueIndex[mj] = -1
                 continue;
             }
-        }
 
-        const atomId = label_atom_id.value(aI);
-        let compId = label_comp_id.value(rI);
+            const residueType = getElementMoleculeType(unit, eI)
+            // skip non-polymer groups
+            if (!ctx.nonPolymer && !isPolymer(residueType)) {
+                atomRadiusType[mj] = VdWLookup[0];
+                serialResidueIndex[mj] = -1
+                continue;
+            }
 
-        // handle modified residues
-        const parentId = model.properties.modifiedResidues.parentId.get(compId);
-        if (parentId !== void 0) compId = parentId;
+            const atomId = label_atom_id(l);
+            const compId = label_comp_id(l);
 
-        if (isNucleic(residueType)) {
-             ctx.atomRadius[aI] = determineRadiusNucl(atomId, element, compId);
-        } else if (residueType === MoleculeType.Protein) {
-            ctx.atomRadius[aI] = determineRadiusAmino(atomId, element, compId);
-        } else {
-            ctx.atomRadius[aI] = handleNonStandardCase(element);
+            if (isNucleic(residueType)) {
+                atomRadiusType[mj] = determineRadiusNucl(atomId, element, compId);
+            } else if (residueType === MoleculeType.Protein) {
+                atomRadiusType[mj] = determineRadiusAmino(atomId, element, compId);
+            } else {
+                atomRadiusType[mj] = handleNonStandardCase(element);
+            }
+            serialResidueIndex[mj] = serialResidueIdx
         }
+        m += elements.length
     }
 }
 
@@ -68,41 +82,37 @@ function computeRange(ctx: ShrakeRupleyContext, begin: number, end: number) {
 function determineRadiusAmino(atomId: string, element: ElementSymbol, compId: string): number {
     switch (element) {
         case 'O':
-        return 5;
+            return 5;
         case 'S':
-        return 6;
+            return 6;
         case 'N':
-        return atomId === 'NZ' ? 4 : 3;
+            return atomId === 'NZ' ? 4 : 3;
         case 'C':
-        switch (atomId) {
-            case 'C': case 'CE1': case'CE2': case 'CE3': case 'CH2': case 'CZ': case 'CZ2': case 'CZ3':
-            return 1;
-            case 'CA': case 'CB': case 'CE': case 'CG1': case 'CG2':
-            return 2;
-            default:
-            switch (compId) {
-                case 'PHE': case 'TRP': case 'TYR': case 'HIS': case 'ASP': case 'ASN':
-                return 1;
-                case 'PRO': case 'LYS': case 'ARG': case 'MET': case 'ILE': case 'LEU':
-                return 2;
-                case 'GLU': case 'GLN':
-                return atomId === 'CD' ? 1 : 2;
+            switch (atomId) {
+                case 'C': case 'CE1': case'CE2': case 'CE3': case 'CH2': case 'CZ': case 'CZ2': case 'CZ3':
+                    return 1;
+                case 'CA': case 'CB': case 'CE': case 'CG1': case 'CG2':
+                    return 2;
+                default:
+                    switch (compId) {
+                        case 'PHE': case 'TRP': case 'TYR': case 'HIS': case 'ASP': case 'ASN':
+                            return 1;
+                        case 'PRO': case 'LYS': case 'ARG': case 'MET': case 'ILE': case 'LEU':
+                            return 2;
+                        case 'GLU': case 'GLN':
+                            return atomId === 'CD' ? 1 : 2;
+                    }
             }
-        }
     }
     return handleNonStandardCase(element);
 }
 
 function determineRadiusNucl(atomId: string, element: ElementSymbol, compId: string): number {
     switch (element) {
-        case 'C':
-        return 7;
-        case 'N':
-        return 8;
-        case 'P':
-        return 9;
-        case 'O':
-        return 5;
+        case 'C': return 7;
+        case 'N': return 8;
+        case 'P': return 9;
+        case 'O': return 5;
     }
     return handleNonStandardCase(element);
 }

src/mol-theme/color.ts

@@ -33,6 +33,7 @@ import { ModelIndexColorThemeProvider } from './color/model-index';
 import { OccupancyColorThemeProvider } from './color/occupancy';
 import { OperatorNameColorThemeProvider } from './color/operator-name';
 import { OperatorHklColorThemeProvider } from './color/operator-hkl';
+import { AccessibleSurfaceAreaColorThemeProvider } from './color/accessible-surface-area';
 
 export type LocationColor = (location: Location, isSecondary: boolean) => Color
 
@@ -74,6 +75,7 @@ namespace ColorTheme {
 }
 
 export const BuiltInColorThemes = {
+    'accessible-surface-area': AccessibleSurfaceAreaColorThemeProvider,
     'carbohydrate-symbol': CarbohydrateSymbolColorThemeProvider,
     'chain-id': ChainIdColorThemeProvider,
     'cross-link': CrossLinkColorThemeProvider,

src/mol-theme/color/accessible-surface-area.ts

@@ -15,7 +15,7 @@ import { ComputedAccessibleSurfaceArea } from '../../mol-model-props/computed/ac
 import { ColorListName, ColorListOptionsScale } from '../../mol-util/color/lists';
 import { VdWLookup } from '../../mol-model-props/computed/accessible-surface-area/shrake-rupley/common';
 
-const DefaultColor = Color(0xFFFFFF)
+const DefaultColor = Color(0xFAFAFA)
 const Description = 'Assigns a color based on the relative accessible surface area of a residue.'
 
 export const AccessibleSurfaceAreaColorThemeParams = {
@@ -35,17 +35,20 @@ export function AccessibleSurfaceAreaColorTheme(ctx: ThemeDataContext, props: PD
         domain: [0.0, 1.0]
     })
 
-    const accessibleSurfaceArea = ctx.structure ? ComputedAccessibleSurfaceArea.get(ctx.structure)!.asa : undefined
+    const accessibleSurfaceArea = ctx.structure ? ComputedAccessibleSurfaceArea.get(ctx.structure.root)?.asa : undefined
+
+    if (accessibleSurfaceArea && ctx.structure) {
+        const { structure } = ctx
+        const { getSerialIndex } = structure.root.serialMapping
+        const { relativeAccessibleSurfaceArea, serialResidueIndex } = accessibleSurfaceArea
 
-    if (accessibleSurfaceArea) {
         color = (location: Location): Color => {
             if (StructureElement.Location.is(location)) {
                 if (Unit.isAtomic(location.unit)) {
-                    const value = accessibleSurfaceArea.relativeAccessibleSurfaceArea[location.unit.residueIndex[location.element]];
-                    return value !== VdWLookup[0] /* signals missing value */ ? scale.color(value) : DefaultColor;
+                    const rSI = serialResidueIndex[getSerialIndex(location.unit, location.element)]
+                    return rSI === -1 ? DefaultColor : scale.color(relativeAccessibleSurfaceArea[rSI])
                 }
             }
-
             return DefaultColor
         }
     } else {
@@ -67,5 +70,5 @@ export const AccessibleSurfaceAreaColorThemeProvider: ColorTheme.Provider<Access
     factory: AccessibleSurfaceAreaColorTheme,
     getParams: getAccessibleSurfaceAreaColorThemeParams,
     defaultValues: PD.getDefaultValues(AccessibleSurfaceAreaColorThemeParams),
-    isApplicable: (ctx: ThemeDataContext) => !!ctx.structure
+    isApplicable: (ctx: ThemeDataContext) => !!ctx.structure,
 }

src/tests/browser/render-asa.ts

@@ -1,150 +0,0 @@
-/**
- * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
- *
- * @author Alexander Rose <alexander.rose@weirdbyte.de>
- */
-
-import './index.html'
-import { Canvas3D } from '../../mol-canvas3d/canvas3d';
-import { CifFrame, CIF } from '../../mol-io/reader/cif'
-import { Model, Structure, StructureElement, Unit } from '../../mol-model/structure';
-import { ColorTheme, LocationColor } from '../../mol-theme/color';
-import { SizeTheme } from '../../mol-theme/size';
-import { CartoonRepresentationProvider } from '../../mol-repr/structure/representation/cartoon';
-import { trajectoryFromMmCIF } from '../../mol-model-formats/structure/mmcif';
-import { Color, ColorScale } from '../../mol-util/color';
-import { Location } from '../../mol-model/location';
-import { ThemeDataContext } from '../../mol-theme/theme';
-import { ParamDefinition as PD } from '../../mol-util/param-definition';
-import { AccessibleSurfaceArea } from '../../mol-model-props/computed/accessible-surface-area/shrake-rupley';
-import { VdWLookup } from '../../mol-model-props/computed/accessible-surface-area/shrake-rupley/common';
-import { resizeCanvas } from '../../mol-canvas3d/util';
-import { ColorListName, ColorListOptions } from '../../mol-util/color/lists';
-
-const parent = document.getElementById('app')!
-parent.style.width = '100%'
-parent.style.height = '100%'
-
-const canvas = document.createElement('canvas')
-parent.appendChild(canvas)
-resizeCanvas(canvas, parent)
-
-const canvas3d = Canvas3D.fromCanvas(canvas)
-canvas3d.animate()
-
-
-async function parseCif(data: string|Uint8Array) {
-    const comp = CIF.parse(data);
-    const parsed = await comp.run();
-    if (parsed.isError) throw parsed;
-    return parsed.result;
-}
-
-async function downloadCif(url: string, isBinary: boolean) {
-    const data = await fetch(url);
-    return parseCif(isBinary ? new Uint8Array(await data.arrayBuffer()) : await data.text());
-}
-
-async function downloadFromPdb(pdb: string) {
-    // const parsed = await downloadCif(`https://files.rcsb.org/download/${pdb}.cif`, false);
-    const parsed = await downloadCif(`https://webchem.ncbr.muni.cz/ModelServer/static/bcif/${pdb}`, true);
-    return parsed.blocks[0];
-}
-
-async function getModels(frame: CifFrame) {
-    return await trajectoryFromMmCIF(frame).run();
-}
-
-async function getStructure(model: Model) {
-    return Structure.ofModel(model);
-}
-
-const reprCtx = {
-    colorThemeRegistry: ColorTheme.createRegistry(),
-    sizeThemeRegistry: SizeTheme.createRegistry()
-}
-function getCartoonRepr() {
-    return CartoonRepresentationProvider.factory(reprCtx, CartoonRepresentationProvider.getParams)
-}
-
-let accessibleSurfaceArea: AccessibleSurfaceArea;
-async function init(props = {}) {
-    const cif = await downloadFromPdb(
-        // '3j3q'
-        // '1aon'
-        // '1acj'
-        // '1pga'
-        // '1brr'
-        '1hrc'
-    )
-    const models = await getModels(cif)
-    const structure = await getStructure(models[0])
-
-    // async compute ASA
-    console.time('computeASA')
-    accessibleSurfaceArea = await AccessibleSurfaceArea.compute(structure)
-    console.timeEnd('computeASA')
-
-    const cartoonRepr = getCartoonRepr()
-
-    // create color theme
-    cartoonRepr.setTheme({
-        color: AccessibleSurfaceAreaColorTheme(reprCtx, { ...PD.getDefaultValues(AccessibleSurfaceAreaColorThemeParams), ...props }),
-        size: reprCtx.sizeThemeRegistry.create('uniform', { structure })
-    })
-    await cartoonRepr.createOrUpdate({ ...CartoonRepresentationProvider.defaultValues, quality: 'auto' }, structure).run()
-
-    canvas3d.add(cartoonRepr)
-    canvas3d.resetCamera()
-}
-
-init()
-
-const DefaultColor = Color(0xFFFFFF)
-const Description = 'Assigns a color based on the relative accessible surface area of a residue.'
-
-export const AccessibleSurfaceAreaColorThemeParams = {
-    list: PD.ColorList<ColorListName>('rainbow', ColorListOptions)
-}
-export type AccessibleSurfaceAreaColorThemeParams = typeof AccessibleSurfaceAreaColorThemeParams
-export function getAccessibleSurfaceAreaColorThemeParams(ctx: ThemeDataContext) {
-    return AccessibleSurfaceAreaColorThemeParams // TODO return copy
-}
-
-export function AccessibleSurfaceAreaColorTheme(ctx: ThemeDataContext, props: PD.Values<AccessibleSurfaceAreaColorThemeParams>): ColorTheme<AccessibleSurfaceAreaColorThemeParams> {
-    let color: LocationColor = () => DefaultColor
-    const scale = ColorScale.create({
-        listOrName: props.list,
-        minLabel: '0.0 (buried)',
-        maxLabel: '1.0 (exposed)',
-        domain: [0.0, 1.0]
-    })
-
-    color = (location: Location): Color => {
-        if (StructureElement.Location.is(location)) {
-            if (Unit.isAtomic(location.unit)) {
-                const value = accessibleSurfaceArea.relativeAccessibleSurfaceArea[location.unit.residueIndex[location.element]];
-                return value !== VdWLookup[0] /* signals missing value */ ? scale.color(value) : DefaultColor;
-            }
-        }
-
-        return DefaultColor
-    }
-
-    return {
-        factory: AccessibleSurfaceAreaColorTheme,
-        granularity: 'group',
-        color,
-        props,
-        description: Description,
-        legend: scale ? scale.legend : undefined
-    }
-}
-
-export const AccessibleSurfaceAreaColorThemeProvider: ColorTheme.Provider<AccessibleSurfaceAreaColorThemeParams> = {
-    label: 'Accessible Surface Area',
-    factory: AccessibleSurfaceAreaColorTheme,
-    getParams: getAccessibleSurfaceAreaColorThemeParams,
-    defaultValues: PD.getDefaultValues(AccessibleSurfaceAreaColorThemeParams),
-    isApplicable: (ctx: ThemeDataContext) => !!ctx.structure
-}

src/tests/browser/render-structure.ts

@@ -15,8 +15,15 @@ import { trajectoryFromMmCIF } from '../../mol-model-formats/structure/mmcif';
 import { MolecularSurfaceRepresentationProvider } from '../../mol-repr/structure/representation/molecular-surface';
 import { BallAndStickRepresentationProvider } from '../../mol-repr/structure/representation/ball-and-stick';
 import { GaussianSurfaceRepresentationProvider } from '../../mol-repr/structure/representation/gaussian-surface';
-import { ComputedSecondaryStructure } from '../../mol-model-props/computed/secondary-structure';
+// import { ComputedSecondaryStructure } from '../../mol-model-props/computed/secondary-structure';
 import { resizeCanvas } from '../../mol-canvas3d/util';
+import { Representation } from '../../mol-repr/representation';
+import { throttleTime } from 'rxjs/operators';
+import { MarkerAction } from '../../mol-util/marker-action';
+import { EveryLoci } from '../../mol-model/loci';
+import { lociLabel } from '../../mol-theme/label';
+import { InteractionsRepresentationProvider } from '../../mol-repr/structure/representation/interactions';
+import { ComputedInteractions } from '../../mol-model-props/computed/interactions';
 
 const parent = document.getElementById('app')!
 parent.style.width = '100%'
@@ -29,6 +36,33 @@ resizeCanvas(canvas, parent)
 const canvas3d = Canvas3D.fromCanvas(canvas)
 canvas3d.animate()
 
+const info = document.createElement('div')
+info.style.position = 'absolute'
+info.style.fontFamily = 'sans-serif'
+info.style.fontSize = '16pt'
+info.style.bottom = '20px'
+info.style.right = '20px'
+info.style.color = 'white'
+parent.appendChild(info)
+
+let prevReprLoci = Representation.Loci.Empty
+canvas3d.input.move.pipe(throttleTime(100)).subscribe(({x, y}) => {
+    const pickingId = canvas3d.identify(x, y)
+    let label = ''
+    if (pickingId) {
+        const reprLoci = canvas3d.getLoci(pickingId)
+        label = lociLabel(reprLoci.loci)
+        if (!Representation.Loci.areEqual(prevReprLoci, reprLoci)) {
+            canvas3d.mark(prevReprLoci, MarkerAction.RemoveHighlight)
+            canvas3d.mark(reprLoci, MarkerAction.Highlight)
+            prevReprLoci = reprLoci
+        }
+    } else {
+        canvas3d.mark({ loci: EveryLoci }, MarkerAction.RemoveHighlight)
+        prevReprLoci = Representation.Loci.Empty
+    }
+    info.innerHTML = label
+})
 
 async function parseCif(data: string|Uint8Array) {
     const comp = CIF.parse(data);
@@ -43,8 +77,7 @@ async function downloadCif(url: string, isBinary: boolean) {
 }
 
 async function downloadFromPdb(pdb: string) {
-    // const parsed = await downloadCif(`https://files.rcsb.org/download/${pdb}.cif`, false);
-    const parsed = await downloadCif(`https://webchem.ncbr.muni.cz/ModelServer/static/bcif/${pdb}`, true);
+    const parsed = await downloadCif(`https://models.rcsb.org/${pdb}.bcif`, true);
     return parsed.blocks[0];
 }
 
@@ -57,6 +90,7 @@ async function getStructure(model: Model) {
 }
 
 const reprCtx = {
+    webgl: canvas3d.webgl,
     colorThemeRegistry: ColorTheme.createRegistry(),
     sizeThemeRegistry: SizeTheme.createRegistry()
 }
@@ -64,6 +98,10 @@ function getCartoonRepr() {
     return CartoonRepresentationProvider.factory(reprCtx, CartoonRepresentationProvider.getParams)
 }
 
+function getInteractionRepr() {
+    return InteractionsRepresentationProvider.factory(reprCtx, InteractionsRepresentationProvider.getParams)
+}
+
 function getBallAndStickRepr() {
     return BallAndStickRepresentationProvider.factory(reprCtx, BallAndStickRepresentationProvider.getParams)
 }
@@ -80,34 +118,54 @@ async function init() {
     const cif = await downloadFromPdb('1crn')
     const models = await getModels(cif)
     const structure = await getStructure(models[0])
-    console.time('computeDSSP')
-    await ComputedSecondaryStructure.attachFromCifOrCompute(structure)
-    console.timeEnd('computeDSSP');
+    // console.time('computeDSSP')
+    // await ComputedSecondaryStructure.attachFromCifOrCompute(structure)
+    // console.timeEnd('computeDSSP');
+
+    // console.time('computeValenceModel')
+    // await ComputedValenceModel.attachFromCifOrCompute(structure)
+    // console.timeEnd('computeValenceModel');
+    // console.log(ComputedValenceModel.get(structure))
+
+    console.time('computeInteractions')
+    await ComputedInteractions.attachFromCifOrCompute(structure)
+    console.timeEnd('computeInteractions');
+    console.log(ComputedInteractions.get(structure))
 
     const show = {
-        cartoon: false,
+        cartoon: true,
+        interaction: true,
         ballAndStick: true,
-        molecularSurface: true,
+        molecularSurface: false,
         gaussianSurface: false,
     }
 
     const cartoonRepr = getCartoonRepr()
+    const interactionRepr = getInteractionRepr()
     const ballAndStickRepr = getBallAndStickRepr()
     const molecularSurfaceRepr = getMolecularSurfaceRepr()
     const gaussianSurfaceRepr = getGaussianSurfaceRepr()
 
     if (show.cartoon) {
         cartoonRepr.setTheme({
-            color: reprCtx.colorThemeRegistry.create('secondary-structure', { structure }),
+            color: reprCtx.colorThemeRegistry.create('element-symbol', { structure }),
             size: reprCtx.sizeThemeRegistry.create('uniform', { structure })
         })
         await cartoonRepr.createOrUpdate({ ...CartoonRepresentationProvider.defaultValues, quality: 'auto' }, structure).run()
     }
 
+    if (show.interaction) {
+        interactionRepr.setTheme({
+            color: reprCtx.colorThemeRegistry.create('interaction-type', { structure }),
+            size: reprCtx.sizeThemeRegistry.create('uniform', { structure })
+        })
+        await interactionRepr.createOrUpdate({ ...InteractionsRepresentationProvider.defaultValues, quality: 'auto' }, structure).run()
+    }
+
     if (show.ballAndStick) {
         ballAndStickRepr.setTheme({
-            color: reprCtx.colorThemeRegistry.create('secondary-structure', { structure }),
-            size: reprCtx.sizeThemeRegistry.create('uniform', { structure })
+            color: reprCtx.colorThemeRegistry.create('element-symbol', { structure }),
+            size: reprCtx.sizeThemeRegistry.create('uniform', { structure }, { value: 1 })
         })
         await ballAndStickRepr.createOrUpdate({ ...BallAndStickRepresentationProvider.defaultValues, quality: 'auto' }, structure).run()
     }
@@ -133,6 +191,7 @@ async function init() {
     }
 
     if (show.cartoon) canvas3d.add(cartoonRepr)
+    if (show.interaction) canvas3d.add(interactionRepr)
     if (show.ballAndStick) canvas3d.add(ballAndStickRepr)
     if (show.molecularSurface) canvas3d.add(molecularSurfaceRepr)
     if (show.gaussianSurface) canvas3d.add(gaussianSurfaceRepr)

webpack.config.js

@@ -98,7 +98,6 @@ module.exports = [
     createApp('model-server-query'),
 
     createBrowserTest('font-atlas'),
-    createBrowserTest('render-asa'),
     createBrowserTest('marching-cubes'),
     createBrowserTest('render-lines'),
     createBrowserTest('render-mesh'),