diff --git a/src/mol-math/geometry/lookup3d/grid.ts b/src/mol-math/geometry/lookup3d/grid.ts index aa9393220b348ac1b0a206fc464bf95d38ca21b5..5270a27600e482edf597e63ec2b821116a16ab47 100644 --- a/src/mol-math/geometry/lookup3d/grid.ts +++ b/src/mol-math/geometry/lookup3d/grid.ts @@ -17,8 +17,8 @@ interface GridLookup3D<T = number> extends Lookup3D<T> { readonly buckets: { readonly offset: ArrayLike<number>, readonly count: ArrayLike<number>, readonly array: ArrayLike<number> } } -function GridLookup3D(data: PositionData, cellSize?: Vec3): GridLookup3D { - return new GridLookup3DImpl(data, cellSize); +function GridLookup3D(data: PositionData, cellSizeOrCount?: Vec3 | number): GridLookup3D { + return new GridLookup3DImpl(data, cellSizeOrCount); } export { GridLookup3D } @@ -47,8 +47,8 @@ class GridLookup3DImpl implements GridLookup3D<number> { return query(this.ctx); } - constructor(data: PositionData, cellSize?: Vec3) { - const structure = build(data, cellSize); + constructor(data: PositionData, cellSizeOrCount?: Vec3 | number) { + const structure = build(data, cellSizeOrCount); this.ctx = createContext(structure); this.boundary = { box: structure.boundingBox, sphere: structure.boundingSphere }; this.buckets = { offset: structure.bucketOffset, count: structure.bucketCounts, array: structure.bucketArray }; @@ -184,7 +184,7 @@ function getBoundary(data: PositionData) { return { boundingBox: boundaryHelper.getBox(), boundingSphere: boundaryHelper.getSphere() }; } -function build(data: PositionData, cellSize?: Vec3) { +function build(data: PositionData, cellSizeOrCount?: Vec3 | number) { const { boundingBox, boundingSphere } = getBoundary(data); // need to expand the grid bounds to avoid rounding errors const expandedBox = Box3D.expand(Box3D.empty(), boundingBox, Vec3.create(0.5, 0.5, 0.5)); @@ -195,13 +195,16 @@ function build(data: PositionData, cellSize?: Vec3) { const elementCount = OrderedSet.size(indices); + const cellCount = typeof cellSizeOrCount === 'number' ? cellSizeOrCount : 32 + const cellSize = Array.isArray(cellSizeOrCount) && cellSizeOrCount + if (cellSize) { size = [Math.ceil(S[0] / cellSize[0]), Math.ceil(S[1] / cellSize[1]), Math.ceil(S[2] / cellSize[2])]; delta = cellSize; } else if (elementCount > 0) { // size of the box - // required "grid volume" so that each cell contains on average 32 elements. - const V = Math.ceil(elementCount / 32); + // required "grid volume" so that each cell contains on average 'cellCount' elements. + const V = Math.ceil(elementCount / cellCount); const f = Math.pow(V / (S[0] * S[1] * S[2]), 1 / 3); size = [Math.ceil(S[0] * f), Math.ceil(S[1] * f), Math.ceil(S[2] * f)]; delta = [S[0] / size[0], S[1] / size[1], S[2] / size[2]]; diff --git a/src/mol-model/structure/model/properties/utils/secondary-structure.ts b/src/mol-model/structure/model/properties/utils/secondary-structure.ts index 2ebb150cf6386d583f2f9777f3abddfbe6651c52..5f977bb97b0ee730906568301c1b1a89c188a8a9 100644 --- a/src/mol-model/structure/model/properties/utils/secondary-structure.ts +++ b/src/mol-model/structure/model/properties/utils/secondary-structure.ts @@ -88,6 +88,9 @@ namespace DSSPType { /** max distance between two C-alpha atoms to check for hbond */ const caMaxDist = 7.0; +/** min distance between two C-alpha atoms to check for hbond */ +const caMinDist = 4.0; + function calcAtomicTraceLookup3D(hierarchy: AtomicHierarchy, conformation: AtomicConformation) { const { x, y, z } = conformation; const { moleculeType, traceElementIndex } = hierarchy.derived.residue @@ -99,7 +102,7 @@ function calcAtomicTraceLookup3D(hierarchy: AtomicHierarchy, conformation: Atomi _proteinResidues[_proteinResidues.length] = i } } - const lookup3d = GridLookup3D({ x, y, z, indices: SortedArray.ofSortedArray(indices) }); + const lookup3d = GridLookup3D({ x, y, z, indices: SortedArray.ofSortedArray(indices) }, 4); const proteinResidues = SortedArray.ofSortedArray<ResidueIndex>(_proteinResidues) return { lookup3d, proteinResidues } } @@ -160,6 +163,8 @@ function calcBackboneHbonds(hierarchy: AtomicHierarchy, conformation: AtomicConf const cPosPrev = Vec3.zero() const oPosPrev = Vec3.zero() + const caMinDistSq = caMinDist * caMinDist + for (let i = 0, il = proteinResidues.length; i < il; ++i) { const oPI = i const oRI = proteinResidues[i] @@ -178,9 +183,11 @@ function calcBackboneHbonds(hierarchy: AtomicHierarchy, conformation: AtomicConf position(cAtom, cPos) position(caAtom, caPos) - const { indices, count } = lookup3d.find(caPos[0], caPos[1], caPos[2], caMaxDist) + const { indices, count, squaredDistances } = lookup3d.find(caPos[0], caPos[1], caPos[2], caMaxDist) for (let j = 0; j < count; ++j) { + if (squaredDistances[j] < caMinDistSq) continue + const nPI = indices[j] // ignore bonds within a residue or to prev or next residue, TODO take chain border into account