Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

dxplugin.C

Go to the documentation of this file.
00001 /***************************************************************************
00002  *cr
00003  *cr            (C) Copyright 1995-2016 The Board of Trustees of the
00004  *cr                        University of Illinois
00005  *cr                         All Rights Reserved
00006  *cr
00007  ***************************************************************************/
00008 
00009 /***************************************************************************
00010  * RCS INFORMATION:
00011  *
00012  *      $RCSfile: dxplugin.C,v $
00013  *      $Author: johns $       $Locker:  $             $State: Exp $
00014  *      $Revision: 1.43 $       $Date: 2016/11/28 05:01:53 $
00015  *
00016  ***************************************************************************/
00017 
00018 /* 
00019  * DX potential maps
00020  *
00021  * Format of the file is:
00022  * # Comments
00023  * .
00024  * .
00025  * .
00026  * object 1 class gridpositions counts xn yn zn
00027  * origin xorg yorg zorg
00028  * delta xdel 0 0
00029  * delta 0 ydel 0
00030  * delta 0 0 zdel
00031  * object 2 class gridconnections counts xn yn zn
00032  * object 3 class array type double rank 0 items { xn*yn*zn } [binary] data follows
00033  * f1 f2 f3
00034  * f4 f5 f6 f7 f8 f9
00035  * .
00036  * .
00037  * .
00038  * object "Dataset name" class field
00039  
00040  * Where xn, yn, and zn are the number of data points along each axis;
00041  * xorg, yorg, and zorg is the origin of the grid, assumed to be in angstroms;
00042  * xdel, ydel, and zdel are the scaling factors to convert grid units to
00043  * angstroms.
00044  *
00045  * Grid data follows, with a single or multiple values per line (maximum 
00046  * allowed linelength is hardcoded into the plugin with ~2000 chars), 
00047  * ordered z fast, y medium, and x slow.
00048  *
00049  * Note that the ordering of grid data in VMD's VolumetricData class
00050  * is transposed, i.e. x changes fastest and z slowest! 
00051  * The file reading and writing routines perform the transpose.
00052  *
00053  * If the optional keyword 'binary' is present, the data is expected to
00054  * be in binary, native endian, single precision IEEE-754 floating point format.
00055  *
00056  * Note: this plugin was written to read the OpenDX files created by the
00057  * APBS program, and thus supports only files that are writting in this style.
00058  * the full OpenDX data format is extremely powerful, complex, and flexible.
00059  *
00060  */
00061 
00062 #include <stdlib.h>
00063 #include <stdio.h>
00064 #include <ctype.h>
00065 #include <math.h>
00066 #include <string.h>
00067 
00068 #if defined(_AIX)
00069 #include <strings.h>
00070 #endif
00071 
00072 #if defined(WIN32) || defined(WIN64)
00073 #define strcasecmp  stricmp
00074 #define strncasecmp strnicmp
00075 #endif
00076 
00077 #include "molfile_plugin.h"
00078 #include "largefiles.h"
00079 
00080 #define THISPLUGIN plugin
00081 #include "vmdconio.h"
00082 
00083 #define LINESIZE 2040
00084 
00085 typedef struct {
00086   FILE *fd;
00087   int nsets;
00088   molfile_volumetric_t *vol;
00089   int isBinary; 
00090 } dx_t;
00091 
00092 
00093 // Get a string from a stream, printing any errors that occur
00094 static char *dxgets(char *s, int n, FILE *stream) {
00095   char *returnVal;
00096 
00097   if (feof(stream)) {
00098     vmdcon_printf(VMDCON_ERROR, "dxplugin) Unexpected end-of-file.\n");
00099     return NULL;
00100   } else if (ferror(stream)) {
00101     vmdcon_printf(VMDCON_ERROR, "dxplugin) Error reading file.\n");
00102     return NULL;
00103   } else {
00104     returnVal = fgets(s, n, stream);
00105     if (returnVal == NULL) {
00106       vmdcon_printf(VMDCON_ERROR, "dxplugin) Error reading line.\n");
00107     }
00108   }
00109 
00110   return returnVal;
00111 }
00112 
00113 
00114 static void *open_dx_read(const char *filepath, const char *filetype,
00115     int *natoms) {
00116   FILE *fd;
00117   dx_t *dx;
00118   char inbuf[LINESIZE];
00119   int xsize, ysize, zsize;
00120   float orig[3], xdelta[3], ydelta[3], zdelta[3];
00121   int isBinary = 0;
00122   
00123   fd = fopen(filepath, "rb");
00124   if (!fd) {
00125     vmdcon_printf(VMDCON_ERROR, "dxplugin) Error opening file.\n");
00126     return NULL;
00127   }
00128 
00129   /* skip comments */
00130   do {
00131     if (dxgets(inbuf, LINESIZE, fd) == NULL) 
00132       return NULL;
00133   } while (inbuf[0] == '#');
00134 
00135   /* get the number of grid points */
00136   if (sscanf(inbuf, "object 1 class gridpositions counts %d %d %d", &xsize, &ysize, &zsize) != 3) {
00137     vmdcon_printf(VMDCON_ERROR, "dxplugin) Error reading grid dimensions.\n");
00138     return NULL;
00139   }
00140 
00141   /* get the cell origin */
00142   if (dxgets(inbuf, LINESIZE, fd) == NULL) {
00143     return NULL;
00144   }
00145   if (sscanf(inbuf, "origin %e %e %e", orig, orig+1, orig+2) != 3) {
00146     vmdcon_printf(VMDCON_ERROR, "dxplugin) Error reading grid origin.\n");
00147     return NULL;
00148   }
00149 
00150   /* get the cell dimensions */
00151   if (dxgets(inbuf, LINESIZE, fd) == NULL) {
00152     return NULL;
00153   }
00154   if (sscanf(inbuf, "delta %e %e %e", xdelta, xdelta+1, xdelta+2) != 3) {
00155     vmdcon_printf(VMDCON_ERROR, "dxplugin) Error reading cell x-dimension.\n");
00156     return NULL;
00157   }
00158 
00159   if (dxgets(inbuf, LINESIZE, fd) == NULL) {
00160     return NULL;
00161   }
00162   if (sscanf(inbuf, "delta %e %e %e", ydelta, ydelta+1, ydelta+2) != 3) {
00163     vmdcon_printf(VMDCON_ERROR, "dxplugin) Error reading cell y-dimension.\n");
00164     return NULL;
00165   }
00166 
00167   if (dxgets(inbuf, LINESIZE, fd) == NULL) {
00168     return NULL;
00169   }
00170   if (sscanf(inbuf, "delta %e %e %e", zdelta, zdelta+1, zdelta+2) != 3) {
00171     vmdcon_printf(VMDCON_ERROR, "dxplugin) Error reading cell z-dimension.\n");
00172     return NULL;
00173   }
00174 
00175   /* skip the next line of the header (described at the beginning of
00176    * the code), which aren't utilized by APBS-produced DX files.  */
00177   if (dxgets(inbuf, LINESIZE, fd) == NULL) 
00178     return NULL;
00179   /* The next line tells us whether to expect ascii or binary format.
00180    * We scan for the word 'binary' somewhere in the line, and if it's found
00181    * we assume binary.
00182    */
00183   if (dxgets(inbuf, LINESIZE, fd) == NULL)
00184     return NULL;
00185   if (strstr(inbuf, "binary")) {
00186       isBinary = 1;
00187   }
00188  
00189   /* allocate and initialize the dx structure */
00190   dx = new dx_t;
00191   dx->fd = fd;
00192   dx->vol = NULL;
00193   dx->isBinary = isBinary;
00194   *natoms = MOLFILE_NUMATOMS_NONE;
00195   dx->nsets = 1; /* this file contains only one data set */
00196 
00197   dx->vol = new molfile_volumetric_t[1];
00198   memset(dx->vol, 0, sizeof(molfile_volumetric_t));
00199   strcpy(dx->vol[0].dataname, "DX map");
00200 
00201   /* Set the unit cell origin and basis vectors */
00202   for (int i=0; i<3; i++) {
00203     dx->vol[0].origin[i] = orig[i];
00204     dx->vol[0].xaxis[i] = xdelta[i] * ((xsize-1 > 0) ? (xsize-1) : 1);
00205     dx->vol[0].yaxis[i] = ydelta[i] * ((ysize-1 > 0) ? (ysize-1) : 1);
00206     dx->vol[0].zaxis[i] = zdelta[i] * ((zsize-1 > 0) ? (zsize-1) : 1);
00207   }
00208 
00209   dx->vol[0].xsize = xsize;
00210   dx->vol[0].ysize = ysize;
00211   dx->vol[0].zsize = zsize;
00212 
00213   /* DX files contain no color information. Taken from edmplugin.C */
00214   dx->vol[0].has_color = 0;
00215 
00216   return dx;
00217 }
00218 
00219 static int read_dx_metadata(void *v, int *nsets, 
00220   molfile_volumetric_t **metadata) {
00221   dx_t *dx = (dx_t *)v;
00222   *nsets = dx->nsets; 
00223   *metadata = dx->vol;  
00224 
00225   return MOLFILE_SUCCESS;
00226 }
00227 
00228 static int read_binary_dx_data(dx_t *dx, int set, float *datablock) {
00229     
00230   int i, j, k;
00231   int xsize = dx->vol[0].xsize;
00232   int ysize = dx->vol[0].ysize;
00233   int zsize = dx->vol[0].zsize;
00234   int xysize = xsize * ysize;
00235   size_t total = xysize * zsize;
00236   float *tmp = (float *)malloc(total*sizeof(float));
00237   if (fread(tmp, sizeof(float), total, dx->fd) != total) {
00238     vmdcon_printf(VMDCON_ERROR, "dxplugin) Failed to read %d binary floats\n", total);
00239     free(tmp);
00240     return MOLFILE_ERROR;
00241   }
00242   // take the transpose - nasty
00243   int ind = 0;
00244   for (i=0; i<xsize; i++) {
00245       for (j=0; j<ysize; j++) {
00246           for (k=0; k<zsize; k++) {
00247               datablock[k * xysize + j*xsize + i] = tmp[ind++];
00248           }
00249       }
00250   }
00251   free( tmp );
00252   return MOLFILE_SUCCESS;
00253 }
00254 
00255 static int read_dx_data(void *v, int set, float *datablock,
00256                          float *colorblock) {
00257   dx_t *dx = (dx_t *)v;
00258   FILE *fd = dx->fd;
00259   char inbuf[LINESIZE];
00260   char *p;
00261   float grid;
00262   int x, y, z, xsize, ysize, zsize, xysize, count, total, i, line;
00263 
00264   if (dx->isBinary)
00265     return read_binary_dx_data(dx, set, datablock);
00266 
00267   xsize = dx->vol[0].xsize;
00268   ysize = dx->vol[0].ysize;
00269   zsize = dx->vol[0].zsize;
00270   xysize = xsize * ysize;
00271   total = xysize * zsize;
00272 
00273   /* Read the values from the file */
00274   x = y = z = line = 0;
00275   for (count = 0; count < total;) {
00276     ++line;
00277     p=dxgets(inbuf, LINESIZE, fd);
00278     if (p == NULL) {
00279       vmdcon_printf(VMDCON_ERROR, "dxplugin) Error reading grid data.\n");
00280       vmdcon_printf(VMDCON_ERROR, "dxplugin) line: %d. item: %d/%d. last data: %s\n", 
00281               line, count, total, inbuf);
00282       return MOLFILE_ERROR;
00283     }
00284 
00285     // chop line into whitespace separated tokens and parse them one by one.
00286     while (*p != '\n' && *p != '\0') {
00287 
00288       // skip over whitespace and try to parse non-blank as number
00289       while (*p != '\0' && (*p == ' ' || *p == '\t' || *p == '\n')) ++p;
00290       i = sscanf(p, "%e", &grid);
00291       if (i < 0) break; // end of line/string. get a new one.
00292       
00293       // a 0 return value means non-parsable as number.
00294       if (i == 0) {
00295         vmdcon_printf(VMDCON_ERROR, "dxplugin) Error parsing grid data.\n");
00296         vmdcon_printf(VMDCON_ERROR, "dxplugin) line: %d. item: %d/%d. data %s\n", 
00297                 line, count, total, p);
00298         return MOLFILE_ERROR;
00299       }
00300 
00301       // success! add to dataset.
00302       if (i == 1) {
00303         ++count;
00304         datablock[x + y*xsize + z*xysize] = grid;
00305         z++;
00306         if (z >= zsize) {
00307           z = 0; y++;
00308           if (y >= ysize) {
00309             y = 0; x++;
00310           }
00311         }
00312       }
00313 
00314       // skip over the parsed text and search for next blank.
00315       while (*p != '\0' && *p != ' ' && *p != '\t' && *p != '\n') ++p;
00316     }
00317   }
00318   
00319   char dxname[256];
00320   while (dxgets(inbuf, LINESIZE, dx->fd)) {
00321     if (sscanf(inbuf, "object \"%[^\"]\" class field", dxname) == 1) {
00322       // a dataset name has been found; override the default
00323       strcpy(dx->vol[0].dataname, dxname);
00324       break;
00325     }
00326   }
00327 
00328   return MOLFILE_SUCCESS;
00329 }
00330 
00331 static void close_dx_read(void *v) {
00332   dx_t *dx = (dx_t *)v;
00333   
00334   fclose(dx->fd);
00335   if (dx->vol != NULL)
00336     delete [] dx->vol; 
00337   delete dx;
00338 }
00339 
00340 static void *open_dx_write(const char *filepath, const char *filetype,
00341         int natoms) {
00342 
00343     FILE *fd = fopen(filepath, "wb");
00344     if (!fd) {
00345         vmdcon_printf(VMDCON_ERROR, 
00346                       "dxplugin) Could not open path '%s' for writing.\n",
00347                       filepath);
00348     }
00349     return fd;
00350 }
00351 
00352 static void close_dx_write(void *v) {
00353     if (v) {
00354         fclose((FILE *)v);
00355     }
00356 }
00357 
00358 /*
00359  *
00360 # Data from APBS
00361 # 
00362 # POTENTIAL (kT/e)
00363 # 
00364 object 1 class gridpositions counts 129 129 129
00365 origin -3.075250e+01 -3.848600e+01 -2.908250e+01
00366 delta 4.687500e-01 0.000000e+00 0.000000e+00
00367 delta 0.000000e+00 6.250000e-01 0.000000e+00
00368 delta 0.000000e+00 0.000000e+00 4.687500e-01
00369 object 2 class gridconnections counts 129 129 129
00370 object 3 class array type double rank 0 items 2146689 data follows
00371 
00372 */
00373 static int write_dx_data(void *v, molfile_volumetric_t *metadata,
00374         float *datablock, float *colorblock) {
00375 
00376     int i, j, k, count;
00377     FILE *fd = (FILE *)v;
00378     const int xsize = metadata->xsize;
00379     const int ysize = metadata->ysize;
00380     const int zsize = metadata->zsize;
00381     const int xysize = xsize * ysize;
00382     const int total = xysize * zsize;
00383 
00384     double xdelta[3], ydelta[3], zdelta[3];
00385     for (i=0; i<3; i++) {
00386       xdelta[i] = metadata->xaxis[i]/(xsize - 1);
00387       ydelta[i] = metadata->yaxis[i]/(ysize - 1);
00388       zdelta[i] = metadata->zaxis[i]/(zsize - 1);
00389     }
00390 
00391     fprintf(fd, "# Data from VMD\n");
00392     fprintf(fd, "# %s\n", metadata->dataname);
00393     fprintf(fd, "object 1 class gridpositions counts %d %d %d\n",
00394             xsize, ysize, zsize);
00395     fprintf(fd, "origin %g %g %g\n", 
00396             metadata->origin[0], metadata->origin[1], metadata->origin[2]);
00397     fprintf(fd, "delta %g %g %g\n", 
00398             xdelta[0], xdelta[1], xdelta[2]);
00399     fprintf(fd, "delta %g %g %g\n", 
00400             ydelta[0], ydelta[1], ydelta[2]);
00401     fprintf(fd, "delta %g %g %g\n", 
00402             zdelta[0], zdelta[1], zdelta[2]);
00403     fprintf(fd, "object 2 class gridconnections counts %d %d %d\n",
00404             xsize, ysize, zsize);
00405 
00406     int useBinary = (getenv("VMDBINARYDX") != NULL);
00407     fprintf(fd, "object 3 class array type double rank 0 items %d %sdata follows\n", total, useBinary ? "binary " : "");
00408     count = 0;
00409     for (i=0; i<xsize; i++) {
00410         for (j=0; j<ysize; j++) {
00411             for (k=0; k<zsize; k++) {
00412                 if (useBinary) {
00413                     fwrite(datablock + k*xysize + j*xsize + i, sizeof(float),
00414                             1, fd);
00415                 } else {
00416                     fprintf(fd, "%g ", datablock[k*xysize + j*xsize + i]);
00417                     if (++count == 3) {
00418                         fprintf(fd, "\n");
00419                         count = 0;
00420                     }
00421                 }
00422             }
00423         }
00424     }
00425     if (!useBinary && count) 
00426       fprintf(fd, "\n");
00427 
00428     // Replace any double quotes (") by single quotes (') in the 
00429     // dataname string to make sure that we don't prematurely
00430     // terminate the string in the dx file.
00431     char *squotes = new char[strlen(metadata->dataname)+1];
00432     strcpy(squotes, metadata->dataname);
00433     char *s = squotes;
00434     // while(s=strchr(s, '"')) *s = '\'';
00435     while(1) { 
00436       s=strchr(s, '"'); 
00437       if (s) { 
00438         *s = '\'';
00439       } else { 
00440         break;
00441       }
00442     }
00443 
00444     // Print last line
00445     fprintf(fd, "object \"%s\" class field\n", squotes);
00446     delete [] squotes;
00447 
00448     fflush(fd);
00449     return MOLFILE_SUCCESS;
00450 }
00451 
00452 /*
00453  * Initialization stuff here
00454  */
00455 static molfile_plugin_t dxplugin;
00456 
00457 VMDPLUGIN_API int VMDPLUGIN_init(void) { 
00458   memset(&dxplugin, 0, sizeof(molfile_plugin_t));
00459   dxplugin.abiversion = vmdplugin_ABIVERSION;
00460   dxplugin.type = MOLFILE_PLUGIN_TYPE;
00461   dxplugin.name = "dx";
00462   dxplugin.prettyname = "DX";
00463   dxplugin.author = "Eamon Caddigan, Justin Gullingsrud, John Stone, Leonardo Trabuco";
00464   dxplugin.majorv = 2;
00465   dxplugin.minorv = 0;
00466   dxplugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
00467   dxplugin.filename_extension = "dx";
00468   dxplugin.open_file_read = open_dx_read;
00469   dxplugin.read_volumetric_metadata = read_dx_metadata;
00470   dxplugin.read_volumetric_data = read_dx_data;
00471   dxplugin.close_file_read = close_dx_read;
00472 #if vmdplugin_ABIVERSION > 9
00473   dxplugin.open_file_write = open_dx_write;
00474   dxplugin.write_volumetric_data = write_dx_data;
00475   dxplugin.close_file_write = close_dx_write;
00476 #endif
00477   return VMDPLUGIN_SUCCESS; 
00478 }
00479 
00480 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
00481   (*cb)(v, (vmdplugin_t *)&dxplugin);
00482   return VMDPLUGIN_SUCCESS;
00483 }
00484 
00485 VMDPLUGIN_API int VMDPLUGIN_fini(void) { return VMDPLUGIN_SUCCESS; }
00486 

Generated on Wed Nov 11 03:06:22 2020 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002