00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include <stdio.h>
00019 #include <stdlib.h>
00020 #include <math.h>
00021 #include <string.h>
00022
00023 #include "molfile_plugin.h"
00024
00025 #ifndef NAN //not a number
00026 const float NAN = sqrtf(-1.f);
00027 #endif
00028
00029 #define TOLERANCE 1e-4
00030
00031 typedef struct {
00032 FILE *fd;
00033 int nsets;
00034 molfile_volumetric_t *vol;
00035 } edm_t;
00036
00037 static void eatline(FILE * fd) {
00038 char readbuf[1025];
00039 fgets(readbuf, 1024, fd);
00040 }
00041
00042 static void *open_edm_read(const char *filepath, const char *filetype,
00043 int *natoms) {
00044 FILE *fd;
00045 edm_t *edm;
00046 int ntitle, na, nb, nc, xsize, ysize, zsize;
00047 int amin, amax, bmin, bmax, cmin, cmax;
00048 float a, b, c, alpha, beta, gamma;
00049 float xdelta, ydelta, zdelta;
00050 float alpha1, beta1, gamma1;
00051 float xaxis[3], yaxis[3], zaxis[3], z1, z2, z3;
00052 int i, convcnt;
00053
00054 fd = fopen(filepath, "rb");
00055 if (!fd)
00056 return NULL;
00057
00058 edm = new edm_t;
00059 edm->fd = fd;
00060 edm->vol = NULL;
00061 *natoms = MOLFILE_NUMATOMS_NONE;
00062
00063 edm->vol = new molfile_volumetric_t[1];
00064
00065 edm->nsets = 1;
00066
00067
00068 eatline(edm->fd);
00069
00070 convcnt = fscanf(edm->fd, "%d", &ntitle);
00071 if (convcnt != 1) {
00072 printf("edmplugin) failed to read in title line count\n");
00073 fclose(edm->fd);
00074 delete [] edm->vol;
00075 delete edm;
00076 return NULL;
00077 }
00078
00079 eatline(edm->fd);
00080
00081
00082 for (i=0; i<ntitle; i++) {
00083 eatline(edm->fd);
00084 }
00085
00086
00087 convcnt = fscanf(edm->fd, "%d %d %d %d %d %d %d %d %d",
00088 &na, &amin, &amax, &nb, &bmin, &bmax, &nc, &cmin, &cmax);
00089 if (convcnt != 9) {
00090 printf("edmplugin) failed to read in box dimensions\n");
00091 fclose(edm->fd);
00092 delete [] edm->vol;
00093 delete edm;
00094 return NULL;
00095 }
00096
00097 eatline(edm->fd);
00098
00099
00100 xsize = amax - amin + 1;
00101 ysize = bmax - bmin + 1;
00102 zsize = cmax - cmin + 1;
00103 edm->vol[0].xsize = xsize;
00104 edm->vol[0].ysize = ysize;
00105 edm->vol[0].zsize = zsize;
00106 edm->vol[0].has_color = 0;
00107
00108
00109 convcnt = fscanf(edm->fd, "%f %f %f %f %f %f",
00110 &a, &b, &c, &alpha, &beta, &gamma);
00111 if (convcnt != 6) {
00112 printf("edmplugin) failed to read in box lengths and angles\n");
00113 fclose(edm->fd);
00114 delete [] edm->vol;
00115 delete edm;
00116 return NULL;
00117 }
00118 eatline(edm->fd);
00119
00120
00121 xdelta = a / (float) na;
00122 ydelta = b / (float) nb;
00123 zdelta = c / (float) nc;
00124
00125 strcpy(edm->vol[0].dataname, "X-PLOR Electron Density Map");
00126
00127
00128 alpha1 = 3.14159265358979323846 * alpha / 180.0;
00129 beta1 = 3.14159265358979323846 * beta / 180.0;
00130 gamma1 = 3.14159265358979323846 * gamma / 180.0;
00131
00132
00133 xaxis[0] = xdelta;
00134 xaxis[1] = 0;
00135 xaxis[2] = 0;
00136
00137 yaxis[0] = cos(gamma1) * ydelta;
00138 yaxis[1] = sin(gamma1) * ydelta;
00139 yaxis[2] = 0;
00140
00141 z1 = cos(beta1);
00142 z2 = (cos(alpha1) - cos(beta1)*cos(gamma1)) / sin(gamma1);
00143 z3 = sqrt(1.0 - z1*z1 - z2*z2);
00144 zaxis[0] = z1 * zdelta;
00145 zaxis[1] = z2 * zdelta;
00146 zaxis[2] = z3 * zdelta;
00147
00148 edm->vol[0].origin[0] = xaxis[0] * amin + yaxis[0] * bmin + zaxis[0] * cmin;
00149 edm->vol[0].origin[1] = yaxis[1] * bmin + zaxis[1] * cmin;
00150 edm->vol[0].origin[2] = zaxis[2] * cmin;
00151
00152 edm->vol[0].xaxis[0] = xaxis[0] * (xsize-1);
00153 edm->vol[0].xaxis[1] = 0;
00154 edm->vol[0].xaxis[2] = 0;
00155
00156 edm->vol[0].yaxis[0] = yaxis[0] * (ysize-1);
00157 edm->vol[0].yaxis[1] = yaxis[1] * (ysize-1);
00158 edm->vol[0].yaxis[2] = 0;
00159
00160 edm->vol[0].zaxis[0] = zaxis[0] * (zsize-1);
00161 edm->vol[0].zaxis[1] = zaxis[1] * (zsize-1);
00162 edm->vol[0].zaxis[2] = zaxis[2] * (zsize-1);
00163
00164
00165
00166 char planeorder[4];
00167 memset(planeorder, 0, sizeof(planeorder));
00168 convcnt = fscanf(edm->fd, "%3s", planeorder);
00169 if (convcnt != 1) {
00170 printf("edmplugin) failed to read in plane order\n");
00171 fclose(edm->fd);
00172 delete [] edm->vol;
00173 delete edm;
00174 return NULL;
00175 }
00176
00177 if (strcmp(planeorder, "ZYX")) {
00178 printf("edmplugin) unsupported plane ordering %s\n", planeorder);
00179 fclose(edm->fd);
00180 delete [] edm->vol;
00181 delete edm;
00182 return NULL;
00183 }
00184 eatline(edm->fd);
00185
00186 return edm;
00187 }
00188
00189 static int read_edm_metadata(void *v, int *nsets,
00190 molfile_volumetric_t **metadata) {
00191 edm_t *edm = (edm_t *)v;
00192 *nsets = edm->nsets;
00193 *metadata = edm->vol;
00194
00195 return MOLFILE_SUCCESS;
00196 }
00197
00198 static int read_edm_data(void *v, int set, float *datablock,
00199 float *colorblock) {
00200 edm_t *edm = (edm_t *)v;
00201 float * cell = datablock;
00202 int z, sentinel, convcnt;
00203 char readbuf[16];
00204
00205 int xsize = edm->vol[0].xsize;
00206 int ysize = edm->vol[0].ysize;
00207 int zsize = edm->vol[0].zsize;
00208
00209
00210 int nperslice = xsize * ysize;
00211 int nlines = (int)(nperslice / 6.0);
00212 if ((nlines * 6) < nperslice) {
00213 nlines++;
00214 }
00215
00216 for (z=0; z<zsize; z++) {
00217 int c;
00218 eatline(edm->fd);
00219
00220
00221 for (c=0; c<nperslice; c++) {
00222 convcnt = fscanf(edm->fd, "%f", cell);
00223 if (convcnt != 1) {
00224 printf("edmplugin) failed reading cell data\n");
00225 printf("edmplugin) cell %d of %d, slice %d\n", c, nperslice, z);
00226 return MOLFILE_ERROR;
00227 }
00228 cell++;
00229 }
00230 eatline(edm->fd);
00231 }
00232
00233
00234 sentinel = 0;
00235 fgets(readbuf, 13, edm->fd);
00236 sscanf(readbuf, "%d", &sentinel);
00237 if (sentinel != -9999) {
00238 printf("edmplugin) EOF sentinel != -9999\n");
00239
00240 }
00241
00242 return MOLFILE_SUCCESS;
00243 }
00244
00245 static void close_edm_read(void *v) {
00246 edm_t *edm = (edm_t *)v;
00247
00248 fclose(edm->fd);
00249 delete [] edm->vol;
00250 delete edm;
00251 }
00252
00253 static void *open_edm_write(const char *filepath, const char *filetype, int natoms) {
00254
00255 FILE *fd = fopen(filepath, "w");
00256 if (!fd) {
00257 fprintf(stderr, "edmplugin) Could not open path '%s' for writing.\n",
00258 filepath);
00259 }
00260 return fd;
00261 }
00262
00263 static void close_edm_write(void *v) {
00264 if (v) {
00265 fclose((FILE *)v);
00266 }
00267 }
00268
00269
00270
00271
00272
00273
00274
00275
00277 float edm_voxel_value_safe(int x, int y, int z, const int xsize, const int ysize, const int zsize, const float *data) {
00278 int xx, yy, zz;
00279 xx = (x > 0) ? ((x < xsize) ? x : xsize-1) : 0;
00280 yy = (y > 0) ? ((y < ysize) ? y : ysize-1) : 0;
00281 zz = (z > 0) ? ((z < zsize) ? z : zsize-1) : 0;
00282 int index = zz*xsize*ysize + yy*xsize + xx;
00283 return data[index];
00284 }
00285
00287 float edm_voxel_value_interpolate(float xv, float yv, float zv, const int xsize, const int ysize, const int zsize, const float *data) {
00288 int x = (int) xv;
00289 int y = (int) yv;
00290 int z = (int) zv;
00291
00292 float xf = xv - x;
00293 float yf = yv - y;
00294 float zf = zv - z;
00295 float xlerps[4];
00296 float ylerps[2];
00297 float tmp;
00298
00299 tmp = edm_voxel_value_safe(x, y, z, xsize, ysize, zsize, data);
00300 xlerps[0] = tmp + xf*(edm_voxel_value_safe(x+1, y, z, xsize, ysize, zsize, data) - tmp);
00301
00302 tmp = edm_voxel_value_safe(x, y+1, z, xsize, ysize, zsize, data);
00303 xlerps[1] = tmp + xf*(edm_voxel_value_safe(x+1, y+1, z, xsize, ysize, zsize, data) - tmp);
00304
00305 tmp = edm_voxel_value_safe(x, y, z+1, xsize, ysize, zsize, data);
00306 xlerps[2] = tmp + xf*(edm_voxel_value_safe(x+1, y, z+1, xsize, ysize, zsize, data) - tmp);
00307
00308 tmp = edm_voxel_value_safe(x, y+1, z+1, xsize, ysize, zsize, data);
00309 xlerps[3] = tmp + xf*(edm_voxel_value_safe(x+1, y+1, z+1, xsize, ysize, zsize, data) - tmp);
00310
00311 ylerps[0] = xlerps[0] + yf*(xlerps[1] - xlerps[0]);
00312 ylerps[1] = xlerps[2] + yf*(xlerps[3] - xlerps[2]);
00313
00314 return ylerps[0] + zf*(ylerps[1] - ylerps[0]);
00315 }
00316
00319 float edm_voxel_value_interpolate_from_coord(float xpos, float ypos, float zpos, const float *origin, const float *xdelta, const float *ydelta, const float *zdelta, const int xsize, const int ysize, const int zsize, float *data) {
00320 xpos = (xpos-origin[0])/xdelta[0];
00321 ypos = (ypos-origin[1])/ydelta[1];
00322 zpos = (zpos-origin[2])/zdelta[2];
00323 int gx = (int) xpos;
00324 int gy = (int) ypos;
00325 int gz = (int) zpos;
00326
00327
00328
00329
00330
00331
00332 if (gx < 0 || gx >= xsize) return 0;
00333 if (gy < 0 || gy >= ysize) return 0;
00334 if (gz < 0 || gz >= zsize) return 0;
00335
00336 return edm_voxel_value_interpolate(xpos, ypos, zpos, xsize, ysize, zsize, data);
00337
00338 }
00339
00340 static int write_edm_data(void *v, molfile_volumetric_t *metadata, float *datablock, float *colorblock) {
00341
00342 FILE *fd = (FILE *)v;
00343 const int xsize = metadata->xsize;
00344 const int ysize = metadata->ysize;
00345 const int zsize = metadata->zsize;
00346
00347 float xaxis[3], yaxis[3], zaxis[3];
00348 float xdelta[3], ydelta[3], zdelta[3];
00349 float origin[3], porigin[3];
00350
00351 int i, j, k;
00352
00353 int na, amin, amax, nb, bmin, bmax, nc, cmin, cmax;
00354 float a, b, c, alpha, beta, gamma;
00355
00356 for (i=0; i<3; i++) {
00357 origin[i] = metadata->origin[i];
00358 xaxis[i] = metadata->xaxis[i];
00359 yaxis[i] = metadata->yaxis[i];
00360 zaxis[i] = metadata->zaxis[i];
00361 xdelta[i] = xaxis[i]/(xsize - 1);
00362 ydelta[i] = yaxis[i]/(ysize - 1);
00363 zdelta[i] = zaxis[i]/(zsize - 1);
00364 }
00365
00366
00367
00368
00369
00370
00371 if (fabs(xaxis[1]) > TOLERANCE || fabs(xaxis[2]) > TOLERANCE ||
00372 fabs(yaxis[0]) > TOLERANCE || fabs(yaxis[2]) > TOLERANCE ||
00373 fabs(zaxis[0]) > TOLERANCE || fabs(zaxis[1]) > TOLERANCE) {
00374 fprintf(stderr, "edmplugin) Could not write X-PLOR file: only orthogonal cells are currently supported.\n");
00375 return MOLFILE_ERROR;
00376 }
00377
00378 amin = (int) floorf(origin[0]/xdelta[0]);
00379 bmin = (int) floorf(origin[1]/ydelta[1]);
00380 cmin = (int) floorf(origin[2]/zdelta[2]);
00381
00382 porigin[0] = amin * xdelta[0];
00383 porigin[1] = bmin * ydelta[1];
00384 porigin[2] = cmin * zdelta[2];
00385
00386 amax = (int) ceilf((xaxis[0]+origin[0])/xdelta[0]);
00387 bmax = (int) ceilf((yaxis[1]+origin[1])/ydelta[1]);
00388 cmax = (int) ceilf((zaxis[2]+origin[2])/zdelta[2]);
00389
00390 na = amax - amin + 1;
00391 nb = bmax - bmin + 1;
00392 nc = cmax - cmin + 1;
00393
00394
00395 a = na*xdelta[0];
00396 b = nb*ydelta[1];
00397 c = nc*zdelta[2];
00398
00399
00400 alpha = beta = gamma = 90;
00401
00402
00403 fprintf(fd,"\n 2 !NTITLE\n");
00404 fprintf(fd,"REMARKS FILENAME=\"\"\n");
00405 fprintf(fd,"REMARKS created by VMD\n");
00406
00407
00408 fprintf(fd,"%d %d %d %d %d %d %d %d %d\n", na, amin, amax, nb, bmin, bmax,
00409 nc, cmin, cmax);
00410 fprintf(fd,"%g %g %g %g %g %g\n", a, b, c, alpha, beta, gamma);
00411
00412
00413 fprintf(fd,"ZYX\n");
00414
00415
00416 int psize = na*nb*nc;
00417 int pxysize = na*nb;
00418
00419
00420 float *pdata = (float*) malloc(psize*sizeof(float));
00421 for (i=0; i<na; i++) {
00422 float xpos = porigin[0] + i*xdelta[0];
00423 for (j=0; j<nb; j++) {
00424 float ypos = porigin[1] + j*ydelta[1];
00425 for (k=0; k<nc; k++) {
00426 float zpos = porigin[2] + k*zdelta[2];
00427 pdata[i + j*na + k*pxysize] = edm_voxel_value_interpolate_from_coord(xpos, ypos, zpos, origin, xdelta, ydelta, zdelta, xsize, ysize, zsize, datablock);
00428 }
00429 }
00430 }
00431
00432
00433 int count = 0;
00434 for (k=0; k<nc; k++) {
00435 if (count % 6) fprintf(fd, "\n");
00436 fprintf(fd, "%8d\n", k);
00437 count=0;
00438 for (j=0; j<nb; j++) {
00439 for (i=0; i<na; i++) {
00440 fprintf(fd, "%12.5e", pdata[k*pxysize + j*na + i]);
00441 if (++count % 6 == 0)
00442 fprintf(fd, "\n");
00443 }
00444 }
00445 }
00446 if (count % 6) fprintf(fd, "\n");
00447
00448
00449 int sentinel = -9999;
00450 fprintf(fd, "%8d\n", sentinel);
00451
00452
00453 double avg = 0;
00454 double stddev = 0;
00455 double sum = 0;
00456 double sum2 = 0;
00457 for (i=0; i<psize; i++) {
00458 sum += pdata[i];
00459 sum2 += pdata[i]*pdata[i];
00460 }
00461 avg = sum/(double)psize;
00462 stddev = psize/(psize-1) * sqrt(sum2/psize - avg*avg);
00463 fprintf(fd, "%g %g\n", avg, stddev);
00464
00465 free(pdata);
00466
00467 fflush(fd);
00468 return MOLFILE_SUCCESS;
00469 }
00470
00471
00472
00473
00474 static molfile_plugin_t plugin;
00475
00476 VMDPLUGIN_API int VMDPLUGIN_init(void) {
00477 memset(&plugin, 0, sizeof(molfile_plugin_t));
00478 plugin.abiversion = vmdplugin_ABIVERSION;
00479 plugin.type = MOLFILE_PLUGIN_TYPE;
00480 plugin.name = "edm";
00481 plugin.prettyname = "XPLOR Electron Density Map";
00482 plugin.author = "John Stone, Leonardo Trabuco";
00483 plugin.majorv = 0;
00484 plugin.minorv = 9;
00485 plugin.is_reentrant = VMDPLUGIN_THREADSAFE;
00486 plugin.filename_extension = "cns,edm,xplor";
00487 plugin.open_file_read = open_edm_read;
00488 plugin.read_volumetric_metadata = read_edm_metadata;
00489 plugin.read_volumetric_data = read_edm_data;
00490 plugin.close_file_read = close_edm_read;
00491 #if vmdplugin_ABIVERSION > 9
00492 plugin.open_file_write = open_edm_write;
00493 plugin.write_volumetric_data = write_edm_data;
00494 plugin.close_file_write = close_edm_write;
00495 #endif
00496 return VMDPLUGIN_SUCCESS;
00497 }
00498
00499 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
00500 (*cb)(v, (vmdplugin_t *)&plugin);
00501 return VMDPLUGIN_SUCCESS;
00502 }
00503
00504 VMDPLUGIN_API int VMDPLUGIN_fini(void) { return VMDPLUGIN_SUCCESS; }
00505
00506
00507
00508 #ifdef TEST_EDMPLUGIN
00509
00510 int main(int argc, char *argv[]) {
00511 int natoms;
00512 void *v;
00513 int i, nsets, set;
00514 molfile_volumetric_t * meta;
00515
00516 while (--argc) {
00517 ++argv;
00518 v = open_edm_read(*argv, "edm", &natoms);
00519 if (!v) {
00520 fprintf(stderr, "open_edm_read failed for file %s\n", *argv);
00521 return 1;
00522 }
00523
00524
00525 if (read_edm_metadata(v, &nsets, &meta)) {
00526 return 1;
00527 }
00528
00529 for (set=0; set<nsets; set++) {
00530 printf("Loading volume set: %d\n", set);
00531
00532 int elements = meta[set].xsize * meta[set].ysize * meta[set].zsize;
00533 printf(" Grid Elements: %d\n", elements);
00534 printf(" Grid dimensions: X: %d Y: %d Z: %d\n",
00535 meta[set].xsize, meta[set].ysize, meta[set].zsize);
00536
00537 float * voldata = (float *) malloc(sizeof(float) * elements);
00538 float * coldata = NULL;
00539
00540 if (meta[set].has_color) {
00541 coldata = (float *) malloc(sizeof(float) * elements * 3);
00542 }
00543
00544
00545 if (read_edm_data(v, set, voldata, coldata)) {
00546 return 1;
00547 }
00548
00549 printf("First 6 elements:\n ");
00550 for (i=0; i<6; i++) {
00551 printf("%g, ", voldata[i]);
00552 }
00553 printf("\n");
00554
00555 printf("Last 6 elements:\n ");
00556 for (i=elements - 6; i<elements; i++) {
00557 printf("%g, ", voldata[i]);
00558 }
00559 printf("\n");
00560 }
00561
00562 close_edm_read(v);
00563 }
00564 return 0;
00565 }
00566
00567 #endif
00568
00569
00570
00571