trajectory cell handling improvements

src/mol-io/reader/dcd/parser.ts

@@ -6,9 +6,7 @@
 
 import { ReaderResult as Result } from '../result';
 import { Task } from '../../../mol-task';
-import { Cell } from '../../../mol-math/geometry/spacegroup/cell';
-import { Vec3 } from '../../../mol-math/linear-algebra';
-import { Mutable } from '../../../mol-util/type-helpers';
+import { Mutable, FiniteArray } from '../../../mol-util/type-helpers';
 import { uint8ToString } from '../../common/binary';
 
 export interface DcdHeader {
@@ -30,7 +28,7 @@ export interface DcdFrame {
     readonly z: ArrayLike<number>
 
     // optional cell
-    readonly cell?: Cell
+    readonly cell?: FiniteArray<number, 6>
 }
 
 export interface DcdFile {
@@ -158,19 +156,14 @@ export function _parseDcd(data: Uint8Array): DcdFile {
 
         if (extraBlock) {
             nextPos += 4; // block start
-            // TODO this is not standardized and we need to add heuristics to handle more variants
-            // cell: A, alpha, B, beta, gamma, C (doubles)
-            const size = Vec3.create(
+            frame.cell = [
                 dv.getFloat64(nextPos, ef),
-                dv.getFloat64(nextPos + 2 * 8, ef),
-                dv.getFloat64(nextPos + 5 * 8, ef)
-            );
-            const anglesInRadians = Vec3.create(
                 dv.getFloat64(nextPos + 1, ef),
+                dv.getFloat64(nextPos + 2 * 8, ef),
                 dv.getFloat64(nextPos + 3 * 8, ef),
-                dv.getFloat64(nextPos + 4 * 8, ef)
-            );
-            frame.cell = Cell.create(size, anglesInRadians);
+                dv.getFloat64(nextPos + 4 * 8, ef),
+                dv.getFloat64(nextPos + 5 * 8, ef)
+            ] as const;
             nextPos += 48;
             nextPos += 4; // block end
         }

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

@@ -1,5 +1,5 @@
 /**
- * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2019-2020 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
  */
@@ -18,6 +18,27 @@ function Cell() {
 }
 
 namespace Cell {
-    export function create(size: Vec3, anglesInRadians: Vec3): Cell { return { size, anglesInRadians }; }
-    export function empty(): Cell { return { size: Vec3(), anglesInRadians: Vec3() }; }
+    export function create(size: Vec3, anglesInRadians: Vec3): Cell {
+        return { size, anglesInRadians };
+    }
+
+    export function empty(): Cell {
+        return create(Vec3(), Vec3());
+    }
+
+    export function fromBasis(x: Vec3, y: Vec3, z: Vec3) {
+        const a = Vec3.magnitude(x);
+        const b = Vec3.magnitude(y);
+        const c = Vec3.magnitude(z);
+
+        const alpha = Math.acos(Vec3.dot(y, z) / (b * c));
+        const beta = Math.acos(Vec3.dot(x, z) / (a * c));
+        const gamma = Math.acos(Vec3.dot(x, y) / (a * b));
+
+        if (a <= 0 || b <= 0 || c <= 0 || alpha >= Math.PI || beta >= Math.PI || gamma >= Math.PI) {
+            return empty();
+        } else {
+            return create(Vec3.create(a, b, c), Vec3.create(alpha, beta, gamma));
+        }
+    }
 }

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

@@ -1,5 +1,5 @@
 /**
- * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2019-2020 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
  */
@@ -7,6 +7,10 @@
 import { Task } from '../../mol-task';
 import { DcdFile } from '../../mol-io/reader/dcd/parser';
 import { Coordinates, Frame, Time } from '../../mol-model/structure/coordinates';
+import { Vec3 } from '../../mol-math/linear-algebra';
+import { degToRad, halfPI } from '../../mol-math/misc';
+import { Cell } from '../../mol-math/geometry/spacegroup/cell';
+import { Mutable } from '../../mol-util/type-helpers';
 
 const charmmTimeUnitFactor = 20.45482949774598;
 
@@ -16,22 +20,55 @@ export function coordinatesFromDcd(dcdFile: DcdFile): Task<Coordinates> {
 
         const { header } = dcdFile;
 
-        let deltaTime = Time(1, 'step');
-        if (header.DELTA) {
-            deltaTime = Time(header.DELTA * charmmTimeUnitFactor, 'ps');
-        }
+        const deltaTime = header.DELTA
+            ? Time(header.DELTA * charmmTimeUnitFactor, 'ps')
+            : Time(1, 'step');
 
-        let offsetTime = Time(0, deltaTime.unit);
-        if (header.ISTART >= 1) {
-            offsetTime = Time((header.ISTART - 1) * deltaTime.value, deltaTime.unit);
-        }
+        const offsetTime = header.ISTART >= 1
+            ? Time((header.ISTART - 1) * deltaTime.value, deltaTime.unit)
+            : Time(0, deltaTime.unit);
 
         const frames: Frame[] = [];
         for (let i = 0, il = dcdFile.frames.length; i < il; ++i) {
-            frames.push({
-                ...dcdFile.frames[i],
-                time: Time(offsetTime.value + deltaTime.value * i, deltaTime.unit)
-            });
+            const dcdFrame = dcdFile.frames[i];
+            const frame: Mutable<Frame> = {
+                elementCount: dcdFrame.elementCount,
+                time: Time(offsetTime.value + deltaTime.value * i, deltaTime.unit),
+
+                x: dcdFrame.x,
+                y: dcdFrame.y,
+                z: dcdFrame.z,
+            };
+
+            if (dcdFrame.cell) {
+                // this is not standardized, using heuristics to handle variants
+                const c = dcdFrame.cell;
+                if (c[1] >= -1 && c[1] <= 1 && c[3] >= -1 && c[3] <= 1 && c[4] >= -1 && c[4] <= 1) {
+                    frame.cell = Cell.create(
+                        Vec3.create(c[0], c[2], c[5]),
+                        Vec3.create(
+                            degToRad(90 - Math.asin(c[1]) * 90 / halfPI),
+                            degToRad(90 - Math.asin(c[3]) * 90 / halfPI),
+                            degToRad(90 - Math.asin(c[4]) * 90 / halfPI)
+                        )
+                    );
+                } else if (
+                    c[0] < 0 || c[1] < 0 || c[2] < 0 || c[3] < 0 || c[4] < 0 || c[5] < 0 ||
+                    c[3] > 180 || c[4] > 180 || c[5] > 180
+                ) {
+                    frame.cell = Cell.fromBasis(
+                        Vec3.create(c[0], c[1], c[3]),
+                        Vec3.create(c[1], c[2], c[4]),
+                        Vec3.create(c[3], c[4], c[5])
+                    );
+                } else {
+                    frame.cell = Cell.create(
+                        Vec3.create(c[0], c[2], c[5]),
+                        Vec3.create(degToRad(c[1]), degToRad(c[3]), degToRad(c[4]))
+                    );
+                }
+            }
+            frames.push(frame);
         }
 
         return Coordinates.create(frames, deltaTime, offsetTime);

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

@@ -39,6 +39,12 @@ namespace ModelSymmetry {
         const isNonStandardCrystalFrame = checkNonStandardCrystalFrame(data.atom_sites, spacegroup);
         return { assemblies, spacegroup, isNonStandardCrystalFrame, ncsOperators: getNcsOperators(data.struct_ncs_oper) };
     }
+
+    export function fromCell(size: Vec3, anglesInRadians: Vec3): Symmetry {
+        const spaceCell = SpacegroupCell.create('P 1', size, anglesInRadians);
+        const spacegroup = Spacegroup.create(spaceCell);
+        return { assemblies: [], spacegroup, isNonStandardCrystalFrame: false };
+    }
 }
 
 function checkNonStandardCrystalFrame(atom_sites: Table<mmCIF_Schema['atom_sites']>, spacegroup: Spacegroup) {

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

@@ -2,12 +2,14 @@
  * Copyright (c) 2020 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 { Task } from '../../mol-task';
 import { XtcFile } from '../../mol-io/reader/xtc/parser';
 import { Coordinates, Frame, Time } from '../../mol-model/structure/coordinates';
 import { Cell } from '../../mol-math/geometry/spacegroup/cell';
+import { Vec3 } from '../../mol-math/linear-algebra';
 
 export function coordinatesFromXtc(file: XtcFile): Task<Coordinates> {
     return Task.create('Parse XTC', async ctx => {
@@ -18,10 +20,13 @@ export function coordinatesFromXtc(file: XtcFile): Task<Coordinates> {
 
         const frames: Frame[] = [];
         for (let i = 0, il = file.frames.length; i < il; ++i) {
+            const box = file.boxes[i];
+            const x = Vec3.fromArray(Vec3(), box, 0);
+            const y = Vec3.fromArray(Vec3(), box, 3);
+            const z = Vec3.fromArray(Vec3(), box, 6);
             frames.push({
                 elementCount: file.frames[i].count,
-                // TODO:
-                cell: Cell.empty(),
+                cell: Cell.fromBasis(x, y, z),
                 x: file.frames[i].x,
                 y: file.frames[i].y,
                 z: file.frames[i].z,

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

@@ -111,6 +111,12 @@ export namespace Model {
                 _staticPropertyData: Object.create(null),
                 _dynamicPropertyData: Object.create(null)
             };
+
+            if (f.cell) {
+                const symmetry = ModelSymmetry.fromCell(f.cell.size, f.cell.anglesInRadians);
+                ModelSymmetry.Provider.set(m, symmetry);
+            }
+
             trajectory.push(m);
         }
         return { trajectory, srcIndexArray };

src/mol-repr/shape/model/unitcell.ts

@@ -34,11 +34,17 @@ const CellRef = {
     model: 'Model'
 };
 
+const CellAttachement = {
+    corner: 'Corner',
+    center: 'Center'
+};
+
 const CellParams = {
     ...Mesh.Params,
     cellColor: PD.Color(ColorNames.orange),
     cellScale: PD.Numeric(2, { min: 0.1, max: 5, step: 0.1 }),
-    ref: PD.Select('model', PD.objectToOptions(CellRef), { isEssential: true })
+    ref: PD.Select('model', PD.objectToOptions(CellRef), { isEssential: true }),
+    attachement: PD.Select('corner', PD.objectToOptions(CellAttachement), { isEssential: true }),
 };
 type MeshParams = typeof CellParams
 
@@ -57,7 +63,13 @@ function getUnitcellMesh(data: UnitcellData, props: UnitcellProps, mesh?: Mesh)
 
     const { fromFractional } = data.symmetry.spacegroup.cell;
 
-    Vec3.floor(tmpRef, data.ref);
+    Vec3.copy(tmpRef, data.ref);
+    if (props.attachement === 'center') {
+        Vec3.trunc(tmpRef, tmpRef);
+        Vec3.subScalar(tmpRef, tmpRef, 0.5);
+    } else {
+        Vec3.floor(tmpRef, tmpRef);
+    }
     Mat4.fromTranslation(tmpTranslate, tmpRef);
     const cellCage = transformCage(cloneCage(unitCage), tmpTranslate);