added dcd format parser

src/mol-io/common/binary.ts

@@ -1,5 +1,5 @@
 /**
- * Copyright (c) 2017-2018 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 Alexander Rose <alexander.rose@weirdbyte.de>
  * @author David Sehnal <david.sehnal@gmail.com>
@@ -16,4 +16,17 @@ export function flipByteOrder(data: Uint8Array, bytes: number) {
         }
     }
     return buffer;
+}
+
+const ChunkSize = 0x7000
+export function uint8ToString(array: Uint8Array) {
+    if (array.length > ChunkSize) {
+        const c = []
+        for (let i = 0; i < array.length; i += ChunkSize) {
+            c.push(String.fromCharCode.apply(null, array.subarray(i, i + ChunkSize)))
+        }
+        return c.join('')
+    } else {
+        return String.fromCharCode.apply(null, array)
+    }
 }

src/mol-io/reader/_spec/dcd.spec.ts

@@ -0,0 +1,75 @@
+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import { parseDcd } from '../dcd/parser';
+
+function createDcdData() {
+    const data = new Uint8Array(4 * 128)
+
+    const dv = new DataView(data.buffer)
+    // const intView = new Int32Array(data.buffer, 0, 23)
+
+    // set little endian
+    dv.setInt32(0, 84)
+
+    // set format string
+    dv.setUint8(4, 'D'.charCodeAt(0))
+    dv.setUint8(5, 'R'.charCodeAt(0))
+    dv.setUint8(6, 'O'.charCodeAt(0))
+    dv.setUint8(7, 'C'.charCodeAt(0))
+
+    dv.setInt32(8, 1) // NSET
+
+    // header end
+    dv.setInt32(22 * 4, 84)
+
+    // title
+    const titleEnd = 164
+    const titleStart = 23 * 4 + 1
+    dv.setInt32(23 * 4, titleEnd)
+    dv.setInt32(titleStart + titleEnd + 4 - 1, titleEnd)
+
+    // atoms
+    const atomStart = 23 * 4 + titleEnd + 8
+    dv.setInt32(atomStart, 4)
+    dv.setInt32(atomStart + 4, 1) // one atom
+    dv.setInt32(atomStart + 8, 4)
+
+    // coords
+    const coordsStart = atomStart + 12
+    dv.setInt32(coordsStart, 4)
+    dv.setFloat32(coordsStart + 4, 0.1)
+    dv.setInt32(coordsStart + 8, 4)
+    dv.setInt32(coordsStart + 12, 4)
+    dv.setFloat32(coordsStart + 16, 0.2)
+    dv.setInt32(coordsStart + 20, 4)
+    dv.setInt32(coordsStart + 24, 4)
+    dv.setFloat32(coordsStart + 28, 0.3)
+    dv.setInt32(coordsStart + 32, 4)
+
+    return data
+}
+
+describe('dcd reader', () => {
+    it('basic', async () => {
+        const data = createDcdData()
+        const parsed = await parseDcd(data).run();
+
+        if (parsed.isError) {
+            throw new Error(parsed.message)
+        }
+
+        const dcdFile = parsed.result;
+        const { header, frames } = dcdFile;
+
+        expect(header.NSET).toBe(1)
+        expect(header.NATOM).toBe(1)
+
+        expect(frames[0].x[0]).toBeCloseTo(0.1, 0.0001);
+        expect(frames[0].y[0]).toBeCloseTo(0.2, 0.0001);
+        expect(frames[0].z[0]).toBeCloseTo(0.3, 0.0001);
+    });
+});

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

@@ -0,0 +1,215 @@
+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+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 { uint8ToString } from '../../common/binary'
+
+export interface DcdHeader {
+    readonly NSET: number,
+    readonly ISTART: number,
+    readonly NSAVC: number,
+    readonly NAMNF: number,
+    readonly DELTA: number,
+    readonly TITLE: string,
+    readonly NATOM: number
+}
+
+export interface DcdFrame {
+    readonly cell: Cell
+    readonly elementCount: number
+
+    // positions
+    readonly x: ArrayLike<number>
+    readonly y: ArrayLike<number>
+    readonly z: ArrayLike<number>
+}
+
+export interface DcdFile {
+    readonly header: DcdHeader
+    readonly frames: DcdFrame[]
+}
+
+export function _parseDcd(data: Uint8Array): DcdFile {
+    // http://www.ks.uiuc.edu/Research/vmd/plugins/molfile/dcdplugin.html
+
+    // The DCD format is structured as follows
+    //   (FORTRAN UNFORMATTED, with Fortran data type descriptions):
+    // HDR     NSET    ISTRT   NSAVC   5-ZEROS NATOM-NFREAT    DELTA   9-ZEROS
+    // `CORD'  #files  step 1  step    zeroes  (zero)          timestep  (zeroes)
+    //                         interval
+    // C*4     INT     INT     INT     5INT    INT             DOUBLE  9INT
+    // ==========================================================================
+    // NTITLE          TITLE
+    // INT (=2)        C*MAXTITL
+    //                 (=32)
+    // ==========================================================================
+    // NATOM
+    // #atoms
+    // INT
+    // ==========================================================================
+    // X(I), I=1,NATOM         (DOUBLE)
+    // Y(I), I=1,NATOM
+    // Z(I), I=1,NATOM
+    // ==========================================================================
+
+    const dv = new DataView(data.buffer)
+
+    const header: Mutable<DcdHeader> = Object.create(null)
+    const frames: DcdFrame[] = []
+
+    let nextPos = 0
+
+    // header block
+
+    const intView = new Int32Array(data.buffer, 0, 23)
+    const ef = intView[0] !== dv.getInt32(0) // endianess flag
+    // swap byte order when big endian (84 indicates little endian)
+    if (intView[0] !== 84) {
+        const n = data.byteLength
+        for (let i = 0; i < n; i += 4) {
+            dv.setFloat32(i, dv.getFloat32(i), true)
+        }
+    }
+    if (intView[0] !== 84) {
+        throw new Error('dcd bad format, header block start')
+    }
+
+    // format indicator, should read 'CORD'
+    const formatString = String.fromCharCode(
+        dv.getUint8(4), dv.getUint8(5),
+        dv.getUint8(6), dv.getUint8(7)
+    )
+    if (formatString !== 'CORD') {
+        throw new Error('dcd bad format, format string')
+    }
+    let isCharmm = false
+    let extraBlock = false
+    let fourDims = false
+    // version field in charmm, unused in X-PLOR
+    if (intView[22] !== 0) {
+        isCharmm = true
+        if (intView[12] !== 0) extraBlock = true
+        if (intView[13] === 1) fourDims = true
+    }
+    header.NSET = intView[2]
+    header.ISTART = intView[3]
+    header.NSAVC = intView[4]
+    header.NAMNF = intView[10]
+
+    if (isCharmm) {
+        header.DELTA = dv.getFloat32(44, ef)
+    } else {
+        header.DELTA = dv.getFloat64(44, ef)
+    }
+
+    if (intView[22] !== 84) {
+        throw new Error('dcd bad format, header block end')
+    }
+    nextPos = nextPos + 21 * 4 + 8
+
+    // title block
+
+    const titleEnd = dv.getInt32(nextPos, ef)
+    const titleStart = nextPos + 1
+    if ((titleEnd - 4) % 80 !== 0) {
+        throw new Error('dcd bad format, title block start')
+    }
+    header.TITLE = uint8ToString(data.subarray(titleStart, titleEnd))
+    if (dv.getInt32(titleStart + titleEnd + 4 - 1, ef) !== titleEnd) {
+        throw new Error('dcd bad format, title block end')
+    }
+
+    nextPos = nextPos + titleEnd + 8
+
+    // natom block
+
+    if (dv.getInt32(nextPos, ef) !== 4) {
+        throw new Error('dcd bad format, natom block start')
+    }
+    header.NATOM = dv.getInt32(nextPos + 4, ef)
+    if (dv.getInt32(nextPos + 8, ef) !== 4) {
+        throw new Error('dcd bad format, natom block end')
+    }
+    nextPos = nextPos + 4 + 8
+
+    // fixed atoms block
+
+    if (header.NAMNF > 0) {
+        // TODO read coordinates and indices of fixed atoms
+        throw new Error('dcd format with fixed atoms unsupported, aborting')
+    }
+
+    // frames
+
+    const natom = header.NATOM
+    const natom4 = natom * 4
+
+    for (let i = 0, n = header.NSET; i < n; ++i) {
+        const frame: Mutable<DcdFrame> = Object.create({
+            elementCount: natom
+        })
+
+        if (extraBlock) {
+            nextPos += 4 // block start
+            // cell: A, alpha, B, beta, gamma, C (doubles)
+            const size = Vec3.create(
+                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 + 3 * 8, ef),
+                dv.getFloat64(nextPos + 4 * 8, ef)
+            )
+            frame.cell = Cell.create(size, anglesInRadians)
+            nextPos += 48
+            nextPos += 4 // block end
+        }
+
+        // xyz coordinates
+        for (let j = 0; j < 3; ++j) {
+            if (dv.getInt32(nextPos, ef) !== natom4) {
+                throw new Error(`dcd bad format, coord block start: ${i}, ${j}`)
+            }
+            nextPos += 4 // block start
+            const c = new Float32Array(data.buffer, nextPos, natom)
+            if (j === 0) frame.x = c
+            else if (j === 1) frame.y = c
+            else frame.z = c
+
+            nextPos += natom4
+            if (dv.getInt32(nextPos, ef) !== natom4) {
+                throw new Error(`dcd bad format, coord block end: ${i}, ${j}`)
+            }
+            nextPos += 4 // block end
+        }
+
+        if (fourDims) {
+            const bytes = dv.getInt32(nextPos, ef)
+            nextPos += 4 + bytes + 4 // block start + skip + block end
+        }
+
+        frames.push(frame)
+    }
+
+    return { header, frames }
+}
+
+export function parseDcd(data: Uint8Array) {
+    return Task.create<Result<DcdFile>>('Parse DCD', async ctx => {
+        try {
+            const dcdFile = _parseDcd(data)
+            return Result.success(dcdFile)
+        } catch (e) {
+            return Result.error(e)
+        }
+    })
+}

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

@@ -0,0 +1,40 @@
+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import { Task } from '../../mol-task';
+import { DcdFile } from '../../mol-io/reader/dcd/parser';
+import { Coordinates, Frame, Time } from '../../mol-model/structure/coordinates';
+
+const charmmTimeUnitFactor = 20.45482949774598
+
+export function coordinatesFromDcd(dcdFile: DcdFile): Task<Coordinates> {
+    console.log('coordinatesFromDcd', dcdFile)
+    return Task.create('Parse DCD', async ctx => {
+        await ctx.update('Converting to coordinates');
+
+        const { header } = dcdFile
+
+        let deltaTime = Time(1, 'step')
+        if (header.DELTA) {
+            deltaTime = Time(header.DELTA * charmmTimeUnitFactor, 'ps')
+        }
+
+        let offsetTime = Time(0, deltaTime.unit)
+        if (header.ISTART >= 1) {
+            offsetTime = Time((header.ISTART - 1) * deltaTime.value, 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)
+            })
+        }
+
+        return Coordinates.create(frames, deltaTime, offsetTime)
+    })
+}