Merge pull request #193 from JonStargaryen/anvil-fixes

ANVIL fixes

src/extensions/anvil/algorithm.ts

@@ -5,10 +5,10 @@
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
  */
 
-import { Structure, StructureElement, StructureProperties } from '../../mol-model/structure';
+import { Structure, StructureElement, StructureProperties, Unit } from '../../mol-model/structure';
 import { Task, RuntimeContext } from '../../mol-task';
 import { CentroidHelper } from '../../mol-math/geometry/centroid-helper';
-import { AccessibleSurfaceAreaProvider } from '../../mol-model-props/computed/accessible-surface-area';
+import { AccessibleSurfaceAreaParams } from '../../mol-model-props/computed/accessible-surface-area';
 import { Vec3 } from '../../mol-math/linear-algebra';
 import { getElementMoleculeType } from '../../mol-model/structure/util';
 import { MoleculeType } from '../../mol-model/structure/model/types';
@@ -16,6 +16,10 @@ import { AccessibleSurfaceArea } from '../../mol-model-props/computed/accessible
 import { ParamDefinition as PD } from '../../mol-util/param-definition';
 import { MembraneOrientation } from './prop';
 
+const LARGE_CA_THRESHOLD = 5000;
+const DEFAULT_UPDATE_INTERVAL = 10;
+const LARGE_CA_UPDATE_INTERVAL = 1;
+
 interface ANVILContext {
     structure: Structure,
 
@@ -24,19 +28,23 @@ interface ANVILContext {
     minThickness: number,
     maxThickness: number,
     asaCutoff: number,
+    adjust: number,
 
     offsets: ArrayLike<number>,
-    exposed: ArrayLike<boolean>,
+    exposed: ArrayLike<number>,
+    hydrophobic: ArrayLike<boolean>,
     centroid: Vec3,
-    extent: number
+    extent: number,
+    large: boolean
 };
 
 export const ANVILParams = {
-    numberOfSpherePoints: PD.Numeric(120, { min: 35, max: 700, step: 1 }, { description: 'Number of spheres/directions to test for membrane placement. Original value is 350.' }),
+    numberOfSpherePoints: PD.Numeric(140, { min: 35, max: 700, step: 1 }, { description: 'Number of spheres/directions to test for membrane placement. Original value is 350.' }),
     stepSize: PD.Numeric(1, { min: 0.25, max: 4, step: 0.25 }, { description: 'Thickness of membrane slices that will be tested' }),
     minThickness: PD.Numeric(20, { min: 10, max: 30, step: 1}, { description: 'Minimum membrane thickness used during refinement' }),
     maxThickness: PD.Numeric(40, { min: 30, max: 50, step: 1}, { description: 'Maximum membrane thickness used during refinement' }),
-    asaCutoff: PD.Numeric(40, { min: 10, max: 100, step: 1 }, { description: 'Absolute ASA cutoff above which residues will be considered' })
+    asaCutoff: PD.Numeric(40, { min: 10, max: 100, step: 1 }, { description: 'Relative ASA cutoff above which residues will be considered' }),
+    adjust: PD.Numeric(14, { min: 0, max: 30, step: 1 }, { description: 'Minimum length of membrane-spanning regions (original values: 14 for alpha-helices and 5 for beta sheets). Set to 0 to not optimize membrane thickness.' })
 };
 export type ANVILParams = typeof ANVILParams
 export type ANVILProps = PD.Values<ANVILParams>
@@ -54,22 +62,33 @@ export function computeANVIL(structure: Structure, props: ANVILProps) {
     });
 }
 
+// avoiding namespace lookup improved performance in Chrome (Aug 2020)
+const v3add = Vec3.add;
+const v3clone = Vec3.clone;
+const v3create = Vec3.create;
+const v3distance = Vec3.distance;
+const v3dot = Vec3.dot;
+const v3magnitude = Vec3.magnitude;
+const v3normalize = Vec3.normalize;
+const v3scale = Vec3.scale;
+const v3scaleAndAdd = Vec3.scaleAndAdd;
+const v3set = Vec3.set;
+const v3squaredDistance = Vec3.squaredDistance;
+const v3sub = Vec3.sub;
+const v3zero = Vec3.zero;
 
 const centroidHelper = new CentroidHelper();
