From a04ed8135db57ffe2e507cd77e35b5ca12b6eb3d Mon Sep 17 00:00:00 2001 From: David Sehnal <david.sehnal@gmail.com> Date: Thu, 19 Apr 2018 16:22:27 +0200 Subject: [PATCH] refactoring grid lookup (step 3) --- src/mol-math/geometry.ts | 12 + src/mol-math/geometry/_spec/lookup3d.spec.ts | 51 ++ src/mol-math/geometry/common.ts | 15 +- src/mol-math/geometry/grid-lookup.ts | 582 +++++++++--------- src/mol-math/geometry/lookup.ts | 0 src/mol-math/geometry/lookup3d/common.ts | 3 +- src/mol-math/geometry/lookup3d/grid.ts | 118 ++-- src/mol-math/geometry/primitives/box3d.ts | 30 +- src/mol-math/geometry/primitives/sphere3d.ts | 58 +- src/mol-math/geometry/symmetry-operator.ts | 2 +- .../structure/model/formats/mmcif/assembly.ts | 2 +- src/mol-model/structure/model/model.ts | 31 +- .../structure/model/properties/symmetry.ts | 2 +- .../structure/model/utils/compute-bonds.ts | 506 +++++++-------- .../structure/structure/element/group.ts | 16 +- .../structure/structure/structure.ts | 2 +- src/mol-model/structure/structure/unit.ts | 2 +- src/perf-tests/lookup3d.ts | 12 +- 18 files changed, 763 insertions(+), 681 deletions(-) create mode 100644 src/mol-math/geometry.ts create mode 100644 src/mol-math/geometry/_spec/lookup3d.spec.ts delete mode 100644 src/mol-math/geometry/lookup.ts diff --git a/src/mol-math/geometry.ts b/src/mol-math/geometry.ts new file mode 100644 index 000000000..f8134ddd0 --- /dev/null +++ b/src/mol-math/geometry.ts @@ -0,0 +1,12 @@ +/* + * Copyright (c) 2018 mol* contributors, licensed under MIT, See LICENSE file for more info. + * + * @author David Sehnal <david.sehnal@gmail.com> + */ + +export * from './geometry/common' +export * from './geometry/symmetry-operator' +export * from './geometry/lookup3d/common' +export * from './geometry/lookup3d/grid' +export * from './geometry/primitives/box3d' +export * from './geometry/primitives/sphere3d' \ No newline at end of file diff --git a/src/mol-math/geometry/_spec/lookup3d.spec.ts b/src/mol-math/geometry/_spec/lookup3d.spec.ts new file mode 100644 index 000000000..af2c30a46 --- /dev/null +++ b/src/mol-math/geometry/_spec/lookup3d.spec.ts @@ -0,0 +1,51 @@ +/* +* Copyright (c) 2018 mol* contributors, licensed under MIT, See LICENSE file for more info. +* +* @author David Sehnal <david.sehnal@gmail.com> +*/ + +import { GridLookup3D } from '../../geometry'; +import { sortArray } from 'mol-data/util'; +import { OrderedSet } from 'mol-data/int'; + +const xs = [0, 0, 1]; +const ys = [0, 1, 0]; +const zs = [0, 0, 0]; +const rs = [0, 0.5, 1/3]; + +describe('GridLookup3d', () => { + it('basic', () => { + const grid = GridLookup3D({ x: xs, y: ys, z: zs, indices: OrderedSet.ofBounds(0, 3) }); + + let r = grid.find(0, 0, 0, 0); + expect(r.count).toBe(1); + expect(r.indices[0]).toBe(0); + + r = grid.find(0, 0, 0, 1); + expect(r.count).toBe(3); + expect(sortArray(r.indices)).toEqual([0, 1, 2]); + }); + + it('radius', () => { + const grid = GridLookup3D({ x: xs, y: ys, z: zs, radius: [0, 0.5, 1 / 3], indices: OrderedSet.ofBounds(0, 3) }); + + let r = grid.find(0, 0, 0, 0); + expect(r.count).toBe(1); + expect(r.indices[0]).toBe(0); + + r = grid.find(0, 0, 0, 0.5); + expect(r.count).toBe(2); + expect(sortArray(r.indices)).toEqual([0, 1]); + }); + + it('indexed', () => { + const grid = GridLookup3D({ x: xs, y: ys, z: zs, indices: OrderedSet.ofSingleton(1), radius: rs }); + + let r = grid.find(0, 0, 0, 0); + expect(r.count).toBe(0); + + r = grid.find(0, 0, 0, 0.5); + expect(r.count).toBe(1); + expect(sortArray(r.indices)).toEqual([0]); + }); +}); \ No newline at end of file diff --git a/src/mol-math/geometry/common.ts b/src/mol-math/geometry/common.ts index 54d5b3b84..195224635 100644 --- a/src/mol-math/geometry/common.ts +++ b/src/mol-math/geometry/common.ts @@ -1,14 +1,17 @@ /* - * Copyright (c) 2018 mol* contributors, licensed under MIT, See LICENSE file for more info. - * - * @author David Sehnal <david.sehnal@gmail.com> - */ +* Copyright (c) 2018 mol* contributors, licensed under MIT, See LICENSE file for more info. +* +* @author David Sehnal <david.sehnal@gmail.com> +*/ + +import { OrderedSet } from 'mol-data/int' export interface PositionData { x: ArrayLike<number>, y: ArrayLike<number>, z: ArrayLike<number>, - radius?: ArrayLike<number>, // subset indices into the x/y/z/radius arrays - indices?: ArrayLike<number>, + indices: OrderedSet, + // optional element radius + radius?: ArrayLike<number> } diff --git a/src/mol-math/geometry/grid-lookup.ts b/src/mol-math/geometry/grid-lookup.ts index fa703c242..46660d95f 100644 --- a/src/mol-math/geometry/grid-lookup.ts +++ b/src/mol-math/geometry/grid-lookup.ts @@ -1,291 +1,291 @@ -/* - * Copyright (c) 2017 MolQL contributors, licensed under MIT, See LICENSE file for more info. - * - * @author David Sehnal <david.sehnal@gmail.com> - * @author Alexander Rose <alexander.rose@weirdbyte.de> - */ - -// TODO: 3d grid lookup (use single cell for small sets), make bounding sphere part of the structure -// TODO: assign radius to points? - -/** - * A "masked" 3D spatial lookup structure. - */ - -import Mask from 'mol-util/mask' - -export type FindFunc = (x: number, y: number, z: number, radius: number) => Result -export type CheckFunc = (x: number, y: number, z: number, radius: number) => boolean - -export interface Result { - readonly count: number, - readonly indices: number[], - readonly squaredDistances: number[] -} - -export interface Positions { x: ArrayLike<number>, y: ArrayLike<number>, z: ArrayLike<number> } - -interface GridLookup { find: (mask: Mask) => FindFunc, check: (mask: Mask) => CheckFunc } - -function GridLookup(positions: Positions): GridLookup { - const tree = build(createInputData(positions)); - return { - find(mask) { - const ctx = QueryContext.create(tree, mask, false); - return function (x: number, y: number, z: number, radius: number) { - QueryContext.update(ctx, x, y, z, radius); - nearest(ctx); - return ctx.buffer; - } - }, - check(mask) { - const ctx = QueryContext.create(tree, mask, true); - return function (x: number, y: number, z: number, radius: number) { - QueryContext.update(ctx, x, y, z, radius); - return nearest(ctx); - } - } - } -} - -interface InputData { - bounds: Box3D, - count: number, - positions: Positions -} - -interface Box3D { - min: number[], - max: number[] -} - -namespace Box3D { - export function createInfinite(): Box3D { - return { - min: [Number.MAX_VALUE, Number.MAX_VALUE, Number.MAX_VALUE], - max: [-Number.MAX_VALUE, -Number.MAX_VALUE, -Number.MAX_VALUE] - } - } -} - -/** -* Query context. Handles the actual querying. -*/ -interface QueryContext<T> { - structure: T, - pivot: number[], - radius: number, - radiusSq: number, - buffer: QueryContext.Buffer, - mask: Mask, - isCheck: boolean -} - -namespace QueryContext { - export interface Buffer extends Result { - count: number; - indices: any[]; - squaredDistances: number[]; - } - - export function add<T>(ctx: QueryContext<T>, distSq: number, index: number) { - const buffer = ctx.buffer; - buffer.squaredDistances[buffer.count] = distSq; - buffer.indices[buffer.count++] = index; - } - - function resetBuffer(buffer: Buffer) { buffer.count = 0; } - - function createBuffer(): Buffer { - return { - indices: [], - count: 0, - squaredDistances: [] - } - } - - /** - * Query the tree and store the result to this.buffer. Overwrites the old result. - */ - export function update<T>(ctx: QueryContext<T>, x: number, y: number, z: number, radius: number) { - ctx.pivot[0] = x; - ctx.pivot[1] = y; - ctx.pivot[2] = z; - ctx.radius = radius; - ctx.radiusSq = radius * radius; - resetBuffer(ctx.buffer); - } - - export function create<T>(structure: T, mask: Mask, isCheck: boolean): QueryContext<T> { - return { - structure, - buffer: createBuffer(), - pivot: [0.1, 0.1, 0.1], - radius: 1.1, - radiusSq: 1.1 * 1.1, - mask, - isCheck - } - } -} - -function createInputData(positions: { x: ArrayLike<number>, y: ArrayLike<number>, z: ArrayLike<number> }): InputData { - const { x, y, z } = positions; - const bounds = Box3D.createInfinite(); - const count = x.length; - const { min, max } = bounds; - - for (let i = 0; i < count; i++) { - min[0] = Math.min(x[i], min[0]); - min[1] = Math.min(y[i], min[1]); - min[2] = Math.min(z[i], min[2]); - max[0] = Math.max(x[i], max[0]); - max[1] = Math.max(y[i], max[1]); - max[2] = Math.max(z[i], max[2]); - } - - return { positions, bounds, count }; -} - -/** - * Adapted from https://github.com/arose/ngl - * MIT License Copyright (C) 2014+ Alexander Rose - */ - -interface SpatialHash { - size: number[], - min: number[], - grid: Uint32Array, - bucketOffset: Uint32Array, - bucketCounts: Int32Array, - bucketArray: Int32Array, - positions: Positions -} - -interface State { - size: number[], - positions: Positions, - bounds: Box3D, - count: number -} - -const enum Constants { Exp = 3 } - -function nearest(ctx: QueryContext<SpatialHash>): boolean { - const { min: [minX, minY, minZ], size: [sX, sY, sZ], bucketOffset, bucketCounts, bucketArray, grid, positions: { x: px, y: py, z: pz } } = ctx.structure; - const { radius: r, radiusSq: rSq, pivot: [x, y, z], isCheck, mask } = ctx; - - const loX = Math.max(0, (x - r - minX) >> Constants.Exp); - const loY = Math.max(0, (y - r - minY) >> Constants.Exp); - const loZ = Math.max(0, (z - r - minZ) >> Constants.Exp); - - const hiX = Math.min(sX, (x + r - minX) >> Constants.Exp); - const hiY = Math.min(sY, (y + r - minY) >> Constants.Exp); - const hiZ = Math.min(sZ, (z + r - minZ) >> Constants.Exp); - - for (let ix = loX; ix <= hiX; ix++) { - for (let iy = loY; iy <= hiY; iy++) { - for (let iz = loZ; iz <= hiZ; iz++) { - const bucketIdx = grid[(((ix * sY) + iy) * sZ) + iz]; - - if (bucketIdx > 0) { - const k = bucketIdx - 1; - const offset = bucketOffset[k]; - const count = bucketCounts[k]; - const end = offset + count; - - for (let i = offset; i < end; i++) { - const idx = bucketArray[i]; - if (!mask.has(idx)) continue; - - const dx = px[idx] - x; - const dy = py[idx] - y; - const dz = pz[idx] - z; - const distSq = dx * dx + dy * dy + dz * dz; - - if (distSq <= rSq) { - if (isCheck) return true; - QueryContext.add(ctx, distSq, idx) - } - } - } - } - } - } - return ctx.buffer.count > 0; -} - -function _build(state: State): SpatialHash { - const { bounds, size: [sX, sY, sZ], positions: { x: px, y: py, z: pz }, count } = state; - const n = sX * sY * sZ; - const { min: [minX, minY, minZ] } = bounds; - - let bucketCount = 0; - const grid = new Uint32Array(n); - const bucketIndex = new Int32Array(count); - for (let i = 0; i < count; i++) { - const x = (px[i] - minX) >> Constants.Exp; - const y = (py[i] - minY) >> Constants.Exp; - const z = (pz[i] - minZ) >> Constants.Exp; - const idx = (((x * sY) + y) * sZ) + z; - if ((grid[idx] += 1) === 1) { - bucketCount += 1 - } - bucketIndex[i] = idx; - } - - const bucketCounts = new Int32Array(bucketCount); - for (let i = 0, j = 0; i < n; i++) { - const c = grid[i]; - if (c > 0) { - grid[i] = j + 1; - bucketCounts[j] = c; - j += 1; - } - } - - const bucketOffset = new Uint32Array(count); - for (let i = 1; i < count; ++i) { - bucketOffset[i] += bucketOffset[i - 1] + bucketCounts[i - 1]; - } - - const bucketFill = new Int32Array(bucketCount); - const bucketArray = new Int32Array(count); - for (let i = 0; i < count; i++) { - const bucketIdx = grid[bucketIndex[i]] - if (bucketIdx > 0) { - const k = bucketIdx - 1; - bucketArray[bucketOffset[k] + bucketFill[k]] = i; - bucketFill[k] += 1; - } - } - - return { - size: state.size, - bucketArray, - bucketCounts, - bucketOffset, - grid, - min: state.bounds.min, - positions: state.positions - } -} - -function build({ positions, bounds, count }: InputData): SpatialHash { - const size = [ - ((bounds.max[0] - bounds.min[0]) >> Constants.Exp) + 1, - ((bounds.max[1] - bounds.min[1]) >> Constants.Exp) + 1, - ((bounds.max[2] - bounds.min[2]) >> Constants.Exp) + 1 - ]; - - const state: State = { - size, - positions, - bounds, - count - } - - return _build(state); -} - -export default GridLookup; \ No newline at end of file +// /* +// * Copyright (c) 2017 MolQL contributors, licensed under MIT, See LICENSE file for more info. +// * +// * @author David Sehnal <david.sehnal@gmail.com> +// * @author Alexander Rose <alexander.rose@weirdbyte.de> +// */ + +// // TODO: 3d grid lookup (use single cell for small sets), make bounding sphere part of the structure +// // TODO: assign radius to points? + +// /** +// * A "masked" 3D spatial lookup structure. +// */ + +// import Mask from 'mol-util/mask' + +// export type FindFunc = (x: number, y: number, z: number, radius: number) => Result +// export type CheckFunc = (x: number, y: number, z: number, radius: number) => boolean + +// export interface Result { +// readonly count: number, +// readonly indices: number[], +// readonly squaredDistances: number[] +// } + +// export interface Positions { x: ArrayLike<number>, y: ArrayLike<number>, z: ArrayLike<number> } + +// interface GridLookup { find: (mask: Mask) => FindFunc, check: (mask: Mask) => CheckFunc } + +// function GridLookup(positions: Positions): GridLookup { +// const tree = build(createInputData(positions)); +// return { +// find(mask) { +// const ctx = QueryContext.create(tree, mask, false); +// return function (x: number, y: number, z: number, radius: number) { +// QueryContext.update(ctx, x, y, z, radius); +// nearest(ctx); +// return ctx.buffer; +// } +// }, +// check(mask) { +// const ctx = QueryContext.create(tree, mask, true); +// return function (x: number, y: number, z: number, radius: number) { +// QueryContext.update(ctx, x, y, z, radius); +// return nearest(ctx); +// } +// } +// } +// } + +// interface InputData { +// bounds: Box3D, +// count: number, +// positions: Positions +// } + +// interface Box3D { +// min: number[], +// max: number[] +// } + +// namespace Box3D { +// export function createInfinite(): Box3D { +// return { +// min: [Number.MAX_VALUE, Number.MAX_VALUE, Number.MAX_VALUE], +// max: [-Number.MAX_VALUE, -Number.MAX_VALUE, -Number.MAX_VALUE] +// } +// } +// } + +// /** +// * Query context. Handles the actual querying. +// */ +// interface QueryContext<T> { +// structure: T, +// pivot: number[], +// radius: number, +// radiusSq: number, +// buffer: QueryContext.Buffer, +// mask: Mask, +// isCheck: boolean +// } + +// namespace QueryContext { +// export interface Buffer extends Result { +// count: number; +// indices: any[]; +// squaredDistances: number[]; +// } + +// export function add<T>(ctx: QueryContext<T>, distSq: number, index: number) { +// const buffer = ctx.buffer; +// buffer.squaredDistances[buffer.count] = distSq; +// buffer.indices[buffer.count++] = index; +// } + +// function resetBuffer(buffer: Buffer) { buffer.count = 0; } + +// function createBuffer(): Buffer { +// return { +// indices: [], +// count: 0, +// squaredDistances: [] +// } +// } + +// /** +// * Query the tree and store the result to this.buffer. Overwrites the old result. +// */ +// export function update<T>(ctx: QueryContext<T>, x: number, y: number, z: number, radius: number) { +// ctx.pivot[0] = x; +// ctx.pivot[1] = y; +// ctx.pivot[2] = z; +// ctx.radius = radius; +// ctx.radiusSq = radius * radius; +// resetBuffer(ctx.buffer); +// } + +// export function create<T>(structure: T, mask: Mask, isCheck: boolean): QueryContext<T> { +// return { +// structure, +// buffer: createBuffer(), +// pivot: [0.1, 0.1, 0.1], +// radius: 1.1, +// radiusSq: 1.1 * 1.1, +// mask, +// isCheck +// } +// } +// } + +// function createInputData(positions: { x: ArrayLike<number>, y: ArrayLike<number>, z: ArrayLike<number> }): InputData { +// const { x, y, z } = positions; +// const bounds = Box3D.createInfinite(); +// const count = x.length; +// const { min, max } = bounds; + +// for (let i = 0; i < count; i++) { +// min[0] = Math.min(x[i], min[0]); +// min[1] = Math.min(y[i], min[1]); +// min[2] = Math.min(z[i], min[2]); +// max[0] = Math.max(x[i], max[0]); +// max[1] = Math.max(y[i], max[1]); +// max[2] = Math.max(z[i], max[2]); +// } + +// return { positions, bounds, count }; +// } + +// /** +// * Adapted from https://github.com/arose/ngl +// * MIT License Copyright (C) 2014+ Alexander Rose +// */ + +// interface SpatialHash { +// size: number[], +// min: number[], +// grid: Uint32Array, +// bucketOffset: Uint32Array, +// bucketCounts: Int32Array, +// bucketArray: Int32Array, +// positions: Positions +// } + +// interface State { +// size: number[], +// positions: Positions, +// bounds: Box3D, +// count: number +// } + +// const enum Constants { Exp = 3 } + +// function nearest(ctx: QueryContext<SpatialHash>): boolean { +// const { min: [minX, minY, minZ], size: [sX, sY, sZ], bucketOffset, bucketCounts, bucketArray, grid, positions: { x: px, y: py, z: pz } } = ctx.structure; +// const { radius: r, radiusSq: rSq, pivot: [x, y, z], isCheck, mask } = ctx; + +// const loX = Math.max(0, (x - r - minX) >> Constants.Exp); +// const loY = Math.max(0, (y - r - minY) >> Constants.Exp); +// const loZ = Math.max(0, (z - r - minZ) >> Constants.Exp); + +// const hiX = Math.min(sX, (x + r - minX) >> Constants.Exp); +// const hiY = Math.min(sY, (y + r - minY) >> Constants.Exp); +// const hiZ = Math.min(sZ, (z + r - minZ) >> Constants.Exp); + +// for (let ix = loX; ix <= hiX; ix++) { +// for (let iy = loY; iy <= hiY; iy++) { +// for (let iz = loZ; iz <= hiZ; iz++) { +// const bucketIdx = grid[(((ix * sY) + iy) * sZ) + iz]; + +// if (bucketIdx > 0) { +// const k = bucketIdx - 1; +// const offset = bucketOffset[k]; +// const count = bucketCounts[k]; +// const end = offset + count; + +// for (let i = offset; i < end; i++) { +// const idx = bucketArray[i]; +// if (!mask.has(idx)) continue; + +// const dx = px[idx] - x; +// const dy = py[idx] - y; +// const dz = pz[idx] - z; +// const distSq = dx * dx + dy * dy + dz * dz; + +// if (distSq <= rSq) { +// if (isCheck) return true; +// QueryContext.add(ctx, distSq, idx) +// } +// } +// } +// } +// } +// } +// return ctx.buffer.count > 0; +// } + +// function _build(state: State): SpatialHash { +// const { bounds, size: [sX, sY, sZ], positions: { x: px, y: py, z: pz }, count } = state; +// const n = sX * sY * sZ; +// const { min: [minX, minY, minZ] } = bounds; + +// let bucketCount = 0; +// const grid = new Uint32Array(n); +// const bucketIndex = new Int32Array(count); +// for (let i = 0; i < count; i++) { +// const x = (px[i] - minX) >> Constants.Exp; +// const y = (py[i] - minY) >> Constants.Exp; +// const z = (pz[i] - minZ) >> Constants.Exp; +// const idx = (((x * sY) + y) * sZ) + z; +// if ((grid[idx] += 1) === 1) { +// bucketCount += 1 +// } +// bucketIndex[i] = idx; +// } + +// const bucketCounts = new Int32Array(bucketCount); +// for (let i = 0, j = 0; i < n; i++) { +// const c = grid[i]; +// if (c > 0) { +// grid[i] = j + 1; +// bucketCounts[j] = c; +// j += 1; +// } +// } + +// const bucketOffset = new Uint32Array(count); +// for (let i = 1; i < count; ++i) { +// bucketOffset[i] += bucketOffset[i - 1] + bucketCounts[i - 1]; +// } + +// const bucketFill = new Int32Array(bucketCount); +// const bucketArray = new Int32Array(count); +// for (let i = 0; i < count; i++) { +// const bucketIdx = grid[bucketIndex[i]] +// if (bucketIdx > 0) { +// const k = bucketIdx - 1; +// bucketArray[bucketOffset[k] + bucketFill[k]] = i; +// bucketFill[k] += 1; +// } +// } + +// return { +// size: state.size, +// bucketArray, +// bucketCounts, +// bucketOffset, +// grid, +// min: state.bounds.min, +// positions: state.positions +// } +// } + +// function build({ positions, bounds, count }: InputData): SpatialHash { +// const size = [ +// ((bounds.max[0] - bounds.min[0]) >> Constants.Exp) + 1, +// ((bounds.max[1] - bounds.min[1]) >> Constants.Exp) + 1, +// ((bounds.max[2] - bounds.min[2]) >> Constants.Exp) + 1 +// ]; + +// const state: State = { +// size, +// positions, +// bounds, +// count +// } + +// return _build(state); +// } + +// export default GridLookup; \ No newline at end of file diff --git a/src/mol-math/geometry/lookup.ts b/src/mol-math/geometry/lookup.ts deleted file mode 100644 index e69de29bb..000000000 diff --git a/src/mol-math/geometry/lookup3d/common.ts b/src/mol-math/geometry/lookup3d/common.ts index 1b4212be1..a5a3e8d58 100644 --- a/src/mol-math/geometry/lookup3d/common.ts +++ b/src/mol-math/geometry/lookup3d/common.ts @@ -32,6 +32,5 @@ export interface Lookup3D { // The result is mutated with each call to find. find(x: number, y: number, z: number, radius: number): Result, check(x: number, y: number, z: number, radius: number): boolean, - boundingBox: Box3D, - boundingSphere: Sphere3D + boundary: { box: Box3D, sphere: Sphere3D } } \ No newline at end of file diff --git a/src/mol-math/geometry/lookup3d/grid.ts b/src/mol-math/geometry/lookup3d/grid.ts index b676904d0..b83bb4272 100644 --- a/src/mol-math/geometry/lookup3d/grid.ts +++ b/src/mol-math/geometry/lookup3d/grid.ts @@ -5,26 +5,47 @@ * @author Alexander Rose <alexander.rose@weirdbyte.de> */ -import { Result/*, Lookup3D*/ } from './common' +import { Result, Lookup3D } from './common' import { Box3D } from '../primitives/box3d'; import { Sphere3D } from '../primitives/sphere3d'; import { PositionData } from '../common'; import { Vec3 } from '../../linear-algebra'; +import { OrderedSet } from 'mol-data/int'; -// class GridLookup3D implements Lookup3D { -// private result = Result.create(); -// private pivot = [0.1, 0.1, 0.1]; -// private radiusSq = 0.1; - -// find(x: number, y: number, z: number, radius: number): Result { -// throw new Error("Method not implemented."); -// } -// check(x: number, y: number, z: number, radius: number): boolean { -// throw new Error("Method not implemented."); -// } -// boundingBox: Box3D; -// boundingSphere: Sphere3D; -// } +function GridLookup3D(data: PositionData): Lookup3D { + return new GridLookup3DImpl(data); +} + +export { GridLookup3D } + +class GridLookup3DImpl implements Lookup3D { + private ctx: QueryContext; + boundary: Lookup3D['boundary']; + + find(x: number, y: number, z: number, radius: number): Result { + this.ctx.x = x; + this.ctx.y = y; + this.ctx.z = z; + this.ctx.radius = radius; + this.ctx.isCheck = false; + query(this.ctx); + return this.ctx.result; + } + check(x: number, y: number, z: number, radius: number): boolean { + this.ctx.x = x; + this.ctx.y = y; + this.ctx.z = z; + this.ctx.radius = radius; + this.ctx.isCheck = true; + return query(this.ctx); + } + + constructor(data: PositionData) { + const structure = build(data); + this.ctx = createContext(structure); + this.boundary = { box: structure.boundingBox, sphere: structure.boundingSphere }; + } +} /** * Adapted from https://github.com/arose/ngl @@ -35,11 +56,11 @@ interface InputData { x: ArrayLike<number>, y: ArrayLike<number>, z: ArrayLike<number>, + indices: OrderedSet, radius?: ArrayLike<number>, - indices: ArrayLike<number> } -interface SpatialHash { +interface Grid3D { size: number[], min: number[], grid: Uint32Array, @@ -64,7 +85,7 @@ interface BuildState { count: number } -function _build(state: BuildState): SpatialHash { +function _build(state: BuildState): Grid3D { const { expandedBox, size: [sX, sY, sZ], data: { x: px, y: py, z: pz, radius, indices }, count, delta } = state; const n = sX * sY * sZ; const { min: [minX, minY, minZ] } = expandedBox; @@ -73,8 +94,8 @@ function _build(state: BuildState): SpatialHash { let bucketCount = 0; const grid = new Uint32Array(n); const bucketIndex = new Int32Array(count); - for (let t = 0, _t = indices.length; t < _t; t++) { - const i = indices[t]; + for (let t = 0, _t = OrderedSet.size(indices); t < _t; t++) { + const i = OrderedSet.getAt(indices, t); const x = Math.floor((px[i] - minX) / delta[0]); const y = Math.floor((py[i] - minY) / delta[1]); const z = Math.floor((pz[i] - minZ) / delta[2]); @@ -87,8 +108,8 @@ function _build(state: BuildState): SpatialHash { } if (radius) { - for (let t = 0, _t = indices.length; t < _t; t++) { - const i = indices[t]; + for (let t = 0, _t = OrderedSet.size(indices); t < _t; t++) { + const i = OrderedSet.getAt(indices, t); if (radius[i] > maxRadius) maxRadius = radius[i]; } } @@ -135,24 +156,22 @@ function _build(state: BuildState): SpatialHash { } } -export function build(data: PositionData) { +function build(data: PositionData) { const boundingBox = Box3D.computeBounding(data); + // need to expand the grid bounds to avoid rounding errors const expandedBox = Box3D.expand(boundingBox, Vec3.create(0.5, 0.5, 0.5)); const boundingSphere = Sphere3D.computeBounding(data); - - // TODO: specialized code that does not require the indices annotation? - const indices = data.indices ? data.indices as number[] : new Int32Array(data.x.length) as any as number[]; - if (!data.indices) { - for (let i = 0, _i = data.x.length; i < _i; i++) indices[i] = i; - } + const { indices } = data; const S = Vec3.sub(Vec3.zero(), expandedBox.max, expandedBox.min); let delta, size; - if (indices.length > 0) { + const elementCount = OrderedSet.size(indices); + + if (elementCount > 0) { // size of the box // required "grid volume" so that each cell contains on average 32 elements. - const V = Math.ceil(indices.length / 32); + const V = Math.ceil(elementCount / 32); 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]]; @@ -175,19 +194,35 @@ export function build(data: PositionData) { expandedBox, boundingBox, boundingSphere, - count: indices.length, + count: elementCount, delta } return _build(state); } -export function inSphere(structure: SpatialHash, x: number, y: number, z: number, radius: number): Result { - const { min, size: [sX, sY, sZ], bucketOffset, bucketCounts, bucketArray, grid, data: { x: px, y: py, z: pz, indices }, delta } = structure; - //const { radius: r, radiusSq: rSq, pivot: [x, y, z], isCheck, mask } = ctx; - const r = radius, rSq = r * r; +interface QueryContext { + structure: Grid3D, + x: number, + y: number, + z: number, + radius: number, + result: Result, + isCheck: boolean +} + +function createContext(structure: Grid3D): QueryContext { + return { structure, x: 0.1, y: 0.1, z: 0.1, radius: 0.1, result: Result.create(), isCheck: false } +} + +function query(ctx: QueryContext): boolean { + const { min, size: [sX, sY, sZ], bucketOffset, bucketCounts, bucketArray, grid, data: { x: px, y: py, z: pz, indices, radius }, delta, maxRadius } = ctx.structure; + const { radius: inputRadius, isCheck, x, y, z, result } = ctx; + + const r = inputRadius + maxRadius; + const rSq = r * r; - const result = Result.create(); + Result.reset(result); const loX = Math.max(0, Math.floor((x - r - min[0]) / delta[0])); const loY = Math.max(0, Math.floor((y - r - min[1]) / delta[1])); @@ -197,7 +232,7 @@ export function inSphere(structure: SpatialHash, x: number, y: number, z: number const hiY = Math.min(sY - 1, Math.floor((y + r - min[1]) / delta[1])); const hiZ = Math.min(sZ - 1, Math.floor((z + r - min[2]) / delta[2])); - if (loX > hiX || loY > hiY || loZ > hiZ) return result; + if (loX > hiX || loY > hiY || loZ > hiZ) return false; for (let ix = loX; ix <= hiX; ix++) { for (let iy = loY; iy <= hiY; iy++) { @@ -211,7 +246,7 @@ export function inSphere(structure: SpatialHash, x: number, y: number, z: number const end = offset + count; for (let i = offset; i < end; i++) { - const idx = indices[bucketArray[i]]; + const idx = OrderedSet.getAt(indices, bucketArray[i]); const dx = px[idx] - x; const dy = py[idx] - y; @@ -219,12 +254,13 @@ export function inSphere(structure: SpatialHash, x: number, y: number, z: number const distSq = dx * dx + dy * dy + dz * dz; if (distSq <= rSq) { - // if (isCheck) return true; + if (maxRadius > 0 && Math.sqrt(distSq) - radius![idx] > inputRadius) continue; + if (isCheck) return true; Result.add(result, bucketArray[i], distSq); } } } } } - return result; + return result.count > 0; } \ No newline at end of file diff --git a/src/mol-math/geometry/primitives/box3d.ts b/src/mol-math/geometry/primitives/box3d.ts index 0bdd9bd9d..21a579516 100644 --- a/src/mol-math/geometry/primitives/box3d.ts +++ b/src/mol-math/geometry/primitives/box3d.ts @@ -6,6 +6,7 @@ import { Vec3 } from '../../linear-algebra' import { PositionData } from '../common' +import { OrderedSet } from 'mol-data/int'; interface Box3D { min: Vec3, max: Vec3 } @@ -15,28 +16,15 @@ namespace Box3D { const max = [-Number.MAX_VALUE, -Number.MAX_VALUE, -Number.MAX_VALUE]; const { x, y, z, indices } = data; - - if (indices) { - for (let t = 0, _t = indices.length; t < _t; t++) { - const i = indices[t]; - min[0] = Math.min(x[i], min[0]); - min[1] = Math.min(y[i], min[1]); - min[2] = Math.min(z[i], min[2]); - max[0] = Math.max(x[i], max[0]); - max[1] = Math.max(y[i], max[1]); - max[2] = Math.max(z[i], max[2]); - } - } else { - for (let i = 0, _i = x.length; i < _i; i++) { - min[0] = Math.min(x[i], min[0]); - min[1] = Math.min(y[i], min[1]); - min[2] = Math.min(z[i], min[2]); - max[0] = Math.max(x[i], max[0]); - max[1] = Math.max(y[i], max[1]); - max[2] = Math.max(z[i], max[2]); - } + for (let t = 0, _t = OrderedSet.size(indices); t < _t; t++) { + const i = OrderedSet.getAt(indices, t); + min[0] = Math.min(x[i], min[0]); + min[1] = Math.min(y[i], min[1]); + min[2] = Math.min(z[i], min[2]); + max[0] = Math.max(x[i], max[0]); + max[1] = Math.max(y[i], max[1]); + max[2] = Math.max(z[i], max[2]); } - return { min: Vec3.create(min[0], min[1], min[2]), max: Vec3.create(max[0], max[1], max[2]) } } diff --git a/src/mol-math/geometry/primitives/sphere3d.ts b/src/mol-math/geometry/primitives/sphere3d.ts index 701daf098..bd74a09a2 100644 --- a/src/mol-math/geometry/primitives/sphere3d.ts +++ b/src/mol-math/geometry/primitives/sphere3d.ts @@ -6,6 +6,7 @@ import { Vec3 } from '../../linear-algebra' import { PositionData } from '../common' +import { OrderedSet } from 'mol-data/int'; interface Sphere3D { center: Vec3, radius: number } @@ -15,44 +16,25 @@ namespace Sphere3D { let cx = 0, cy = 0, cz = 0; let radiusSq = 0; - if (indices) { - for (let t = 0, _t = indices.length; t < _t; t++) { - const i = indices[t]; - cx += x[i]; - cy += y[i]; - cz += z[i]; - } - - if (indices.length > 0) { - cx /= indices.length; - cy /= indices.length; - cz /= indices.length; - } - - for (let t = 0, _t = indices.length; t < _t; t++) { - const i = indices[t]; - const dx = x[i] - cx, dy = y[i] - cy, dz = z[i] - cz; - const d = dx * dx + dy * dy + dz * dz; - if (d > radiusSq) radiusSq = d; - } - } else { - for (let i = 0, _i = x.length; i < _i; i++) { - cx += x[i]; - cy += y[i]; - cz += z[i]; - } - - if (x.length > 0) { - cx /= x.length; - cy /= x.length; - cz /= x.length; - } - - for (let i = 0, _i = x.length; i < _i; i++) { - const dx = x[i] - cx, dy = y[i] - cy, dz = z[i] - cz; - const d = dx * dx + dy * dy + dz * dz; - if (d > radiusSq) radiusSq = d; - } + const size = OrderedSet.size(indices); + for (let t = 0; t < size; t++) { + const i = OrderedSet.getAt(indices, t); + cx += x[i]; + cy += y[i]; + cz += z[i]; + } + + if (size > 0) { + cx /= size; + cy /= size; + cz /= size; + } + + for (let t = 0; t < size; t++) { + const i = OrderedSet.getAt(indices, t); + const dx = x[i] - cx, dy = y[i] - cy, dz = z[i] - cz; + const d = dx * dx + dy * dy + dz * dz; + if (d > radiusSq) radiusSq = d; } return { center: Vec3.create(cx, cy, cz), radius: Math.sqrt(radiusSq) }; diff --git a/src/mol-math/geometry/symmetry-operator.ts b/src/mol-math/geometry/symmetry-operator.ts index dcc83ade5..6b65fafd0 100644 --- a/src/mol-math/geometry/symmetry-operator.ts +++ b/src/mol-math/geometry/symmetry-operator.ts @@ -61,7 +61,7 @@ namespace SymmetryOperator { } } -export default SymmetryOperator +export { SymmetryOperator } interface Projections { x(index: number): number, y(index: number): number, z(index: number): number } diff --git a/src/mol-model/structure/model/formats/mmcif/assembly.ts b/src/mol-model/structure/model/formats/mmcif/assembly.ts index 5d83ca871..a03ff07ae 100644 --- a/src/mol-model/structure/model/formats/mmcif/assembly.ts +++ b/src/mol-model/structure/model/formats/mmcif/assembly.ts @@ -5,7 +5,7 @@ */ import { Mat4, Tensor } from 'mol-math/linear-algebra' -import SymmetryOperator from 'mol-math/geometry/symmetry-operator' +import { SymmetryOperator } from 'mol-math/geometry/symmetry-operator' import Format from '../../format' import { Assembly, OperatorGroup, OperatorGroups } from '../../properties/symmetry' import { Queries as Q } from '../../../query' diff --git a/src/mol-model/structure/model/model.ts b/src/mol-model/structure/model/model.ts index e6ed614f4..99a3877aa 100644 --- a/src/mol-model/structure/model/model.ts +++ b/src/mol-model/structure/model/model.ts @@ -5,16 +5,12 @@ */ import UUID from 'mol-util/uuid' -import GridLookup from 'mol-math/geometry/grid-lookup' import Format from './format' import Hierarchy from './properties/hierarchy' import Conformation from './properties/conformation' import Symmetry from './properties/symmetry' -import Bonds from './properties/bonds' import CoarseGrained from './properties/coarse-grained' -import computeBonds from './utils/compute-bonds' - import from_gro from './formats/gro' import from_mmCIF from './formats/mmcif' @@ -38,8 +34,7 @@ interface Model extends Readonly<{ atomCount: number, }> { - '@spatialLookup'?: GridLookup, - '@bonds'?: Bonds + } { } namespace Model { @@ -49,18 +44,18 @@ namespace Model { case 'mmCIF': return from_mmCIF(format); } } - export function spatialLookup(model: Model): GridLookup { - if (model['@spatialLookup']) return model['@spatialLookup']!; - const lookup = GridLookup(model.conformation); - model['@spatialLookup'] = lookup; - return lookup; - } - export function bonds(model: Model): Bonds { - if (model['@bonds']) return model['@bonds']!; - const bonds = computeBonds(model); - model['@bonds'] = bonds; - return bonds; - } + // export function spatialLookup(model: Model): GridLookup { + // if (model['@spatialLookup']) return model['@spatialLookup']!; + // const lookup = GridLookup(model.conformation); + // model['@spatialLookup'] = lookup; + // return lookup; + // } + // export function bonds(model: Model): Bonds { + // if (model['@bonds']) return model['@bonds']!; + // const bonds = computeBonds(model); + // model['@bonds'] = bonds; + // return bonds; + // } } export default Model \ No newline at end of file diff --git a/src/mol-model/structure/model/properties/symmetry.ts b/src/mol-model/structure/model/properties/symmetry.ts index dfe599c2a..d459241b4 100644 --- a/src/mol-model/structure/model/properties/symmetry.ts +++ b/src/mol-model/structure/model/properties/symmetry.ts @@ -4,7 +4,7 @@ * @author David Sehnal <david.sehnal@gmail.com> */ -import SymmetryOperator from 'mol-math/geometry/symmetry-operator' +import { SymmetryOperator } from 'mol-math/geometry/symmetry-operator' import { arrayFind } from 'mol-data/util' import { Query } from '../../query' import { Model } from '../../model' diff --git a/src/mol-model/structure/model/utils/compute-bonds.ts b/src/mol-model/structure/model/utils/compute-bonds.ts index 9b9170574..6db9047fd 100644 --- a/src/mol-model/structure/model/utils/compute-bonds.ts +++ b/src/mol-model/structure/model/utils/compute-bonds.ts @@ -1,253 +1,253 @@ -/** - * Copyright (c) 2017 MolQL contributors, licensed under MIT, See LICENSE file for more info. - * - * @author David Sehnal <david.sehnal@gmail.com> - */ - -import Mask from 'mol-util/mask' -import Model from '../model' -import { BondType, ElementSymbol } from '../types' -import Bonds from '../properties/bonds' -import { StructConn, ComponentBondInfo } from '../formats/mmcif/bonds' - -export interface BondComputationParameters { - maxHbondLength: number, - forceCompute: boolean -} - -// H,D,T are all mapped to H -const __ElementIndex: { [e: string]: number | undefined } = { 'H': 0, 'h': 0, 'D': 0, 'd': 0, 'T': 0, 't': 0, 'He': 2, 'HE': 2, 'he': 2, 'Li': 3, 'LI': 3, 'li': 3, 'Be': 4, 'BE': 4, 'be': 4, 'B': 5, 'b': 5, 'C': 6, 'c': 6, 'N': 7, 'n': 7, 'O': 8, 'o': 8, 'F': 9, 'f': 9, 'Ne': 10, 'NE': 10, 'ne': 10, 'Na': 11, 'NA': 11, 'na': 11, 'Mg': 12, 'MG': 12, 'mg': 12, 'Al': 13, 'AL': 13, 'al': 13, 'Si': 14, 'SI': 14, 'si': 14, 'P': 15, 'p': 15, 'S': 16, 's': 16, 'Cl': 17, 'CL': 17, 'cl': 17, 'Ar': 18, 'AR': 18, 'ar': 18, 'K': 19, 'k': 19, 'Ca': 20, 'CA': 20, 'ca': 20, 'Sc': 21, 'SC': 21, 'sc': 21, 'Ti': 22, 'TI': 22, 'ti': 22, 'V': 23, 'v': 23, 'Cr': 24, 'CR': 24, 'cr': 24, 'Mn': 25, 'MN': 25, 'mn': 25, 'Fe': 26, 'FE': 26, 'fe': 26, 'Co': 27, 'CO': 27, 'co': 27, 'Ni': 28, 'NI': 28, 'ni': 28, 'Cu': 29, 'CU': 29, 'cu': 29, 'Zn': 30, 'ZN': 30, 'zn': 30, 'Ga': 31, 'GA': 31, 'ga': 31, 'Ge': 32, 'GE': 32, 'ge': 32, 'As': 33, 'AS': 33, 'as': 33, 'Se': 34, 'SE': 34, 'se': 34, 'Br': 35, 'BR': 35, 'br': 35, 'Kr': 36, 'KR': 36, 'kr': 36, 'Rb': 37, 'RB': 37, 'rb': 37, 'Sr': 38, 'SR': 38, 'sr': 38, 'Y': 39, 'y': 39, 'Zr': 40, 'ZR': 40, 'zr': 40, 'Nb': 41, 'NB': 41, 'nb': 41, 'Mo': 42, 'MO': 42, 'mo': 42, 'Tc': 43, 'TC': 43, 'tc': 43, 'Ru': 44, 'RU': 44, 'ru': 44, 'Rh': 45, 'RH': 45, 'rh': 45, 'Pd': 46, 'PD': 46, 'pd': 46, 'Ag': 47, 'AG': 47, 'ag': 47, 'Cd': 48, 'CD': 48, 'cd': 48, 'In': 49, 'IN': 49, 'in': 49, 'Sn': 50, 'SN': 50, 'sn': 50, 'Sb': 51, 'SB': 51, 'sb': 51, 'Te': 52, 'TE': 52, 'te': 52, 'I': 53, 'i': 53, 'Xe': 54, 'XE': 54, 'xe': 54, 'Cs': 55, 'CS': 55, 'cs': 55, 'Ba': 56, 'BA': 56, 'ba': 56, 'La': 57, 'LA': 57, 'la': 57, 'Ce': 58, 'CE': 58, 'ce': 58, 'Pr': 59, 'PR': 59, 'pr': 59, 'Nd': 60, 'ND': 60, 'nd': 60, 'Pm': 61, 'PM': 61, 'pm': 61, 'Sm': 62, 'SM': 62, 'sm': 62, 'Eu': 63, 'EU': 63, 'eu': 63, 'Gd': 64, 'GD': 64, 'gd': 64, 'Tb': 65, 'TB': 65, 'tb': 65, 'Dy': 66, 'DY': 66, 'dy': 66, 'Ho': 67, 'HO': 67, 'ho': 67, 'Er': 68, 'ER': 68, 'er': 68, 'Tm': 69, 'TM': 69, 'tm': 69, 'Yb': 70, 'YB': 70, 'yb': 70, 'Lu': 71, 'LU': 71, 'lu': 71, 'Hf': 72, 'HF': 72, 'hf': 72, 'Ta': 73, 'TA': 73, 'ta': 73, 'W': 74, 'w': 74, 'Re': 75, 'RE': 75, 're': 75, 'Os': 76, 'OS': 76, 'os': 76, 'Ir': 77, 'IR': 77, 'ir': 77, 'Pt': 78, 'PT': 78, 'pt': 78, 'Au': 79, 'AU': 79, 'au': 79, 'Hg': 80, 'HG': 80, 'hg': 80, 'Tl': 81, 'TL': 81, 'tl': 81, 'Pb': 82, 'PB': 82, 'pb': 82, 'Bi': 83, 'BI': 83, 'bi': 83, 'Po': 84, 'PO': 84, 'po': 84, 'At': 85, 'AT': 85, 'at': 85, 'Rn': 86, 'RN': 86, 'rn': 86, 'Fr': 87, 'FR': 87, 'fr': 87, 'Ra': 88, 'RA': 88, 'ra': 88, 'Ac': 89, 'AC': 89, 'ac': 89, 'Th': 90, 'TH': 90, 'th': 90, 'Pa': 91, 'PA': 91, 'pa': 91, 'U': 92, 'u': 92, 'Np': 93, 'NP': 93, 'np': 93, 'Pu': 94, 'PU': 94, 'pu': 94, 'Am': 95, 'AM': 95, 'am': 95, 'Cm': 96, 'CM': 96, 'cm': 96, 'Bk': 97, 'BK': 97, 'bk': 97, 'Cf': 98, 'CF': 98, 'cf': 98, 'Es': 99, 'ES': 99, 'es': 99, 'Fm': 100, 'FM': 100, 'fm': 100, 'Md': 101, 'MD': 101, 'md': 101, 'No': 102, 'NO': 102, 'no': 102, 'Lr': 103, 'LR': 103, 'lr': 103, 'Rf': 104, 'RF': 104, 'rf': 104, 'Db': 105, 'DB': 105, 'db': 105, 'Sg': 106, 'SG': 106, 'sg': 106, 'Bh': 107, 'BH': 107, 'bh': 107, 'Hs': 108, 'HS': 108, 'hs': 108, 'Mt': 109, 'MT': 109, 'mt': 109 }; - -const __ElementBondThresholds: { [e: number]: number | undefined } = { 0: 1.42, 1: 1.42, 3: 2.7, 4: 2.7, 6: 1.75, 7: 1.6, 8: 1.52, 11: 2.7, 12: 2.7, 13: 2.7, 14: 1.9, 15: 1.9, 16: 1.9, 17: 1.8, 19: 2.7, 20: 2.7, 21: 2.7, 22: 2.7, 23: 2.7, 24: 2.7, 25: 2.7, 26: 2.7, 27: 2.7, 28: 2.7, 29: 2.7, 30: 2.7, 31: 2.7, 33: 2.68, 37: 2.7, 38: 2.7, 39: 2.7, 40: 2.7, 41: 2.7, 42: 2.7, 43: 2.7, 44: 2.7, 45: 2.7, 46: 2.7, 47: 2.7, 48: 2.7, 49: 2.7, 50: 2.7, 55: 2.7, 56: 2.7, 57: 2.7, 58: 2.7, 59: 2.7, 60: 2.7, 61: 2.7, 62: 2.7, 63: 2.7, 64: 2.7, 65: 2.7, 66: 2.7, 67: 2.7, 68: 2.7, 69: 2.7, 70: 2.7, 71: 2.7, 72: 2.7, 73: 2.7, 74: 2.7, 75: 2.7, 76: 2.7, 77: 2.7, 78: 2.7, 79: 2.7, 80: 2.7, 81: 2.7, 82: 2.7, 83: 2.7, 87: 2.7, 88: 2.7, 89: 2.7, 90: 2.7, 91: 2.7, 92: 2.7, 93: 2.7, 94: 2.7, 95: 2.7, 96: 2.7, 97: 2.7, 98: 2.7, 99: 2.7, 100: 2.7, 101: 2.7, 102: 2.7, 103: 2.7, 104: 2.7, 105: 2.7, 106: 2.7, 107: 2.7, 108: 2.7, 109: 2.88 }; - -const __ElementPairThresholds: { [e: number]: number | undefined } = { 0: 0.8, 20: 1.31, 27: 1.3, 35: 1.3, 44: 1.05, 54: 1, 60: 1.84, 72: 1.88, 84: 1.75, 85: 1.56, 86: 1.76, 98: 1.6, 99: 1.68, 100: 1.63, 112: 1.55, 113: 1.59, 114: 1.36, 129: 1.45, 144: 1.6, 170: 1.4, 180: 1.55, 202: 2.4, 222: 2.24, 224: 1.91, 225: 1.98, 243: 2.02, 269: 2, 293: 1.9, 480: 2.3, 512: 2.3, 544: 2.3, 612: 2.1, 629: 1.54, 665: 1, 813: 2.6, 854: 2.27, 894: 1.93, 896: 2.1, 937: 2.05, 938: 2.06, 981: 1.62, 1258: 2.68, 1309: 2.33, 1484: 1, 1763: 2.14, 1823: 2.48, 1882: 2.1, 1944: 1.72, 2380: 2.34, 3367: 2.44, 3733: 2.11, 3819: 2.6, 3821: 2.36, 4736: 2.75, 5724: 2.73, 5959: 2.63, 6519: 2.84, 6750: 2.87, 8991: 2.81 }; - -const __DefaultBondingRadius = 2.001; - -const MetalsSet = (function () { - const metals = ['LI', 'NA', 'K', 'RB', 'CS', 'FR', 'BE', 'MG', 'CA', 'SR', 'BA', 'RA', 'AL', 'GA', 'IN', 'SN', 'TL', 'PB', 'BI', 'SC', 'TI', 'V', 'CR', 'MN', 'FE', 'CO', 'NI', 'CU', 'ZN', 'Y', 'ZR', 'NB', 'MO', 'TC', 'RU', 'RH', 'PD', 'AG', 'CD', 'LA', 'HF', 'TA', 'W', 'RE', 'OS', 'IR', 'PT', 'AU', 'HG', 'AC', 'RF', 'DB', 'SG', 'BH', 'HS', 'MT', 'CE', 'PR', 'ND', 'PM', 'SM', 'EU', 'GD', 'TB', 'DY', 'HO', 'ER', 'TM', 'YB', 'LU', 'TH', 'PA', 'U', 'NP', 'PU', 'AM', 'CM', 'BK', 'CF', 'ES', 'FM', 'MD', 'NO', 'LR']; - const set = new Set<number>(); - for (const m of metals) { - set.add(__ElementIndex[m]!); - } - return set; -})(); - -function pair(a: number, b: number) { - if (a < b) return (a + b) * (a + b + 1) / 2 + b; - else return (a + b) * (a + b + 1) / 2 + a; -} - -function idx(e: ElementSymbol) { - const i = __ElementIndex[e as any as string]; - if (i === void 0) return -1; - return i; -} - -function pairThreshold(i: number, j: number) { - if (i < 0 || j < 0) return -1; - const r = __ElementPairThresholds[pair(i, j)]; - if (r === void 0) return -1; - return r; -} - -function threshold(i: number) { - if (i < 0) return __DefaultBondingRadius; - const r = __ElementBondThresholds[i]; - if (r === void 0) return __DefaultBondingRadius; - return r; -} - -const H_ID = __ElementIndex['H']!; -function isHydrogen(i: number) { - return i === H_ID; -} - -function computePerAtomBonds(atomA: number[], atomB: number[], _order: number[], _flags: number[], atomCount: number) { - const bucketSizes = new Int32Array(atomCount); - const bucketOffsets = new Int32Array(atomCount + 1) as any as number[]; - const bucketFill = new Int32Array(atomCount); - - for (const i of atomA) bucketSizes[i]++; - for (const i of atomB) bucketSizes[i]++; - - let offset = 0; - for (let i = 0; i < atomCount; i++) { - bucketOffsets[i] = offset; - offset += bucketSizes[i]; - } - bucketOffsets[atomCount] = offset; - - const neighbor = new Int32Array(offset) as any as number[]; - const flags = new Uint16Array(offset) as any as number[]; - const order = new Int8Array(offset) as any as number[]; - - for (let i = 0, _i = atomA.length; i < _i; i++) { - const a = atomA[i], b = atomB[i], f = _flags[i], o = _order[i]; - - const oa = bucketOffsets[a] + bucketFill[a]; - const ob = bucketOffsets[b] + bucketFill[b]; - - neighbor[oa] = b; - flags[oa] = f; - order[oa] = o; - bucketFill[a]++; - - neighbor[ob] = a; - flags[ob] = f; - order[ob] = o; - bucketFill[b]++; - } - - return { - offsets: bucketOffsets, - neighbor, - flags, - order - }; -} - -function _computeBonds(model: Model, params: BondComputationParameters): Bonds { - const MAX_RADIUS = 3; - - const { x, y, z } = model.conformation; - const { atomCount } = model; - const { residueKey } = model.hierarchy - const { type_symbol, label_atom_id, label_alt_id } = model.hierarchy.atoms; - const { label_comp_id } = model.hierarchy.residues; - const query3d = Model.spatialLookup(model).find(Mask.always(atomCount)); - - const structConn = model.sourceData.kind === 'mmCIF' ? StructConn.create(model) : void 0 - const component = model.sourceData.kind === 'mmCIF' ? ComponentBondInfo.create(model) : void 0 - - const atomA: number[] = []; - const atomB: number[] = []; - const flags: number[] = []; - const order: number[] = []; - - let lastResidue = -1; - let componentMap: Map<string, Map<string, { flags: number, order: number }>> | undefined = void 0; - - for (let aI = 0; aI < atomCount; aI++) { - const raI = residueKey.value(aI); - // const rowA = dataIndex[aI]; // TODO - const rowA = aI; - - if (!params.forceCompute && raI !== lastResidue) { - const resn = label_comp_id.value(rowA)!; - if (!!component && component.entries.has(resn)) { - componentMap = component.entries.get(resn)!.map; - } else { - componentMap = void 0; - } - } - lastResidue = raI; - - const componentPairs = componentMap ? componentMap.get(label_atom_id.value(rowA)) : void 0; - - const aeI = idx(type_symbol.value(rowA)!); - - const { indices, count, squaredDistances } = query3d(x[aI], y[aI], z[aI], MAX_RADIUS); - const isHa = isHydrogen(aeI); - const thresholdA = threshold(aeI); - const altA = label_alt_id.value(rowA); - const metalA = MetalsSet.has(aeI); - const structConnEntries = params.forceCompute ? void 0 : structConn && structConn.getAtomEntries(aI); - - for (let ni = 0; ni < count; ni++) { - const bI = indices[ni]; - if (bI <= aI) continue; - - // const rowB = dataIndex[bI]; // TODO - const rowB = bI; - - const altB = label_alt_id.value(rowB); - if (altA && altB && altA !== altB) continue; - - const beI = idx(type_symbol.value(rowB)!); - const isMetal = metalA || MetalsSet.has(beI); - - const rbI = residueKey.value(bI); - // handle "component dictionary" bonds. - if (raI === rbI && componentPairs) { - const e = componentPairs.get(label_atom_id.value(rowB)!); - if (e) { - atomA[atomA.length] = aI; - atomB[atomB.length] = bI; - order[order.length] = e.order; - let flag = e.flags; - if (isMetal) { - if (flag | BondType.Flag.Covalent) flag ^= BondType.Flag.Covalent; - flag |= BondType.Flag.MetallicCoordination; - } - flags[flags.length] = flag; - } - continue; - } - - const isHb = isHydrogen(beI); - if (isHa && isHb) continue; - - const dist = Math.sqrt(squaredDistances[ni]); - if (dist === 0) continue; - - // handle "struct conn" bonds. - if (structConnEntries && structConnEntries.length) { - let added = false; - for (const se of structConnEntries) { - for (const p of se.partners) { - if (p.atomIndex === bI) { - atomA[atomA.length] = aI; - atomB[atomB.length] = bI; - flags[flags.length] = se.flags; - order[order.length] = se.order; - added = true; - break; - } - } - if (added) break; - } - if (added) continue; - } - - if (isHa || isHb) { - if (dist < params.maxHbondLength) { - atomA[atomA.length] = aI; - atomB[atomB.length] = bI; - order[order.length] = 1; - flags[flags.length] = BondType.Flag.Covalent | BondType.Flag.Computed; // TODO: check if correct - } - continue; - } - - const thresholdAB = pairThreshold(aeI, beI); - const pairingThreshold = thresholdAB > 0 - ? thresholdAB - : beI < 0 ? thresholdA : Math.max(thresholdA, threshold(beI)); - - - if (dist <= pairingThreshold) { - atomA[atomA.length] = aI; - atomB[atomB.length] = bI; - order[order.length] = 1; - flags[flags.length] = (isMetal ? BondType.Flag.MetallicCoordination : BondType.Flag.Covalent) | BondType.Flag.Computed; - } - } - } - - const bonds = computePerAtomBonds(atomA, atomB, order, flags, atomCount); - return { - offset: bonds.offsets, - neighbor: bonds.neighbor, - flags: bonds.flags, - order: bonds.order, - count: atomA.length - }; -} - -export default function computeBonds(model: Model, params?: Partial<BondComputationParameters>) { - return _computeBonds(model, { - maxHbondLength: (params && params.maxHbondLength) || 1.15, - forceCompute: !!(params && params.forceCompute), - }); -} \ No newline at end of file +// /** +// * Copyright (c) 2017 MolQL contributors, licensed under MIT, See LICENSE file for more info. +// * +// * @author David Sehnal <david.sehnal@gmail.com> +// */ + +// import Mask from 'mol-util/mask' +// import Model from '../model' +// import { BondType, ElementSymbol } from '../types' +// import Bonds from '../properties/bonds' +// import { StructConn, ComponentBondInfo } from '../formats/mmcif/bonds' + +// export interface BondComputationParameters { +// maxHbondLength: number, +// forceCompute: boolean +// } + +// // H,D,T are all mapped to H +// const __ElementIndex: { [e: string]: number | undefined } = { 'H': 0, 'h': 0, 'D': 0, 'd': 0, 'T': 0, 't': 0, 'He': 2, 'HE': 2, 'he': 2, 'Li': 3, 'LI': 3, 'li': 3, 'Be': 4, 'BE': 4, 'be': 4, 'B': 5, 'b': 5, 'C': 6, 'c': 6, 'N': 7, 'n': 7, 'O': 8, 'o': 8, 'F': 9, 'f': 9, 'Ne': 10, 'NE': 10, 'ne': 10, 'Na': 11, 'NA': 11, 'na': 11, 'Mg': 12, 'MG': 12, 'mg': 12, 'Al': 13, 'AL': 13, 'al': 13, 'Si': 14, 'SI': 14, 'si': 14, 'P': 15, 'p': 15, 'S': 16, 's': 16, 'Cl': 17, 'CL': 17, 'cl': 17, 'Ar': 18, 'AR': 18, 'ar': 18, 'K': 19, 'k': 19, 'Ca': 20, 'CA': 20, 'ca': 20, 'Sc': 21, 'SC': 21, 'sc': 21, 'Ti': 22, 'TI': 22, 'ti': 22, 'V': 23, 'v': 23, 'Cr': 24, 'CR': 24, 'cr': 24, 'Mn': 25, 'MN': 25, 'mn': 25, 'Fe': 26, 'FE': 26, 'fe': 26, 'Co': 27, 'CO': 27, 'co': 27, 'Ni': 28, 'NI': 28, 'ni': 28, 'Cu': 29, 'CU': 29, 'cu': 29, 'Zn': 30, 'ZN': 30, 'zn': 30, 'Ga': 31, 'GA': 31, 'ga': 31, 'Ge': 32, 'GE': 32, 'ge': 32, 'As': 33, 'AS': 33, 'as': 33, 'Se': 34, 'SE': 34, 'se': 34, 'Br': 35, 'BR': 35, 'br': 35, 'Kr': 36, 'KR': 36, 'kr': 36, 'Rb': 37, 'RB': 37, 'rb': 37, 'Sr': 38, 'SR': 38, 'sr': 38, 'Y': 39, 'y': 39, 'Zr': 40, 'ZR': 40, 'zr': 40, 'Nb': 41, 'NB': 41, 'nb': 41, 'Mo': 42, 'MO': 42, 'mo': 42, 'Tc': 43, 'TC': 43, 'tc': 43, 'Ru': 44, 'RU': 44, 'ru': 44, 'Rh': 45, 'RH': 45, 'rh': 45, 'Pd': 46, 'PD': 46, 'pd': 46, 'Ag': 47, 'AG': 47, 'ag': 47, 'Cd': 48, 'CD': 48, 'cd': 48, 'In': 49, 'IN': 49, 'in': 49, 'Sn': 50, 'SN': 50, 'sn': 50, 'Sb': 51, 'SB': 51, 'sb': 51, 'Te': 52, 'TE': 52, 'te': 52, 'I': 53, 'i': 53, 'Xe': 54, 'XE': 54, 'xe': 54, 'Cs': 55, 'CS': 55, 'cs': 55, 'Ba': 56, 'BA': 56, 'ba': 56, 'La': 57, 'LA': 57, 'la': 57, 'Ce': 58, 'CE': 58, 'ce': 58, 'Pr': 59, 'PR': 59, 'pr': 59, 'Nd': 60, 'ND': 60, 'nd': 60, 'Pm': 61, 'PM': 61, 'pm': 61, 'Sm': 62, 'SM': 62, 'sm': 62, 'Eu': 63, 'EU': 63, 'eu': 63, 'Gd': 64, 'GD': 64, 'gd': 64, 'Tb': 65, 'TB': 65, 'tb': 65, 'Dy': 66, 'DY': 66, 'dy': 66, 'Ho': 67, 'HO': 67, 'ho': 67, 'Er': 68, 'ER': 68, 'er': 68, 'Tm': 69, 'TM': 69, 'tm': 69, 'Yb': 70, 'YB': 70, 'yb': 70, 'Lu': 71, 'LU': 71, 'lu': 71, 'Hf': 72, 'HF': 72, 'hf': 72, 'Ta': 73, 'TA': 73, 'ta': 73, 'W': 74, 'w': 74, 'Re': 75, 'RE': 75, 're': 75, 'Os': 76, 'OS': 76, 'os': 76, 'Ir': 77, 'IR': 77, 'ir': 77, 'Pt': 78, 'PT': 78, 'pt': 78, 'Au': 79, 'AU': 79, 'au': 79, 'Hg': 80, 'HG': 80, 'hg': 80, 'Tl': 81, 'TL': 81, 'tl': 81, 'Pb': 82, 'PB': 82, 'pb': 82, 'Bi': 83, 'BI': 83, 'bi': 83, 'Po': 84, 'PO': 84, 'po': 84, 'At': 85, 'AT': 85, 'at': 85, 'Rn': 86, 'RN': 86, 'rn': 86, 'Fr': 87, 'FR': 87, 'fr': 87, 'Ra': 88, 'RA': 88, 'ra': 88, 'Ac': 89, 'AC': 89, 'ac': 89, 'Th': 90, 'TH': 90, 'th': 90, 'Pa': 91, 'PA': 91, 'pa': 91, 'U': 92, 'u': 92, 'Np': 93, 'NP': 93, 'np': 93, 'Pu': 94, 'PU': 94, 'pu': 94, 'Am': 95, 'AM': 95, 'am': 95, 'Cm': 96, 'CM': 96, 'cm': 96, 'Bk': 97, 'BK': 97, 'bk': 97, 'Cf': 98, 'CF': 98, 'cf': 98, 'Es': 99, 'ES': 99, 'es': 99, 'Fm': 100, 'FM': 100, 'fm': 100, 'Md': 101, 'MD': 101, 'md': 101, 'No': 102, 'NO': 102, 'no': 102, 'Lr': 103, 'LR': 103, 'lr': 103, 'Rf': 104, 'RF': 104, 'rf': 104, 'Db': 105, 'DB': 105, 'db': 105, 'Sg': 106, 'SG': 106, 'sg': 106, 'Bh': 107, 'BH': 107, 'bh': 107, 'Hs': 108, 'HS': 108, 'hs': 108, 'Mt': 109, 'MT': 109, 'mt': 109 }; + +// const __ElementBondThresholds: { [e: number]: number | undefined } = { 0: 1.42, 1: 1.42, 3: 2.7, 4: 2.7, 6: 1.75, 7: 1.6, 8: 1.52, 11: 2.7, 12: 2.7, 13: 2.7, 14: 1.9, 15: 1.9, 16: 1.9, 17: 1.8, 19: 2.7, 20: 2.7, 21: 2.7, 22: 2.7, 23: 2.7, 24: 2.7, 25: 2.7, 26: 2.7, 27: 2.7, 28: 2.7, 29: 2.7, 30: 2.7, 31: 2.7, 33: 2.68, 37: 2.7, 38: 2.7, 39: 2.7, 40: 2.7, 41: 2.7, 42: 2.7, 43: 2.7, 44: 2.7, 45: 2.7, 46: 2.7, 47: 2.7, 48: 2.7, 49: 2.7, 50: 2.7, 55: 2.7, 56: 2.7, 57: 2.7, 58: 2.7, 59: 2.7, 60: 2.7, 61: 2.7, 62: 2.7, 63: 2.7, 64: 2.7, 65: 2.7, 66: 2.7, 67: 2.7, 68: 2.7, 69: 2.7, 70: 2.7, 71: 2.7, 72: 2.7, 73: 2.7, 74: 2.7, 75: 2.7, 76: 2.7, 77: 2.7, 78: 2.7, 79: 2.7, 80: 2.7, 81: 2.7, 82: 2.7, 83: 2.7, 87: 2.7, 88: 2.7, 89: 2.7, 90: 2.7, 91: 2.7, 92: 2.7, 93: 2.7, 94: 2.7, 95: 2.7, 96: 2.7, 97: 2.7, 98: 2.7, 99: 2.7, 100: 2.7, 101: 2.7, 102: 2.7, 103: 2.7, 104: 2.7, 105: 2.7, 106: 2.7, 107: 2.7, 108: 2.7, 109: 2.88 }; + +// const __ElementPairThresholds: { [e: number]: number | undefined } = { 0: 0.8, 20: 1.31, 27: 1.3, 35: 1.3, 44: 1.05, 54: 1, 60: 1.84, 72: 1.88, 84: 1.75, 85: 1.56, 86: 1.76, 98: 1.6, 99: 1.68, 100: 1.63, 112: 1.55, 113: 1.59, 114: 1.36, 129: 1.45, 144: 1.6, 170: 1.4, 180: 1.55, 202: 2.4, 222: 2.24, 224: 1.91, 225: 1.98, 243: 2.02, 269: 2, 293: 1.9, 480: 2.3, 512: 2.3, 544: 2.3, 612: 2.1, 629: 1.54, 665: 1, 813: 2.6, 854: 2.27, 894: 1.93, 896: 2.1, 937: 2.05, 938: 2.06, 981: 1.62, 1258: 2.68, 1309: 2.33, 1484: 1, 1763: 2.14, 1823: 2.48, 1882: 2.1, 1944: 1.72, 2380: 2.34, 3367: 2.44, 3733: 2.11, 3819: 2.6, 3821: 2.36, 4736: 2.75, 5724: 2.73, 5959: 2.63, 6519: 2.84, 6750: 2.87, 8991: 2.81 }; + +// const __DefaultBondingRadius = 2.001; + +// const MetalsSet = (function () { +// const metals = ['LI', 'NA', 'K', 'RB', 'CS', 'FR', 'BE', 'MG', 'CA', 'SR', 'BA', 'RA', 'AL', 'GA', 'IN', 'SN', 'TL', 'PB', 'BI', 'SC', 'TI', 'V', 'CR', 'MN', 'FE', 'CO', 'NI', 'CU', 'ZN', 'Y', 'ZR', 'NB', 'MO', 'TC', 'RU', 'RH', 'PD', 'AG', 'CD', 'LA', 'HF', 'TA', 'W', 'RE', 'OS', 'IR', 'PT', 'AU', 'HG', 'AC', 'RF', 'DB', 'SG', 'BH', 'HS', 'MT', 'CE', 'PR', 'ND', 'PM', 'SM', 'EU', 'GD', 'TB', 'DY', 'HO', 'ER', 'TM', 'YB', 'LU', 'TH', 'PA', 'U', 'NP', 'PU', 'AM', 'CM', 'BK', 'CF', 'ES', 'FM', 'MD', 'NO', 'LR']; +// const set = new Set<number>(); +// for (const m of metals) { +// set.add(__ElementIndex[m]!); +// } +// return set; +// })(); + +// function pair(a: number, b: number) { +// if (a < b) return (a + b) * (a + b + 1) / 2 + b; +// else return (a + b) * (a + b + 1) / 2 + a; +// } + +// function idx(e: ElementSymbol) { +// const i = __ElementIndex[e as any as string]; +// if (i === void 0) return -1; +// return i; +// } + +// function pairThreshold(i: number, j: number) { +// if (i < 0 || j < 0) return -1; +// const r = __ElementPairThresholds[pair(i, j)]; +// if (r === void 0) return -1; +// return r; +// } + +// function threshold(i: number) { +// if (i < 0) return __DefaultBondingRadius; +// const r = __ElementBondThresholds[i]; +// if (r === void 0) return __DefaultBondingRadius; +// return r; +// } + +// const H_ID = __ElementIndex['H']!; +// function isHydrogen(i: number) { +// return i === H_ID; +// } + +// function computePerAtomBonds(atomA: number[], atomB: number[], _order: number[], _flags: number[], atomCount: number) { +// const bucketSizes = new Int32Array(atomCount); +// const bucketOffsets = new Int32Array(atomCount + 1) as any as number[]; +// const bucketFill = new Int32Array(atomCount); + +// for (const i of atomA) bucketSizes[i]++; +// for (const i of atomB) bucketSizes[i]++; + +// let offset = 0; +// for (let i = 0; i < atomCount; i++) { +// bucketOffsets[i] = offset; +// offset += bucketSizes[i]; +// } +// bucketOffsets[atomCount] = offset; + +// const neighbor = new Int32Array(offset) as any as number[]; +// const flags = new Uint16Array(offset) as any as number[]; +// const order = new Int8Array(offset) as any as number[]; + +// for (let i = 0, _i = atomA.length; i < _i; i++) { +// const a = atomA[i], b = atomB[i], f = _flags[i], o = _order[i]; + +// const oa = bucketOffsets[a] + bucketFill[a]; +// const ob = bucketOffsets[b] + bucketFill[b]; + +// neighbor[oa] = b; +// flags[oa] = f; +// order[oa] = o; +// bucketFill[a]++; + +// neighbor[ob] = a; +// flags[ob] = f; +// order[ob] = o; +// bucketFill[b]++; +// } + +// return { +// offsets: bucketOffsets, +// neighbor, +// flags, +// order +// }; +// } + +// function _computeBonds(model: Model, params: BondComputationParameters): Bonds { +// const MAX_RADIUS = 3; + +// const { x, y, z } = model.conformation; +// const { atomCount } = model; +// const { residueKey } = model.hierarchy +// const { type_symbol, label_atom_id, label_alt_id } = model.hierarchy.atoms; +// const { label_comp_id } = model.hierarchy.residues; +// const query3d = Model.spatialLookup(model).find(Mask.always(atomCount)); + +// const structConn = model.sourceData.kind === 'mmCIF' ? StructConn.create(model) : void 0 +// const component = model.sourceData.kind === 'mmCIF' ? ComponentBondInfo.create(model) : void 0 + +// const atomA: number[] = []; +// const atomB: number[] = []; +// const flags: number[] = []; +// const order: number[] = []; + +// let lastResidue = -1; +// let componentMap: Map<string, Map<string, { flags: number, order: number }>> | undefined = void 0; + +// for (let aI = 0; aI < atomCount; aI++) { +// const raI = residueKey.value(aI); +// // const rowA = dataIndex[aI]; // TODO +// const rowA = aI; + +// if (!params.forceCompute && raI !== lastResidue) { +// const resn = label_comp_id.value(rowA)!; +// if (!!component && component.entries.has(resn)) { +// componentMap = component.entries.get(resn)!.map; +// } else { +// componentMap = void 0; +// } +// } +// lastResidue = raI; + +// const componentPairs = componentMap ? componentMap.get(label_atom_id.value(rowA)) : void 0; + +// const aeI = idx(type_symbol.value(rowA)!); + +// const { indices, count, squaredDistances } = query3d(x[aI], y[aI], z[aI], MAX_RADIUS); +// const isHa = isHydrogen(aeI); +// const thresholdA = threshold(aeI); +// const altA = label_alt_id.value(rowA); +// const metalA = MetalsSet.has(aeI); +// const structConnEntries = params.forceCompute ? void 0 : structConn && structConn.getAtomEntries(aI); + +// for (let ni = 0; ni < count; ni++) { +// const bI = indices[ni]; +// if (bI <= aI) continue; + +// // const rowB = dataIndex[bI]; // TODO +// const rowB = bI; + +// const altB = label_alt_id.value(rowB); +// if (altA && altB && altA !== altB) continue; + +// const beI = idx(type_symbol.value(rowB)!); +// const isMetal = metalA || MetalsSet.has(beI); + +// const rbI = residueKey.value(bI); +// // handle "component dictionary" bonds. +// if (raI === rbI && componentPairs) { +// const e = componentPairs.get(label_atom_id.value(rowB)!); +// if (e) { +// atomA[atomA.length] = aI; +// atomB[atomB.length] = bI; +// order[order.length] = e.order; +// let flag = e.flags; +// if (isMetal) { +// if (flag | BondType.Flag.Covalent) flag ^= BondType.Flag.Covalent; +// flag |= BondType.Flag.MetallicCoordination; +// } +// flags[flags.length] = flag; +// } +// continue; +// } + +// const isHb = isHydrogen(beI); +// if (isHa && isHb) continue; + +// const dist = Math.sqrt(squaredDistances[ni]); +// if (dist === 0) continue; + +// // handle "struct conn" bonds. +// if (structConnEntries && structConnEntries.length) { +// let added = false; +// for (const se of structConnEntries) { +// for (const p of se.partners) { +// if (p.atomIndex === bI) { +// atomA[atomA.length] = aI; +// atomB[atomB.length] = bI; +// flags[flags.length] = se.flags; +// order[order.length] = se.order; +// added = true; +// break; +// } +// } +// if (added) break; +// } +// if (added) continue; +// } + +// if (isHa || isHb) { +// if (dist < params.maxHbondLength) { +// atomA[atomA.length] = aI; +// atomB[atomB.length] = bI; +// order[order.length] = 1; +// flags[flags.length] = BondType.Flag.Covalent | BondType.Flag.Computed; // TODO: check if correct +// } +// continue; +// } + +// const thresholdAB = pairThreshold(aeI, beI); +// const pairingThreshold = thresholdAB > 0 +// ? thresholdAB +// : beI < 0 ? thresholdA : Math.max(thresholdA, threshold(beI)); + + +// if (dist <= pairingThreshold) { +// atomA[atomA.length] = aI; +// atomB[atomB.length] = bI; +// order[order.length] = 1; +// flags[flags.length] = (isMetal ? BondType.Flag.MetallicCoordination : BondType.Flag.Covalent) | BondType.Flag.Computed; +// } +// } +// } + +// const bonds = computePerAtomBonds(atomA, atomB, order, flags, atomCount); +// return { +// offset: bonds.offsets, +// neighbor: bonds.neighbor, +// flags: bonds.flags, +// order: bonds.order, +// count: atomA.length +// }; +// } + +// export default function computeBonds(model: Model, params?: Partial<BondComputationParameters>) { +// return _computeBonds(model, { +// maxHbondLength: (params && params.maxHbondLength) || 1.15, +// forceCompute: !!(params && params.forceCompute), +// }); +// } \ No newline at end of file diff --git a/src/mol-model/structure/structure/element/group.ts b/src/mol-model/structure/structure/element/group.ts index ff542066d..ec7544470 100644 --- a/src/mol-model/structure/structure/element/group.ts +++ b/src/mol-model/structure/structure/element/group.ts @@ -5,12 +5,15 @@ */ import { OrderedSet } from 'mol-data/int' +import { Lookup3D, GridLookup3D } from 'mol-math/geometry'; import Unit from '../unit' interface ElementGroup { elements: OrderedSet, // Unique identifier of the group, usable as partial key for various "caches". - key: number + key: number, + + __lookup3d__?: Lookup3D } namespace ElementGroup { @@ -61,6 +64,17 @@ namespace ElementGroup { return createNew(set); } + export function getLookup3d(unit: Unit, group: ElementGroup) { + if (group.__lookup3d__) return group.__lookup3d__; + if (Unit.isAtomic(unit)) { + const { x, y, z } = unit.model.conformation; + group.__lookup3d__ = GridLookup3D({ x, y, z, indices: group.elements }); + return group.__lookup3d__; + } + + throw 'not implemented'; + } + let _id = 0; function nextKey() { const ret = _id; diff --git a/src/mol-model/structure/structure/structure.ts b/src/mol-model/structure/structure/structure.ts index b79ab2905..5f02951ce 100644 --- a/src/mol-model/structure/structure/structure.ts +++ b/src/mol-model/structure/structure/structure.ts @@ -6,7 +6,7 @@ import { OrderedSet, Iterator } from 'mol-data/int' import { UniqueArray } from 'mol-data/generic' -import SymmetryOperator from 'mol-math/geometry/symmetry-operator' +import { SymmetryOperator } from 'mol-math/geometry/symmetry-operator' import { Model, Format } from '../model' import Unit from './unit' import ElementSet from './element/set' diff --git a/src/mol-model/structure/structure/unit.ts b/src/mol-model/structure/structure/unit.ts index 73a503a2a..76accfeaf 100644 --- a/src/mol-model/structure/structure/unit.ts +++ b/src/mol-model/structure/structure/unit.ts @@ -4,7 +4,7 @@ * @author David Sehnal <david.sehnal@gmail.com> */ -import SymmetryOperator from 'mol-math/geometry/symmetry-operator' +import { SymmetryOperator } from 'mol-math/geometry/symmetry-operator' import ElementGroup from './element/group' import { Model } from '../model' diff --git a/src/perf-tests/lookup3d.ts b/src/perf-tests/lookup3d.ts index 548fa35a1..3ccd9483f 100644 --- a/src/perf-tests/lookup3d.ts +++ b/src/perf-tests/lookup3d.ts @@ -5,8 +5,9 @@ import CIF from 'mol-io/reader/cif' import { Structure, Model } from 'mol-model/structure' import { Run } from 'mol-task'; -import { build, inSphere } from 'mol-math/geometry/lookup3d/grid'; +import { GridLookup3D } from 'mol-math/geometry'; import { sortArray } from 'mol-data/util'; +import { OrderedSet } from 'mol-data/int'; require('util.promisify').shim(); const readFileAsync = util.promisify(fs.readFile); @@ -42,13 +43,14 @@ export async function readCIF(path: string) { export async function test() { const { mmcif } = await readCIF('e:/test/quick/1tqn_updated.cif'); - const data = build({ x: mmcif.atom_site.Cartn_x.toArray(), y: mmcif.atom_site.Cartn_y.toArray(), z: mmcif.atom_site.Cartn_z.toArray(), - //indices: [0, 1, 2, 3] + const lookup = GridLookup3D({ x: mmcif.atom_site.Cartn_x.toArray(), y: mmcif.atom_site.Cartn_y.toArray(), z: mmcif.atom_site.Cartn_z.toArray(), + indices: OrderedSet.ofBounds(0, 4), + radius: [1, 1, 1, 1] //indices: [1] }); - console.log(data.boundingBox, data.boundingSphere); + console.log(lookup.boundary.box, lookup.boundary.sphere); - const result = inSphere(data, -30.07, 8.178, -13.897, 10); + const result = lookup.find(-30.07, 8.178, -13.897, 1); console.log(result.count, sortArray(result.indices)); } -- GitLab