From 5221e929dbacbbfa87ca6949475664528ac3ee87 Mon Sep 17 00:00:00 2001
From: Alexander Rose <alexander.rose@weirdbyte.de>
Date: Sat, 16 Feb 2019 14:12:20 -0800
Subject: [PATCH] better support for ncs operators

---
 .../geometry/spacegroup/construction.ts       |  2 +-
 src/mol-math/geometry/symmetry-operator.ts    | 22 +++++++++++--------
 .../structure/mmcif/parser.ts                 |  3 ++-
 .../structure/structure/structure.ts          |  2 +-
 src/mol-model/structure/structure/symmetry.ts | 17 ++++++++++++--
 5 files changed, 32 insertions(+), 14 deletions(-)

diff --git a/src/mol-math/geometry/spacegroup/construction.ts b/src/mol-math/geometry/spacegroup/construction.ts
index 076a849f6..a347c60e7 100644
--- a/src/mol-math/geometry/spacegroup/construction.ts
+++ b/src/mol-math/geometry/spacegroup/construction.ts
@@ -90,7 +90,7 @@ namespace Spacegroup {
 
     export function getSymmetryOperator(spacegroup: Spacegroup, index: number, i: number, j: number, k: number): SymmetryOperator {
         const operator = updateOperatorMatrix(spacegroup, index, i, j, k, Mat4.zero());
-        return SymmetryOperator.create(`${index + 1}_${5 + i}${5 + j}${5 + k}`, operator, { id: '', operList: [] }, Vec3.create(i, j, k));
+        return SymmetryOperator.create(`${index + 1}_${5 + i}${5 + j}${5 + k}`, operator, { id: '', operList: [] }, '', Vec3.create(i, j, k));
     }
 
     function getOperatorMatrix(ids: number[]) {
diff --git a/src/mol-math/geometry/symmetry-operator.ts b/src/mol-math/geometry/symmetry-operator.ts
index 0c8297c8c..5064c7359 100644
--- a/src/mol-math/geometry/symmetry-operator.ts
+++ b/src/mol-math/geometry/symmetry-operator.ts
@@ -10,12 +10,15 @@ interface SymmetryOperator {
     readonly name: string,
 
     readonly assembly: {
-        /** pointer to `pdbx_struct_assembly.id` */
+        /** pointer to `pdbx_struct_assembly.id` or empty string */
         readonly id: string
-        /** pointers to `pdbx_struct_oper_list_id` */
+        /** pointers to `pdbx_struct_oper_list.id` or empty list */
         readonly operList: string[]
     }
 
+    /** pointer to `struct_ncs_oper.id` or empty string */
+    readonly ncsId: string,
+
     readonly hkl: Vec3,
 
     readonly matrix: Mat4,
@@ -31,11 +34,12 @@ namespace SymmetryOperator {
 
     const RotationEpsilon = 0.0001;
 
-    export function create(name: string, matrix: Mat4, assembly: SymmetryOperator['assembly'], hkl?: Vec3): SymmetryOperator {
+    export function create(name: string, matrix: Mat4, assembly: SymmetryOperator['assembly'], ncsId?: string, hkl?: Vec3): SymmetryOperator {
         const _hkl = hkl ? Vec3.clone(hkl) : Vec3.zero();
-        if (Mat4.isIdentity(matrix)) return { name, assembly, matrix, inverse: Mat4.identity(), isIdentity: true, hkl: _hkl };
+        ncsId = ncsId || ''
+        if (Mat4.isIdentity(matrix)) return { name, assembly, matrix, inverse: Mat4.identity(), isIdentity: true, hkl: _hkl, ncsId };
         if (!Mat4.isRotationAndTranslation(matrix, RotationEpsilon)) throw new Error(`Symmetry operator (${name}) must be a composition of rotation and translation.`);
-        return { name, assembly, matrix, inverse: Mat4.invert(Mat4.zero(), matrix), isIdentity: false, hkl: _hkl };
+        return { name, assembly, matrix, inverse: Mat4.invert(Mat4.zero(), matrix), isIdentity: false, hkl: _hkl, ncsId };
     }
 
     export function checkIfRotationAndTranslation(rot: Mat3, offset: Vec3) {
@@ -49,7 +53,7 @@ namespace SymmetryOperator {
         return Mat4.isRotationAndTranslation(matrix, RotationEpsilon);
     }
 
-    export function ofRotationAndOffset(name: string, rot: Mat3, offset: Vec3) {
+    export function ofRotationAndOffset(name: string, rot: Mat3, offset: Vec3, ncsId?: string) {
         const t = Mat4.identity();
         for (let i = 0; i < 3; i++) {
             for (let j = 0; j < 3; j++) {
@@ -57,16 +61,16 @@ namespace SymmetryOperator {
             }
         }
         Mat4.setTranslation(t, offset);
-        return create(name, t, { id: '', operList: [] });
+        return create(name, t, { id: '', operList: [] }, ncsId);
     }
 
     /**
      * Apply the 1st and then 2nd operator. ( = second.matrix * first.matrix).
-     * Keep `name`, `assembly` and `hkl` properties from second.
+     * Keep `name`, `assembly`, `ncsId` and `hkl` properties from second.
      */
     export function compose(first: SymmetryOperator, second: SymmetryOperator) {
         const matrix = Mat4.mul(Mat4.zero(), second.matrix, first.matrix);
-        return create(second.name, matrix, second.assembly, second.hkl);
+        return create(second.name, matrix, second.assembly, second.ncsId, second.hkl);
     }
 
     export interface CoordinateMapper<T extends number> { (index: T, slot: Vec3): Vec3 }
diff --git a/src/mol-model-formats/structure/mmcif/parser.ts b/src/mol-model-formats/structure/mmcif/parser.ts
index 4a16920ba..934de9467 100644
--- a/src/mol-model-formats/structure/mmcif/parser.ts
+++ b/src/mol-model-formats/structure/mmcif/parser.ts
@@ -75,7 +75,8 @@ function getNcsOperators(format: mmCIF_Format) {
         const m = Tensor.toMat3(matrixSpace, matrix.value(i));
         const v = Tensor.toVec3(vectorSpace, vector.value(i));
         if (!SymmetryOperator.checkIfRotationAndTranslation(m, v)) continue;
-        opers[opers.length] = SymmetryOperator.ofRotationAndOffset(`ncs_${id.value(i)}`, m, v);
+        const ncsId = id.value(i)
+        opers[opers.length] = SymmetryOperator.ofRotationAndOffset(`ncs_${ncsId}`, m, v, ncsId);
     }
     return opers;
 }
diff --git a/src/mol-model/structure/structure/structure.ts b/src/mol-model/structure/structure/structure.ts
index b1c520601..b86668743 100644
--- a/src/mol-model/structure/structure/structure.ts
+++ b/src/mol-model/structure/structure/structure.ts
@@ -381,7 +381,7 @@ namespace Structure {
         const units: Unit[] = [];
         for (const u of s.units) {
             const old = u.conformation.operator;
-            const op = SymmetryOperator.create(old.name, transform, { id: '', operList: [] }, old.hkl);
+            const op = SymmetryOperator.create(old.name, transform, { id: '', operList: [] }, old.ncsId, old.hkl);
             units.push(u.applyOperator(u.id, op));
         }
 
diff --git a/src/mol-model/structure/structure/symmetry.ts b/src/mol-model/structure/structure/symmetry.ts
index 162d3151c..1b9170414 100644
--- a/src/mol-model/structure/structure/symmetry.ts
+++ b/src/mol-model/structure/structure/symmetry.ts
@@ -8,7 +8,7 @@
 import { SortedArray } from 'mol-data/int';
 import { EquivalenceClasses } from 'mol-data/util';
 import { Spacegroup, SpacegroupCell, SymmetryOperator } from 'mol-math/geometry';
-import { Vec3 } from 'mol-math/linear-algebra';
+import { Vec3, Mat4 } from 'mol-math/linear-algebra';
 import { RuntimeContext, Task } from 'mol-task';
 import { ModelSymmetry } from '../model';
 import { QueryContext, StructureSelection } from '../query';
@@ -98,13 +98,26 @@ function getOperators(symmetry: ModelSymmetry, ijkMin: Vec3, ijkMax: Vec3) {
         operators[0] = Spacegroup.getSymmetryOperator(spacegroup, 0, 0, 0, 0)
     }
 
+    const { ncsOperators } = symmetry
+    const ncsCount = (ncsOperators && ncsOperators.length) || 0
+
     for (let op = 0; op < spacegroup.operators.length; op++) {
         for (let i = ijkMin[0]; i <= ijkMax[0]; i++) {
             for (let j = ijkMin[1]; j <= ijkMax[1]; j++) {
                 for (let k = ijkMin[2]; k <= ijkMax[2]; k++) {
                     // we have added identity as the 1st operator.
                     if (op === 0 && i === 0 && j === 0 && k === 0) continue;
-                    operators[operators.length] = Spacegroup.getSymmetryOperator(spacegroup, op, i, j, k);
+                    const symOp = Spacegroup.getSymmetryOperator(spacegroup, op, i, j, k);
+                    if (ncsCount) {
+                        for (let u = 0; u < ncsCount; ++u) {
+                            const ncsOp = ncsOperators![u]
+                            const matrix = Mat4.mul(Mat4.zero(), symOp.matrix, ncsOp.matrix)
+                            const operator = SymmetryOperator.create(`${symOp.name} ${ncsOp.name}`, matrix, symOp.assembly, ncsOp.ncsId, symOp.hkl);
+                            operators[operators.length] = operator;
+                        }
+                    } else {
+                        operators[operators.length] = symOp;
+                    }
                 }
             }
         }
-- 
GitLab