diff --git a/src/extensions/anvil/algorithm.ts b/src/extensions/anvil/algorithm.ts
index 63928b2ab7f2d6522ab113561af57e24842ee295..1465ed681fdbf72c1934efd495f7826cebad124f 100644
--- a/src/extensions/anvil/algorithm.ts
+++ b/src/extensions/anvil/algorithm.ts
@@ -16,6 +16,8 @@ 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;
+
interface ANVILContext {
structure: Structure,
@@ -27,7 +29,7 @@ interface ANVILContext {
adjust: number,
offsets: ArrayLike<number>,
- exposed: ArrayLike<boolean>,
+ exposed: ArrayLike<number>,
centroid: Vec3,
extent: number
};
@@ -56,7 +58,6 @@ export function computeANVIL(structure: Structure, props: ANVILProps) {
});
}
-
const centroidHelper = new CentroidHelper();
function initialize(structure: Structure, props: ANVILProps): ANVILContext {
const l = StructureElement.Location.create(structure);
@@ -65,7 +66,7 @@ function initialize(structure: Structure, props: ANVILProps): ANVILContext {
centroidHelper.reset();
let offsets = new Array<number>(elementCount);
- let exposed = new Array<boolean>(elementCount);
+ let exposed = new Array<number>(elementCount);
const accessibleSurfaceArea = structure && AccessibleSurfaceAreaProvider.get(structure);
const asa = accessibleSurfaceArea.value!;
@@ -73,6 +74,7 @@ function initialize(structure: Structure, props: ANVILProps): ANVILContext {
const vec = Vec3();
let m = 0;
+ let n = 0;
for (let i = 0, il = structure.units.length; i < il; ++i) {
const unit = structure.units[i];
const { elements } = unit;
@@ -102,15 +104,16 @@ function initialize(structure: Structure, props: ANVILProps): ANVILContext {
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) / MaxAsa[label_comp_id(l)] > asaCutoff;
- m++;
+ offsets[m++] = structure.serialMapping.getSerialIndex(l.unit, l.element);
+ if (AccessibleSurfaceArea.getValue(l, asa) / MaxAsa[label_comp_id(l)] > asaCutoff) {
+ exposed[n++] = structure.serialMapping.getSerialIndex(l.unit, l.element);
+ }
}
}
// omit potentially empty tail
offsets = offsets.slice(0, m);
- exposed = exposed.slice(0, m);
+ exposed = exposed.slice(0, n);
// calculate centroid and extent
centroidHelper.finishedIncludeStep();
@@ -121,7 +124,6 @@ function initialize(structure: Structure, props: ANVILProps): ANVILContext {
centroidHelper.radiusStep(vec);
}
const extent = Math.sqrt(centroidHelper.radiusSq);
- console.log(`center: ${centroid}, radius: ${Math.sqrt(centroidHelper.radiusSq)}`);
return {
...props,
@@ -139,19 +141,19 @@ export async function calculate(runtime: RuntimeContext, structure: Structure, p
const initialHphobHphil = HphobHphil.filtered(ctx);
if (runtime.shouldUpdate) {
- await runtime.update({ message: 'Placing initial membrane...' });
+ await runtime.update({ message: 'Placing initial membrane...', current: 1, max: 3 });
}
const initialMembrane = findMembrane(ctx, generateSpherePoints(ctx, ctx.numberOfSpherePoints), initialHphobHphil)!;
if (runtime.shouldUpdate) {
- await runtime.update({ message: 'Refining membrane placement...' });
+ await runtime.update({ message: 'Refining membrane placement...', current: 2, max: 3 });
}
const refinedMembrane = findMembrane(ctx, findProximateAxes(ctx, initialMembrane), initialHphobHphil)!;
let membrane = initialMembrane.qmax! > refinedMembrane.qmax! ? initialMembrane : refinedMembrane;
- if (ctx.adjust) {
+ if (ctx.adjust && ctx.offsets.length < LARGE_CA_THRESHOLD) {
if (runtime.shouldUpdate) {
- await runtime.update({ message: 'Adjusting membrane thickness...' });
+ await runtime.update({ message: 'Adjusting membrane thickness...', current: 3, max: 3 });
}
membrane = adjustThickness(ctx, membrane, initialHphobHphil);
}
@@ -195,7 +197,7 @@ namespace MembraneCandidate {
return {
planePoint1: c1,
planePoint2: c2,
- stats: stats
+ stats
};
}
@@ -240,19 +242,14 @@ function findMembrane(ctx: ANVILContext, spherePoints: Vec3[], initialStats: Hph
const dMax = Math.max(d1, d2);
// evaluate how well this membrane slice embeddeds the peculiar residues
- const stats = HphobHphil.filtered(ctx, (testPoint: Vec3) => _isInMembranePlane(testPoint, diam, dMin, dMax));
+ const stats = HphobHphil.filtered(ctx, { normalVector: diam, dMin, dMax });
qvartemp.push(MembraneCandidate.initial(c1, c2, stats));
}
let jmax = Math.floor((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;
-
+ for (let i = 0, il = qvartemp.length - 1 - jmax; i < il; i++) {
let hphob = 0;
let hphil = 0;
for (let j = 0; j < jmax; j++) {
@@ -271,7 +268,7 @@ function findMembrane(ctx: ANVILContext, spherePoints: Vec3[], initialStats: Hph
const qvaltest = qValue(stats, initialStats);
if (qvaltest >= qmax) {
qmax = qvaltest;
- membrane = MembraneCandidate.scored(spherePoint, c1, c2, stats, qmax, centroid);
+ membrane = MembraneCandidate.scored(spherePoint, qvartemp[i].planePoint1, qvartemp[i + jmax].planePoint2, stats, qmax, centroid);
}
}
}
@@ -525,25 +522,20 @@ namespace HphobHphil {
}
const testPoint = Vec3();
- export function filtered(ctx: ANVILContext, filter?: (test: Vec3) => boolean): HphobHphil {
- const { offsets, exposed, structure } = ctx;
+ 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 { x, y, z } = StructureProperties.atom;
let hphob = 0;
let hphil = 0;
- for (let k = 0, kl = offsets.length; k < kl; k++) {
- // ignore buried residues
- if (!exposed[k]) {
- continue;
- }
-
- setLocation(l, structure, offsets[k]);
+ for (let k = 0, kl = exposed.length; k < kl; k++) {
+ setLocation(l, structure, exposed[k]);
// testPoints have to be in putative membrane layer
if (filter) {
Vec3.set(testPoint, x(l), y(l), z(l));
- if (!filter(testPoint)) {
+ if (!_isInMembranePlane(testPoint, filter.normalVector, filter.dMin, filter.dMax)) {
continue;
}
}