Skip to content
Snippets Groups Projects
Commit 5d022ffe authored by David Sehnal's avatar David Sehnal
Browse files

Fixes to group bond computation

parent 837f5eae
No related branches found
No related tags found
No related merge requests found
...@@ -10,7 +10,7 @@ require('util.promisify').shim(); ...@@ -10,7 +10,7 @@ require('util.promisify').shim();
// import { Table } from 'mol-data/db' // import { Table } from 'mol-data/db'
import CIF from 'mol-io/reader/cif' import CIF from 'mol-io/reader/cif'
import { Model, Structure, ElementSet, Unit } from 'mol-model/structure' import { Model, Structure, ElementSet, Unit, ElementGroup } from 'mol-model/structure'
import { Run, Progress } from 'mol-task' import { Run, Progress } from 'mol-task'
import { OrderedSet } from 'mol-data/int'; import { OrderedSet } from 'mol-data/int';
...@@ -22,7 +22,8 @@ async function parseCif(data: string|Uint8Array) { ...@@ -22,7 +22,8 @@ async function parseCif(data: string|Uint8Array) {
} }
async function getPdb(pdb: string) { async function getPdb(pdb: string) {
const data = await fetch(`https://files.rcsb.org/download/${pdb}.cif`) //const data = await fetch(`https://files.rcsb.org/download/${pdb}.cif`)
const data = await fetch(`http://www.ebi.ac.uk/pdbe/static/entry/${pdb}_updated.cif`);
const parsed = await parseCif(await data.text()) const parsed = await parseCif(await data.text())
return CIF.schema.mmCIF(parsed.result.blocks[0]) return CIF.schema.mmCIF(parsed.result.blocks[0])
} }
...@@ -39,7 +40,6 @@ export function atomLabel(model: Model, aI: number) { ...@@ -39,7 +40,6 @@ export function atomLabel(model: Model, aI: number) {
function printBonds(structure: Structure) { function printBonds(structure: Structure) {
const { units, elements } = structure; const { units, elements } = structure;
const unitIds = ElementSet.unitIndices(elements); const unitIds = ElementSet.unitIndices(elements);
...@@ -50,11 +50,18 @@ function printBonds(structure: Structure) { ...@@ -50,11 +50,18 @@ function printBonds(structure: Structure) {
const { count, offset, neighbor } = Unit.getGroupBonds(unit, group); const { count, offset, neighbor } = Unit.getGroupBonds(unit, group);
const { model } = unit; const { model } = unit;
for (let j = 0; j < count; ++j) { if (!count) continue;
for (let j = 0; j < offset.length - 1; ++j) {
const start = offset[j]; const start = offset[j];
const end = offset[j + 1]; const end = offset[j + 1];
for (let bI = start; bI < end; bI++) {
console.log(`${atomLabel(model, j)} -- ${atomLabel(model, neighbor[bI])}`) if (end <= start) continue;
const aI = ElementGroup.getAt(group, j);
for (let _bI = start; _bI < end; _bI++) {
const bI = ElementGroup.getAt(group, neighbor[_bI])
console.log(`${atomLabel(model, aI)} -- ${atomLabel(model, bI)}`);
} }
} }
} }
......
...@@ -256,7 +256,7 @@ function query(ctx: QueryContext): boolean { ...@@ -256,7 +256,7 @@ function query(ctx: QueryContext): boolean {
if (distSq <= rSq) { if (distSq <= rSq) {
if (maxRadius > 0 && Math.sqrt(distSq) - radius![idx] > inputRadius) continue; if (maxRadius > 0 && Math.sqrt(distSq) - radius![idx] > inputRadius) continue;
if (isCheck) return true; if (isCheck) return true;
Result.add(result, idx, distSq); Result.add(result, bucketArray[i], distSq);
} }
} }
} }
......
...@@ -112,7 +112,7 @@ function _computeBonds(unit: Unit.Atomic, atoms: ElementGroup, params: BondCompu ...@@ -112,7 +112,7 @@ function _computeBonds(unit: Unit.Atomic, atoms: ElementGroup, params: BondCompu
const { x, y, z } = unit.model.conformation; const { x, y, z } = unit.model.conformation;
const atomCount = ElementGroup.size(atoms); const atomCount = ElementGroup.size(atoms);
const { residueKey } = unit.model.hierarchy; const { residueIndex } = unit;
const { type_symbol, label_atom_id, label_alt_id } = unit.model.hierarchy.atoms; const { type_symbol, label_atom_id, label_alt_id } = unit.model.hierarchy.atoms;
const { label_comp_id } = unit.model.hierarchy.residues; const { label_comp_id } = unit.model.hierarchy.residues;
const query3d = Unit.getLookup3d(unit, atoms); const query3d = Unit.getLookup3d(unit, atoms);
...@@ -130,12 +130,10 @@ function _computeBonds(unit: Unit.Atomic, atoms: ElementGroup, params: BondCompu ...@@ -130,12 +130,10 @@ function _computeBonds(unit: Unit.Atomic, atoms: ElementGroup, params: BondCompu
for (let _aI = 0; _aI < atomCount; _aI++) { for (let _aI = 0; _aI < atomCount; _aI++) {
const aI = ElementGroup.getAt(atoms, _aI); const aI = ElementGroup.getAt(atoms, _aI);
const raI = residueKey.value(aI); const raI = residueIndex[aI];
// const rowA = dataIndex[aI]; // TODO
const rowA = aI;
if (!params.forceCompute && raI !== lastResidue) { if (!params.forceCompute && raI !== lastResidue) {
const resn = label_comp_id.value(rowA)!; const resn = label_comp_id.value(raI)!;
if (!!component && component.entries.has(resn)) { if (!!component && component.entries.has(resn)) {
componentMap = component.entries.get(resn)!.map; componentMap = component.entries.get(resn)!.map;
} else { } else {
...@@ -144,37 +142,35 @@ function _computeBonds(unit: Unit.Atomic, atoms: ElementGroup, params: BondCompu ...@@ -144,37 +142,35 @@ function _computeBonds(unit: Unit.Atomic, atoms: ElementGroup, params: BondCompu
} }
lastResidue = raI; lastResidue = raI;
const componentPairs = componentMap ? componentMap.get(label_atom_id.value(rowA)) : void 0; const componentPairs = componentMap ? componentMap.get(label_atom_id.value(aI)) : void 0;
const aeI = idx(type_symbol.value(rowA)!); const aeI = idx(type_symbol.value(aI)!);
const { indices, count, squaredDistances } = query3d.find(x[aI], y[aI], z[aI], MAX_RADIUS); const { indices, count, squaredDistances } = query3d.find(x[aI], y[aI], z[aI], MAX_RADIUS);
const isHa = isHydrogen(aeI); const isHa = isHydrogen(aeI);
const thresholdA = threshold(aeI); const thresholdA = threshold(aeI);
const altA = label_alt_id.value(rowA); const altA = label_alt_id.value(aI);
const metalA = MetalsSet.has(aeI); const metalA = MetalsSet.has(aeI);
const structConnEntries = params.forceCompute ? void 0 : structConn && structConn.getAtomEntries(aI); const structConnEntries = params.forceCompute ? void 0 : structConn && structConn.getAtomEntries(aI);
for (let ni = 0; ni < count; ni++) { for (let ni = 0; ni < count; ni++) {
const bI = indices[ni]; const _bI = indices[ni];
const bI = ElementGroup.getAt(atoms, _bI);
if (bI <= aI) continue; if (bI <= aI) continue;
// const rowB = dataIndex[bI]; // TODO const altB = label_alt_id.value(bI);
const rowB = bI;
const altB = label_alt_id.value(rowB);
if (altA && altB && altA !== altB) continue; if (altA && altB && altA !== altB) continue;
const beI = idx(type_symbol.value(rowB)!); const beI = idx(type_symbol.value(bI)!);
const isMetal = metalA || MetalsSet.has(beI); const isMetal = metalA || MetalsSet.has(beI);
const rbI = residueKey.value(bI); const rbI = residueIndex[bI];
// handle "component dictionary" bonds. // handle "component dictionary" bonds.
if (raI === rbI && componentPairs) { if (raI === rbI && componentPairs) {
const e = componentPairs.get(label_atom_id.value(rowB)!); const e = componentPairs.get(label_atom_id.value(bI)!);
if (e) { if (e) {
atomA[atomA.length] = aI; atomA[atomA.length] = _aI;
atomB[atomB.length] = bI; atomB[atomB.length] = _bI;
order[order.length] = e.order; order[order.length] = e.order;
let flag = e.flags; let flag = e.flags;
if (isMetal) { if (isMetal) {
...@@ -198,8 +194,8 @@ function _computeBonds(unit: Unit.Atomic, atoms: ElementGroup, params: BondCompu ...@@ -198,8 +194,8 @@ function _computeBonds(unit: Unit.Atomic, atoms: ElementGroup, params: BondCompu
for (const se of structConnEntries) { for (const se of structConnEntries) {
for (const p of se.partners) { for (const p of se.partners) {
if (p.atomIndex === bI) { if (p.atomIndex === bI) {
atomA[atomA.length] = aI; atomA[atomA.length] = _aI;
atomB[atomB.length] = bI; atomB[atomB.length] = _bI;
flags[flags.length] = se.flags; flags[flags.length] = se.flags;
order[order.length] = se.order; order[order.length] = se.order;
added = true; added = true;
...@@ -213,8 +209,8 @@ function _computeBonds(unit: Unit.Atomic, atoms: ElementGroup, params: BondCompu ...@@ -213,8 +209,8 @@ function _computeBonds(unit: Unit.Atomic, atoms: ElementGroup, params: BondCompu
if (isHa || isHb) { if (isHa || isHb) {
if (dist < params.maxHbondLength) { if (dist < params.maxHbondLength) {
atomA[atomA.length] = aI; atomA[atomA.length] = _aI;
atomB[atomB.length] = bI; atomB[atomB.length] = _bI;
order[order.length] = 1; order[order.length] = 1;
flags[flags.length] = BondType.Flag.Covalent | BondType.Flag.Computed; // TODO: check if correct flags[flags.length] = BondType.Flag.Covalent | BondType.Flag.Computed; // TODO: check if correct
} }
...@@ -228,8 +224,8 @@ function _computeBonds(unit: Unit.Atomic, atoms: ElementGroup, params: BondCompu ...@@ -228,8 +224,8 @@ function _computeBonds(unit: Unit.Atomic, atoms: ElementGroup, params: BondCompu
if (dist <= pairingThreshold) { if (dist <= pairingThreshold) {
atomA[atomA.length] = aI; atomA[atomA.length] = _aI;
atomB[atomB.length] = bI; atomB[atomB.length] = _bI;
order[order.length] = 1; order[order.length] = 1;
flags[flags.length] = (isMetal ? BondType.Flag.MetallicCoordination : BondType.Flag.Covalent) | BondType.Flag.Computed; flags[flags.length] = (isMetal ? BondType.Flag.MetallicCoordination : BondType.Flag.Covalent) | BondType.Flag.Computed;
} }
...@@ -237,6 +233,7 @@ function _computeBonds(unit: Unit.Atomic, atoms: ElementGroup, params: BondCompu ...@@ -237,6 +233,7 @@ function _computeBonds(unit: Unit.Atomic, atoms: ElementGroup, params: BondCompu
} }
const bonds = computePerAtomBonds(atomA, atomB, order, flags, atomCount); const bonds = computePerAtomBonds(atomA, atomB, order, flags, atomCount);
return { return {
offset: bonds.offsets, offset: bonds.offsets,
neighbor: bonds.neighbor, neighbor: bonds.neighbor,
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment