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

uhbdplugin.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: uhbdplugin.C,v $
00013  *      $Author: johns $       $Locker:  $             $State: Exp $
00014  *      $Revision: 1.14 $       $Date: 2016/11/28 05:01:54 $
00015  *
00016  ***************************************************************************/
00017 
00018 /*
00019  * Plugin by Alexander Spaar for reading UHBD grid files
00020  * UHBD related docs:
00021  *   http://adrik.bchs.uh.edu/uhbd.html
00022  *   http://adrik.bchs.uh.edu/uhbd/
00023  *   http://mccammon.ucsd.edu/uhbd.html
00024  */
00025 
00026 #include <stdlib.h>
00027 #include <stdio.h>
00028 #include <ctype.h>
00029 #include <math.h>
00030 #include <string.h>
00031 
00032 #if defined(_AIX)
00033 #include <strings.h>
00034 #endif
00035 
00036 #if defined(WIN32) || defined(WIN64)
00037 #define strcasecmp  stricmp
00038 #define strncasecmp strnicmp
00039 #endif
00040 
00041 #include "molfile_plugin.h"
00042 #include "endianswap.h"
00043 
00044 #define LINESIZE 85
00045 
00046 typedef struct {
00047   FILE *fd;
00048   int nsets;
00049   molfile_volumetric_t *vol;
00050   float scale;  // 0 for ascii files, otherwise signals binary and provides
00051                 // scale factor for written data.
00052   int doswap;   // true if byteswapping needed
00053 } uhbd_t;
00054 
00055 
00056 // Get a string from a stream, printing any errors that occur
00057 static char *uhbdgets(char *s, int n, FILE *stream, const char *msg) {
00058   char *returnVal;
00059 
00060   if (feof(stream)) {
00061     printf("%s", msg);
00062     printf("uhbdplugin) Unexpected end-of-file.\n");
00063     return NULL;
00064   } else if (ferror(stream)) {
00065     printf("%s", msg);
00066     printf("uhbdplugin) Error reading file.\n");
00067     return NULL;
00068   } else {
00069     returnVal = fgets(s, n, stream);
00070     if (returnVal == NULL) {
00071       printf("%s", msg);
00072       printf("uhbdplugin) Encountered EOF or error reading line.\n");
00073     }
00074   }
00075 
00076   return returnVal;
00077 }
00078 
00079 static void *open_uhbd_read(const char *filepath, const char *filetype,
00080     int *natoms) {
00081   FILE *fd;
00082   uhbd_t *uhbd;
00083   char inbuf[LINESIZE];
00084   int xsize, ysize, zsize;
00085   float orig[3], ra, o[3], s[3]; //xdelta[3], ydelta[3], zdelta[3];
00086   int headersize = 0; // for binary format
00087   int doswap = 0; // true if byteswapping needed
00088   float scale=0;    // scale parameter in binary uhbd files
00089   
00090   if ((fd = fopen(filepath, "rb")) == NULL) {
00091     printf("uhbdplugin) Error opening file.\n");
00092     return NULL;
00093   }
00094 
00095   // check if the first part of file is a binary 160, which would indicate
00096   // that this is a binary rather than ascii UHBD file.
00097   fread(&headersize, sizeof(int), 1, fd);
00098   if (headersize == 160) {
00099     printf("uhbdplugin) Detected binary .grd file in native endian\n");
00100     doswap = 0;
00101   } else {
00102     swap4_unaligned(&headersize, 1);
00103     if (headersize == 160) {
00104       printf("uhbdplugin) Detected binary .grd file in opposite endian\n");
00105       doswap = 1;
00106     } else {
00107       headersize = 0;
00108     }
00109   }
00110   if (headersize == 160) {
00111     int iparams[8];
00112     float fparams[4];
00113     // we're doing binary
00114     // read in the entire header, plus the trailing header size
00115     char buf[164];
00116     if (fread(buf, 1, 160, fd) != 160) {
00117       fprintf(stderr, "uhbdplugin) Error: incomplete header in .grd file.\n");
00118       fclose(fd);
00119       return NULL;
00120     }
00121     // format of header is 72 character title, followed by:
00122     // scale dum2 grdflag, idum2 km one km im jm km h ox oy oz
00123     // The first two and last four parameters are floats, the rest ints.
00124     memcpy(&scale,  buf + 72, sizeof(float));
00125     memcpy(iparams, buf + 72 + 8, 32);
00126     memcpy(fparams, buf + 72 + 40, 16);
00127     if (doswap) {
00128       swap4_unaligned(&scale, 1);
00129       swap4_unaligned(iparams, 8);
00130       swap4_unaligned(fparams, 4);
00131     }
00132     xsize = iparams[5];
00133     ysize = iparams[6];
00134     zsize = iparams[7];
00135     ra = fparams[0];
00136     orig[0] = fparams[1];
00137     orig[1] = fparams[2];
00138     orig[2] = fparams[3];
00139 
00140   } else {
00141     rewind(fd);
00142     // Read the header
00143     if (uhbdgets(inbuf, LINESIZE, fd, 
00144         "uhbdplugin) error while skipping header lines\n") == NULL) 
00145       return NULL;
00146     if (uhbdgets(inbuf, LINESIZE, fd,
00147         "uhbdplugin) error while skipping header lines\n") == NULL) 
00148       return NULL;
00149   
00150     /* get grid dimensions, spacing and origin */
00151     if (uhbdgets(inbuf, LINESIZE, fd,
00152         "uhbdplugin) error while getting grid dimensions\n") == NULL) {
00153       return NULL;
00154     }
00155     if (sscanf(inbuf, "%d %d %d %e %e %e %e", &xsize, &ysize, &zsize, &ra, orig, orig+1, orig+2)  != 7) {
00156       printf("uhbdplugin) Error reading grid dimensions, spacing and origin.\n");
00157       return NULL;
00158     }
00159     if (uhbdgets(inbuf, LINESIZE, fd,
00160         "uhbdplugin) error while skipping header lines\n") == NULL) 
00161       return NULL;
00162     if (uhbdgets(inbuf, LINESIZE, fd,
00163         "uhbdplugin) error while skipping header lines\n") == NULL) 
00164       return NULL;
00165   }
00166 
00167   /* allocate and initialize the uhbd structure */
00168   uhbd = new uhbd_t;
00169   uhbd->fd = fd;
00170   uhbd->vol = NULL;
00171   *natoms = MOLFILE_NUMATOMS_NONE;
00172   uhbd->nsets = 1; /* this file contains only one data set */
00173   uhbd->scale = scale; // set by binary files only
00174   uhbd->doswap = doswap;
00175 
00176   uhbd->vol = new molfile_volumetric_t[1];
00177   strcpy(uhbd->vol[0].dataname, 
00178       headersize ? "UHBD binary Electron Density Map"
00179                  : "UHBD ascii Electron Density Map");
00180 
00181   /* Set the unit cell origin and basis vectors */
00182   for (int i=0; i<3; i++) {
00183     uhbd->vol[0].origin[i] = orig[i] + ra;
00184     o[i] = uhbd->vol[0].origin[i];
00185   }
00186 
00187   uhbd->vol[0].xaxis[0] = ra * (xsize-1);
00188   uhbd->vol[0].xaxis[1] = 0;
00189   uhbd->vol[0].xaxis[2] = 0;
00190 
00191   uhbd->vol[0].yaxis[0] = 0;
00192   uhbd->vol[0].yaxis[1] = ra * (ysize-1);
00193   uhbd->vol[0].yaxis[2] = 0;
00194 
00195   uhbd->vol[0].zaxis[0] = 0;
00196   uhbd->vol[0].zaxis[1] = 0;
00197   uhbd->vol[0].zaxis[2] = ra * (zsize-1);
00198 
00199   s[0] = uhbd->vol[0].xaxis[0];
00200   s[1] = uhbd->vol[0].yaxis[1];
00201   s[2] = uhbd->vol[0].zaxis[2];
00202 
00203   uhbd->vol[0].xsize = xsize;
00204   uhbd->vol[0].ysize = ysize;
00205   uhbd->vol[0].zsize = zsize;
00206 
00207   /* UHBD files contain no color information. */
00208   uhbd->vol[0].has_color = 0;
00209 
00210   return uhbd;
00211 }
00212 
00213 static int read_uhbd_metadata(void *v, int *nsets, 
00214   molfile_volumetric_t **metadata) {
00215   uhbd_t *uhbd = (uhbd_t *)v;
00216   *nsets = uhbd->nsets; 
00217   *metadata = uhbd->vol;  
00218 
00219   return MOLFILE_SUCCESS;
00220 }
00221 
00222 static int read_uhbd_data(void *v, int set, float *datablock,
00223                          float *colorblock) {
00224   uhbd_t *uhbd = (uhbd_t *)v;
00225   FILE *fd = uhbd->fd;
00226   char inbuf[LINESIZE];
00227   float grid[6];
00228   int z, xsize, ysize, zsize, xysize, count, count2, total, i;
00229 
00230   xsize = uhbd->vol[0].xsize;
00231   ysize = uhbd->vol[0].ysize;
00232   zsize = uhbd->vol[0].zsize;
00233   xysize = xsize * ysize;
00234   total = xysize * zsize;
00235 
00236   if (uhbd->scale) {
00237     int headerblock[6];
00238     // The binary data is written one z-slice at a time.  First, a 
00239     // block of three ints indicating which slice is next, of the form
00240     // kk, im, jm, where kk is the on-based slice index, and im and jm
00241     // are the x and y dimensions.  We just read past them and also read the
00242     // start of the next block, which is the actual data.
00243     for (z=0; z<zsize; z++) {
00244       if (fread(headerblock, sizeof(int), 6, fd) != 6) {
00245         fprintf(stderr, "uhbdplugin) Error reading header block in binary uhbd file\n");
00246         return MOLFILE_ERROR;
00247       }
00248       // TODO: some sanity checks on the header block.
00249       if (fread(datablock + z*xysize, sizeof(float), xysize, fd) != xysize) {
00250         fprintf(stderr, "uhbdplugin) Error reading data block in binary uhbd file\n");
00251         return MOLFILE_ERROR;
00252       }
00253       // read the trailing block delimiter
00254       fseek(fd, 4, SEEK_CUR);
00255     }
00256     if (uhbd->doswap) {
00257       swap4_aligned(datablock, total);
00258     }
00259     return MOLFILE_SUCCESS;
00260   }
00261 
00262   /* Read the values from the file */
00263   for (z = 0; z < zsize; z++) {
00264     // read header
00265     if (uhbdgets(inbuf, LINESIZE, fd, 
00266         "uhbdplugin) error while getting density plane indices\n") == NULL)
00267       return MOLFILE_ERROR;
00268 
00269     // read data
00270     for (count = 0; count < xysize/6; count++) {
00271       if (uhbdgets(inbuf, LINESIZE, fd,
00272           "uhbdplugin) error while getting density values\n") == NULL)
00273         return MOLFILE_ERROR;
00274 
00275       if (sscanf(inbuf, "%e %e %e %e %e %e", &grid[0], &grid[1], &grid[2], &grid[3], &grid[4], &grid[5]) != 6) {
00276         printf("uhbdplugin) Error reading grid data.\n");
00277         return MOLFILE_ERROR;
00278       }
00279     
00280       for (i = 0; i < 6; i++) { 
00281         datablock[i + count*6 + z*xysize] = grid[i];
00282       }
00283     }
00284 
00285     if ((xysize%6) != 0) {
00286       if (uhbdgets(inbuf, LINESIZE, fd, 
00287           "uhbdplugin) error reading data elements modulo 6\n") == NULL)
00288         return MOLFILE_ERROR;
00289 
00290       count2 = sscanf(inbuf, "%e %e %e %e %e %e", &grid[0], &grid[1], &grid[2], &grid[3], &grid[4], &grid[5]);
00291       if (count2 != (xysize%6)) {
00292         printf("uhbdplugin) Error: incorrect number of data points.\n");
00293         return MOLFILE_ERROR;
00294       }
00295 
00296       for (i = 0; i < count2; i++) {
00297         datablock[i + (count+1)*6 + z*xysize] = grid[i];
00298       }
00299     }
00300   }
00301 
00302   return MOLFILE_SUCCESS;
00303 }
00304 
00305 
00306 static void close_uhbd_read(void *v) {
00307   uhbd_t *uhbd = (uhbd_t *)v;
00308   
00309   fclose(uhbd->fd);
00310   if (uhbd->vol != NULL)
00311     delete [] uhbd->vol; 
00312   delete uhbd;
00313 }
00314 
00315 
00316 /*
00317  * Initialization stuff here
00318  */
00319 static molfile_plugin_t plugin;
00320 
00321 VMDPLUGIN_API int VMDPLUGIN_init(void) {
00322   memset(&plugin, 0, sizeof(molfile_plugin_t));
00323   plugin.abiversion = vmdplugin_ABIVERSION;
00324   plugin.type = MOLFILE_PLUGIN_TYPE;
00325   plugin.name = "uhbd";
00326   plugin.prettyname = "UHBD Grid";
00327   plugin.author = "Alexander Spaar, Justin Gullingsrud";
00328   plugin.majorv = 0;
00329   plugin.minorv = 5;
00330   plugin.is_reentrant = VMDPLUGIN_THREADSAFE;
00331   plugin.filename_extension = "uhbdgrd,grd";
00332   plugin.open_file_read = open_uhbd_read;
00333   plugin.read_volumetric_metadata = read_uhbd_metadata;
00334   plugin.read_volumetric_data = read_uhbd_data;
00335   plugin.close_file_read = close_uhbd_read;
00336   return VMDPLUGIN_SUCCESS;
00337 }
00338 
00339 
00340 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
00341   (*cb)(v, (vmdplugin_t *)&plugin);
00342   return VMDPLUGIN_SUCCESS;
00343 }
00344 
00345 VMDPLUGIN_API int VMDPLUGIN_fini(void) { return VMDPLUGIN_SUCCESS; }
00346 

Generated on Thu Nov 14 02:15:37 2024 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002