dssp tweak, caMinDist and grid cellCount

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]];

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