lookup3d at structure level, color provider

src/mol-model/structure/structure/accessible-surface-area.ts

@@ -0,0 +1,321 @@
+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Sebastian Bittrich <sebastian.bittrich@rcsb.org>
+ */
+
+import Structure from './structure';
+import { Task, RuntimeContext } from 'mol-task';
+import { BitFlags } from 'mol-util';
+import { ParamDefinition as PD } from 'mol-util/param-definition'
+import { Vec3 } from 'mol-math/linear-algebra';
+import { isPolymer, ElementSymbol, isNucleic, MoleculeType } from '../model/types';
+import { VdwRadius } from '../model/properties/atomic';
+import { isHydrogen, getElementIdx } from './unit/links/common';
+
+namespace AccessibleSurfaceArea {
+    // Chothia's amino acid atoms vdw radii
+    const trigonalCarbonVdw = 1.76;
+    const tetrahedralCarbonVdw = 1.87;
+    const trigonalNitrogenVdw = 1.65;
+    const tetrahedralNitrogenVdw = 1.50;
+    // deviating radii from definition in types.ts
+    const oxygenVdw = 1.40;
+    const sulfurVdw = 1.85;
+    // Chothia's nucleotide atom vdw radii
+    const nucCarbonVdw = 1.80;
+    const nucNitrogenVdw = 1.60;
+    const nucPhosphorusVdw = 1.40;
+    export const missingValue = -1.0;
+
+    /**
+     * 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<AccessibleSurfaceAreaComputationParams>> = {}) {
+        params = { ...PD.getDefaultValues(AccessibleSurfaceAreaComputationParams), ...params };
+        return Task.create('Compute Accessible Surface Area', async rtctx => {
+            return await _compute(rtctx, structure, params);
+        }).run();
+    }
+
+    async function _compute(rtctx: RuntimeContext, structure: Structure, params: Partial<PD.Values<AccessibleSurfaceAreaComputationParams>> = {}): Promise<AccessibleSurfaceArea> {
+        const ctx = initialize(rtctx, structure, params);
+
+        assignRadiusForHeavyAtoms(ctx);
+        computePerResidue(ctx);
+        normalizeAccessibleSurfaceArea(ctx);
+
+        return {
+            atomRadius: ctx.atomRadius!,
+            accessibleSurfaceArea: ctx.accessibleSurfaceArea!,
+            relativeAccessibleSurfaceArea: ctx.relativeAccessibleSurfaceArea!,
+            buried: (index: number) => ctx.relativeAccessibleSurfaceArea![index] < 0.16 // TODO this doesnt seem super elegant
+        };
+    }
+
+    interface AccessibleSurfaceAreaContext {
+        rtctx: RuntimeContext,
+        structure: Structure,
+        params: Partial<PD.Values<AccessibleSurfaceAreaComputationParams>>,
+        spherePoints: Vec3[],
+        cons: number,
+        maxLookupRadius: number,
+        atomRadius?: Float32Array, // TODO there are only 5-10 unique values in this array - rather than storing values, a int pointing to a dictionary will be far more memory efficient
+        accessibleSurfaceArea?: Float32Array,
+        relativeAccessibleSurfaceArea?: Float32Array
+    }
+
+    function normalizeAccessibleSurfaceArea(ctx: AccessibleSurfaceAreaContext) {
+        const { accessibleSurfaceArea, relativeAccessibleSurfaceArea, structure } = ctx;
+        const { residues, derived } = structure.model.atomicHierarchy;
+
+        for (let i = 0; i < residues.label_comp_id.rowCount; ++i) {
+            // skip entities not part of a polymer chain
+            if (!ctx.params.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;
+        }
+    }
+
+    async function computePerResidue(ctx: AccessibleSurfaceAreaContext) {
+        const { structure, atomRadius, accessibleSurfaceArea, spherePoints, cons, params, maxLookupRadius } = ctx;
+        const { probeSize } = params;
+        const { model, elementCount: atomCount } = 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]);
+        const aPos = Vec3.zero();
+        const bPos = Vec3.zero();
+
+        for (let aI = 0; aI < atomCount; ++aI) {
+            if (aI % 10000 === 0) {
+                console.log(`calculating accessible surface area, current: ${ aI }, max: ${ atomCount }`);
+                if (ctx.rtctx.shouldUpdate) {
+                    await ctx.rtctx.update({ message: 'calculating accessible surface area', current: aI, max: atomCount });
+                }
+            }
+
+            const radius1 = atomRadius![aI];
+            if (radius1 === missingValue) continue;
+
+            // pre-filter by lookup3d
+            // 36275 ms - lookup ~3000 ms
+            const { count, units, indices, squaredDistances } = lookup3d.find(x[aI], y[aI], z[aI], maxLookupRadius);
+
+            // 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 = atomRadius![bI];
+                if (aI === bI || radius2 === missingValue) 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];
+                const 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 cutoff3 = (atomRadius![nI] + probeSize!) * (atomRadius![nI] + probeSize!);
+                    if (Vec3.squaredDistance(testPoint, bPos) < cutoff3) {
+                        accessible = false;
+                        break;
+                    }
+                }
+
+                if (accessible) ++accessiblePointCount;
+            }
+
+            accessibleSurfaceArea![residueAtomSegments.index[aI]] += cons * accessiblePointCount * scalar * scalar;
+        }
+    }
+
+    function assignRadiusForHeavyAtoms(ctx: AccessibleSurfaceAreaContext) {
+        const { structure } = ctx;
+        const { model, elementCount: atomCount } = 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;
+        const residueCount = moleculeType.length;
+
+        // with atom and residue count at hand: initialize arrays
+        ctx.atomRadius = new Float32Array(atomCount - 1);
+        ctx.accessibleSurfaceArea = new Float32Array(residueCount - 1);
+        ctx.relativeAccessibleSurfaceArea = new Float32Array(residueCount - 1);
+
+        for (let aI = 0; aI < atomCount; ++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] = missingValue;
+                continue;
+            }
+
+            const residueType = moleculeType[rI];
+            // skip non-polymer groups
+            if (!ctx.params.nonPolymer) {
+                if (!isPolymer(residueType)) {
+                    ctx.atomRadius[aI] = missingValue;
+                    continue;
+                }
+            }
+
+            const atomId = label_atom_id.value(aI);
+            const residueId = label_comp_id.value(rI);
+            if (isNucleic(residueType)) {
+                ctx.atomRadius[aI] = determineRadiusNucl(atomId, element, residueId);
+            } else if (residueType === MoleculeType.protein) {
+                ctx.atomRadius[aI] = determineRadiusAmino(atomId, element, residueId);
+            } else {
+                ctx.atomRadius[aI] = VdwRadius(element);
+            }
+        }
+    }
+
+    /**
+     * Gets the van der Waals radius of the given atom following the values defined by Chothia (1976)
+     * J.Mol.Biol.105,1-14. NOTE: the vdw values defined by the paper assume no Hydrogens and thus "inflates" slightly
+     * the heavy atoms to account for Hydrogens.
+     */
+    function determineRadiusAmino(atomId: string, element: ElementSymbol, compId: string): number {
+        switch (element) {
+            case 'O':
+            return oxygenVdw;
+            case 'S':
+            return sulfurVdw;
+            case 'N':
+            return atomId === 'NZ' ? tetrahedralNitrogenVdw : trigonalNitrogenVdw;
+            case 'C':
+            switch (atomId) {
+                case 'C': case 'CE1': case'CE2': case 'CE3': case 'CH2': case 'CZ': case 'CZ2': case 'CZ3':
+                return trigonalCarbonVdw;
+                case 'CA': case 'CB': case 'CE': case 'CG1': case 'CG2':
+                return tetrahedralCarbonVdw;
+                default:
+                switch (compId) {
+                    case 'PHE': case 'TRP': case 'TYR': case 'HIS': case 'ASP': case 'ASN':
+                    return trigonalCarbonVdw;
+                    case 'PRO': case 'LYS': case 'ARG': case 'MET': case 'ILE': case 'LEU':
+                    return tetrahedralCarbonVdw;
+                    case 'GLU': case 'GLN':
+                    return atomId === 'CD' ? trigonalCarbonVdw : tetrahedralCarbonVdw;
+                }
+            }
+        }
+        return VdwRadius(element);
+    }
+
+    function determineRadiusNucl(atomId: string, element: ElementSymbol, compId: string): number {
+        switch (element) {
+            case 'C':
+            return nucCarbonVdw;
+            case 'N':
+            return nucNitrogenVdw;
+            case 'P':
+            return nucPhosphorusVdw;
+            case 'O':
+            return oxygenVdw;
+        }
+        return VdwRadius(element);
+    }
+
+    function initialize(rtctx: RuntimeContext, structure: Structure, params: Partial<PD.Values<AccessibleSurfaceAreaComputationParams>>): AccessibleSurfaceAreaContext {
+        return {
+            rtctx: rtctx,
+            structure: structure,
+            params: params,
+            spherePoints: generateSpherePoints(params.numberOfSpherePoints!),
+            cons: 4.0 * Math.PI / params.numberOfSpherePoints!,
+            maxLookupRadius: 2 * params.probeSize! + 2 * tetrahedralCarbonVdw // 2x probe size + 2x largest VdW
+        }
+    }
+
+    /** Creates a collection of points on a sphere by the Golden Section Spiral algorithm. */
+    function generateSpherePoints(numberOfSpherePoints: number): Vec3[] {
+        const points: Vec3[] = [];
+        const inc = Math.PI * (3.0 - Math.sqrt(5.0));
+        const offset = 2.0 / numberOfSpherePoints;
+        for (let k = 0; k < numberOfSpherePoints; ++k) {
+            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;
+        }
+        return points;
+    }
+
+    /** Maximum accessible surface area observed for amino acids. Taken from: http://dx.doi.org/10.1371/journal.pone.0080635 */
+    export const MaxAsa = {
+        'ALA': 121.0,
+        'ARG': 265.0,
+        'ASN': 187.0,
+        'ASP': 187.0,
+        'CYS': 148.0,
+        'GLU': 214.0,
+        'GLN': 214.0,
+        'GLY': 97.0,
+        'HIS': 216.0,
+        'ILE': 195.0,
+        'LEU': 191.0,
+        'LYS': 230.0,
+        'MET': 203.0,
+        'PHE': 228.0,
+        'PRO': 154.0,
+        'SER': 143.0,
+        'THR': 163.0,
+        'TRP': 264.0,
+        'TYR': 255.0,
+        'VAL': 165.0
+    }
+    export const DefaultMaxAsa = 121.0
+
+    export const AccessibleSurfaceAreaComputationParams = {
+        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.' })
+    }
+    export type AccessibleSurfaceAreaComputationParams = typeof AccessibleSurfaceAreaComputationParams
+
+    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
+        }
+    }
+}
+
+interface AccessibleSurfaceArea {
+    readonly atomRadius?: ArrayLike<number>
+    readonly accessibleSurfaceArea?: ArrayLike<number>
+    readonly relativeAccessibleSurfaceArea?: ArrayLike<number>
+    buried(index: number): boolean
+}
+
+export { AccessibleSurfaceArea }

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

