calculation of atom ASA

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

@@ -27,25 +27,42 @@ export function computeModelASA(hierarchy: AtomicHierarchy, conformation: Atomic
     calculateASA(ctx);
 }
 
+const valueForIgnoredAtom = -1.0;
+
 function calculateASA(ctx: ASAContext) {
+    const { probeSize, spherePoints, cons } = ctx;
     const { atoms, residues, residueAtomSegments, derived } = ctx.hierarchy;
     const { type_symbol } = atoms;
+    const { x, y, z } = ctx.conformation;
+
     const radii: number[] = [];
-    const asa: number[] = [];
+    const atomAsa: number[] = [];
+    const residueAsa: number[] = [];
+
+    console.log(ctx.hierarchy);
+    console.log(ctx.conformation);
+
+    const position = (i: number, v: Vec3) => Vec3.set(v, x[i], y[i], z[i]);
+    const a1Pos = Vec3.zero();
+    const a2Pos = Vec3.zero();
 
     // 1. extract all heavy atoms
     // TODO can be more elegant by obtaining residue info once and then using offset to navigate to next residue
     for (let i = 0; i < type_symbol.rowCount; ++i) {
         // skip hydrogen atoms
-        if (type_symbol.value(i) === 'H' || type_symbol.value(i) === 'D' || type_symbol.value(i) === 'T')
+        if (type_symbol.value(i) === 'H' || type_symbol.value(i) === 'D' || type_symbol.value(i) === 'T') {
+            radii[radii.length] = valueForIgnoredAtom; // -1 is used to signal atoms not to be considered downstream
             continue;
+        }
 
         // determine group this atom belongs to
         const groupIndex = residueAtomSegments.index[i];
 
         // skip entities not part of a peptide chain
-        if (derived.residue.moleculeType[groupIndex] !== 4)
+        if (derived.residue.moleculeType[groupIndex] !== 4) {
+            radii[radii.length] = valueForIgnoredAtom;
             continue;
+        }
 
         const atomId = atoms.label_atom_id.value(i);
         const compId = residues.label_comp_id.value(groupIndex);
@@ -54,11 +71,66 @@ function calculateASA(ctx: ASAContext) {
         radii[radii.length] = determineRadius(atomId, element, compId);
     }
 
+    // 3. calculate the individual ASA of each atom
+    // TODO again might be more elegant to use offsets
+    // TODO distance is symmetric, omit redudant calcuations
+    for (let i = 0; i < radii.length; ++i) {
+        const radius = radii[i];
+
+        // skip invalid entries
+        if (radius === valueForIgnoredAtom) {
+            atomAsa[atomAsa.length] = valueForIgnoredAtom;
+            continue;
+        }
+
+        position(i, a1Pos);
+
+        // collect all neighboring atoms
+        const cutoff = probeSize + probeSize + radius;
+        const neighborIndices: number[] = [];
+        for (let k = 0; k < radii.length; ++k) {
+            const radius2 = radii[k];
+            if (i === k || radius2 === valueForIgnoredAtom)
+                continue;
+
+            position(k, a2Pos);
+
+            if (Vec3.distance(a1Pos, a2Pos) < cutoff + radius2) {
+                neighborIndices[neighborIndices.length] = k;
+            }
+        }
+
+        const r = probeSize + radius;
+        let accessiblePoints = 0;
+
+        for (let k = 0; k < spherePoints.length; ++k) {
+            const point = spherePoints[k];
+            let accessible = true;
+            const testPoint = [point[0] * r + a1Pos[0], point[1] * r + a1Pos[1], point[2] * r + a1Pos[2]] as Vec3;
+
+            for (let j = 0; j < neighborIndices.length; ++j) {
+                const naI = neighborIndices[j];
+                // store position of neighboring atom in a2Pos
+                position(naI, a2Pos);
+                const neighboringAtomRadius = radii[naI] + probeSize;
+                if (Vec3.squaredDistance(testPoint, a2Pos) < neighboringAtomRadius * neighboringAtomRadius) {
+                    accessible = false;
+                    break;
+                }
+            }
+
+            if (accessible) {
+                ++accessiblePoints;
+            }
+        }
+
+        atomAsa[atomAsa.length] = cons * accessiblePoints * r * r;
+    }
+
     // 3. for each residue
-    // 3a. calculate the individual ASA of each atom
     // 3b. sum up
     // 3c. (optionally) normalize by maximum value expected for a particular amino acid - needs lookup of max values
-
+    console.log(atomAsa);
 }
 
 /**