import { base64ToArray } from "./utilities";
import { Vector3, Matrix4 } from "./WebGL/math";
import { VASP } from "./parsers/VASP";
import { CUBE } from "./parsers/CUBE";
import { inflate } from "pako";
interface VolumeDataOptions {
negate?: boolean;
normalize?: boolean;
};
/**
* $3Dmol.VolumeData stores volumetric data. This includes file parsing
* functionality.
*
* @class
* @param {string} str - volumetric data
* @param {string} format - format of supplied data (cube, dx, vasp); append .gz if compressed
* @param {Object} options - normalize (zero mean, unit variance), negate
*/
export class VolumeData {
unit = {
x: 1,
y: 1,
z: 1
}; // scale of each voxel
origin = {
x: 0,
y: 0,
z: 0
}; // origin (bottom "left", not center)
size = {
x: 0,
y: 0,
z: 0
}; // number of voxels in each direction
data = new Float32Array([]); // actual floating point data, arranged
// x->y->z
matrix: any = null; //if set must transform data
inversematrix: Matrix4|null = null;
dimensionorder: any;
isbinary = new Set<string>(['ccp4','CCP4']);
constructor(str: any, format: string, options?: VolumeDataOptions) {
format = format.toLowerCase();
if (/\.gz$/.test(format)) {
//unzip gzipped files
format = format.replace(/\.gz$/, '');
try {
if ((this as any)[format] && this.isbinary.has(format)) {
if (typeof (str) == "string") {
//assume base64 encoded
str = base64ToArray(str);
}
str = inflate(str);
}
else {
str = new TextDecoder("utf-8").decode(inflate(str));
}
} catch (err) {
console.error(err);
}
}
if ((this as any)[format]) {
if (this.isbinary.has(format) && typeof (str) == "string") {
str = base64ToArray(str);
}
(this as any)[format](str);
}
if (options) {
if (options.negate) {
for (let i = 0, n = this.data.length; i < n; i++) {
this.data[i] = -this.data[i];
}
}
if (options.normalize) {
var total = 0.0;
for (let i = 0, n = this.data.length; i < n; i++) {
total += this.data[i];
}
var mean = total / this.data.length;
total = 0;
for (let i = 0, n = this.data.length; i < n; i++) {
var diff = this.data[i] - mean;
total += diff * diff; //variance is ave of squared difference with mean
}
var variance = total / this.data.length;
//now normalize
for (let i = 0, n = this.data.length; i < n; i++) {
this.data[i] = (this.data[i] - mean) / variance;
}
}
}
}
/**
* @function $3Dmol.VolumeData.getIndex
* @param {number} x,y,z - the coordinates
* @returns - index into flat array closest to provided coordinate; -1 if invalid
*/
getIndex(x: number, y: number, z: number) {
if (this.matrix) {
//all transformation is done through matrix multiply
if (this.inversematrix == null) {
this.inversematrix = new Matrix4().getInverse(this.matrix);
}
var pt = new Vector3(x, y, z);
pt = pt.applyMatrix4(this.inversematrix as Matrix4);
x = pt.x;
y = pt.y;
z = pt.z;
} else { //use simple origin/unit transform
x -= this.origin.x;
y -= this.origin.y;
z -= this.origin.z;
x /= this.unit.x;
y /= this.unit.y;
z /= this.unit.z;
}
x = Math.round(x);
y = Math.round(y);
z = Math.round(z);
if (x < 0 || x >= this.size.x) return -1;
if (y < 0 || y >= this.size.y) return -1;
if (z < 0 || z >= this.size.z) return -1;
return x * this.size.y * this.size.z + y * this.size.z + z;
};
/**
* @function $3Dmol.VolumeData.getVal
* @param {number} x,y,z - the coordinates
* @returns - value closest to provided coordinate; zero if coordinate invalid
*/
getVal(x: number, y: number, z: number) {
let i = this.getIndex(x, y, z);
if (i < 0) return 0;
return this.data[i];
};
getCoordinates = function (index: number) {
var x = index / (this.size.y * this.size.z);
var y = index % (this.size.y * this.size.z);
var z = index % this.size.z;
x *= this.unit.x;
y *= this.unit.y;
z *= this.unit.z;
x += this.origin.x;
y += this.origin.y;
z += this.origin.z;
return { x: x, y: y, z: z };
};
/*
* parse vasp data
* Essentially this parser converts the CHGCAR data into
* cube data. It has been adapted from 'chg2cube.pl' found in
* http://theory.cm.utexas.edu/vtsttools/
*/
vasp = function (str: string) {
var lines = str.replace(/^\s+/, "").split(/[\n\r]/);
var atomicData = VASP(str)[0];
var natoms = atomicData.length;
if (natoms == 0) {
console.warn("No good formating of CHG or CHGCAR file, not atomic information provided in the file.");
this.data = [];
return;
}
// Assume atomic units
// var unittype = "bohr/hartree";
var l_units = 1.889725992;
var e_units = 0.036749309;
// copied from $3Dmol.Parsers.vasp
var convFactor = parseFloat(lines[1]);
// This is how Vasp reads in the basis We need the l_units in order to
// compute the volume of the cell. Afterwards to obtain the axis for the
// voxels we have to remove this unit and divide by the number of voxels in
// each dimension
var v: string[];
v = lines[2].replace(/^\s+/, "").split(/\s+/);
var xVec = new Vector3(parseFloat(v[0]), parseFloat(v[1]), parseFloat(v[2])).multiplyScalar(convFactor * l_units);
v = lines[3].replace(/^\s+/, "").split(/\s+/);
var yVec = new Vector3(parseFloat(v[0]), parseFloat(v[1]), parseFloat(v[2])).multiplyScalar(convFactor * l_units);
v = lines[4].replace(/^\s+/, "").split(/\s+/);
var zVec = new Vector3(parseFloat(v[0]), parseFloat(v[1]), parseFloat(v[2])).multiplyScalar(convFactor * l_units);
// correct volume for non-orthognal box (expansion by minors)
var vol = xVec.x * (yVec.y * zVec.z - zVec.y * yVec.z) - yVec.x * (xVec.y * zVec.z - zVec.y * xVec.z) + zVec.x * (xVec.y * yVec.z - yVec.y * xVec.z);
vol = Math.abs(vol) / (Math.pow(l_units, 3));
var vol_scale = 1.0 / (vol); //This Only for CHGCAR files
// We splice the structure information
// 2 (header) + 3 (vectors) + 2 (atoms) + 1 (vaspMode) + natoms (coords) + 1 (blank line)
lines.splice(0, 2 + 3 + 2 + 1 + natoms + 1);
var lineArr = lines[0].replace(/^\s+/, "").replace(/\s+/g, " ").split(" ");
var nX = Math.abs(parseFloat(lineArr[0]));
var nY = Math.abs(parseFloat(lineArr[1]));
var nZ = Math.abs(parseFloat(lineArr[2]));
var origin = this.origin = new Vector3(0, 0, 0);
this.size = { x: nX, y: nY, z: nZ };
this.unit = new Vector3(xVec.x, yVec.y, zVec.z);
// resize the vectors accordingly
xVec = xVec.multiplyScalar(1 / (l_units * nX));
yVec = yVec.multiplyScalar(1 / (l_units * nY));
zVec = zVec.multiplyScalar(1 / (l_units * nZ));
if (xVec.y != 0 || xVec.z != 0 || yVec.x != 0 || yVec.z != 0 || zVec.x != 0
|| zVec.y != 0) {
//need a transformation matrix
this.matrix = new Matrix4(xVec.x, yVec.x, zVec.x, 0, xVec.y, yVec.y, zVec.y, 0, xVec.z, yVec.z, zVec.z, 0, 0, 0, 0, 1);
//include translation in matrix
this.matrix = this.matrix.multiplyMatrices(this.matrix,
new Matrix4().makeTranslation(origin.x, origin.y, origin.z));
//all translation and scaling done by matrix, so reset origin and unit
this.origin = new Vector3(0, 0, 0);
this.unit = new Vector3(1, 1, 1);
}
lines.splice(0, 1); // Remove the dimension line
var raw = lines.join(" ");
raw = raw.replace(/^\s+/, '');
var rawArray = raw.split(/[\s\r]+/);
rawArray.splice(nX * nY * nZ + 1);
var preConvertedData = Float32Array.from(rawArray, parseFloat); //We still have to format it to get the density
for (var i = 0; i < preConvertedData.length; i++) {
preConvertedData[i] = preConvertedData[i] * vol_scale * e_units;
}
this.data = preConvertedData;
};
// parse dx data - does not support all features of the file format
dx = function (str: string) {
var lines = str.split(/[\n\r]+/);
var m: string[];
var recounts = /gridpositions\s+counts\s+(\d+)\s+(\d+)\s+(\d+)/;
var reorig = /^origin\s+(\S+)\s+(\S+)\s+(\S+)/;
var redelta = /^delta\s+(\S+)\s+(\S+)\s+(\S+)/;
var follows = /data follows/;
var i = 0;
for (i = 0; i < lines.length; i++) {
var line = lines[i];
if ((m = recounts.exec(line))) {
var nX = parseInt(m[1]);
var nY = parseInt(m[2]);
var nZ = parseInt(m[3]);
this.size = { x: nX, y: nY, z: nZ };
}
else if ((m = redelta.exec(line))) {
var xunit = parseFloat(m[1]);
if (parseFloat(m[2]) != 0 || parseFloat(m[3]) != 0) {
console.warn("Non-orthogonal delta matrix not currently supported in dx format");
}
i += 1;
line = lines[i];
m = redelta.exec(line);
if (m == null) {
console.error("Parse error in dx delta matrix");
return;
}
var yunit = parseFloat(m[2]);
if (parseFloat(m[1]) != 0 || parseFloat(m[3]) != 0) {
console.warn("Non-orthogonal delta matrix not currently supported in dx format");
}
i += 1;
line = lines[i];
m = redelta.exec(line);
if (m == null) {
console.error("Parse error in dx delta matrix");
return;
}
var zunit = parseFloat(m[3]);
if (parseFloat(m[1]) != 0 || parseFloat(m[2]) != 0) {
console.warn("Non-orthogonal delta matrix not currently supported in dx format");
}
this.unit = new Vector3(xunit, yunit, zunit);
}
else if ((m = reorig.exec(line))) {
var xorig = parseFloat(m[1]);
var yorig = parseFloat(m[2]);
var zorig = parseFloat(m[3]);
this.origin = new Vector3(xorig, yorig, zorig);
} else if ((m = follows.exec(line))) {
break;
}
}
i += 1;
if (!this.size || !this.origin || !this.unit || !this.size) {
console.error("Error parsing dx format");
return;
}
var raw = lines.splice(i).join(" ");
var rawArray = raw.split(/[\s\r]+/);
this.data = Float32Array.from(rawArray, parseFloat);
};
// parse cube data
cube(str: string) {
var lines = str.split(/\r?\n/);
if (lines.length < 6)
return;
var cryst = CUBE(str, {}).modelData[0].cryst;
var lineArr = lines[2].replace(/^\s+/, "").replace(/\s+/g, " ").split(" ");
var atomsnum = parseFloat(lineArr[0]); //includes sign, which indicates presence of oribital line in header
var natoms = Math.abs(atomsnum);
this.origin = cryst.origin;
this.size = cryst.size;
this.unit = cryst.unit;
this.matrix = cryst.matrix4;
var headerlines = 6;
if (atomsnum < 0) headerlines++; //see: http://www.ks.uiuc.edu/Research/vmd/plugins/molfile/cubeplugin.html
var raw = lines.splice(natoms + headerlines).join(" ");
raw = raw.replace(/^\s+/, '');
var rawArray = raw.split(/[\s\r]+/);
this.data = Float32Array.from(rawArray, parseFloat);
};
//parse cp4 files
ccp4(bin: Int8Array) {
// http://www.ccp4.ac.uk/html/maplib.html#description
//code from ngl: https://github.com/arose/ngl/blob/master/js/ngl/parser.js
var header:any = {};
bin = new Int8Array(bin);
var intView = new Int32Array(bin.buffer, 0, 56);
var floatView = new Float32Array(bin.buffer, 0, 56);
var dv = new DataView(bin.buffer);
// 53 MAP Character string 'MAP ' to identify file type
header.MAP = String.fromCharCode(
dv.getUint8(52 * 4), dv.getUint8(52 * 4 + 1),
dv.getUint8(52 * 4 + 2), dv.getUint8(52 * 4 + 3)
);
// 54 MACHST Machine stamp indicating machine type which wrote file
// 17 and 17 for big-endian or 68 and 65 for little-endian
header.MACHST = [dv.getUint8(53 * 4), dv.getUint8(53 * 4 + 1)];
// swap byte order when big endian
if (header.MACHST[0] === 17 && header.MACHST[1] === 17) {
var n = bin.byteLength;
for (var i = 0; i < n; i += 4) {
dv.setFloat32(i, dv.getFloat32(i), true);
}
}
header.NX = intView[0]; // NC - columns (fastest changing)
header.NY = intView[1]; // NR - rows
header.NZ = intView[2]; // NS - sections (slowest changing)
// mode
// 0 image : signed 8-bit bytes range -128 to 127
// 1 image : 16-bit halfwords
// 2 image : 32-bit reals
// 3 transform : complex 16-bit integers
// 4 transform : complex 32-bit reals
// 6 image : unsigned 16-bit range 0 to 65535
// 16 image: unsigned char * 3 (for rgb data, non-standard)
//
// Note: Mode 2 is the normal mode used in the CCP4 programs.
// Other modes than 2 and 0 may NOT WORK
header.MODE = intView[3];
// start
header.NXSTART = intView[4]; // NCSTART - first column
header.NYSTART = intView[5]; // NRSTART - first row
header.NZSTART = intView[6]; // NSSTART - first section
// intervals
header.MX = intView[7]; // intervals along x
header.MY = intView[8]; // intervals along y
header.MZ = intView[9]; // intervals along z
// cell length (Angstroms in CCP4)
header.xlen = floatView[10];
header.ylen = floatView[11];
header.zlen = floatView[12];
// cell angle (Degrees)
header.alpha = floatView[13];
header.beta = floatView[14];
header.gamma = floatView[15];
// axis correspondence (1,2,3 for X,Y,Z)
header.MAPC = intView[16]; // column
header.MAPR = intView[17]; // row
header.MAPS = intView[18]; // section
// density statistics
header.DMIN = floatView[19];
header.DMAX = floatView[20];
header.DMEAN = floatView[21];
// space group number 0 or 1 (default=0)
header.ISPG = intView[22];
// number of bytes used for symmetry data (0 or 80)
header.NSYMBT = intView[23];
// Flag for skew transformation, =0 none, =1 if foll
header.LSKFLG = intView[24];
// 26-34 SKWMAT Skew matrix S (in order S11, S12, S13, S21 etc) if
// LSKFLG .ne. 0.
// 35-37 SKWTRN Skew translation t if LSKFLG != 0.
// Skew transformation is from standard orthogonal
// coordinate frame (as used for atoms) to orthogonal
// map frame, as Xo(map) = S * (Xo(atoms) - t)
// 38 future use (some of these are used by the MSUBSX routines
// . " in MAPBRICK, MAPCONT and FRODO)
// . " (all set to zero by default)
// . "
// 52 "
// 50-52 origin in X,Y,Z used for transforms
header.originX = floatView[49];
header.originY = floatView[50];
header.originZ = floatView[51];
// 53 MAP Character string 'MAP ' to identify file type
// => see top of this parser
// 54 MACHST Machine stamp indicating machine type which wrote file
// => see top of this parser
// Rms deviation of map from mean density
header.ARMS = floatView[54];
// 56 NLABL Number of labels being used
// 57-256 LABEL(20,10) 10 80 character text labels (ie. A4 format)
//create transformation matrix, code mostly copied from ngl
var h = header;
var basisX: Array<any> = [
h.xlen,
0,
0
];
var basisY: Array<any> = [
h.ylen * Math.cos(Math.PI / 180.0 * h.gamma),
h.ylen * Math.sin(Math.PI / 180.0 * h.gamma),
0
];
var basisZ: Array<any> = [
h.zlen * Math.cos(Math.PI / 180.0 * h.beta),
h.zlen * (
Math.cos(Math.PI / 180.0 * h.alpha)
- Math.cos(Math.PI / 180.0 * h.gamma)
* Math.cos(Math.PI / 180.0 * h.beta)
) / Math.sin(Math.PI / 180.0 * h.gamma),
0
];
basisZ[2] = Math.sqrt(
h.zlen * h.zlen * Math.sin(Math.PI / 180.0 * h.beta) *
Math.sin(Math.PI / 180.0 * h.beta) - basisZ[1] * basisZ[1]
);
var basis: Array<any> = [0, basisX, basisY, basisZ];
var nxyz: Array<any> = [0, h.MX, h.MY, h.MZ];
var mapcrs: Array<any> = [0, h.MAPC, h.MAPR, h.MAPS];
this.matrix = new Matrix4();
this.matrix.set(
basis[mapcrs[1]][0] / nxyz[mapcrs[1]],
basis[mapcrs[2]][0] / nxyz[mapcrs[2]],
basis[mapcrs[3]][0] / nxyz[mapcrs[3]],
0,
basis[mapcrs[1]][1] / nxyz[mapcrs[1]],
basis[mapcrs[2]][1] / nxyz[mapcrs[2]],
basis[mapcrs[3]][1] / nxyz[mapcrs[3]],
0,
basis[mapcrs[1]][2] / nxyz[mapcrs[1]],
basis[mapcrs[2]][2] / nxyz[mapcrs[2]],
basis[mapcrs[3]][2] / nxyz[mapcrs[3]],
0,
0, 0, 0, 1
);
//include translation in matrix, NXSTART etc are an offset in grid space
this.matrix = this.matrix.multiplyMatrices(
this.matrix,
new Matrix4().makeTranslation(
h.NXSTART + h.originX,
h.NYSTART + h.originY,
h.NZSTART + h.originZ)
);
//all translation and scaling done by matrix, so reset origin and unit
this.origin = new Vector3(0, 0, 0);
this.unit = new Vector3(1, 1, 1);
this.size = { x: header.NX, y: header.NY, z: header.NZ };
this.dimensionorder = [header.MAPC, header.MAPR, header.MAPS];
var data = new Float32Array(bin.buffer, 1024 + header.NSYMBT);
//data must by (slowest changing) x,y,z (fastest changing)
var NX = header.NX, NY = header.NY, NZ = header.NZ;
this.data = new Float32Array(NX * NY * NZ);
for (let i = 0; i < NX; i++) {
for (let j = 0; j < NY; j++) {
for (let k = 0; k < NZ; k++) {
//should I be concerned that I'm not using mapc?
this.data[((i * NY) + j) * NZ + k] = data[((k * NY) + j) * NX + i];
}
}
}
};
};