@@ -18,8 +18,6 @@ import { IntMap, SortedArray } from 'mol-data/int';
 import { hash2, hashFnv32a } from 'mol-data/util';
 import { getAtomicPolymerElements, getCoarsePolymerElements, getAtomicGapElements, getCoarseGapElements } from './util/polymer';
 import { getNucleotideElements } from './util/nucleotide';
-import { computeAccessibleSurfaceArea } from './unit/accessible-surface-area/compute';
-import { AccessibleSurfaceArea } from './unit/accessible-surface-area/data';
 
 /**
  * A building block of a structure that corresponds to an atomic or
@@ -193,12 +191,6 @@ namespace Unit {
             return this.props.nucleotideElements.ref;
         }
 
-        get accessibleSurfaceArea() {
-            if (this.props.accessibleSurfaceArea.ref) return this.props.accessibleSurfaceArea.ref;
-            this.props.accessibleSurfaceArea.ref = computeAccessibleSurfaceArea(this);
-            return this.props.accessibleSurfaceArea.ref;
-        }
-
         getResidueIndex(elementIndex: StructureElement.UnitIndex) {
             return this.model.atomicHierarchy.residueAtomSegments.index[this.elements[elementIndex]];
         }
@@ -223,7 +215,6 @@ namespace Unit {
         polymerElements: ValueRef<SortedArray<ElementIndex> | undefined>
         gapElements: ValueRef<SortedArray<ElementIndex> | undefined>
         nucleotideElements: ValueRef<SortedArray<ElementIndex> | undefined>
-        accessibleSurfaceArea: ValueRef<AccessibleSurfaceArea | undefined>
     }
 
     function AtomicProperties(): AtomicProperties {
@@ -233,8 +224,7 @@ namespace Unit {
             rings: ValueRef.create(void 0),
             polymerElements: ValueRef.create(void 0),
             gapElements: ValueRef.create(void 0),
-            nucleotideElements: ValueRef.create(void 0),
-            accessibleSurfaceArea: ValueRef.create(void 0)
+            nucleotideElements: ValueRef.create(void 0)
         };
     }
 

src/mol-model/structure/structure/unit/accessible-surface-area/compute.ts

@@ -1,265 +0,0 @@
-/**
- * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
- *
- * @author Sebastian Bittrich <sebastian.bittrich@rcsb.org>
- */
-
-import Unit from '../../unit'
-import { Vec3 } from 'mol-math/linear-algebra';
-import { AccessibleSurfaceAreaComputationParams, AccessibleSurfaceArea, SolventAccessibility } from './data';
-import { isHydrogen, getElementIdx } from '../links/common'; // TODO these functions are relevant for many tasks: move them somewhere actually common
-import { ElementSymbol, MaxAsa, DefaultMaxAsa, isPolymer, isNucleic, MoleculeType } from 'mol-model/structure/model/types';
-import { VdwRadius } from 'mol-model/structure/model/properties/atomic/measures';
-import { ParamDefinition as PD } from 'mol-util/param-definition'
-
-// Chothia's amino acid atoms vdw radii
-const trigonalCarbonVdw = 1.76;
-const tetrahedralCarbonVdw = 1.87;
-const trigonalNitrogenVdw = 1.65;
-const tetrahedralNitrogenVdw = 1.50;
-/** deviating radii from definition in types.ts */
-const oxygenVdw = 1.40;
-const sulfurVdw = 1.85;
-// Chothia's nucleotide atom vdw radii
-const nucCarbonVdw = 1.80;
-const nucNitrogenVdw = 1.60;
-const nucPhosphorusVdw = 1.40;
-const missingAccessibleSurfaceAreaValue = -1.0;
-
-interface AccessibleSurfaceAreaContext {
-    unit: Unit.Atomic,
-    params: PD.Values<AccessibleSurfaceAreaComputationParams>,
-    spherePoints: Vec3[],
-    cons: number,
-    maxLookupRadius: number,
-    atomRadius?: Float32Array,
-    accessibleSurfaceArea?: Float32Array,
-    relativeAccessibleSurfaceArea?: Float32Array,
-    buried?: Uint8Array
-}
-
-/**
- * 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).
- */
-function computeAccessibleSurfaceArea(unit: Unit.Atomic, params?: PD.Values<AccessibleSurfaceAreaComputationParams>): AccessibleSurfaceArea {
-    if (!params) params = PD.getDefaultValues(AccessibleSurfaceAreaComputationParams);
-
-    // TODO non-polymer flag is currently useless as hetatms are located in different units - aim is not to color them, but to compute values correctly - relates to correct ASA computation for inter-chain contacts
-    console.log(`computing accessible surface area for unit #${ unit.id + 1 } - ${ params.numberOfSpherePoints } points, ${ params.probeSize } probe size, ${ params.nonPolymer ? 'honoring' : 'ignoring'} non-polymer atoms`);
-
-    const ctx = initialize(unit, params);
-    assignRadiusForHeavyAtoms(ctx);
-    computePerResidue(ctx);
-    normalizeAccessibleSurfaceArea(ctx);
-
-    return {
-        atomRadius: ctx.atomRadius!,
-        accessibleSurfaceArea: ctx.accessibleSurfaceArea!,
-        relativeAccessibleSurfaceArea: ctx.relativeAccessibleSurfaceArea!,
-        buried: ctx.buried!
-    };
-}
-
-function normalizeAccessibleSurfaceArea(ctx: AccessibleSurfaceAreaContext) {
-    const { residues, derived } = ctx.unit.model.atomicHierarchy;
-    const { accessibleSurfaceArea, relativeAccessibleSurfaceArea } = ctx;
-
-    for (let i = 0; i < residues.label_comp_id.rowCount; ++i) {
-        // skip entities not part of a polymer chain
-        if (!ctx.params.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;
-        ctx.buried![i] |= (rasa < ctx.params.buriedRasaThreshold ? SolventAccessibility.Flag.BURIED : SolventAccessibility.Flag.ACCESSIBLE)
-    }
-}
-
-/**
- * notes on performance - scenario: compute for first 10 units of 3j3q @ 960 sphere points
- * lookup3d + refinement: ~5000ms
- * naive approach: ~5600ms - higher variance
- */
-function computePerResidue(ctx: AccessibleSurfaceAreaContext) { // runs at roughly 5000 ms
-    const { atomRadius, accessibleSurfaceArea, maxLookupRadius, spherePoints, cons } = ctx;
-    const { probeSize } = ctx.params;
-    const { elements: atoms, residueIndex } = ctx.unit;
-    const { x, y, z } = ctx.unit.model.atomicConformation;
-    const atomCount = ctx.unit.elements.length;
-    const { lookup3d } = ctx.unit;
-
-    const position = (i: number, v: Vec3) => Vec3.set(v, x[i], y[i], z[i]);
-    const a1Pos = Vec3.zero();
-    const a2Pos = Vec3.zero();
-
-    for (let _aI = 0; _aI < atomCount; ++_aI) {
-        // map the atom index of this unit to the actual 'global' atom index
-        const aI = atoms[_aI];
-        const radii1 = atomRadius![aI];
-        if (radii1 === missingAccessibleSurfaceAreaValue) continue;
-
-        // find suitable neighbor candidates by lookup
-        const { indices, count } = lookup3d.find(x[aI], y[aI], z[aI], maxLookupRadius);
-        position(aI, a1Pos);
-
-        // refine list by actual criterion
-        const cutoff = probeSize + probeSize + radii1;
-        const filteredIndicies = []; // TODO might be better to use IntArray here and reuse it - how to find safe upper limit of possible neighborhood count - BioJava mentions 60 as relatively safe upper bound
-        for (let ni = 0; ni < count; ni++) {
-            const _bI = indices[ni];
-            const bI = atoms[_bI];
-            const radii2 = atomRadius![bI];
-            if (bI === aI || radii2 === missingAccessibleSurfaceAreaValue) continue;
-
-            const cutoff2 = (cutoff + radii2) * (cutoff + radii2);
-            // accurately check for neighborhood
-            position(bI, a2Pos);
-            if (Vec3.squaredDistance(a1Pos, a2Pos) < cutoff2) {
-                filteredIndicies[filteredIndicies.length] = bI;
-            }
-        }
-
-        // test sphere points
-        const r = probeSize + radii1;
-        let accessiblePointCount = 0;
-        for (let si = 0; si < spherePoints.length; ++si) {
-            const spherePoint = spherePoints[si];
-            const testPoint = [spherePoint[0] * r + a1Pos[0], spherePoint[1] * r + a1Pos[1], spherePoint[2] * r + a1Pos[2]] as Vec3;
-            let accessible = true;
-
-            for (let ni = 0; ni < filteredIndicies.length; ++ni) {
-                const naI = filteredIndicies[ni];
-                position(naI, a2Pos);
-                const cutoff3 = (atomRadius![naI] + probeSize) * (atomRadius![naI] + probeSize);
-                if (Vec3.squaredDistance(testPoint, a2Pos) < cutoff3) {
-                    accessible = false;
-                    break;
-                }
-            }
-
-            if (accessible) ++accessiblePointCount;
-        }
-
-        const value = cons * accessiblePointCount * r * r;
-        accessibleSurfaceArea![residueIndex[aI]] += value;
-        // +30% computation by normalizing partial solutions
-        // relativeAccessibleSurfaceArea[residueIndex[aI]] += value * (NormalizationFactors as any)[residueIndex[aI]];
-    }
-}
-
-function assignRadiusForHeavyAtoms(ctx: AccessibleSurfaceAreaContext) {
-    const atomCount = ctx.unit.elements.length;
-    const { elements: atoms, residueIndex } = ctx.unit;
-    const { residues } = ctx.unit.model.atomicHierarchy;
-    const { moleculeType } = ctx.unit.model.atomicHierarchy.derived.residue;
-    const { type_symbol, label_atom_id } = ctx.unit.model.atomicHierarchy.atoms;
-    const { label_comp_id } = ctx.unit.model.atomicHierarchy.residues;
-
-    const residueCount = residues.label_comp_id.rowCount;
-    ctx.atomRadius = new Float32Array(atomCount - 1);
-    ctx.accessibleSurfaceArea = new Float32Array(residueCount - 1);
-    ctx.relativeAccessibleSurfaceArea = new Float32Array(residueCount - 1);
-    ctx.buried = new Uint8Array(residueCount - 1);
-
-    for (let _aI = 0; _aI < atomCount; ++_aI) {
-        const aI =  atoms[_aI];
-        const raI = residueIndex[aI];
-        const aeI = getElementIdx(type_symbol.value(aI)!);
-
-        // skip hydrogen atoms
-        if (isHydrogen(aeI)) {
-            ctx.atomRadius[aI] = missingAccessibleSurfaceAreaValue;
-            continue;
-        }
-
-        // skip non-polymer groups
-        if (!ctx.params.nonPolymer) {
-            if (!isPolymer(moleculeType[raI])) {
-                ctx.atomRadius[aI] = missingAccessibleSurfaceAreaValue;
-                continue;
-            }
-        }
-
-        const atomId = label_atom_id.value(aI);
-        const element = type_symbol.value(aI);
-        const resn = label_comp_id.value(raI)!;
-
-        ctx.atomRadius[aI] = isNucleic(moleculeType[raI]) ? determineRadiusNucl(atomId, element, resn) : moleculeType[raI] === MoleculeType.protein ? determineRadiusAmino(atomId, element, resn) : VdwRadius(element);
-    }
-}
-
-/**
- * Gets the van der Waals radius of the given atom following the values defined by Chothia (1976)
- * J.Mol.Biol.105,1-14. NOTE: the vdw values defined by the paper assume no Hydrogens and thus "inflates" slightly
- * the heavy atoms to account for Hydrogens.
- */
-function determineRadiusAmino(atomId: string, element: ElementSymbol, compId: string): number {
-    switch (element) {
-        case 'O':
-        return oxygenVdw;
-        case 'S':
-        return sulfurVdw;
-        case 'N':
-        return atomId === 'NZ' ? tetrahedralNitrogenVdw : trigonalNitrogenVdw;
-        case 'C':
-        switch (atomId) {
-            case 'C': case 'CE1': case'CE2': case 'CE3': case 'CH2': case 'CZ': case 'CZ2': case 'CZ3':
-            return trigonalCarbonVdw;
-            case 'CA': case 'CB': case 'CE': case 'CG1': case 'CG2':
-            return tetrahedralCarbonVdw;
-            default:
-            switch (compId) {
-                case 'PHE': case 'TRP': case 'TYR': case 'HIS': case 'ASP': case 'ASN':
-                return trigonalCarbonVdw;
-                case 'PRO': case 'LYS': case 'ARG': case 'MET': case 'ILE': case 'LEU':
-                return tetrahedralCarbonVdw;
-                case 'GLU': case 'GLN':
-                return atomId === 'CD' ? trigonalCarbonVdw : tetrahedralCarbonVdw;
-            }
-        }
-    }
-    return VdwRadius(element);
-}
-
-function determineRadiusNucl(atomId: string, element: ElementSymbol, compId: string): number {
-    switch (element) {
-        case 'C':
-        return nucCarbonVdw;
-        case 'N':
-        return nucNitrogenVdw;
-        case 'P':
-        return nucPhosphorusVdw;
-        case 'O':
-        return oxygenVdw;
-    }
-    return VdwRadius(element);
-}
-
-function initialize(unit: Unit.Atomic, params: PD.Values<AccessibleSurfaceAreaComputationParams>): AccessibleSurfaceAreaContext {
-    return {
-        unit: unit,
-        params: params,
-        spherePoints: generateSpherePoints(params.numberOfSpherePoints),
-        cons: 4.0 * Math.PI / params.numberOfSpherePoints,
-        maxLookupRadius: 1.4 + 1.4 + 1.87 + 1.87
-    }
-}
-
-/** Creates a collection of points on a sphere by the Golden Section Spiral algorithm. */
-function generateSpherePoints(numberOfSpherePoints: number): Vec3[] {
-    const points: Vec3[] = [];
-    const inc = Math.PI * (3.0 - Math.sqrt(5.0));
-    const offset = 2.0 / numberOfSpherePoints;
-    for (let k = 0; k < numberOfSpherePoints; ++k) {
-        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;
-    }
-    return points;
-}
-
-export { computeAccessibleSurfaceArea, missingAccessibleSurfaceAreaValue }

src/mol-model/structure/structure/unit/accessible-surface-area/data.ts

@@ -1,35 +0,0 @@
-/**
- * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
- *
- * @author Sebastian Bittrich <sebastian.bittrich@rcsb.org>
- */
-
-import { ParamDefinition as PD } from 'mol-util/param-definition'
-import { BitFlags } from 'mol-util';
-
-export interface AccessibleSurfaceArea {
-    readonly atomRadius: ArrayLike<number>,
-    readonly accessibleSurfaceArea: ArrayLike<number>,
-    readonly relativeAccessibleSurfaceArea: ArrayLike<number>,
-    readonly buried: Uint8Array
-}
-
-export const AccessibleSurfaceAreaComputationParams = {
-    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.' })
-}
-
-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 type AccessibleSurfaceAreaComputationParams = typeof AccessibleSurfaceAreaComputationParams

src/mol-model/structure/structure/util/lookup3d.ts

@@ -5,12 +5,30 @@
  */
 
 import Structure from '../structure'
-import { Lookup3D, GridLookup3D, Result, Box3D, Sphere3D } from 'mol-math/geometry';
+import { Lookup3D, GridLookup3D, Box3D, Sphere3D, Result } from 'mol-math/geometry';
 import { Vec3 } from 'mol-math/linear-algebra';
 import { computeStructureBoundary } from './boundary';
 import { OrderedSet } from 'mol-data/int';
 import { StructureUniqueSubsetBuilder } from './unique-subset-builder';
 import StructureElement from '../element';
+import Unit from '../unit';
+
+export interface StructureResult extends Result<number> {
+    units: Unit[]
+}
+
+export namespace StructureResult {
+    export function add(result: StructureResult, unit: Unit, index: number, distSq: number) {
+        result.indices[result.count] = index;
+        result.units[result.count] = unit;
+        result.squaredDistances[result.count] = distSq;
+        result.count++;
+    }
+
+    export function create(): StructureResult {
+        return { count: 0, indices: [], units: [], squaredDistances: [] };
+    }
+}
 
 export class StructureLookup3D {
     private unitLookup: Lookup3D;
@@ -20,28 +38,27 @@ export class StructureLookup3D {
         return this.unitLookup.find(x, y, z, radius);
     }
 
-    // TODO: find another efficient way how to implement this instead of using "tuple".
-    // find(x: number, y: number, z: number, radius: number): Result<Element.Packed> {
-    //     Result.reset(this.result);
-    //     const { units } = this.structure;
-    //     const closeUnits = this.unitLookup.find(x, y, z, radius);
-    //     if (closeUnits.count === 0) return this.result;
-
-    //     for (let t = 0, _t = closeUnits.count; t < _t; t++) {
-    //         const unit = units[closeUnits.indices[t]];
-    //         Vec3.set(this.pivot, x, y, z);
-    //         if (!unit.conformation.operator.isIdentity) {
-    //             Vec3.transformMat4(this.pivot, this.pivot, unit.conformation.operator.inverse);
-    //         }
-    //         const unitLookup = unit.lookup3d;
-    //         const groupResult = unitLookup.find(this.pivot[0], this.pivot[1], this.pivot[2], radius);
-    //         for (let j = 0, _j = groupResult.count; j < _j; j++) {
-    //             Result.add(this.result, Element.Packed.create(unit.id, groupResult.indices[j]), groupResult.squaredDistances[j]);
-    //         }
-    //     }
-
-    //     return this.result;
-    // }
+    private result: StructureResult = StructureResult.create();
+    find(x: number, y: number, z: number, radius: number): StructureResult {
+        Result.reset(this.result);
+        const { units } = this.structure;
+        const closeUnits = this.unitLookup.find(x, y, z, radius);
+        if (closeUnits.count === 0) return this.result;
+
+        for (let t = 0, _t = closeUnits.count; t < _t; t++) {
+            const unit = units[closeUnits.indices[t]];
+            Vec3.set(this.pivot, x, y, z);
+            if (!unit.conformation.operator.isIdentity) {
+                Vec3.transformMat4(this.pivot, this.pivot, unit.conformation.operator.inverse);
+            }
+            const unitLookup = unit.lookup3d;
+            const groupResult = unitLookup.find(this.pivot[0], this.pivot[1], this.pivot[2], radius);
+            for (let j = 0, _j = groupResult.count; j < _j; j++) {
+                StructureResult.add(this.result, unit, groupResult.indices[j], groupResult.squaredDistances[j]);
+            }
+        }
+        return this.result;
+    }
 
     findIntoBuilder(x: number, y: number, z: number, radius: number, builder: StructureUniqueSubsetBuilder) {
         const { units } = this.structure;
@@ -98,8 +115,6 @@ export class StructureLookup3D {
         }
     }
 
-
-
     check(x: number, y: number, z: number, radius: number): boolean {
         const { units } = this.structure;
         const closeUnits = this.unitLookup.find(x, y, z, radius);

src/mol-theme/color.ts

@@ -7,7 +7,6 @@
 import { Color } from 'mol-util/color';
 import { Location } from 'mol-model/location';
 import { ColorType } from 'mol-geo/geometry/color-data';
-import { AccessibleSurfaceAreaColorThemeProvider } from './color/accessible-surface-area'
 import { CarbohydrateSymbolColorThemeProvider } from './color/carbohydrate-symbol';
 import { UniformColorThemeProvider } from './color/uniform';
 import { deepEqual } from 'mol-util';
@@ -63,7 +62,6 @@ 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

@@ -1,69 +0,0 @@
-/**
- * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
- *
- * @author Sebastian Bittrich <sebastian.bittrich@rcsb.org>
- */
-
-import { Color } from 'mol-util/color';
-import { Location } from 'mol-model/location';
-import { ColorTheme, LocationColor } from '../color';
-import { ParamDefinition as PD } from 'mol-util/param-definition'
-import { ThemeDataContext } from '../theme';
-import { ColorListOptions, ColorListName, ColorScale } from 'mol-util/color/scale';
-import { StructureElement, Unit } from 'mol-model/structure';
-import { missingAccessibleSurfaceAreaValue } from 'mol-model/structure/structure/unit/accessible-surface-area/compute';
-
-const DefaultColor = Color(0xFFFFFF)
-const Description = 'Assigns a color based on the relative accessible surface area of a residue.'
-
-export const AccessibleSurfaceAreaColorThemeParams = {
-    list: PD.ColorScale<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
-    let scale: ColorScale | undefined = undefined
-
-    if (ctx.structure) {
-        scale = ColorScale.create({
-            listOrName: props.list,
-            minLabel: 'Start',
-            maxLabel: 'End'
-        })
-        const scaleColor = scale.color
-
-        color = (location: Location): Color => {
-            if (StructureElement.isLocation(location)) {
-                if (Unit.isAtomic(location.unit)) {
-                    const value = location.unit.accessibleSurfaceArea.relativeAccessibleSurfaceArea[location.unit.residueIndex[location.element]];
-                    return value !== missingAccessibleSurfaceAreaValue ? scaleColor(value) : DefaultColor;
-                }
-            }
-
-            return DefaultColor
-        }
-    } else {
-        color = () => 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

@@ -7,11 +7,17 @@
 import './index.html'
 import { Canvas3D } from 'mol-canvas3d/canvas3d';
 import CIF, { CifFrame } from 'mol-io/reader/cif'
-import { Model, Structure } from 'mol-model/structure';
-import { ColorTheme } from 'mol-theme/color';
+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 { AccessibleSurfaceArea } from 'mol-model/structure/structure/accessible-surface-area';
+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 { ColorListName, ColorListOptions } from 'mol-util/color/scale';
 
 const parent = document.getElementById('app')!
 parent.style.width = '100%'
@@ -60,23 +66,78 @@ function getCartoonRepr() {
     return CartoonRepresentationProvider.factory(reprCtx, CartoonRepresentationProvider.getParams)
 }
 
-async function init() {
-    const cif = await downloadFromPdb(/*'3j3q'*/'1hrc')
+let accessibleSurfaceArea: AccessibleSurfaceArea;
+async function init(props = {}) {
+    const cif = await downloadFromPdb(
+        // '3j3q'
+        '1aon'
+        // '1acj'
+        )
     const models = await getModels(cif)
-
     const structure = await getStructure(models[0])
+
+    // async compute ASA
+    accessibleSurfaceArea = await AccessibleSurfaceArea.compute(structure)
+
     const cartoonRepr = getCartoonRepr()
 
-    console.time('ASA');
+    // create color theme
     cartoonRepr.setTheme({
-        color: reprCtx.colorThemeRegistry.create('accessible-surface-area', { structure }),
+        color: AccessibleSurfaceAreaColorTheme(reprCtx, { ...PD.getDefaultValues(AccessibleSurfaceAreaColorThemeParams), ...props }),
         size: reprCtx.sizeThemeRegistry.create('uniform', { structure })
     })
     await cartoonRepr.createOrUpdate({ ...CartoonRepresentationProvider.defaultValues, quality: 'auto' }, structure).run()
-    console.timeEnd('ASA');
 
     canvas3d.add(cartoonRepr)
     canvas3d.resetCamera()
 }
 
-init()
+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.ColorScale<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.isLocation(location)) {
+            if (Unit.isAtomic(location.unit)) {
+                const value = accessibleSurfaceArea.relativeAccessibleSurfaceArea![location.unit.residueIndex[location.element]];
+                return value !== AccessibleSurfaceArea.missingValue ? 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
+}