src/parser/dcd-parser.js
/**
* @file Dcd Parser
* @author Alexander Rose <alexander.rose@weirdbyte.de>
* @private
*/
import { Debug, Log, ParserRegistry } from '../globals.js'
import { uint8ToString, ensureBuffer } from '../utils.js'
import TrajectoryParser from './trajectory-parser.js'
const charmmTimeUnitFactor = 20.45482949774598
class DcdParser extends TrajectoryParser {
get type () { return 'dcd' }
get isBinary () { return true }
_parse () {
// 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
// ==========================================================================
if (Debug) Log.time('DcdParser._parse ' + this.name)
const bin = ensureBuffer(this.streamer.data)
const dv = new DataView(bin)
const f = this.frames
const coordinates = f.coordinates
const boxes = f.boxes
const header = {}
let nextPos = 0
// header block
const intView = new Int32Array(bin, 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 = bin.byteLength
for (let i = 0; i < n; i += 4) {
dv.setFloat32(i, dv.getFloat32(i), true)
}
}
if (intView[ 0 ] !== 84) {
Log.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') {
Log.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) {
Log.error('dcd bad format, header block end')
}
nextPos = nextPos + 21 * 4 + 8
// title block
const titleLength = dv.getInt32(nextPos, ef)
const titlePos = nextPos + 1
if ((titleLength - 4) % 80 !== 0) {
Log.error('dcd bad format, title block start')
}
header.TITLE = uint8ToString(
new Uint8Array(bin, titlePos, titleLength)
)
if (dv.getInt32(titlePos + titleLength + 4 - 1, ef) !== titleLength) {
Log.error('dcd bad format, title block end')
}
nextPos = nextPos + titleLength + 8
// natom block
if (dv.getInt32(nextPos, ef) !== 4) {
Log.error('dcd bad format, natom block start')
}
header.NATOM = dv.getInt32(nextPos + 4, ef)
if (dv.getInt32(nextPos + 8, ef) !== 4) {
Log.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
Log.error('dcd format with fixed atoms unsupported, aborting')
return
}
// frames
const natom = header.NATOM
const natom4 = natom * 4
for (let i = 0, n = header.NSET; i < n; ++i) {
if (extraBlock) {
nextPos += 4 // block start
// unitcell: A, alpha, B, beta, gamma, C (doubles)
const box = new Float32Array(9)
box[ 0 ] = dv.getFloat64(nextPos, ef)
box[ 4 ] = dv.getFloat64(nextPos + 2 * 8, ef)
box[ 8 ] = dv.getFloat64(nextPos + 5 * 8, ef)
boxes.push(box)
nextPos += 48
nextPos += 4 // block end
}
// xyz coordinates
const coord = new Float32Array(natom * 3)
for (let j = 0; j < 3; ++j) {
if (dv.getInt32(nextPos, ef) !== natom4) {
Log.error('dcd bad format, coord block start', i, j)
}
nextPos += 4 // block start
const c = new Float32Array(bin, nextPos, natom)
for (let k = 0; k < natom; ++k) {
coord[ 3 * k + j ] = c[ k ]
}
nextPos += natom4
if (dv.getInt32(nextPos, ef) !== natom4) {
Log.error('dcd bad format, coord block end', i, j)
}
nextPos += 4 // block end
}
coordinates.push(coord)
if (fourDims) {
const bytes = dv.getInt32(nextPos, ef)
nextPos += 4 + bytes + 4 // block start + skip + block end
}
}
if (header.DELTA) {
f.deltaTime = header.DELTA * charmmTimeUnitFactor
}
if (header.ISTART >= 1) {
f.timeOffset = (header.ISTART - 1) * f.deltaTime
}
// console.log(header)
// console.log(header.TITLE)
// console.log('isCharmm', isCharmm, 'extraBlock', extraBlock, 'fourDims, fourDims)
if (Debug) Log.timeEnd('DcdParser._parse ' + this.name)
}
}
ParserRegistry.add('dcd', DcdParser)
export default DcdParser