Parse and generate NCS operators

data/mmcif-field-names.csv

@@ -144,6 +144,12 @@ struct_keywords.entry_id
 struct_keywords.pdbx_keywords
 struct_keywords.text
 
+struct_ncs_oper.id
+struct_ncs_oper.code
+struct_ncs_oper.matrix
+struct_ncs_oper.vector
+struct_ncs_oper.details
+
 struct_sheet_range.sheet_id
 struct_sheet_range.id
 struct_sheet_range.beg_label_comp_id

src/apps/structure-info/model.ts

@@ -18,6 +18,7 @@ import { openCif, downloadCif } from './helpers';
 import { BitFlags } from 'mol-util';
 import { SecondaryStructureType } from 'mol-model/structure/model/types';
 import { UnitRings } from 'mol-model/structure/structure/unit/rings';
+import { Vec3 } from 'mol-math/linear-algebra';
 
 
 async function downloadFromPdb(pdb: string) {
@@ -182,6 +183,16 @@ export function printUnits(structure: Structure) {
     }
 }
 
+export function printSymmetryInfo(model: Model) {
+    console.log('\nSymmetry Info\n=============');
+    const { symmetry } = model;
+    const { size, anglesInRadians } = symmetry.spacegroup.cell;
+    console.log(`Spacegroup: ${symmetry.spacegroup.name} size: ${Vec3.toString(size)} angles: ${Vec3.toString(anglesInRadians)}`);
+    console.log(`Assembly names: ${symmetry.assemblies.map(a => a.id).join(', ')}`);
+    // NCS example: 1auy
+    console.log(`NCS operators: ${symmetry.ncsOperators && symmetry.ncsOperators.map(a => a.name).join(', ')}`);
+}
+
 export function printIHMModels(model: Model) {
     if (!model.coarseHierarchy.isDefined) return false;
     console.log('\nIHM Models\n=============');
@@ -194,6 +205,7 @@ async function run(frame: CifFrame) {
     //printSequence(models[0]);
     //printIHMModels(models[0]);
     printUnits(structure);
+    printSymmetryInfo(models[0]);
     //printRings(structure);
     //printLinks(structure, true, true);
     //printModRes(models[0]);

src/mol-io/reader/cif/schema/mmcif.ts

@@ -409,6 +409,13 @@ export const mmCIF_Schema = {
         experiment_type: Aliased<'Fraction of bulk' | 'Single molecule'>(str),
         details: str,
     },
+    struct_ncs_oper: {
+        id: str,
+        code: str,
+        matrix: Matrix(3, 3),
+        vector: Vector(3),
+        details: str
+    },
     ihm_modeling_post_process: {
         id: int,
         protocol_id: int,

src/mol-math/geometry/symmetry-operator.ts

@@ -4,7 +4,7 @@
  * @author David Sehnal <david.sehnal@gmail.com>
  */
 
-import { Vec3, Mat4 } from '../linear-algebra/3d'
+import { Vec3, Mat4, Mat3 } from '../linear-algebra/3d'
 
 interface SymmetryOperator {
     readonly name: string,
@@ -30,6 +30,17 @@ namespace SymmetryOperator {
         return { name, matrix, inverse: Mat4.invert(Mat4.zero(), matrix), isIdentity: false, hkl: _hkl };
     }
 
+    export function ofRotationAndOffset(name: string, rot: Mat3, offset: Vec3) {
+        const t = Mat4.identity();
+        for (let i = 0; i < 3; i++) {
+            for (let j = 0; j < 3; j++) {
+                Mat4.setValue(t, i, j, Mat3.getValue(rot, i, j));
+            }
+        }
+        Mat4.setTranslation(t, offset);
+        return create(name, t);
+    }
+
     // Apply the 1st and then 2nd operator. ( = second.matrix * first.matrix)
     export function compose(first: SymmetryOperator, second: SymmetryOperator) {
         const matrix = Mat4.mul(Mat4.zero(), second.matrix, first.matrix);

src/mol-math/linear-algebra/3d/mat3.ts

@@ -108,6 +108,10 @@ namespace Mat3 {
         a[3 * j + i] = value;
     }
 
+    export function getValue(a: Mat3, i: number, j: number) {
+        return a[3 * j + i];
+    }
+
     /**
      * Copy the values from one Mat3 to another
      */

src/mol-math/linear-algebra/3d/vec3.ts

@@ -367,6 +367,10 @@ namespace Vec3 {
     export function isZero(v: Vec3) {
         return v[0] === 0 && v[1] === 0 && v[2] === 0
     }
+
+    export function toString(a: Vec3) {
+        return `[${a[0]} ${a[1]} ${a[2]}]`;
+    }
 }
 
 export default Vec3

src/mol-model/structure/model/formats/mmcif.ts

@@ -6,8 +6,8 @@
 
 import { Column, Table } from 'mol-data/db';
 import { Interval, Segmentation } from 'mol-data/int';
-import { Spacegroup, SpacegroupCell } from 'mol-math/geometry';
-import { Vec3 } from 'mol-math/linear-algebra';
+import { Spacegroup, SpacegroupCell, SymmetryOperator } from 'mol-math/geometry';
+import { Vec3, Tensor, Mat4 } from 'mol-math/linear-algebra';
 import { Task } from 'mol-task';
 import UUID from 'mol-util/uuid';
 import Format from '../format';
@@ -22,7 +22,7 @@ import { getIHMCoarse } from './mmcif/ihm';
 import { getSecondaryStructureMmCif } from './mmcif/secondary-structure';
 import { getSequence } from './mmcif/sequence';
 import { sortAtomSite } from './mmcif/sort';
-import { mmCIF_Database } from 'mol-io/reader/cif/schema/mmcif';
+import { mmCIF_Database, mmCIF_Schema } from 'mol-io/reader/cif/schema/mmcif';
 
 import mmCIF_Format = Format.mmCIF
 type AtomSite = mmCIF_Database['atom_site']
@@ -90,7 +90,7 @@ function getSymmetry(format: mmCIF_Format): ModelSymmetry {
     const assemblies = createAssemblies(format);
     const spacegroup = getSpacegroup(format);
     const isNonStandardCrytalFrame = checkNonStandardCrystalFrame(format, spacegroup);
-    return { assemblies, spacegroup, isNonStandardCrytalFrame };
+    return { assemblies, spacegroup, isNonStandardCrytalFrame, ncsOperators: getNcsOperators(format) };
 }
 
 function checkNonStandardCrystalFrame(format: mmCIF_Format, spacegroup: Spacegroup) {
@@ -111,6 +111,22 @@ function getSpacegroup(format: mmCIF_Format): Spacegroup {
     return Spacegroup.create(spaceCell);
 }
 
+function getNcsOperators(format: mmCIF_Format) {
+    const { struct_ncs_oper } = format.data;
+    if (struct_ncs_oper._rowCount === 0) return void 0;
+    const { id, matrix, vector } = struct_ncs_oper;
+
+    const matrixSpace = mmCIF_Schema.struct_ncs_oper.matrix.space, vectorSpace = mmCIF_Schema.struct_ncs_oper.vector.space;
+
+    const opers: SymmetryOperator[] = [];
+    for (let i = 0; i < struct_ncs_oper._rowCount; i++) {
+        const m = Tensor.toMat3(matrixSpace, matrix.value(i));
+        const v = Tensor.toVec3(vectorSpace, vector.value(i));
+        opers[i] = SymmetryOperator.ofRotationAndOffset(`ncs_${id.value(i)}`, m, v);
+    }
+    return opers;
+}
+
 function isHierarchyDataEqual(a: AtomicData, b: AtomicData) {
     // need to cast because of how TS handles type resolution for interfaces https://github.com/Microsoft/TypeScript/issues/15300
     return Table.areEqual(a.chains as Table<ChainsSchema>, b.chains as Table<ChainsSchema>)

src/mol-model/structure/model/properties/symmetry.ts

@@ -45,6 +45,7 @@ interface ModelSymmetry {
     readonly assemblies: ReadonlyArray<Assembly>,
     readonly spacegroup: Spacegroup,
     readonly isNonStandardCrytalFrame: boolean,
+    readonly ncsOperators?: ReadonlyArray<SymmetryOperator>,
 
     // optionally cached operators from [-3, -3, -3] to [3, 3, 3]
     _operators_333?: SymmetryOperator[]

src/mol-model/structure/structure/symmetry.ts

@@ -51,6 +51,11 @@ namespace StructureSymmetry {
         return Task.create('Build Symmetry', ctx => findSymmetryRange(ctx, structure, ijkMin, ijkMax));
     }
 
+    /** Builds NCS structure, returns the original if NCS operators are not present. */
+    export function buildNcs(structure: Structure) {
+        return Task.create('Build NCS', ctx => _buildNCS(ctx, structure));
+    }
+
     function hashUnit(u: Unit) {
         return hash2(u.invariantId, SortedArray.hashCode(u.elements));
     }
@@ -98,27 +103,37 @@ function getOperators(symmetry: ModelSymmetry, ijkMin: Vec3, ijkMax: Vec3) {
     return operators;
 }
 
-async function findSymmetryRange(ctx: RuntimeContext, structure: Structure, ijkMin: Vec3, ijkMax: Vec3) {
-    const models = Structure.getModels(structure);
-    if (models.length !== 1) throw new Error('Can only build symmetries from structures based on 1 model.');
-
-    const { spacegroup } = models[0].symmetry;
-    if (SpacegroupCell.isZero(spacegroup.cell)) return structure;
-
-    const operators = getOperators(models[0].symmetry, ijkMin, ijkMax);
+function assembleOperators(structure: Structure, operators: ReadonlyArray<SymmetryOperator>) {
     const assembler = Structure.Builder();
-
     const { units } = structure;
     for (const oper of operators) {
         for (const unit of units) {
             assembler.addWithOperator(unit, oper);
         }
-        if (ctx.shouldUpdate) await ctx.update('Building symmetry...');
     }
-
     return assembler.getStructure();
 }
 
+async function _buildNCS(ctx: RuntimeContext, structure: Structure) {
+    const models = Structure.getModels(structure);
+    if (models.length !== 1) throw new Error('Can only build NCS from structures based on 1 model.');
+
+    const operators = models[0].symmetry.ncsOperators;
+    if (!operators || !operators.length) return structure;
+    return assembleOperators(structure, operators);
+}
+
+async function findSymmetryRange(ctx: RuntimeContext, structure: Structure, ijkMin: Vec3, ijkMax: Vec3) {
+    const models = Structure.getModels(structure);
+    if (models.length !== 1) throw new Error('Can only build symmetries from structures based on 1 model.');
+
+    const { spacegroup } = models[0].symmetry;
+    if (SpacegroupCell.isZero(spacegroup.cell)) return structure;
+
+    const operators = getOperators(models[0].symmetry, ijkMin, ijkMax);
+    return assembleOperators(structure, operators);
+}
+
 async function findMatesRadius(ctx: RuntimeContext, structure: Structure, radius: number) {
     const models = Structure.getModels(structure);
     if (models.length !== 1) throw new Error('Can only build symmetries from structures based on 1 model.');