-function initialize(structure: Structure, props: ANVILProps): ANVILContext {
+async function initialize(structure: Structure, props: ANVILProps, accessibleSurfaceArea: AccessibleSurfaceArea): Promise<ANVILContext> {
     const l = StructureElement.Location.create(structure);
-    const { label_atom_id, x, y, z } = StructureProperties.atom;
-    const elementCount = structure.polymerResidueCount;
+    const { label_atom_id, label_comp_id, x, y, z } = StructureProperties.atom;
+    const asaCutoff = props.asaCutoff / 100;
     centroidHelper.reset();
 
-    let offsets = new Int32Array(elementCount);
-    let exposed = new Array<boolean>(elementCount);
-
-    const accessibleSurfaceArea = structure && AccessibleSurfaceAreaProvider.get(structure);
-    const asa = accessibleSurfaceArea.value!;
+    const offsets = new Array<number>();
+    const exposed = new Array<number>();
+    const hydrophobic = new Array<boolean>();
 
-    const vec = Vec3();
-    let m = 0;
+    const vec = v3zero();
     for (let i = 0, il = structure.units.length; i < il; ++i) {
         const unit = structure.units[i];
         const { elements } = unit;
@@ -85,64 +104,82 @@ function initialize(structure: Structure, props: ANVILProps): ANVILContext {
             }
 
             // only CA is considered for downstream operations
-            if (label_atom_id(l) !== 'CA') {
+            if (label_atom_id(l) !== 'CA' && label_atom_id(l) !== 'BB') {
+                continue;
+            }
+
+            // original ANVIL only considers canonical amino acids
+            if (!MaxAsa[label_comp_id(l)]) {
                 continue;
             }
 
             // while iterating use first pass to compute centroid
-            Vec3.set(vec, x(l), y(l), z(l));
+            v3set(vec, x(l), y(l), z(l));
             centroidHelper.includeStep(vec);
 
             // keep track of offsets and exposed state to reuse
-            offsets[m] = structure.serialMapping.getSerialIndex(l.unit, l.element);
-            exposed[m] = AccessibleSurfaceArea.getValue(l, asa) > props.asaCutoff;
-
-            m++;
+            offsets.push(structure.serialMapping.getSerialIndex(l.unit, l.element));
+            if (AccessibleSurfaceArea.getValue(l, accessibleSurfaceArea) / MaxAsa[label_comp_id(l)] > asaCutoff) {
+                exposed.push(structure.serialMapping.getSerialIndex(l.unit, l.element));
+                hydrophobic.push(isHydrophobic(label_comp_id(l)));
+            }
         }
     }
 
-    // omit potentially empty tail1
-    offsets = offsets.slice(0, m);
-    exposed = exposed.slice(0, m);
-
     // calculate centroid and extent
     centroidHelper.finishedIncludeStep();
-    const centroid = centroidHelper.center;
+    const centroid = v3clone(centroidHelper.center);
     for (let k = 0, kl = offsets.length; k < kl; k++) {
         setLocation(l, structure, offsets[k]);
-        Vec3.set(vec, x(l), y(l), z(l));
+        v3set(vec, x(l), y(l), z(l));
         centroidHelper.radiusStep(vec);
     }
     const extent = 1.2 * Math.sqrt(centroidHelper.radiusSq);
 
     return {
         ...props,
-        structure: structure,
-
-        offsets: offsets,
-        exposed: exposed,
-        centroid: centroid,
-        extent: extent
+        structure,
+
+        offsets,
+        exposed,
+        hydrophobic,
+        centroid,
+        extent,
+        large: offsets.length > LARGE_CA_THRESHOLD
     };
 }
 
 export async function calculate(runtime: RuntimeContext, structure: Structure, params: ANVILProps): Promise<MembraneOrientation> {
-    const { label_comp_id } = StructureProperties.atom;
+    // can't get away with the default 92 points here
+    const asaProps = { ...PD.getDefaultValues(AccessibleSurfaceAreaParams), probeSize: 4.0, traceOnly: true, numberOfSpherePoints: 184 };
+    const accessibleSurfaceArea = await AccessibleSurfaceArea.compute(structure, asaProps).runInContext(runtime);
 
-    const ctx = initialize(structure, params);
-    const initialHphobHphil = HphobHphil.filtered(ctx, label_comp_id);
+    const ctx = await initialize(structure, params, accessibleSurfaceArea);
+    const initialHphobHphil = HphobHphil.initial(ctx);
 
-    const initialMembrane = findMembrane(ctx, generateSpherePoints(ctx, ctx.numberOfSpherePoints), initialHphobHphil, label_comp_id);
-    const alternativeMembrane = findMembrane(ctx, findProximateAxes(ctx, initialMembrane), initialHphobHphil, label_comp_id);
+    const initialMembrane = (await findMembrane(runtime, 'Placing initial membrane...', ctx, generateSpherePoints(ctx, ctx.numberOfSpherePoints), initialHphobHphil))!;
+    const refinedMembrane = (await findMembrane(runtime, 'Refining membrane placement...', ctx, findProximateAxes(ctx, initialMembrane), initialHphobHphil))!;
+    let membrane = initialMembrane.qmax! > refinedMembrane.qmax! ? initialMembrane : refinedMembrane;
+
+    if (ctx.adjust && !ctx.large) {
+        membrane = await adjustThickness(runtime, 'Adjusting membrane thickness...', ctx, membrane, initialHphobHphil);
+    }
 
-    const membrane = initialMembrane.qmax! > alternativeMembrane.qmax! ? initialMembrane : alternativeMembrane;
+    const normalVector = v3zero();
+    const center =  v3zero();
+    v3sub(normalVector, membrane.planePoint1, membrane.planePoint2);
+    v3normalize(normalVector, normalVector);
+
+    v3add(center, membrane.planePoint1, membrane.planePoint2);
+    v3scale(center, center, 0.5);
+    const extent = adjustExtent(ctx, membrane, center);
 
     return {
         planePoint1: membrane.planePoint1,
         planePoint2: membrane.planePoint2,
-        normalVector: membrane.normalVector!,
-        radius: ctx.extent,
-        centroid: ctx.centroid
+        normalVector,
+        centroid: center,
+        radius: extent
     };
 }
 
@@ -160,82 +197,79 @@ namespace MembraneCandidate {
         return {
             planePoint1: c1,
             planePoint2: c2,
-            stats: stats
+            stats
         };
     }
 
-    export function scored(spherePoint: Vec3, c1: Vec3, c2: Vec3, stats: HphobHphil, qmax: number, centroid: Vec3): MembraneCandidate {
-        const diam_vect = Vec3();
-        Vec3.sub(diam_vect, centroid, spherePoint);
+    export function scored(spherePoint: Vec3, planePoint1: Vec3, planePoint2: Vec3, stats: HphobHphil, qmax: number, centroid: Vec3): MembraneCandidate {
+        const normalVector = v3zero();
+        v3sub(normalVector, centroid, spherePoint);
         return {
-            planePoint1: c1,
-            planePoint2: c2,
-            stats: stats,
-            normalVector: diam_vect,
-            spherePoint: spherePoint,
-            qmax: qmax
+            planePoint1,
+            planePoint2,
+            stats,
+            normalVector,
+            spherePoint,
+            qmax
         };
     }
 }
 
-function findMembrane(ctx: ANVILContext, spherePoints: Vec3[], initialStats: HphobHphil, label_comp_id: StructureElement.Property<string>): MembraneCandidate {
-    const { centroid, stepSize, minThickness, maxThickness } = ctx;
+async function findMembrane(runtime: RuntimeContext, message: string | undefined, ctx: ANVILContext, spherePoints: Vec3[], initialStats: HphobHphil): Promise<MembraneCandidate | undefined> {
+    const { centroid, stepSize, minThickness, maxThickness, large } = ctx;
     // best performing membrane
-    let membrane: MembraneCandidate;
+    let membrane: MembraneCandidate | undefined;
     // score of the best performing membrane
     let qmax = 0;
 
     // construct slices of thickness 1.0 along the axis connecting the centroid and the spherePoint
-    const diam = Vec3();
-    for (let i = 0, il = spherePoints.length; i < il; i++) {
-        const spherePoint = spherePoints[i];
-        Vec3.sub(diam, centroid, spherePoint);
-        Vec3.scale(diam, diam, 2);
-        const diamNorm = Vec3.magnitude(diam);
-        const qvartemp = [];
+    const diam = v3zero();
+    for (let n = 0, nl = spherePoints.length; n < nl; n++) {
+        if (runtime.shouldUpdate && message && (n + 1) % (large ? LARGE_CA_UPDATE_INTERVAL : DEFAULT_UPDATE_INTERVAL) === 0) {
+            await runtime.update({ message, current: (n + 1), max: nl });
+        }
 
+        const spherePoint = spherePoints[n];
+        v3sub(diam, centroid, spherePoint);
+        v3scale(diam, diam, 2);
+        const diamNorm = v3magnitude(diam);
+
+        const sliceStats = HphobHphil.sliced(ctx, stepSize, spherePoint, diam, diamNorm);
+        const qvartemp = [];
         for (let i = 0, il = diamNorm - stepSize; i < il; i += stepSize) {
-            const c1 = Vec3();
-            const c2 = Vec3();
-            Vec3.scaleAndAdd(c1, spherePoint, diam, i / diamNorm);
-            Vec3.scaleAndAdd(c2, spherePoint, diam, (i + stepSize) / diamNorm);
+            const c1 = v3zero();
+            const c2 = v3zero();
+            v3scaleAndAdd(c1, spherePoint, diam, i / diamNorm);
+            v3scaleAndAdd(c2, spherePoint, diam, (i + stepSize) / diamNorm);
 
             // evaluate how well this membrane slice embeddeds the peculiar residues
-            const stats = HphobHphil.filtered(ctx, label_comp_id, (testPoint: Vec3) => isInMembranePlane(testPoint, diam, c1, c2));
+            const stats = sliceStats[Math.round(i / stepSize)];
             qvartemp.push(MembraneCandidate.initial(c1, c2, stats));
         }
 
-        let jmax = (minThickness / stepSize) - 1;
-
-        for (let width = 0, widthl = maxThickness; width < widthl;) {
-            const imax = qvartemp.length - 1 - jmax;
-
-            for (let i = 0, il = imax; i < il; i++) {
-                const c1 = qvartemp[i].planePoint1;
-                const c2 = qvartemp[i + jmax].planePoint2;
+        let jmax = Math.floor((minThickness / stepSize) - 1);
 
+        for (let width = 0, widthl = maxThickness; width <= widthl;) {
+            for (let i = 0, il = qvartemp.length - 1 - jmax; i < il; i++) {
                 let hphob = 0;
                 let hphil = 0;
-                let total = 0;
                 for (let j = 0; j < jmax; j++) {
                     const ij = qvartemp[i + j];
                     if (j === 0 || j === jmax - 1) {
-                        hphob += 0.5 * ij.stats.hphob;
+                        hphob += Math.floor(0.5 * ij.stats.hphob);
                         hphil += 0.5 * ij.stats.hphil;
                     } else {
                         hphob += ij.stats.hphob;
                         hphil += ij.stats.hphil;
                     }
-                    total += ij.stats.total;
                 }
 
-                const stats = HphobHphil.of(hphob, hphil, total);
-
                 if (hphob !== 0) {
+                    const stats = { hphob, hphil };
                     const qvaltest = qValue(stats, initialStats);
-                    if (qvaltest > qmax) {
+                    if (qvaltest >= qmax) {
                         qmax = qvaltest;
-                        membrane = MembraneCandidate.scored(spherePoint, c1, c2, HphobHphil.of(hphob, hphil, total), qmax, centroid);
+                        membrane = MembraneCandidate.scored(spherePoint, qvartemp[i].planePoint1, qvartemp[i + jmax].planePoint2, stats, qmax, centroid);
                     }
                 }
             }
@@ -244,7 +278,180 @@ function findMembrane(ctx: ANVILContext, spherePoints: Vec3[], initialStats: Hph
         }
     }
 
-    return membrane!;
+    return membrane;
+}
+
+/** Adjust membrane thickness by maximizing the number of membrane segments. */
+async function adjustThickness(runtime: RuntimeContext, message: string | undefined, ctx: ANVILContext, membrane: MembraneCandidate, initialHphobHphil: HphobHphil): Promise<MembraneCandidate> {
+    const { minThickness, large } = ctx;
+    const step = 0.3;
+    let maxThickness = v3distance(membrane.planePoint1, membrane.planePoint2);
+
+    let maxNos = membraneSegments(ctx, membrane).length;
+    let optimalThickness = membrane;
+
+    let n = 0;
+    const nl = Math.ceil((maxThickness - minThickness) / step);
+    while (maxThickness > minThickness) {
+        n++;
+        if (runtime.shouldUpdate && message && n % (large ? LARGE_CA_UPDATE_INTERVAL : DEFAULT_UPDATE_INTERVAL) === 0) {
+            await runtime.update({ message, current: n, max: nl });
+        }
+
+        const p = {
+            ...ctx,
+            maxThickness,
+            stepSize: step
+        };
+        const temp = await findMembrane(runtime, void 0, p, [membrane.spherePoint!], initialHphobHphil);
+        if (temp) {
+            const nos = membraneSegments(ctx, temp).length;
+            if (nos > maxNos) {
+                maxNos = nos;
+                optimalThickness = temp;
+            }
+        }
+        maxThickness -= step;
+    }
+
+    return optimalThickness;
+}
+
+/** Report auth_seq_ids for all transmembrane segments. Will reject segments that are shorter than the adjust parameter specifies. Missing residues are considered in-membrane. */
+function membraneSegments(ctx: ANVILContext, membrane: MembraneCandidate): ArrayLike<{ start: number, end: number }> {
+    const { offsets, structure, adjust } = ctx;
+    const { normalVector, planePoint1, planePoint2 } = membrane;
+    const { units } = structure;
+    const { elementIndices, unitIndices } = structure.serialMapping;
+    const testPoint = v3zero();
+    const { auth_seq_id } = StructureProperties.residue;
+
+    const d1 = -v3dot(normalVector!, planePoint1);
+    const d2 = -v3dot(normalVector!, planePoint2);
+    const dMin = Math.min(d1, d2);
+    const dMax = Math.max(d1, d2);
+
+    const inMembrane: { [k: string]: Set<number> } = Object.create(null);
+    const outMembrane: { [k: string]: Set<number> } = Object.create(null);
+    const segments: Array<{ start: number, end: number }> = [];
+    let authAsymId;
+    let lastAuthAsymId = null;
+    let authSeqId;
+    let lastAuthSeqId = units[0].model.atomicHierarchy.residues.auth_seq_id.value((units[0] as Unit.Atomic).chainIndex[0]) - 1;
+    let startOffset = 0;
+    let endOffset = 0;
+
+    // collect all residues in membrane layer
+    for (let k = 0, kl = offsets.length; k < kl; k++) {
+        const unit = units[unitIndices[offsets[k]]];
+        if (!Unit.isAtomic(unit)) throw 'Property only available for atomic models.';
+        const elementIndex = elementIndices[offsets[k]];
+
+        authAsymId = unit.model.atomicHierarchy.chains.auth_asym_id.value(unit.chainIndex[elementIndex]);
+        if (authAsymId !== lastAuthAsymId) {
+            if (!inMembrane[authAsymId]) inMembrane[authAsymId] = new Set<number>();
+            if (!outMembrane[authAsymId]) outMembrane[authAsymId] = new Set<number>();
+            lastAuthAsymId = authAsymId;
+        }
+
+        authSeqId = unit.model.atomicHierarchy.residues.auth_seq_id.value(unit.residueIndex[elementIndex]);
+        v3set(testPoint, unit.conformation.x(elementIndex), unit.conformation.y(elementIndex), unit.conformation.z(elementIndex));
+        if (_isInMembranePlane(testPoint, normalVector!, dMin, dMax)) {
+            inMembrane[authAsymId].add(authSeqId);
+        } else {
+            outMembrane[authAsymId].add(authSeqId);
+        }
+    }
+
+    for (let k = 0, kl = offsets.length; k < kl; k++) {
+        const unit = units[unitIndices[offsets[k]]];
+        if (!Unit.isAtomic(unit)) throw 'Property only available for atomic models.';
+        const elementIndex = elementIndices[offsets[k]];
+
+        authAsymId = unit.model.atomicHierarchy.chains.auth_asym_id.value(unit.chainIndex[elementIndex]);
+        authSeqId = unit.model.atomicHierarchy.residues.auth_seq_id.value(unit.residueIndex[elementIndex]);
+        if (inMembrane[authAsymId].has(authSeqId)) {
+            // chain change
+            if (authAsymId !== lastAuthAsymId) {
+                segments.push({ start: startOffset, end: endOffset });
+                lastAuthAsymId = authAsymId;
+                startOffset = k;
+                endOffset = k;
+            }
+
+            // sequence gaps
+            if (authSeqId !== lastAuthSeqId + 1) {
+                if (outMembrane[authAsymId].has(lastAuthSeqId + 1)) {
+                    segments.push({ start: startOffset, end: endOffset });
+                    startOffset = k;
+                }
+                lastAuthSeqId = authSeqId;
+                endOffset = k;
+            } else  {
+                lastAuthSeqId++;
+                endOffset++;
+            }
+        }
+    }
+    segments.push({ start: startOffset, end: endOffset });
+
+    const l = StructureElement.Location.create(structure);
+    let startAuth;
+    let endAuth;
+    const refinedSegments: Array<{ start: number, end: number }> = [];
+    for (let k = 0, kl = segments.length; k < kl; k++) {
+        const { start, end } = segments[k];
+        if (start === 0 || end === offsets.length - 1) continue;
+
+        // evaluate residues 1 pos outside of membrane
+        setLocation(l, structure, offsets[start - 1]);
+        v3set(testPoint, l.unit.conformation.x(l.element), l.unit.conformation.y(l.element), l.unit.conformation.z(l.element));
+        const d3 = -v3dot(normalVector!, testPoint);
+
+        setLocation(l, structure, offsets[end + 1]);
+        v3set(testPoint, l.unit.conformation.x(l.element), l.unit.conformation.y(l.element), l.unit.conformation.z(l.element));
+        const d4 = -v3dot(normalVector!, testPoint);
+
+        if (Math.min(d3, d4) < dMin && Math.max(d3, d4) > dMax) {
+            // reject this refinement
+            setLocation(l, structure, offsets[start]);
+            startAuth = auth_seq_id(l);
+            setLocation(l, structure, offsets[end]);
+            endAuth = auth_seq_id(l);
+            if (Math.abs(startAuth - endAuth) + 1 < adjust) {
+                return [];
+            }
+            refinedSegments.push(segments[k]);
+        }
+    }
+
+    return refinedSegments;
+}
+
+/** Filter for membrane residues and calculate the final extent of the membrane layer */
+function adjustExtent(ctx: ANVILContext, membrane: MembraneCandidate, centroid: Vec3): number {
+    const { offsets, structure } = ctx;
+    const { normalVector, planePoint1, planePoint2 } = membrane;
+    const l = StructureElement.Location.create(structure);
+    const testPoint = v3zero();
+    const { x, y, z } = StructureProperties.atom;
+
+    const d1 = -v3dot(normalVector!, planePoint1);
+    const d2 = -v3dot(normalVector!, planePoint2);
+    const dMin = Math.min(d1, d2);
+    const dMax = Math.max(d1, d2);
+    let extent = 0;
+
+    for (let k = 0, kl = offsets.length; k < kl; k++) {
+        setLocation(l, structure, offsets[k]);
+        v3set(testPoint, x(l), y(l), z(l));
+        if (_isInMembranePlane(testPoint, normalVector!, dMin, dMax)) {
+            const dsq = v3squaredDistance(testPoint, centroid);
+            if (dsq > extent) extent = dsq;
+        }
+    }
+
+    return Math.sqrt(extent);
 }
 
 function qValue(currentStats: HphobHphil, initialStats: HphobHphil): number {
@@ -262,23 +469,27 @@ function qValue(currentStats: HphobHphil, initialStats: HphobHphil): number {
 }
 
 export function isInMembranePlane(testPoint: Vec3, normalVector: Vec3, planePoint1: Vec3, planePoint2: Vec3): boolean {
-    const d1 = -Vec3.dot(normalVector, planePoint1);
-    const d2 = -Vec3.dot(normalVector, planePoint2);
-    const d = -Vec3.dot(normalVector, testPoint);
-    return d > Math.min(d1, d2) && d < Math.max(d1, d2);
+    const d1 = -v3dot(normalVector, planePoint1);
+    const d2 = -v3dot(normalVector, planePoint2);
+    return _isInMembranePlane(testPoint, normalVector, Math.min(d1, d2), Math.max(d1, d2));
 }
 
-// generates a defined number of points on a sphere with radius = extent around the specified centroid
+function _isInMembranePlane(testPoint: Vec3, normalVector: Vec3, min: number, max: number): boolean {
+    const d = -v3dot(normalVector, testPoint);
+    return d > min && d < max;
+}
+
+/** Generates a defined number of points on a sphere with radius = extent around the specified centroid */
 function generateSpherePoints(ctx: ANVILContext, numberOfSpherePoints: number): Vec3[] {
     const { centroid, extent } = ctx;
     const points = [];
     let oldPhi = 0, h, theta, phi;
     for(let k = 1, kl = numberOfSpherePoints + 1; k < kl; k++) {
-        h = -1 + 2 * (k - 1) / (numberOfSpherePoints - 1);
+        h = -1 + 2 * (k - 1) / (2 * numberOfSpherePoints - 1);
         theta = Math.acos(h);
-        phi = (k === 1 || k === numberOfSpherePoints) ? 0 : (oldPhi + 3.6 / Math.sqrt(numberOfSpherePoints * (1 - h * h))) % (2 * Math.PI);
+        phi = (k === 1 || k === numberOfSpherePoints) ? 0 : (oldPhi + 3.6 / Math.sqrt(2 * numberOfSpherePoints * (1 - h * h))) % (2 * Math.PI);
 
-        const point = Vec3.create(
+        const point = v3create(
             extent * Math.sin(phi) * Math.sin(theta) + centroid[0],
             extent * Math.cos(theta) + centroid[1],
             extent * Math.cos(phi) * Math.sin(theta) + centroid[2]
@@ -290,18 +501,18 @@ function generateSpherePoints(ctx: ANVILContext, numberOfSpherePoints: number):
     return points;
 }
 
-// generates sphere points close to that of the initial membrane
+/** Generates sphere points close to that of the initial membrane */
 function findProximateAxes(ctx: ANVILContext, membrane: MembraneCandidate): Vec3[] {
     const { numberOfSpherePoints, extent } = ctx;
     const points = generateSpherePoints(ctx, 30000);
     let j = 4;
     let sphere_pts2: Vec3[] = [];
+    const s = 2 * extent / numberOfSpherePoints;
     while (sphere_pts2.length < numberOfSpherePoints) {
-        const d = 2 * extent / numberOfSpherePoints + j;
-        const dsq = d * d;
+        const dsq = (s + j) * (s + j);
         sphere_pts2 = [];
         for (let i = 0, il = points.length; i < il; i++) {
-            if (Vec3.squaredDistance(points[i], membrane.spherePoint!) < dsq) {
+            if (v3squaredDistance(points[i], membrane.spherePoint!) < dsq) {
                 sphere_pts2.push(points[i]);
             }
         }
@@ -312,56 +523,80 @@ function findProximateAxes(ctx: ANVILContext, membrane: MembraneCandidate): Vec3
 
 interface HphobHphil {
     hphob: number,
-    hphil: number,
-    total: number
+    hphil: number
 }
 
 namespace HphobHphil {
-    export function of(hphob: number, hphil: number, total?: number) {
-        return {
-            hphob: hphob,
-            hphil: hphil,
-            total: !!total ? total : hphob + hphil
-        };
-    }
-
-    const testPoint = Vec3();
-    export function filtered(ctx: ANVILContext, label_comp_id: StructureElement.Property<string>, filter?: (test: Vec3) => boolean): HphobHphil {
-        const { offsets, exposed, structure } = ctx;
-        const l = StructureElement.Location.create(structure);
-        const { x, y, z } = StructureProperties.atom;
+    export function initial(ctx: ANVILContext): HphobHphil {
+        const { exposed, hydrophobic } = ctx;
         let hphob = 0;
         let hphil = 0;
-        for (let k = 0, kl = offsets.length; k < kl; k++) {
-            // ignore buried residues
-            if (!exposed[k]) {
-                continue;
+        for (let k = 0, kl = exposed.length; k < kl; k++) {
+            if (hydrophobic[k]) {
+                hphob++;
+            } else {
+                hphil++;
             }
+        }
+        return { hphob, hphil };
+    }
 
-            setLocation(l, structure, offsets[k]);
-            Vec3.set(testPoint, x(l), y(l), z(l));
-
-            // testPoints have to be in putative membrane layer
-            if (filter && !filter(testPoint)) {
-                continue;
-            }
+    const testPoint = v3zero();
+    export function sliced(ctx: ANVILContext, stepSize: number, spherePoint: Vec3, diam: Vec3, diamNorm: number): HphobHphil[] {
+        const { exposed, hydrophobic, structure } = ctx;
+        const { units, serialMapping } = structure;
+        const { unitIndices, elementIndices } = serialMapping;
+        const sliceStats: HphobHphil[] = [];
+        for (let i = 0, il = diamNorm - stepSize; i < il; i += stepSize) {
+            sliceStats[sliceStats.length] = { hphob: 0, hphil: 0 };
+        }
 
-            if (isHydrophobic(label_comp_id(l))) {
-                hphob++;
+        for (let i = 0, il = exposed.length; i < il; i++) {
+            const unit = units[unitIndices[exposed[i]]];
+            const elementIndex = elementIndices[exposed[i]];
+            v3set(testPoint, unit.conformation.x(elementIndex), unit.conformation.y(elementIndex), unit.conformation.z(elementIndex));
+            v3sub(testPoint, testPoint, spherePoint);
+            if (hydrophobic[i]) {
+                sliceStats[Math.floor(v3dot(testPoint, diam) / diamNorm / stepSize)].hphob++;
             } else {
-                hphil++;
+                sliceStats[Math.floor(v3dot(testPoint, diam) / diamNorm / stepSize)].hphil++;
             }
         }
-        return of(hphob, hphil);
+        return sliceStats;
     }
 }
 
-// ANVIL-specific (not general) definition of membrane-favoring amino acids
-const HYDROPHOBIC_AMINO_ACIDS = new Set(['ALA', 'CYS', 'GLY', 'HIS', 'ILE', 'LEU', 'MET', 'PHE', 'SER', 'THR', 'VAL']);
+/** ANVIL-specific (not general) definition of membrane-favoring amino acids */
+const HYDROPHOBIC_AMINO_ACIDS = new Set(['ALA', 'CYS', 'GLY', 'HIS', 'ILE', 'LEU', 'MET', 'PHE', 'SER', 'TRP', 'VAL']);
+/** Returns true if ANVIL considers this as amino acid that favors being embedded in a membrane */
 export function isHydrophobic(label_comp_id: string): boolean {
     return HYDROPHOBIC_AMINO_ACIDS.has(label_comp_id);
 }
 
+/** Accessible surface area used for normalization. ANVIL uses 'Total-Side REL' values from NACCESS, from: Hubbard, S. J., & Thornton, J. M. (1993). naccess. Computer Program, Department of Biochemistry and Molecular Biology, University College London, 2(1). */
+export const MaxAsa: { [k: string]: number } = {
+    'ALA': 69.41,
+    'ARG': 201.25,
+    'ASN': 106.24,
+    'ASP': 102.69,
+    'CYS': 96.75,
+    'GLU': 134.74,
+    'GLN': 140.99,
+    'GLY': 32.33,
+    'HIS': 147.08,
+    'ILE': 137.96,
+    'LEU': 141.12,
+    'LYS': 163.30,
+    'MET': 156.64,
+    'PHE': 164.11,
+    'PRO': 119.90,
+    'SER': 78.11,
+    'THR': 101.70,
+    'TRP': 211.26,
+    'TYR': 177.38,
+    'VAL': 114.28
+};
+
 function setLocation(l: StructureElement.Location, structure: Structure, serialIndex: number) {
     l.structure = structure;
     l.unit = structure.units[structure.serialMapping.unitIndices[serialIndex]];

src/extensions/anvil/prop.ts

@@ -11,7 +11,6 @@ import { CustomPropertyDescriptor } from '../../mol-model/custom-property';
 import { ANVILParams, ANVILProps, computeANVIL, isInMembranePlane } from './algorithm';
 import { CustomStructureProperty } from '../../mol-model-props/common/custom-structure-property';
 import { CustomProperty } from '../../mol-model-props/common/custom-property';
-import { AccessibleSurfaceAreaProvider } from '../../mol-model-props/computed/accessible-surface-area';
 import { Vec3 } from '../../mol-math/linear-algebra';
 import { QuerySymbolRuntime } from '../../mol-script/runtime/query/base';
 import { CustomPropSymbol } from '../../mol-script/language/symbol';
@@ -76,7 +75,6 @@ export const MembraneOrientationProvider: CustomStructureProperty.Provider<Membr
 });
 
 async function computeAnvil(ctx: CustomProperty.Context, data: Structure, props: Partial<ANVILProps>): Promise<MembraneOrientation> {
-    await AccessibleSurfaceAreaProvider.attach(ctx, data);
     const p = { ...PD.getDefaultValues(ANVILParams), ...props };
     return await computeANVIL(data, p).runInContext(ctx.runtime);
 }

src/extensions/anvil/representation.ts

@@ -9,7 +9,6 @@ import { ParamDefinition as PD } from '../../mol-util/param-definition';
 import { Vec3, Mat4 } from '../../mol-math/linear-algebra';
 import { Representation, RepresentationContext, RepresentationParamsGetter } from '../../mol-repr/representation';
 import { Structure } from '../../mol-model/structure';
-import { Spheres } from '../../mol-geo/geometry/spheres/spheres';
 import { StructureRepresentationProvider, StructureRepresentation, StructureRepresentationStateBuilder } from '../../mol-repr/structure/representation';
 import { MembraneOrientation } from './prop';
 import { ThemeRegistryContext } from '../../mol-theme/theme';
@@ -30,18 +29,9 @@ import { CustomProperty } from '../../mol-model-props/common/custom-property';
 
 const SharedParams = {
     color: PD.Color(ColorNames.lightgrey),
-    radiusFactor: PD.Numeric(0.8333, { min: 0.1, max: 3.0, step: 0.01 }, { description: 'Scale the radius of the membrane layer' })
+    radiusFactor: PD.Numeric(1.2, { min: 0.1, max: 3.0, step: 0.01 }, { description: 'Scale the radius of the membrane layer' })
 };
 
-const BilayerSpheresParams = {
-    ...Spheres.Params,
-    ...SharedParams,
-    sphereSize: PD.Numeric(1, { min: 0.1, max: 10, step: 0.1 }, { description: 'Size of spheres that represent membrane planes' }),
-    density: PD.Numeric(1, { min: 0.25, max: 10, step: 0.25 }, { description: 'Distance between spheres'})
-};
-export type BilayerSpheresParams = typeof BilayerSpheresParams
-export type BilayerSpheresProps = PD.Values<BilayerSpheresParams>
-
 const BilayerPlanesParams = {
     ...Mesh.Params,
     ...SharedParams,
@@ -66,7 +56,6 @@ const MembraneOrientationVisuals = {
 };
 
 export const MembraneOrientationParams = {
-    ...BilayerSpheresParams,
     ...BilayerPlanesParams,
     ...BilayerRimsParams,
     visuals: PD.MultiSelect(['bilayer-planes', 'bilayer-rims'], PD.objectToOptions(MembraneOrientationVisuals)),
@@ -104,16 +93,16 @@ function membraneLabel(data: Structure) {
 }
 
 function getBilayerRims(ctx: RuntimeContext, data: Structure, props: BilayerRimsProps, shape?: Shape<Lines>): Shape<Lines> {
-    const { planePoint1: p1, planePoint2: p2, centroid, normalVector: normal, radius } = MembraneOrientationProvider.get(data).value!;
+    const { planePoint1: p1, planePoint2: p2, centroid, radius } = MembraneOrientationProvider.get(data).value!;
     const scaledRadius = props.radiusFactor * radius;
     const builder = LinesBuilder.create(128, 64, shape?.geometry);
-    getLayerCircle(builder, p1, centroid, normal, scaledRadius, props);
-    getLayerCircle(builder, p2, centroid, normal, scaledRadius, props);
+    getLayerCircle(builder, p1, centroid, scaledRadius, props);
+    getLayerCircle(builder, p2, centroid, scaledRadius, props);
     return Shape.create('Bilayer rims', data, builder.getLines(), () => props.color, () => props.linesSize, () => membraneLabel(data));
 }
 
-function getLayerCircle(builder: LinesBuilder, p: Vec3, centroid: Vec3, normal: Vec3, radius: number, props: BilayerRimsProps, shape?: Shape<Lines>) {
-    const circle = getCircle(p, centroid, normal, radius);
+function getLayerCircle(builder: LinesBuilder, p: Vec3, centroid: Vec3, radius: number, props: BilayerRimsProps, shape?: Shape<Lines>) {
+    const circle = getCircle(p, centroid, radius);
     const { indices, vertices } = circle;
     for (let j = 0, jl = indices.length; j < jl; j += 3) {
         if (props.dashedLines && j % 2 === 1) continue; // draw every other segment to get dashes
@@ -130,8 +119,13 @@ function getLayerCircle(builder: LinesBuilder, p: Vec3, centroid: Vec3, normal:
 }
 
 const tmpMat = Mat4();
-function getCircle(p: Vec3, centroid: Vec3, normal: Vec3, radius: number) {
-    Mat4.targetTo(tmpMat, p, centroid, normal);
+const tmpV = Vec3();
+function getCircle(p: Vec3, centroid: Vec3, radius: number) {
+    if (Vec3.dot(Vec3.unitY, Vec3.sub(tmpV, p, centroid)) === 0) {
+        Mat4.targetTo(tmpMat, p, centroid, Vec3.unitY);
+    } else {
+        Mat4.targetTo(tmpMat, p, centroid, Vec3.unitX);
+    }
     Mat4.setTranslation(tmpMat, p);
     Mat4.mul(tmpMat, tmpMat, Mat4.rotX90);
 
@@ -140,16 +134,16 @@ function getCircle(p: Vec3, centroid: Vec3, normal: Vec3, radius: number) {
 }
 
 function getBilayerPlanes(ctx: RuntimeContext, data: Structure, props: BilayerPlanesProps, shape?: Shape<Mesh>): Shape<Mesh> {
-    const { planePoint1: p1, planePoint2: p2, centroid, normalVector: normal, radius } = MembraneOrientationProvider.get(data).value!;
+    const { planePoint1: p1, planePoint2: p2, centroid, radius } = MembraneOrientationProvider.get(data).value!;
     const state = MeshBuilder.createState(128, 64, shape && shape.geometry);
     const scaledRadius = props.radiusFactor * radius;
-    getLayerPlane(state, p1, centroid, normal, scaledRadius);
-    getLayerPlane(state, p2, centroid, normal, scaledRadius);
+    getLayerPlane(state, p1, centroid, scaledRadius);
+    getLayerPlane(state, p2, centroid, scaledRadius);
     return Shape.create('Bilayer planes', data, MeshBuilder.getMesh(state), () => props.color, () => 1, () => membraneLabel(data));
 }
 
-function getLayerPlane(state: MeshBuilder.State, p: Vec3, centroid: Vec3, normal: Vec3, radius: number) {
-    const circle = getCircle(p, centroid, normal, radius);
+function getLayerPlane(state: MeshBuilder.State, p: Vec3, centroid: Vec3, radius: number) {
+    const circle = getCircle(p, centroid, radius);
     state.currentGroup = 0;
     MeshBuilder.addPrimitive(state, Mat4.id, circle);
     MeshBuilder.addPrimitiveFlipped(state, Mat4.id, circle);

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

@@ -13,12 +13,14 @@ import { Structure, StructureElement, StructureProperties } from '../../../mol-m
 import { assignRadiusForHeavyAtoms } from './shrake-rupley/radii';
 import { ShrakeRupleyContext, VdWLookup, MaxAsa, DefaultMaxAsa } from './shrake-rupley/common';
 import { computeArea } from './shrake-rupley/area';
+import { SortedArray } from '../../../mol-data/int';
 
 export const ShrakeRupleyComputationParams = {
     numberOfSpherePoints: PD.Numeric(92, { min: 12, max: 360, 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.' })
+    nonPolymer: PD.Boolean(false, { description: 'Include non-polymer atoms as occluders.' }),
+    traceOnly: PD.Boolean(false, { description: 'Compute only using alpha-carbons, if true increase probeSize accordingly (e.g., 4 A). Considers only canonical amino acids.' })
 };
 export type ShrakeRupleyComputationParams = typeof ShrakeRupleyComputationParams
 export type ShrakeRupleyComputationProps = PD.Values<ShrakeRupleyComputationParams>
@@ -57,12 +59,13 @@ namespace AccessibleSurfaceArea {
 
     function initialize(structure: Structure, props: ShrakeRupleyComputationProps): ShrakeRupleyContext {
         const { elementCount, atomicResidueCount } = structure;
-        const { probeSize, nonPolymer, numberOfSpherePoints } = props;
+        const { probeSize, nonPolymer, traceOnly, numberOfSpherePoints } = props;
 
         return {
             structure,
             probeSize,
             nonPolymer,
+            traceOnly,
             spherePoints: generateSpherePoints(numberOfSpherePoints),
             scalingConstant: 4.0 * Math.PI / numberOfSpherePoints,
             maxLookupRadius: 2 * props.probeSize + 2 * VdWLookup[2], // 2x probe size + 2x largest VdW
@@ -99,9 +102,8 @@ namespace AccessibleSurfaceArea {
     }
 
     export function getValue(location: StructureElement.Location, accessibleSurfaceArea: AccessibleSurfaceArea) {
-        const { getSerialIndex } = location.structure.root.serialMapping;
         const { area, serialResidueIndex } = accessibleSurfaceArea;
-        const rSI = serialResidueIndex[getSerialIndex(location.unit, location.element)];
+        const rSI = serialResidueIndex[SortedArray.indexOf(SortedArray.ofSortedArray(location.structure.root.serialMapping.elementIndices), location.element)];
         if (rSI === -1) return -1;
         return area[rSI];
     }

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

@@ -6,9 +6,7 @@
  */
 
 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
@@ -28,77 +26,65 @@ export async function computeArea(runtime: RuntimeContext, ctx: ShrakeRupleyCont
     }
 }
 
-const aPos = Vec3();
-const bPos = Vec3();
-const testPoint = Vec3();
-
-function setLocation(l: StructureElement.Location, structure: Structure, serialIndex: number) {
-    l.structure = structure;
-    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, area, spherePoints, scalingConstant, maxLookupRadius, probeSize } = ctx;
-    const aLoc = StructureElement.Location.create(structure);
-    const bLoc = StructureElement.Location.create(structure);
-    const { x, y, z } = StructureProperties.atom;
-    const { lookup3d, serialMapping, unitIndexMap } = structure;
-    const { cumulativeUnitElementCount } = serialMapping;
+    const { lookup3d, serialMapping, unitIndexMap, units } = structure;
+    const { cumulativeUnitElementCount, elementIndices, unitIndices } = serialMapping;
 
     for (let aI = begin; aI < end; ++aI) {
-        const radius1 = VdWLookup[atomRadiusType[aI]];
-        if (radius1 === VdWLookup[0]) continue;
+        const vdw1 = VdWLookup[atomRadiusType[aI]];
+        if (vdw1 === VdWLookup[0]) continue;
 
-        setLocation(aLoc, structure, aI);
-        const aX = x(aLoc);
-        const aY = y(aLoc);
-        const aZ = z(aLoc);
+        const aUnit = units[unitIndices[aI]];
+        const aElementIndex = elementIndices[aI];
+        const aX = aUnit.conformation.x(aElementIndex);
+        const aY = aUnit.conformation.y(aElementIndex);
+        const aZ = aUnit.conformation.z(aElementIndex);
 
         // pre-filter by lookup3d (provides >10x speed-up compared to naive evaluation)
-        const { count, units, indices, squaredDistances } = lookup3d.find(aX, aY, aZ, maxLookupRadius);
+        const { count, units: lUnits, indices, squaredDistances } = lookup3d.find(aX, aY, aZ, maxLookupRadius);
 
+        // see optimizations proposed in Eisenhaber et al., 1995 (https://doi.org/10.1002/jcc.540160303)
         // collect neighbors for each atom
-        const cutoff1 = probeSize + probeSize + radius1;
+        const radius1 = probeSize + vdw1;
+        const cutoff1 = probeSize + radius1;
         const neighbors = []; // TODO reuse
         for (let iI = 0; iI < count; ++iI) {
-            const bUnit = units[iI];
-            StructureElement.Location.set(bLoc, ctx.structure, bUnit, bUnit.elements[indices[iI]]);
+            const bUnit = lUnits[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;
+            const bElementIndex = elementIndices[bI];
+
+            const vdw2 = VdWLookup[atomRadiusType[bI]];
+            if ((aUnit === bUnit && aElementIndex === bElementIndex) || vdw2 === VdWLookup[0]) continue;
+
+            const radius2 = probeSize + vdw2;
+            if (squaredDistances[iI] < (cutoff1 + vdw2) * (cutoff1 + vdw2)) {
+                const bElementIndex = elementIndices[bI];
+                // while here: compute values for later lookup
+                neighbors[neighbors.length] = [squaredDistances[iI],
+                    (squaredDistances[iI] + radius1 * radius1 - radius2 * radius2) / (2 * radius1),
+                    bUnit.conformation.x(bElementIndex) - aX,
+                    bUnit.conformation.y(bElementIndex) - aY,
+                    bUnit.conformation.z(bElementIndex) - aZ];
             }
         }
 
-        // 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;
+        // sort ascendingly by distance for improved downstream performance
+        neighbors.sort((a, b) => a[0] - b[0]);
 
-            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;
+        let accessiblePointCount = 0;
+        sl: for (let sI = 0; sI < spherePoints.length; ++sI) {
+            const [sX, sY, sZ] = spherePoints[sI];
+            for (let nI = 0; nI < neighbors.length; ++nI) {
+                const [, sqRadius, nX, nY, nZ] = neighbors[nI];
+                const dot = sX * nX + sY * nY + sZ * nZ;
+                if (dot > sqRadius) {
+                    continue sl;
                 }
             }
-
-            if (accessible) ++accessiblePointCount;
+            ++accessiblePointCount;
         }
 
-        area[serialResidueIndex[aI]] += scalingConstant * accessiblePointCount * scale * scale;
+        area[serialResidueIndex[aI]] += scalingConstant * accessiblePointCount * radius1 * radius1;
     }
 }

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

@@ -12,6 +12,7 @@ export interface ShrakeRupleyContext {
     spherePoints: Vec3[],
     probeSize: number,
     nonPolymer: boolean,
+    traceOnly: boolean,
     scalingConstant: number,
     maxLookupRadius: number,
     atomRadiusType: Int8Array,

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

@@ -5,7 +5,7 @@
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
  */
 
-import { ShrakeRupleyContext, VdWLookup } from './common';
+import { MaxAsa, ShrakeRupleyContext, VdWLookup } from './common';
 import { getElementIdx, isHydrogen } from '../../../../mol-model/structure/structure/unit/bonds/common';
 import { isPolymer, isNucleic, MoleculeType, ElementSymbol } from '../../../../mol-model/structure/model/types';
 import { VdwRadius } from '../../../../mol-model/structure/model/properties/atomic';
@@ -45,21 +45,26 @@ export function assignRadiusForHeavyAtoms(ctx: ShrakeRupleyContext) {
 
             // skip hydrogen atoms
             if (isHydrogen(elementIdx)) {
-                atomRadiusType[mj] = VdWLookup[0];
+                atomRadiusType[mj] = 0;
                 serialResidueIndex[mj] = -1;
                 continue;
             }
 
+            const atomId = label_atom_id(l);
             const moleculeType = getElementMoleculeType(unit, eI);
             // skip water and optionally non-polymer groups
             if (moleculeType === MoleculeType.Water || (!ctx.nonPolymer && !isPolymer(moleculeType))) {
-                atomRadiusType[mj] = VdWLookup[0];
+                atomRadiusType[mj] = 0;
                 serialResidueIndex[mj] = -1;
                 continue;
             }
 
-            const atomId = label_atom_id(l);
             const compId = label_comp_id(l);
+            if (ctx.traceOnly && ((atomId !== 'CA' && atomId !== 'BB') || !MaxAsa[compId])) {
+                atomRadiusType[mj] = 0;
+                serialResidueIndex[mj] = serialResidueIdx;
+                continue;
+            }
 
             if (isNucleic(moleculeType)) {
                 atomRadiusType[mj] = determineRadiusNucl(atomId, element, compId);