calculate symmetry operators for transform around deposited coordinates

src/mol-math/geometry/spacegroup/construction.ts

@@ -1,7 +1,8 @@
 /**
- * Copyright (c) 2018 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2018-2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author David Sehnal <david.sehnal@gmail.com>
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
  */
 
 import { Vec3, Mat4 } from '../../linear-algebra'
@@ -73,7 +74,7 @@ namespace SpacegroupCell {
 
 
 namespace Spacegroup {
-    // P1 with [1, 1, 1] cell.
+    /** P1 with [1, 1, 1] cell */
     export const ZeroP1 = create(SpacegroupCell.Zero);
 
     export function create(cell: SpacegroupCell): Spacegroup {
@@ -81,18 +82,56 @@ namespace Spacegroup {
         return { name: SpacegroupNames[cell.index], cell, operators };
     }
 
-    const _tempVec = Vec3.zero(), _tempMat = Mat4.zero();
-    export function updateOperatorMatrix(spacegroup: Spacegroup, index: number, i: number, j: number, k: number, target: Mat4) {
-        _tempVec[0] = i;
-        _tempVec[1] = j;
-        _tempVec[2] = k;
-
-        Mat4.fromTranslation(_tempMat, _tempVec);
-        return Mat4.mul(target, Mat4.mul(target, Mat4.mul(target, spacegroup.cell.fromFractional, _tempMat), spacegroup.operators[index]), spacegroup.cell.toFractional);
+    const _ijkVec = Vec3();
+    const _tempMat = Mat4();
+    export function setOperatorMatrix(spacegroup: Spacegroup, index: number, i: number, j: number, k: number, target: Mat4) {
+        Vec3.set(_ijkVec, i, j, k);
+
+        Mat4.fromTranslation(_tempMat, _ijkVec);
+        return Mat4.mul(
+            target,
+            Mat4.mul(
+                target,
+                Mat4.mul(target, spacegroup.cell.fromFractional, _tempMat),
+                spacegroup.operators[index]
+            ),
+            spacegroup.cell.toFractional
+        );
     }
 
     export function getSymmetryOperator(spacegroup: Spacegroup, index: number, i: number, j: number, k: number): SymmetryOperator {
-        const operator = updateOperatorMatrix(spacegroup, index, i, j, k, Mat4.zero());
+        const operator = setOperatorMatrix(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), index);
+    }
+
+    const _translationRef = Vec3()
+    const _translationRefSymop = Vec3()
+    const _translationSymop = Vec3()
+    export function setOperatorMatrixRef(spacegroup: Spacegroup, index: number, i: number, j: number, k: number, ref: Vec3, target: Mat4) {
+        Vec3.set(_ijkVec, i, j, k);
+        Vec3.floor(_translationRef, ref)
+
+        Mat4.copy(target, spacegroup.operators[index])
+
+        Vec3.floor(_translationRefSymop, Vec3.transformMat4(_translationRefSymop, ref, target))
+
+        Mat4.getTranslation(_translationSymop, target)
+        Vec3.sub(_translationSymop, _translationSymop, _translationRefSymop)
+        Vec3.add(_translationSymop, _translationSymop, _translationRef)
+        Vec3.add(_translationSymop, _translationSymop, _ijkVec)
+
+        Mat4.setTranslation(target, _translationSymop)
+        Mat4.mul(target, spacegroup.cell.fromFractional, target)
+        Mat4.mul(target, target, spacegroup.cell.toFractional)
+        return target
+    }
+
+    /**
+     * Get Symmetry operator for transformation around the given
+     * reference point `ref` in fractional coordinates
+     */
+    export function getSymmetryOperatorRef(spacegroup: Spacegroup, index: number, i: number, j: number, k: number, ref: Vec3) {
+        const operator = setOperatorMatrixRef(spacegroup, index, i, j, k, ref, Mat4.zero());
         return SymmetryOperator.create(`${index + 1}_${5 + i}${5 + j}${5 + k}`, operator, { id: '', operList: [] }, '', Vec3.create(i, j, k), index);
     }
 

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

@@ -1,5 +1,5 @@
 /**
- * Copyright (c) 2017 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2017-2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author David Sehnal <david.sehnal@gmail.com>
  */
@@ -9,6 +9,7 @@ import { arrayFind } from '../../../../mol-data/util'
 import { StructureQuery } from '../../query'
 import { Model } from '../../model'
 import { Spacegroup } from '../../../../mol-math/geometry';
+import { Vec3 } from '../../../../mol-math/linear-algebra';
 
 /** Determine an atom set and a list of operators that should be applied to that set  */
 export interface OperatorGroup {
@@ -47,8 +48,14 @@ interface ModelSymmetry {
     readonly isNonStandardCrytalFrame: boolean,
     readonly ncsOperators?: ReadonlyArray<SymmetryOperator>,
 
-    // optionally cached operators from [-3, -3, -3] to [3, 3, 3]
-    _operators_333?: SymmetryOperator[]
+    /**
+     * optionally cached operators from [-3, -3, -3] to [3, 3, 3]
+     * around reference point `ref` in fractional coordinates
+     */
+    _operators_333?: {
+        ref: Vec3,
+        operators: SymmetryOperator[]
+    }
 }
 
 namespace ModelSymmetry {

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

@@ -10,7 +10,7 @@ import { EquivalenceClasses } from '../../../mol-data/util';
 import { Spacegroup, SpacegroupCell, SymmetryOperator } from '../../../mol-math/geometry';
 import { Vec3, Mat4 } from '../../../mol-math/linear-algebra';
 import { RuntimeContext, Task } from '../../../mol-task';
-import { ModelSymmetry } from '../model';
+import { ModelSymmetry, Model } from '../model';
 import { QueryContext, StructureSelection } from '../query';
 import Structure from './structure';
 import Unit from './unit';
@@ -95,7 +95,7 @@ namespace StructureSymmetry {
     }
 }
 
-function getOperators(symmetry: ModelSymmetry, ijkMin: Vec3, ijkMax: Vec3) {
+function getOperators(symmetry: ModelSymmetry, ijkMin: Vec3, ijkMax: Vec3, modelCenter: Vec3) {
     const { spacegroup, ncsOperators } = symmetry;
     const ncsCount = (ncsOperators && ncsOperators.length) || 0
     const operators: SymmetryOperator[] = [];
@@ -107,13 +107,17 @@ function getOperators(symmetry: ModelSymmetry, ijkMin: Vec3, ijkMax: Vec3) {
         operators[0] = Spacegroup.getSymmetryOperator(spacegroup, 0, 0, 0, 0)
     }
 
+    const { toFractional } = spacegroup.cell
+    const ref = Vec3.transformMat4(Vec3(), modelCenter, toFractional)
+
     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++) {
                     // check if we have added identity as the 1st operator.
                     if (!ncsCount && op === 0 && i === 0 && j === 0 && k === 0) continue;
-                    const symOp = Spacegroup.getSymmetryOperator(spacegroup, op, i, j, k);
+
+                    const symOp = Spacegroup.getSymmetryOperatorRef(spacegroup, op, i, j, k, ref)
                     if (ncsCount) {
                         for (let u = 0; u < ncsCount; ++u) {
                             const ncsOp = ncsOperators![u]
@@ -131,10 +135,15 @@ function getOperators(symmetry: ModelSymmetry, ijkMin: Vec3, ijkMax: Vec3) {
     return operators;
 }
 
-function getOperatorsCached333(symmetry: ModelSymmetry) {
-    if (typeof symmetry._operators_333 !== 'undefined') return symmetry._operators_333;
-    symmetry._operators_333 = getOperators(symmetry, Vec3.create(-3, -3, -3), Vec3.create(3, 3, 3));
-    return symmetry._operators_333;
+function getOperatorsCached333(symmetry: ModelSymmetry, ref: Vec3) {
+    if (symmetry._operators_333 && Vec3.equals(ref, symmetry._operators_333.ref)) {
+        return symmetry._operators_333.operators;
+    }
+    symmetry._operators_333 = {
+        ref: Vec3.clone(ref),
+        operators: getOperators(symmetry, Vec3.create(-3, -3, -3), Vec3.create(3, 3, 3), ref)
+    };
+    return symmetry._operators_333.operators;
 }
 
 function assembleOperators(structure: Structure, operators: ReadonlyArray<SymmetryOperator>) {
@@ -164,7 +173,8 @@ async function findSymmetryRange(ctx: RuntimeContext, structure: Structure, ijkM
     const { spacegroup } = models[0].symmetry;
     if (SpacegroupCell.isZero(spacegroup.cell)) return structure;
 
-    const operators = getOperators(models[0].symmetry, ijkMin, ijkMax);
+    const modelCenter = Model.getCenter(models[0])
+    const operators = getOperators(models[0].symmetry, ijkMin, ijkMax, modelCenter);
     return assembleOperators(structure, operators);
 }
 
@@ -177,7 +187,8 @@ async function findMatesRadius(ctx: RuntimeContext, structure: Structure, radius
     if (SpacegroupCell.isZero(spacegroup.cell)) return structure;
 
     if (ctx.shouldUpdate) await ctx.update('Initialing...');
-    const operators = getOperatorsCached333(symmetry);
+    const modelCenter = Model.getCenter(models[0])
+    const operators = getOperatorsCached333(symmetry, modelCenter);
     const lookup = structure.lookup3d;
 
     const assembler = Structure.Builder();