diff --git a/src/mol-model/structure/model/properties/utils/accessible-surface-area.ts b/src/mol-model/structure/model/properties/utils/accessible-surface-area.ts index a21e2489b12e6fedd98086ba7fa1036e1a3397aa..16f9a4c6e31922969ad612774c7e91e87db8ef66 100644 --- a/src/mol-model/structure/model/properties/utils/accessible-surface-area.ts +++ b/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); } /**