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

cubeplugin.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: cubeplugin.C,v $
00013  *      $Author: johns $       $Locker:  $             $State: Exp $
00014  *      $Revision: 1.32 $       $Date: 2016/11/28 05:01:53 $
00015  *
00016  ***************************************************************************/
00017 
00018 //
00019 // Plugin reader for Gaussian "cube" files
00020 //
00021 // Gaussian "cube" file format described here (as of 2013-09-24):
00022 // http://www.gaussian.com/g_tech/g_ur/u_cubegen.htm
00023 //
00024 
00025 #include <stdio.h>
00026 #include <stdlib.h>
00027 #include <math.h>
00028 #include <string.h>
00029 
00030 #include "molfile_plugin.h"
00031 #include "unit_conversion.h"
00032 
00033 #ifndef M_PI_2
00034 #define M_PI_2 1.57079632679489661922
00035 #endif
00036 
00037 #include "periodic_table.h"
00038 
00039 #define THISPLUGIN plugin
00040 #include "vmdconio.h"
00041 
00042 static const float bohr = BOHR_TO_ANGS;
00043 
00044 // A format-independent structure to hold unit cell data
00045 typedef struct {
00046   float A, B, C, alpha, beta, gamma;
00047 } cube_box;
00048 
00049 typedef struct {
00050   FILE *fd;              // file descriptor
00051   int nsets;             // number of volume datasets
00052   int numatoms;          // number of atoms
00053   bool coord;            // has coordinate data
00054   long crdpos, datapos;  // seek offsets for coords and data
00055   char *file_name;       // original filename 
00056   float *datacache;      // temporary cache of orbital data prior to conversion
00057   molfile_volumetric_t *vol; // volume data
00058   float origin[3];       // origin, stored for periodic display hack 
00059   float rotmat[3][3];    // rotation matrix, stored for periodic display hack
00060   cube_box box;          // unit cell dimensions
00061 } cube_t;
00062 
00063 // Converts box basis vectors to A, B, C, alpha, beta, and gamma.  
00064 // Stores values in cube_box struct, which should be allocated before calling
00065 // this function.
00066 static int cube_readbox(cube_box *box, float *x, float *y, float *z) {
00067   float A, B, C;
00068 
00069   if (!box) {
00070     return 1;
00071   }
00072 
00073   // provide defaults
00074   box->A = 10.0;
00075   box->B = 10.0;
00076   box->C = 10.0;
00077   box->alpha = 90.0;
00078   box->beta  = 90.0;
00079   box->gamma = 90.0;
00080 
00081   // A, B, C are the lengths of the x, y, z vectors, respectively
00082   A = sqrt( x[0]*x[0] + x[1]*x[1] + x[2]*x[2] );
00083   B = sqrt( y[0]*y[0] + y[1]*y[1] + y[2]*y[2] );
00084   C = sqrt( z[0]*z[0] + z[1]*z[1] + z[2]*z[2] );
00085   if ((A<=0) || (B<=0) || (C<=0)) {
00086     return 1;
00087   }
00088   box->A = A;
00089   box->B = B;
00090   box->C = C;
00091 
00092   // gamma, beta, alpha are the angles between the x & y, x & z, y & z
00093   // vectors, respectively
00094   box->gamma = acos( (x[0]*y[0]+x[1]*y[1]+x[2]*y[2])/(A*B) ) * 90.0/M_PI_2;
00095   box->beta = acos( (x[0]*z[0]+x[1]*z[1]+x[2]*z[2])/(A*C) ) * 90.0/M_PI_2;
00096   box->alpha = acos( (y[0]*z[0]+y[1]*z[1]+y[2]*z[2])/(B*C) ) * 90.0/M_PI_2; 
00097 
00098   return 0;
00099 }
00100 
00101 // calculate and store origin and rotation matrix to realign everything later.
00102 static void cube_buildrotmat(cube_t *cube, float *o, float *a, float *b)
00103 {
00104   // we rotate first around y and z to align a along the x-axis...
00105   const double len   = sqrt(a[0]*a[0] + a[1]*a[1]);
00106   const double phi   = atan2((double) a[2], (double) len);
00107   const double theta = atan2((double) a[1], (double) a[0]);
00108 
00109   const double cph = cos(phi);
00110   const double cth = cos(theta);
00111   const double sph = sin(phi);
00112   const double sth = sin(theta);
00113 
00114   // ...then we rotate around x to put b into the xy-plane.
00115   const double psi = atan2(-sph*cth*b[0] - sph*sth*b[1] + cph*b[2],-sth*b[0] + cth*b[1]);
00116   const double cps = cos(psi);
00117   const double sps = sin(psi);
00118 
00119   const double r[3][3] = { 
00120     {               cph*cth,                    cph*sth,      sph},
00121     {-sth*cps - sph*cth*sps,      cth*cps - sph*sth*sps,  cph*sps},
00122     { sth*sps - sph*cth*cps,     -cth*sps - sph*sth*cps,  cph*cps}
00123   };
00124 
00125   for (int i=0; i<3; ++i) {
00126     cube->origin[i] = o[i];
00127     for (int j=0; j<3; ++j) {
00128       cube->rotmat[i][j] = r[i][j];
00129     }
00130   }
00131 }
00132 
00133 // prototype.
00134 static void close_cube_read(void *v);
00135 
00136 static void *open_cube_read(const char *filepath, const char *filetype,
00137     int *natoms) {
00138   FILE *fd;
00139   cube_t *cube;
00140   int xsize, ysize, zsize;
00141   float a[3], b[3], c[3];
00142   int i;
00143   char readbuf[1024];
00144  
00145   fd = fopen(filepath, "rb");
00146   if (!fd) 
00147     return NULL;
00148 
00149   cube = new cube_t;
00150   cube->fd = fd;
00151   cube->vol = NULL;
00152   cube->coord = false;
00153   cube->file_name = strdup(filepath);
00154   cube->datacache = NULL;
00155 
00156   // initialize origin and rotmat to sensible defaults.
00157   for (i=0; i<3; ++i) {
00158       for (int j=0; j<3; ++j) {
00159           cube->rotmat[i][j] = 0.0;
00160       }
00161   }
00162   for (i=0; i<3; ++i) {
00163       cube->origin[i] = 0.0;
00164       cube->rotmat[i][i] = 1.0;
00165   }
00166 
00167   molfile_volumetric_t voltmpl; // base information for all data sets.
00168 
00169   // read in cube file header information
00170   fgets(readbuf, 1023, cube->fd);    // go on to next line
00171 
00172   // identify this file, and read title string into dataset info
00173   strcpy(voltmpl.dataname, "Gaussian Cube: ");
00174   strncat(voltmpl.dataname, readbuf, 240);      // 240 is max space left after
00175                                                 //   "Gaussian Cube: "
00176   fgets(readbuf, 1023, cube->fd); // skip second header line 
00177 
00178   // read number of atoms and cube origin
00179   // XXX: As of Gaussian09 revD01 this line contains an
00180   // additional, yet undocumented (integer) number
00181   // we read the entire and parse only the known numbers
00182   if ((fgets(readbuf, 255, cube->fd) == NULL) ||
00183       (sscanf(readbuf, "%d%f%f%f",
00184               &cube->numatoms,
00185               &voltmpl.origin[0],
00186               &voltmpl.origin[1],
00187               &voltmpl.origin[2]) != 4 )) {
00188     close_cube_read(cube);
00189     return NULL;
00190   }
00191 
00192   if (cube->numatoms > 0) {   // density cube file
00193     cube->nsets = 1;          // this cube file contains only one data set
00194   } else {
00195     // cube file with orbitals => multiple densities.
00196     cube->numatoms = - cube->numatoms;
00197     cube->nsets = 0;          // we don't know yet how many data sets.
00198   }
00199   *natoms = cube->numatoms;
00200 
00201 
00202   // read in cube axes and sizes
00203   if ((fgets(readbuf, 255, cube->fd) == NULL) ||
00204       (sscanf(readbuf, "%d%f%f%f",&xsize, &a[0], &a[1], &a[2]) != 4)) {
00205     close_cube_read(cube);
00206     return NULL;
00207   }
00208   if ((fgets(readbuf, 255, cube->fd) == NULL) ||
00209       (sscanf(readbuf, "%d%f%f%f",&ysize, &b[0], &b[1], &b[2]) != 4)) {
00210     close_cube_read(cube);
00211     return NULL;
00212   }
00213   if ((fgets(readbuf, 255, cube->fd) == NULL) ||
00214       (sscanf(readbuf, "%d%f%f%f",&zsize, &c[0], &c[1], &c[2]) != 4)) {
00215     close_cube_read(cube);
00216     return NULL;
00217   }
00218 
00219   // calculate number of samples in each dimension
00220   voltmpl.xsize = xsize;
00221   voltmpl.ysize = ysize;
00222   voltmpl.zsize = zsize;
00223   voltmpl.has_color = 0;
00224 
00225   // to make the periodic display work, we need to rotate
00226   // the cellvectors (and the coordinates) in such a way,
00227   // that the a-vector is collinear with the x-axis and
00228   // the b-vector is in the xy-plane. 
00229   cube_buildrotmat(cube, voltmpl.origin, a, b);
00230   // print warning, if the rotation will be significant:
00231   if ((fabs((double) a[1]) + fabs((double) a[2]) + fabs((double) b[2]))
00232       > 0.001) {
00233     vmdcon_printf(VMDCON_WARN, 
00234             "cubeplugin) Coordinates will be rotated to comply \n");
00235     vmdcon_printf(VMDCON_WARN, 
00236             "cubeplugin) with VMD's conventions for periodic display.\n");
00237   }
00238 
00239   // all dimensional units are always in Bohr
00240   // so scale axes and origin correctly.
00241   // NOTE: the angstroms are only allowed in input.
00242   voltmpl.origin[0] *= bohr; 
00243   voltmpl.origin[1] *= bohr;
00244   voltmpl.origin[2] *= bohr;
00245 
00246   // store aligned axes.
00247   for (i=0; i<3; ++i) {
00248     voltmpl.xaxis[i] = cube->rotmat[i][0] * a[0] 
00249       + cube->rotmat[i][1] * a[1] + cube->rotmat[i][2] * a[2];
00250 
00251     voltmpl.yaxis[i] = cube->rotmat[i][0] * b[0] 
00252       + cube->rotmat[i][1] * b[1] + cube->rotmat[i][2] * b[2];
00253     
00254     voltmpl.zaxis[i] = cube->rotmat[i][0] * c[0] 
00255       + cube->rotmat[i][1] * c[1] + cube->rotmat[i][2] * c[2];
00256   }
00257 
00258   voltmpl.xaxis[0] *= bohr * xsize;
00259   voltmpl.xaxis[1] *= bohr * xsize; 
00260   voltmpl.xaxis[2] *= bohr * xsize; 
00261 
00262   voltmpl.yaxis[0] *= bohr * ysize; 
00263   voltmpl.yaxis[1] *= bohr * ysize; 
00264   voltmpl.yaxis[2] *= bohr * ysize; 
00265 
00266   voltmpl.zaxis[0] *= bohr * zsize; 
00267   voltmpl.zaxis[1] *= bohr * zsize; 
00268   voltmpl.zaxis[2] *= bohr * zsize; 
00269 
00270   /*   As of VMD version 1.8.3, volumetric data points are 
00271    *   expected to represent the center of a grid box. cube format 
00272    *   volumetric data represents the value at the edges of the 
00273    *   grid boxes, so we need to shift the internal origin by half
00274    *   a grid box diagonal to have the data at the correct position.
00275    *   This will need to be changed again when the plugin interface 
00276    *   is updated to explicitly allow point/face-centered data sets.
00277    */
00278   voltmpl.origin[0] -= 0.5 * ( voltmpl.xaxis[0] / (double) xsize
00279                     + voltmpl.yaxis[0] / (double) ysize
00280                     + voltmpl.zaxis[0] / (double) zsize);
00281   voltmpl.origin[1] -= 0.5 * ( voltmpl.xaxis[1] / (double) xsize
00282                     + voltmpl.yaxis[1] / (double) ysize
00283                     + voltmpl.zaxis[1] / (double) zsize);
00284   voltmpl.origin[2] -= 0.5 * ( voltmpl.xaxis[2] / (double) xsize
00285                     + voltmpl.yaxis[2] / (double) ysize
00286                     + voltmpl.zaxis[2] / (double) zsize);
00287 
00288 #if defined(TEST_PLUGIN)
00289   printf("cell before rotation:\n");
00290   printf("a: %12.8f %12.8f %12.8f\n", a[0]*xsize*bohr, a[1]*ysize*bohr, a[2]*zsize*bohr);
00291   printf("b: %12.8f %12.8f %12.8f\n", b[0]*xsize*bohr, b[1]*ysize*bohr, b[2]*zsize*bohr);
00292   printf("c: %12.8f %12.8f %12.8f\n", c[0]*xsize*bohr, c[1]*ysize*bohr, c[2]*zsize*bohr);
00293 
00294   printf("cell after rotation:\n");
00295   printf("x: %12.8f %12.8f %12.8f\n", voltmpl.xaxis[0], voltmpl.xaxis[1], voltmpl.xaxis[2]);
00296   printf("y: %12.8f %12.8f %12.8f\n", voltmpl.yaxis[0], voltmpl.yaxis[1], voltmpl.yaxis[2]);
00297   printf("z: %12.8f %12.8f %12.8f\n", voltmpl.zaxis[0], voltmpl.zaxis[1], voltmpl.zaxis[2]);
00298 #endif
00299 
00300   // store the unit cell information for later perusal.
00301   if (cube_readbox(&(cube->box), voltmpl.xaxis, voltmpl.yaxis, voltmpl.zaxis)) {
00302     vmdcon_printf(VMDCON_WARN, "cubeplugin) Calculation of unit cell "
00303                   "size failed. Continuing anyways...\n");
00304   }
00305 
00306   cube->crdpos = ftell(cube->fd); // and record beginning of coordinates
00307   // XXX fseek()/ftell() are incompatible with 64-bit LFS I/O implementations, 
00308   // hope we don't read any files >= 2GB...
00309 
00310   if (cube->nsets >0) { 
00311     int i;
00312 
00313     // density cube file, copy voltmpl into the cube struct.
00314     cube->vol = new molfile_volumetric_t[1];
00315     memcpy(cube->vol, &voltmpl, sizeof(voltmpl));
00316 
00317     // skip over coordinates to find the start of volumetric data
00318     for (i=0; i < cube->numatoms; i++)
00319       fgets(readbuf, 1023, cube->fd);
00320 
00321     cube->datapos = ftell(cube->fd); // and record beginning of data
00322     // XXX fseek()/ftell() are incompatible with 64-bit LFS I/O, 
00323     // hope we don't read any files >= 2GB...
00324   } else {              
00325     int i;
00326 
00327     // orbital cube file. we now have to read the orbitals section
00328     // skip over coordinates
00329     for (i=0; i < cube->numatoms; i++)
00330       fgets(readbuf, 1023, cube->fd);
00331       
00332     fscanf(cube->fd, "%d", &cube->nsets);
00333     vmdcon_printf(VMDCON_INFO, "cubeplugin) found %d orbitals\n", cube->nsets);
00334     cube->vol = new molfile_volumetric_t[cube->nsets];
00335     for (i=0; i < cube->nsets; ++i) {
00336       int orb;
00337       fscanf(cube->fd, "%d", &orb);
00338       memcpy(&cube->vol[i], &voltmpl, sizeof(voltmpl));
00339       sprintf(cube->vol[i].dataname, "Gaussian Cube: Orbital %d", orb);
00340     }
00341       
00342     fgets(readbuf, 1023, cube->fd);  // gobble up rest of line
00343     cube->datapos = ftell(cube->fd); // and record beginning of data
00344     // XXX fseek()/ftell() are incompatible with 64-bit LFS I/O, 
00345     // hope we don't read any files >= 2GB...
00346   }
00347 
00348   return cube;
00349 }
00350 
00351   
00352 static int read_cube_structure(void *v, int *optflags, molfile_atom_t *atoms) {
00353   int i, j;
00354   char *k;
00355   molfile_atom_t *atom;
00356 
00357   cube_t *cube = (cube_t *)v;
00358 
00359   // go to coordinates
00360   fseek(cube->fd, cube->crdpos, SEEK_SET);
00361   // XXX fseek()/ftell() are incompatible with 64-bit LFS I/O implementations, 
00362   // hope we don't read any files >= 2GB...
00363 
00364   /* set mass and radius from PTE, charge from atoms section of cube file. */
00365   *optflags = MOLFILE_ATOMICNUMBER | MOLFILE_MASS | MOLFILE_RADIUS | MOLFILE_CHARGE; 
00366 
00367   for(i=0;i<cube->numatoms;i++) {
00368     int idx;
00369     float chrg;
00370     char fbuffer[1024];
00371 
00372     atom = atoms + i;
00373 
00374     k = fgets(fbuffer, 1024, cube->fd);
00375     j=sscanf(fbuffer, "%d %f %*f %*f %*f", &idx, &chrg);
00376     if (k == NULL) {
00377       vmdcon_printf(VMDCON_ERROR, "cube structure) missing atom(s) in "
00378                                   "file '%s'\n",cube->file_name);
00379       vmdcon_printf(VMDCON_ERROR, "cube structure) expecting '%d' atoms,"
00380                                   " found only '%d'\n",cube->numatoms,i+1);
00381       return MOLFILE_ERROR;
00382     } else if (j < 2) {
00383       vmdcon_printf(VMDCON_INFO, "cube structure) missing atom data in file"
00384                                  " '%s' for atom '%d'\n",cube->file_name,i+1);
00385       return MOLFILE_ERROR;
00386     }
00387 
00388     // assign atom symbol to number. flag unknown or dummy atoms with X.
00389     atom->atomicnumber = idx;
00390     strncpy(atom->name, get_pte_label(idx), sizeof(atom->name));
00391     strncpy(atom->type, atom->name, sizeof(atom->type));
00392     atom->mass = get_pte_mass(idx);
00393     atom->radius = get_pte_vdw_radius(idx);
00394     atom->resname[0] = '\0';
00395     atom->resid = 1;
00396     atom->chain[0] = '\0';
00397     atom->segid[0] = '\0';
00398     atom->charge = chrg;
00399     /* skip to the end of line */
00400   }
00401 
00402   return MOLFILE_SUCCESS;
00403 }
00404 
00405 
00406 static int read_cube_timestep(void *v, int natoms, molfile_timestep_t *ts) {
00407   int i, j, n;
00408   char fbuffer[1024];
00409   float x, y, z;
00410   char *k;
00411   
00412   cube_t *cube = (cube_t *)v;
00413 
00414   // there is only one set of coordinates
00415   if (cube->coord) return MOLFILE_EOF;
00416   cube->coord = true;
00417 
00418   // jump to coordinate position
00419   fseek(cube->fd, cube->crdpos, SEEK_SET);
00420   // XXX fseek()/ftell() are incompatible with 64-bit LFS I/O implementations, 
00421   // hope we don't read any files >= 2GB...
00422  
00423   /* read the coordinates */
00424   for (i=0; i<cube->numatoms; i++) {
00425     k = fgets(fbuffer, 1024, cube->fd);
00426     j = sscanf(fbuffer, "%*d %*f %f %f %f", &x, &y, &z);
00427     
00428     if (k == NULL) {
00429       return MOLFILE_ERROR;
00430     } else if (j < 3) {
00431       vmdcon_printf(VMDCON_ERROR, "cube timestep) missing type or coordinate(s)"
00432                     " in file '%s' for atom '%d'\n",cube->file_name,i+1);
00433       return MOLFILE_ERROR;
00434     } else if (j>=3) {
00435       if (ts != NULL) { 
00436         // Only save coords if we're given a timestep pointer, 
00437         // otherwise assume that VMD wants us to skip past it.
00438         
00439         // In order to make the periodic display work, we need to
00440         // rotate the coordinates around the origin by the stored
00441         // rotation matrix. All coordinates are in Bohr, so they
00442         // must be converted to Angstrom, too.
00443         x -= cube->origin[0];
00444         y -= cube->origin[1];
00445         z -= cube->origin[2];
00446         
00447         for (n=0; n<3; ++n) {
00448           ts->coords[3*i + n] = bohr*(cube->origin[n] 
00449                                       + cube->rotmat[n][0] * x
00450                                       + cube->rotmat[n][1] * y
00451                                       + cube->rotmat[n][2] * z);
00452         }
00453       }
00454     } else {
00455       break;
00456     }
00457   }
00458   // set unitcell dimensions from cached data.
00459   if (ts != NULL) { 
00460       ts->A = cube->box.A;
00461       ts->B = cube->box.B;
00462       ts->C = cube->box.C;
00463       ts->alpha = cube->box.alpha;
00464       ts->beta  = cube->box.beta;
00465       ts->gamma = cube->box.gamma;
00466   }
00467   
00468   return MOLFILE_SUCCESS;
00469 }
00470 
00471 static int read_cube_metadata(void *v, int *nsets, 
00472   molfile_volumetric_t **metadata) {
00473   cube_t *cube = (cube_t *)v;
00474   *nsets = cube->nsets; 
00475   *metadata = cube->vol;  
00476 
00477   return MOLFILE_SUCCESS;
00478 }
00479 
00480 static int read_cube_data(void *v, int set, float *datablock, float *colorblock) {
00481   cube_t *cube = (cube_t *)v;
00482 
00483   vmdcon_printf(VMDCON_INFO, "cubeplugin) trying to read cube data set %d\n", set);
00484 
00485   int xsize = cube->vol[set].xsize; 
00486   int ysize = cube->vol[set].ysize;
00487   int zsize = cube->vol[set].zsize;
00488   int xysize = xsize*ysize;
00489   int nsize = cube->nsets;
00490   int nzsize = nsize*zsize;
00491   int nyzsize = nsize*zsize*ysize;
00492   int x, y, z;
00493 
00494   // go to data
00495   fseek(cube->fd, cube->datapos, SEEK_SET);
00496   // XXX fseek()/ftell() are incompatible with 64-bit LFS I/O implementations, 
00497   // hope we don't read any files >= 2GB...
00498 
00499   // read the data values in 
00500   if (cube->nsets == 1) { // density cube file
00501     for (x=0; x<xsize; x++) {
00502       for (y=0; y<ysize; y++) {
00503         for (z=0; z<zsize; z++) {
00504           if (fscanf(cube->fd, "%f", &datablock[z*xysize + y*xsize + x]) != 1) {
00505               return MOLFILE_ERROR;
00506           }
00507         }
00508       } 
00509     }
00510   } else {
00511     // XXX we may wish to examine this strategy for alternatives that provide
00512     // the same performance but without the extra copy, but it makes sense 
00513     // for now.  
00514 
00515     // Since the orbital cube file stores the data orb1(a1,b1,c1), orb2(a1,b1,c1), 
00516     // ... orbn(a1,b1,c1), orb1(a1,b1,c2), orb2(a1,a1,c2), ..., orbn(ai,bj,ck)
00517     // we have to cache the whole data set of have any kind of reasonable performance.
00518     // otherwise we would have to read (and parse!) the whole file over and over again.
00519     // this way we have to do it only once at the temporary expense of some memory.
00520     if (cube->datacache == NULL) {
00521       int points = xsize*ysize*zsize * nsize; // number of grid cells * nsets
00522       int i;
00523 
00524       // let people know what is going on.
00525       vmdcon_printf(VMDCON_INFO, "cubeplugin) creating %d MByte cube orbital"
00526                     " cache.\n", (int) (points*sizeof(float)) / 1048576);
00527       cube->datacache = new float[points];
00528             
00529       for (i=0; i < points; ++i) {
00530         if (fscanf(cube->fd, "%f", &cube->datacache[i]) != 1) {
00531           return MOLFILE_ERROR;
00532         }
00533 
00534         // print an ascii progress bar so impatient people do not get scared.
00535         // this does not translate well with vmdcon, so we keep it.
00536         // we'd need some way to flag replacing a line instead of writing it.
00537         if ((i % (1048576/sizeof(float))) == 0) {  // one dot per MB.
00538           fprintf(stderr, "."); 
00539         }
00540       }
00541     }
00542       
00543     for (x=0; x<xsize; x++) {
00544       for (y=0; y<ysize; y++) {
00545         for (z=0; z<zsize; z++) {
00546           datablock[z*xysize + y*xsize + x] = cube->datacache[x*nyzsize + y*nzsize + z*nsize + set];
00547         }
00548       }
00549     }
00550   }
00551 
00552   return MOLFILE_SUCCESS;
00553 }
00554 
00555 static void close_cube_read(void *v) {
00556   cube_t *cube = (cube_t *)v;
00557 
00558   fclose(cube->fd);
00559   if (cube->vol) {
00560     delete[] cube->vol; 
00561   }
00562   free(cube->file_name);
00563   if (cube->datacache) { 
00564     vmdcon_printf(VMDCON_INFO, "cubeplugin) freeing cube orbital cache.\n");
00565     delete[] cube->datacache; 
00566   }  
00567 
00568   delete cube;
00569 }
00570 
00571 /*
00572  * Initialization stuff here
00573  */
00574 
00575 int VMDPLUGIN_init(void) { 
00576   memset(&plugin, 0, sizeof(molfile_plugin_t));
00577   plugin.abiversion = vmdplugin_ABIVERSION;
00578   plugin.type = MOLFILE_PLUGIN_TYPE;
00579   plugin.name = "cube";
00580   plugin.prettyname = "Gaussian Cube";
00581   plugin.author = "Axel Kohlmeyer, John Stone";
00582   plugin.majorv = 1;
00583   plugin.minorv = 2;
00584   plugin.is_reentrant = VMDPLUGIN_THREADSAFE;
00585   plugin.filename_extension = "cub,cube";
00586   plugin.open_file_read = open_cube_read;
00587   plugin.read_structure = read_cube_structure;
00588   plugin.read_next_timestep = read_cube_timestep;
00589   plugin.close_file_read = close_cube_read;
00590   plugin.read_volumetric_metadata = read_cube_metadata;
00591   plugin.read_volumetric_data = read_cube_data;
00592   return VMDPLUGIN_SUCCESS; 
00593 }
00594 
00595 int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
00596   (*cb)(v, (vmdplugin_t *)&plugin);
00597   return VMDPLUGIN_SUCCESS;
00598 }
00599 
00600 int VMDPLUGIN_fini(void) { return VMDPLUGIN_SUCCESS; }
00601 
00602 
00603 
00604 #ifdef TEST_PLUGIN
00605 
00606 int main(int argc, char *argv[]) {
00607   int natoms;
00608   void *v;
00609   int i, nsets, set;
00610   molfile_volumetric_t * meta;
00611 
00612   while (--argc) {
00613     ++argv;
00614     v = open_cube_read(*argv, "cube", &natoms);
00615     if (!v) {
00616       fprintf(stderr, "open_cube_read failed for file %s\n", *argv);
00617       return 1;
00618     }
00619     fprintf(stderr, "open_cube_read succeeded for file %s\n", *argv);
00620 
00621     // try loading the EDM metadata now
00622     if (read_cube_metadata(v, &nsets, &meta)) {
00623       return 1; // failed to load cube file
00624     }
00625     fprintf(stderr, "read_cube_metadata succeeded for file %s\n", *argv);
00626 
00627     for (set=0; set<nsets; set++) {
00628       printf("Loading volume set: %d\n", set);   
00629       
00630       int elements = meta[set].xsize * meta[set].ysize * meta[set].zsize;
00631       printf("   Grid Elements: %d\n", elements);
00632       printf(" Grid dimensions: X: %d Y: %d Z: %d\n", 
00633              meta[set].xsize, meta[set].ysize, meta[set].zsize);
00634 
00635       float * voldata = (float *) malloc(sizeof(float) * elements);
00636       float * coldata = NULL;
00637 
00638       if (meta[set].has_color) {
00639         coldata = (float *) malloc(sizeof(float) * elements * 3);
00640       }
00641 
00642       // try loading the data sets now
00643       if (read_cube_data(v, set, voldata, coldata)) {
00644         return 1; // failed to load cube file
00645       }
00646 
00647       printf("First 6 elements:\n   ");
00648       for (i=0; i<6; i++) {
00649         printf("%g, ", voldata[i]);
00650       }
00651       printf("\n"); 
00652 
00653       printf("Last 6 elements:\n   ");
00654       for (i=elements - 6; i<elements; i++) {
00655         printf("%g, ", voldata[i]);
00656       }
00657       printf("\n"); 
00658     }
00659 
00660     close_cube_read(v);
00661   }
00662   return 0;
00663 }
00664 
00665 #endif

Generated on Sat Dec 14 03:10:17 2024 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002