diff --git a/src/extensions/anvil/algorithm.ts b/src/extensions/anvil/algorithm.ts
index b800a30643a0d57ca2b822c39689f48f048f293b..5425be9ca23e12a3d3ad129b3a9d323f0eb58204 100644
--- a/src/extensions/anvil/algorithm.ts
+++ b/src/extensions/anvil/algorithm.ts
@@ -18,6 +18,7 @@ import { MembraneOrientation } from './prop';
const LARGE_CA_THRESHOLD = 5000;
const UPDATE_INTERVAL = 10;
+const EMPTY_SET = new Set<number>();
interface ANVILContext {
structure: Structure,
@@ -31,6 +32,7 @@ interface ANVILContext {
offsets: ArrayLike<number>,
exposed: ArrayLike<number>,
+ hydrophobic: ArrayLike<boolean>,
centroid: Vec3,
extent: number
};
@@ -65,7 +67,6 @@ 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;
@@ -83,6 +84,7 @@ async function initialize(runtime: RuntimeContext, structure: Structure, props:
let offsets = new Array<number>(elementCount);
let exposed = new Array<number>(elementCount);
+ let hydrophobic = new Array<boolean>(elementCount);
// can't get away with the default 92 points here
const p = { ...PD.getDefaultValues(AccessibleSurfaceAreaParams), ...props, probeSize: 4.0, traceOnly: true, numberOfSpherePoints: 184 };
@@ -123,7 +125,9 @@ async function initialize(runtime: RuntimeContext, structure: Structure, props:
// keep track of offsets and exposed state to reuse
offsets[m++] = structure.serialMapping.getSerialIndex(l.unit, l.element);
if (AccessibleSurfaceArea.getValue(l, accessibleSurfaceArea) / MaxAsa[label_comp_id(l)] > asaCutoff) {
- exposed[n++] = structure.serialMapping.getSerialIndex(l.unit, l.element);
+ exposed[n] = structure.serialMapping.getSerialIndex(l.unit, l.element);
+ hydrophobic[n] = isHydrophobic(label_comp_id(l));
+ n++;
}
}
}
@@ -131,6 +135,7 @@ async function initialize(runtime: RuntimeContext, structure: Structure, props:
// omit potentially empty tail
offsets = offsets.slice(0, m);
exposed = exposed.slice(0, n);
+ hydrophobic = hydrophobic.splice(0, n);
// calculate centroid and extent
centroidHelper.finishedIncludeStep();
@@ -148,6 +153,7 @@ async function initialize(runtime: RuntimeContext, structure: Structure, props:
offsets,
exposed,
+ hydrophobic,
centroid,
extent
};
@@ -216,7 +222,8 @@ namespace MembraneCandidate {
}
async function findMembrane(runtime: RuntimeContext, message: string | undefined, ctx: ANVILContext, spherePoints: Vec3[], initialStats: HphobHphil): Promise<MembraneCandidate | undefined> {
- const { centroid, stepSize, minThickness, maxThickness } = ctx;
+ const { centroid, stepSize, minThickness, maxThickness, exposed, structure } = ctx;
+ const { x, y, z } = StructureProperties.atom;
// best performing membrane
let membrane: MembraneCandidate | undefined;
// score of the best performing membrane
@@ -232,21 +239,33 @@ async function findMembrane(runtime: RuntimeContext, message: string | undefined
const spherePoint = spherePoints[n];
v3sub(diam, centroid, spherePoint);
v3scale(diam, diam, 2);
- const diamNorm = v3magnitude(diam);
+ const diamDot = v3dot(diam, diam);
+ const diamNorm = Math.sqrt(diamDot);
const qvartemp = [];
+ const filter: Map<number, Set<number>> = new Map<number, Set<number>>();
+ const l = StructureElement.Location.create(structure);
+ const testPoint = v3zero();
+ for (let i = 0, il = exposed.length; i < il; i++) {
+ setLocation(l, structure, exposed[i]);
+ v3set(testPoint, x(l), y(l), z(l));
+ v3sub(testPoint, testPoint, spherePoint);
+ const membership = Math.floor(v3dot(testPoint, diam) / diamNorm / stepSize);
+ if (!filter.has(membership)) {
+ filter.set(membership, new Set<number>([i]));
+ } else {
+ filter.get(membership)?.add(i);
+ }
+ }
+
for (let i = 0, il = diamNorm - stepSize; i < il; i += stepSize) {
const c1 = v3zero();
const c2 = v3zero();
v3scaleAndAdd(c1, spherePoint, diam, i / diamNorm);
v3scaleAndAdd(c2, spherePoint, diam, (i + stepSize) / diamNorm);
- const d1 = -v3dot(diam, c1);
- const d2 = -v3dot(diam, c2);
- const dMin = Math.min(d1, d2);
- const dMax = Math.max(d1, d2);
// evaluate how well this membrane slice embeddeds the peculiar residues
- const stats = HphobHphil.filtered(ctx, { normalVector: diam, dMin, dMax });
+ const stats = HphobHphil.filtered(ctx, filter.get(i) || EMPTY_SET);
qvartemp.push(MembraneCandidate.initial(c1, c2, stats));
}
@@ -531,25 +550,31 @@ namespace HphobHphil {
};
}
- const testPoint = v3zero();
- export function filtered(ctx: ANVILContext, filter?: { normalVector: Vec3, dMin: number, dMax: number }): HphobHphil {
- const { exposed, structure } = ctx;
- const { label_comp_id } = StructureProperties.atom;
- const l = StructureElement.Location.create(structure);
+ // const testPoint = v3zero();
+ // export function filtered(ctx: ANVILContext, filter?: { normalVector: Vec3, dMin: number, dMax: number }): HphobHphil {
+ export function filtered(ctx: ANVILContext, filter?: Set<number>): HphobHphil {
+ // const { exposed, structure, hydrophobic } = ctx;
+ const { exposed, hydrophobic } = ctx;
+ // const l = StructureElement.Location.create(structure);
let hphob = 0;
let hphil = 0;
for (let k = 0, kl = exposed.length; k < kl; k++) {
- setLocation(l, structure, exposed[k]);
+ // setLocation(l, structure, exposed[k]);
// testPoints have to be in putative membrane layer
- if (filter) {
- v3set(testPoint, l.unit.conformation.x(l.element), l.unit.conformation.y(l.element), l.unit.conformation.z(l.element));
- if (!_isInMembranePlane(testPoint, filter.normalVector, filter.dMin, filter.dMax)) {
- continue;
- }
+ if (filter && !filter.has(k)) {
+ continue;
}
- if (isHydrophobic(label_comp_id(l))) {
+ // testPoints have to be in putative membrane layer
+ // if (filter) {
+ // v3set(testPoint, l.unit.conformation.x(l.element), l.unit.conformation.y(l.element), l.unit.conformation.z(l.element));
+ // if (!_isInMembranePlane(testPoint, filter.normalVector, filter.dMin, filter.dMax)) {
+ // continue;
+ // }
+ // }
+
+ if (hydrophobic[k]) {
hphob++;
} else {
hphil++;