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

xsfplugin.C

Go to the documentation of this file.
00001 /***************************************************************************
00002  * RCS INFORMATION:
00003  *
00004  *      $RCSfile: xsfplugin.C,v $
00005  *      $Author: johns $       $Locker:  $             $State: Exp $
00006  *      $Revision: 1.20 $       $Date: 2016/11/09 17:53:21 $
00007  *
00008  ***************************************************************************/
00009 
00010 //
00011 // Molefile plugin for xsf/axsf format files as created by the 
00012 // Quantum Espresso software package (http://www.quantum-espresso.org) 
00013 // a.k.a. PWScf (http://www.pwscf.org/) and XCrySDen (http://www.xcrysden.org/
00014 //
00015 // a file format description is available at:
00016 // http://www.xcrysden.org/doc/XSF.html
00017 //
00018 
00019 #include "largefiles.h"   /* platform dependent 64-bit file I/O defines */
00020 
00021 #include <stdio.h>
00022 #include <stdlib.h>
00023 #include <math.h>
00024 #include <string.h>
00025 
00026 #include "molfile_plugin.h"
00027 
00028 #ifndef M_PI_2
00029 #define M_PI_2 1.57079632679489661922
00030 #endif
00031 
00032 #include "periodic_table.h"
00033 
00034 static const char *xsf_symtab[] = {
00035   "(unknown keyword)", "#",
00036   "BEGIN_INFO", "END_INFO",
00037   "BEGIN_BLOCK_DATAGRID_2D", "END_BLOCK_DATAGRID_2D",
00038   "BEGIN_DATAGRID_2D", "END_DATAGRID_2D",
00039   "BEGIN_BLOCK_DATAGRID_3D", "END_BLOCK_DATAGRID_3D",
00040   "BEGIN_DATAGRID_3D", "END_DATAGRID_3D",
00041   "BEGIN_BLOCK_BANDGRID_3D", "END_BLOCK_BANDGRID_3D",
00042   "ATOMS", "ANIMSTEPS", "BAND",
00043   "MOLECULE", "POLYMER", "SLAB", "CRYSTAL",
00044   "PRIMVEC", "CONVVEC", "PRIMCOORD", "CONVCOORD"
00045 };
00046 
00047 
00048 typedef enum {
00049   xsf_UNKNOWN = 0,  xsf_COMMENT,
00050   xsf_BEGINFO,      xsf_ENDINFO, 
00051   xsf_BEG_2D_BLOCK, xsf_END_2D_BLOCK,
00052   xsf_BEG_2D_DATA,  xsf_END_2D_DATA,
00053   xsf_BEG_3D_BLOCK, xsf_END_3D_BLOCK,
00054   xsf_BEG_3D_DATA,  xsf_END_3D_DATA,
00055   xsf_BEG_3D_BAND,  xsf_END_3D_BAND,
00056   xsf_ATOMS, xsf_ANIMSTEPS, xsf_BAND,
00057   xsf_MOLECULE, xsf_POLYMER, xsf_SLAB, xsf_CRYSTAL,
00058   xsf_PRIMVEC, xsf_CONVVEC, xsf_PRIMCOORD, xsf_CONVCOORD,
00059   xsf_NR_KEYWORDS
00060 } xsf_keyword_t;
00061 
00062 
00063 // list of known alternatives to the keywords above
00064 // last entry has to be an xsf_UNKNOWN.
00065 static const struct {
00066   const char *name;
00067   xsf_keyword_t kw;
00068 }  xsf_aliases[] = {
00069   { "DATAGRID_2D", xsf_BEG_2D_DATA },
00070   { "DATAGRID_3D", xsf_BEG_3D_DATA },
00071   { "BEGIN_BLOCK_DATAGRID2D", xsf_BEG_2D_BLOCK },
00072   { "BEGIN_BLOCK_DATAGRID3D", xsf_BEG_3D_BLOCK },
00073   { "END_BLOCK_DATAGRID2D", xsf_END_2D_BLOCK },
00074   { "END_BLOCK_DATAGRID3D", xsf_END_3D_BLOCK },
00075   { NULL,          xsf_UNKNOWN     }
00076 };
00077 
00078 
00079 static xsf_keyword_t lookup_keyword(const char* word) {
00080   int i, j;
00081   
00082   if (word == 0) return xsf_UNKNOWN;
00083   
00084   // find start of word.
00085   j=0;
00086   for (i=0; i < (int)strlen(word); ++i) {
00087     j=i;
00088     if (!isspace(word[i])) 
00089       break;
00090   }
00091   
00092   for (i=1; i < xsf_NR_KEYWORDS; ++i) {
00093     if (0 == strncmp(word + j, xsf_symtab[i], strlen(xsf_symtab[i])))
00094       return (xsf_keyword_t) i;
00095   }
00096 
00097   // check for known aliases/alternatives
00098   for (i=0; xsf_aliases[i].kw != xsf_UNKNOWN; ++i) {
00099     const char *name = xsf_aliases[i].name;
00100 
00101     if (0 == strncmp(word + j, name, strlen(name)))
00102       return xsf_aliases[i].kw;
00103   }
00104 
00105   return xsf_UNKNOWN;
00106 }
00107 
00108 
00109 
00110 // A format-independent structure to hold unit cell data
00111 typedef struct {
00112   float A, B, C, alpha, beta, gamma, cell[3][3];
00113 } xsf_box;
00114 
00115 typedef struct {
00116   FILE *fd;                     // file descriptor
00117   int nvolsets;                 // number of volumetric datasets
00118   int numatoms;                 // number of atoms
00119   int animsteps;                // for comparison.
00120   int numsteps;                 // number of coordinate sets.
00121   bool coord;                   // has coordinate data
00122   char *file_name;              // original filename 
00123   xsf_keyword_t pbctype;        // type of periodicity (none/polymer/slab/crystal)
00124   molfile_volumetric_t *vol;    // volume set metadata 
00125   int numvolmeta;               // number of entries in *vol
00126   float origin[3];              // origin, stored for periodic display hack 
00127   float rotmat[3][3];           // rotation matrix, stored for periodic display hack
00128   float invmat[3][3];           // reciprocal cell matrix (for PBC wrapping).
00129   xsf_box box;                  // unit cell dimensions (for VMD).
00130 } xsf_t;
00131 
00132 
00133 // Converts box basis vectors to A, B, C, alpha, beta, and gamma.  
00134 // Stores values in xsf_box struct, which should be allocated before calling
00135 // this function.
00136 static int xsf_readbox(xsf_box *box, float *x, float *y, float *z) {
00137   float A, B, C;
00138   int i;
00139 
00140   if (!box) {
00141     return 1;
00142   }
00143 
00144   // provide defaults
00145   box->A = 10.0;
00146   box->B = 10.0;
00147   box->C = 10.0;
00148   box->alpha = 90.0;
00149   box->beta  = 90.0;
00150   box->gamma = 90.0;
00151 
00152   // A, B, C are the lengths of the x, y, z vectors, respectively
00153   A = sqrt( x[0]*x[0] + x[1]*x[1] + x[2]*x[2] );
00154   B = sqrt( y[0]*y[0] + y[1]*y[1] + y[2]*y[2] );
00155   C = sqrt( z[0]*z[0] + z[1]*z[1] + z[2]*z[2] );
00156   if ((A<=0) || (B<=0) || (C<=0)) {
00157     return 1;
00158   }
00159   box->A = A;
00160   box->B = B;
00161   box->C = C;
00162 
00163   // gamma, beta, alpha are the angles between the x & y, x & z, y & z
00164   // vectors, respectively
00165   box->gamma = acos( (x[0]*y[0]+x[1]*y[1]+x[2]*y[2])/(A*B) ) * 90.0/M_PI_2;
00166   box->beta = acos( (x[0]*z[0]+x[1]*z[1]+x[2]*z[2])/(A*C) ) * 90.0/M_PI_2;
00167   box->alpha = acos( (y[0]*z[0]+y[1]*z[1]+y[2]*z[2])/(B*C) ) * 90.0/M_PI_2; 
00168 
00169   // copy original cell vectors for PBC wrapping.
00170   for (i=0; i<3; ++i) {
00171     box->cell[i][0] = x[i];
00172     box->cell[i][1] = y[i];
00173     box->cell[i][2] = z[i];
00174   }
00175   
00176   return 0;
00177 }
00178 
00179 
00180 // calculate and store rotation matrix to realign everything later.
00181 static void xsf_buildrotmat(xsf_t *xsf, float *a, float *b) {
00182   // we rotate first around y and z to align a along the x-axis...
00183   const double len   = sqrt(a[0]*a[0] + a[1]*a[1]);
00184   const double phi   = atan2((double) a[2], (double) len);
00185   const double theta = atan2((double) a[1], (double) a[0]);
00186 
00187   const double cph = cos(phi);
00188   const double cth = cos(theta);
00189   const double sph = sin(phi);
00190   const double sth = sin(theta);
00191 
00192   // ...then we rotate around x to put b into the xy-plane.
00193   const double psi = atan2(-sph*cth*b[0] - sph*sth*b[1] + cph*b[2],-sth*b[0] + cth*b[1]);
00194   const double cps = cos(psi);
00195   const double sps = sin(psi);
00196 
00197   const double r[3][3] = { 
00198     {               cph*cth,                cph*sth,      sph},
00199     {-sth*cps - sph*cth*sps,  cth*cps - sph*sth*sps,  cph*sps},
00200     { sth*sps - sph*cth*cps, -cth*sps - sph*sth*cps,  cph*cps}
00201   };
00202 
00203   for (int i=0; i<3; ++i) {
00204     for (int j=0; j<3; ++j) {
00205       xsf->rotmat[i][j] = r[i][j];
00206     }
00207   }
00208 }
00209 
00210 
00211 static void xsf_buildinvmat(xsf_t *xsf, float *a, float *b, float *c) {
00212   float det, id;
00213   
00214   det = a[0]*b[1]*c[2] + b[0]*c[1]*a[2] + c[0]*a[1]*b[2] 
00215         - a[0]*c[1]*b[2] - b[0]*a[1]*c[2] - c[0]*b[1]*a[2];
00216   
00217   id = 1.0 / det;
00218   xsf->invmat[0][0] = id * ( b[1]*c[2]-b[2]*c[1] );
00219   xsf->invmat[1][0] = id * ( a[2]*c[1]-a[1]*c[2] );
00220   xsf->invmat[2][0] = id * ( a[1]*b[2]-a[2]*b[1] );
00221   xsf->invmat[0][1] = id * ( b[2]*c[0]-b[0]*c[2] );
00222   xsf->invmat[1][1] = id * ( a[0]*c[2]-a[2]*c[0] );
00223   xsf->invmat[2][1] = id * ( a[2]*b[0]-a[0]*b[2] );
00224   xsf->invmat[0][2] = id * ( b[0]*c[1]-b[1]*c[0] );
00225   xsf->invmat[1][2] = id * ( a[1]*c[0]-a[0]*c[1] );
00226   xsf->invmat[2][2] = id * ( a[0]*b[1]-a[1]*b[0] );
00227 }
00228 
00229 
00230 // read a line and forget the data
00231 static void eatline(FILE *fd) {
00232   char readbuf[1025];
00233   fgets(readbuf, 1024, fd);    // go on to next line
00234 }  
00235 
00236 
00237 static bool xsf_read_cell(FILE *fd, float *a, float *b, float *c) {
00238   return (9 == fscanf(fd, "%f%f%f%f%f%f%f%f%f", 
00239                       &a[0],&a[1],&a[2],
00240                       &b[0],&b[1],&b[2],
00241                       &c[0],&c[1],&c[2]));
00242 }
00243 
00244 // forward declaration for cleanup handling...
00245 static void close_xsf_read(void *v);
00246 
00247 static void *open_xsf_read(const char *filepath, const char *filetype,
00248                            int *natoms) {
00249   FILE *fd;
00250   xsf_t *xsf;
00251   int i,j;
00252   
00253   fd = fopen(filepath, "rb");
00254   if (!fd) 
00255     return NULL;
00256 
00257   xsf = new xsf_t;
00258   xsf->fd = fd;
00259   xsf->vol = NULL;
00260   xsf->numvolmeta = 0;
00261   xsf->coord = false;
00262   xsf->nvolsets = 0;
00263   xsf->numatoms = 0;
00264   xsf->numsteps = 0;
00265   xsf->file_name = strdup(filepath);
00266   // this will disable pbc wrapping of coordinates by default
00267   xsf->pbctype = xsf_MOLECULE; 
00268 
00269   // initialize origin and rotmat to sensible defaults.
00270   for (i=0; i<3; ++i) {
00271     for (j=0; j<3; ++j) {
00272       xsf->rotmat[i][j] = 0.0;
00273     }
00274   }
00275   for (i=0; i<3; ++i) {
00276     xsf->origin[i] = 0.0;
00277     xsf->rotmat[i][i] = 1.0;
00278   }
00279 
00280   // since there can be all kinds of data in the file, 
00281   // we start by scanning through the whole file and analyse
00282   // the available info.
00283   char readbuf[256]; // line buffer
00284   xsf_keyword_t kw;
00285 
00286   // we loop until can't read anymore.
00287   do {
00288     if (NULL == fgets(readbuf, 256, xsf->fd)) break;
00289 
00290     again:
00291     kw = lookup_keyword(readbuf);
00292 #ifdef TEST_PLUGIN
00293     fprintf(stderr, "keyword: %d / %s", kw, readbuf);
00294 #endif          
00295     
00296     switch (kw) {
00297       case xsf_ANIMSTEPS:
00298 #ifdef TEST_PLUGIN
00299       {
00300         int n;
00301         if (1 == sscanf(readbuf, "%*s%d", &n)) {
00302           fprintf(stderr, "ANIMSTEPS: found %d steps\n", n);
00303         }
00304       }
00305 #endif          
00306       break;
00307         
00308       case xsf_ATOMS: // no specification for the number of atoms, so we
00309                       // try to figure them out
00310         ++xsf->numsteps;
00311         if (xsf->numatoms == 0) { // count atoms only, if we don't know how many
00312           while (fgets(readbuf, 256, xsf->fd)) {
00313             float x,y,z;
00314             // the coordinate lines are <index> <x> <y> <z> with optional forces.
00315             if (3 == sscanf(readbuf, "%*s%f%f%f", &x, &y, &z)) {
00316               ++xsf->numatoms;
00317             } else {
00318               // we've most likely read the next keyword. reparse buffer.
00319               goto again;
00320               break;
00321             }
00322           }
00323 #ifdef TEST_PLUGIN
00324           fprintf(stderr, "ATOMS: found %d atoms\n", xsf->numatoms);
00325 #endif          
00326         } else { // skip over the lines
00327           int n;
00328           for (n=0; n < xsf->numatoms; ++n) eatline(xsf->fd);
00329         }
00330         break;
00331         
00332       case xsf_PRIMCOORD: // number of atoms is in the next line
00333         if (fgets(readbuf, 256, xsf->fd)) {
00334           int mol, mult;
00335           
00336           if (xsf->numatoms == 0) {
00337             if (2 == sscanf(readbuf, "%d%d", &mol, &mult)) {
00338               xsf->numatoms = mol * mult;
00339             } else {
00340               xsf->numatoms = mol;
00341             }
00342           }
00343           // skip over atom coordinates
00344           int n;
00345           for (n=0; n < xsf->numatoms; ++n) eatline(xsf->fd);
00346           ++xsf->numsteps; 
00347 
00348 #ifdef TEST_PLUGIN
00349           fprintf(stderr, "PRIMCOORD: found %d atoms\n", xsf->numatoms);
00350 #endif          
00351         }
00352         break;
00353         
00354       case xsf_CONVCOORD: // number of atoms is in the next line
00355         if (fgets(readbuf, 256, xsf->fd)) {
00356           int mol, mult, num;
00357           
00358           num = 0;
00359           if (2 == sscanf(readbuf, "%d%d", &mol, &mult)) {
00360             num = mol * mult;
00361           }
00362           
00363           // skip over atom coordinates
00364           int n;
00365           for (n=0; n < num; ++n) eatline(xsf->fd);
00366 #ifdef TEST_PLUGIN
00367           fprintf(stderr, "CONVCOORD: ignoring %d atoms\n", num);
00368 #endif          
00369         }
00370         break;
00371 
00372       case xsf_PRIMVEC: // store primitive cell info for rotation of the volumetric data grid vectors
00373       {
00374         float a[3], b[3], c[3];
00375         
00376         if (xsf_read_cell(xsf->fd, a, b, c)) {
00377           xsf_buildrotmat(xsf, a, b);
00378         } else {
00379           fprintf(stderr, "xsfplugin) WARNING: error reading unit cell. ignoring unit cell info.\n");
00380         }
00381       }
00382       break;
00383       
00384       case xsf_CONVVEC: // ignore conventional cells.
00385       {
00386         int n;
00387         for (n=0; n < 3; ++n)
00388           eatline(xsf->fd);
00389       }
00390       break;
00391 
00392       case xsf_BEG_3D_BLOCK: // analyse number of 3d-data sets
00393         // ordinarily the parsing of the metadata would be done in read_xsf_metadata()
00394         // but since we have to move through the whole file to count its contents,
00395         // it is faster to parse it here and just pass the data later.
00396 
00397         if (xsf->vol == NULL) { // initialize the volume set list with 32 entries
00398           xsf->numvolmeta = 32;
00399           xsf->vol = new molfile_volumetric_t[xsf->numvolmeta];
00400         }
00401         
00402         // next line is title, then check for data blocks
00403         fgets(readbuf, 256, xsf->fd);
00404         printf("xsfplugin) found grid data block: %s", readbuf);
00405 
00406         do { // loop until we reach the end of the whole block 
00407              // or run out of data.
00408           if (NULL == fgets(readbuf, 256, xsf->fd)) break;
00409           switch (lookup_keyword(readbuf)) {
00410             case xsf_BEG_3D_DATA: // fallthrough
00411             {
00412               int n;
00413               molfile_volumetric_t *set;
00414               float a[3], b[3], c[3];
00415               float orig[3];
00416               ++xsf->nvolsets;
00417 
00418               // double the size of the cache for metainfo, if needed
00419               if (xsf->nvolsets > xsf->numvolmeta) {
00420                 molfile_volumetric_t *ptr = xsf->vol;
00421                 xsf->vol = new molfile_volumetric_t[2 * xsf->numvolmeta];
00422                 memcpy((void *)xsf->vol, (void *)ptr, xsf->numvolmeta*sizeof(molfile_volumetric_t));
00423                 xsf->numvolmeta *= 2;
00424                 delete[] ptr;
00425               }
00426 
00427               // get a handle to the current volume set meta data
00428               set = &(xsf->vol[xsf->nvolsets - 1]);
00429               set->has_color = 0;
00430 
00431               // the begin mark is also the title of the data set.
00432               // we need the exact name to later find the start of the data set.
00433               strncpy(set->dataname, readbuf, 255);
00434               
00435               // next is the number of grid points, the origin and
00436               // the spanning vectors of the data block
00437               fgets(readbuf, 256, xsf->fd);
00438               sscanf(readbuf, "%d%d%d", &(set->xsize), &(set->ysize), &(set->zsize));
00439               fgets(readbuf, 256, xsf->fd);
00440               sscanf(readbuf, "%f%f%f", &orig[0], &orig[1], &orig[2]);
00441               fgets(readbuf, 256, xsf->fd);
00442               sscanf(readbuf, "%f%f%f", &a[0], &a[1], &a[2]);
00443               fgets(readbuf, 256, xsf->fd);
00444               sscanf(readbuf, "%f%f%f", &b[0], &b[1], &b[2]);
00445               fgets(readbuf, 256, xsf->fd);
00446               sscanf(readbuf, "%f%f%f", &c[0], &c[1], &c[2]);
00447               
00448               // we need to fix up the size of the data points, since xsf file 
00449               // store the data points at the borders on both sides.
00450               -- set->xsize; -- set->ysize; -- set->zsize;
00451 
00452               // store the realigned axes.
00453               for (n=0; n<3; ++n) {
00454                 set->origin[n] = xsf->rotmat[n][0] * orig[0] 
00455                   + xsf->rotmat[n][1] * orig[1] + xsf->rotmat[n][2] * orig[2];
00456 
00457                 set->xaxis[n] = xsf->rotmat[n][0] * a[0] 
00458                   + xsf->rotmat[n][1] * a[1] + xsf->rotmat[n][2] * a[2];
00459                 
00460                 set->yaxis[n] = xsf->rotmat[n][0] * b[0] 
00461                   + xsf->rotmat[n][1] * b[1] + xsf->rotmat[n][2] * b[2];
00462     
00463                 set->zaxis[n] = xsf->rotmat[n][0] * c[0] 
00464                   + xsf->rotmat[n][1] * c[1] + xsf->rotmat[n][2] * c[2];
00465               }
00466 
00467               do { // loop until we reach the end of the data set
00468                 fgets(readbuf, 256, xsf->fd);
00469               } while (xsf_END_3D_DATA != lookup_keyword(readbuf));
00470 
00471               /*   as of VMD version 1.8.3, volumetric data points are 
00472                *   expected to represent the center of a grid box. xsf format 
00473                *   volumetric data represents the value at the edges of the 
00474                *   grid boxes, so we need to shift the internal origin by half 
00475                *   a grid box diagonal to have the data at the correct position 
00476                *   This will need to be changed again when the plugin interface
00477                *   is updated to explicitly allow point/face-centered data sets.
00478                */
00479               set->origin[0] -= 0.5 * ( set->xaxis[0] / (double) set->xsize
00480                                         + set->yaxis[0] / (double) set->ysize
00481                                         + set->zaxis[0] / (double) set->zsize);
00482               set->origin[1] -= 0.5 * ( set->xaxis[1] / (double) set->xsize
00483                                         + set->yaxis[1] / (double) set->ysize
00484                                         + set->zaxis[1] / (double) set->zsize);
00485               set->origin[2] -= 0.5 * ( set->xaxis[2] / (double) set->xsize
00486                                         + set->yaxis[2] / (double) set->ysize
00487                                         + set->zaxis[2] / (double) set->zsize);
00488             }
00489             break;
00490 
00491             default:
00492               break;
00493           }
00494         } while (xsf_END_3D_BLOCK != lookup_keyword(readbuf));
00495         
00496 #ifdef TEST_PLUGIN
00497         fprintf(stderr, "found %d volumetric data sets\n", xsf->nvolsets);
00498 #endif          
00499         break;
00500 
00501       case xsf_BEG_2D_BLOCK:
00502         do { // skip over data
00503           fgets(readbuf, 256, xsf->fd);
00504         } while (xsf_END_2D_BLOCK != lookup_keyword(readbuf));
00505         break;
00506         
00507         // periodicity encoding. needed for coordinate wrapping.
00508       case xsf_MOLECULE: // fallthrough
00509       case xsf_SLAB:     // fallthrough
00510       case xsf_POLYMER:  // fallthrough
00511       case xsf_CRYSTAL:
00512         xsf->pbctype = kw;
00513         break;
00514 
00515       case xsf_COMMENT:  // fallthrough
00516       case xsf_UNKNOWN:  // fallthrough
00517       default:                  // ignore everything unknown
00518         break;
00519         
00520     }
00521   } while (! (feof(xsf->fd) || ferror(xsf->fd)));
00522 #ifdef TEST_PLUGIN
00523   fprintf(stderr, "total of %d coordinate sets\n", xsf->numsteps);
00524 #endif
00525   rewind(xsf->fd);
00526   *natoms = xsf->numatoms;
00527   return xsf;
00528 }
00529 
00530   
00531 static int read_xsf_structure(void *v, int *optflags, molfile_atom_t *atoms) {
00532   int i;
00533   xsf_t *xsf = (xsf_t *)v;
00534 
00535   // return immediately if there is no structure in this file.
00536   if (xsf->numatoms < 1) return MOLFILE_SUCCESS;
00537 
00538   
00539   // go to beginning of file and find first set of coordinates
00540   rewind(xsf->fd);
00541 
00542   // we loop until we have found, read and parsed the first 
00543   // set of atomic coordinates. we only accept ATOMS and PRIMCOORD
00544   // sections. if we happen to find a PRIMVEC section, too, we use
00545   // that to set the default for animations.
00546   char readbuf[1024]; // line buffer
00547   do {
00548     if (NULL == fgets(readbuf, 256, xsf->fd)) break;
00549 
00550     switch (lookup_keyword(readbuf)) {
00551 
00552       case xsf_PRIMCOORD: // number of atoms is in the next line. skip.
00553         eatline(xsf->fd);
00554         // fallthrough
00555       case xsf_ATOMS: // no specification for the number of atoms, 
00556                       // we can use the same parser for both sections.
00557 
00558         /* we set atom mass and VDW radius from the PTE. */
00559         *optflags = MOLFILE_ATOMICNUMBER | MOLFILE_MASS | MOLFILE_RADIUS;
00560 
00561         for(i=0; i<xsf->numatoms; ++i) {
00562           int j;
00563           char *k;
00564           float coord;
00565           molfile_atom_t *atom;
00566 
00567           char buffer[1024];
00568           k = fgets(readbuf, 1024, xsf->fd);
00569           atom = atoms + i;
00570           j=sscanf(readbuf, "%s %f %f %f", buffer, &coord, &coord, &coord);
00571           if (k == NULL) {
00572             fprintf(stderr, "xsfplugin) structure missing atom(s) in file '%s'\n",xsf->file_name);
00573             fprintf(stderr, "xsfplugin) expecting '%d' atoms, found only '%d'\n",xsf->numatoms,i+1);
00574             return MOLFILE_ERROR;
00575           } else if (j < 4) {
00576             fprintf(stderr, "xsfplugin) missing type or coordinate(s) in file '%s' for atom '%d'\n",
00577                     xsf->file_name, i+1);
00578             return MOLFILE_ERROR;
00579           }
00580 
00581           /* handle the case if the first item is an ordinal number 
00582            * from the PTE */
00583           if (isdigit(buffer[0])) {
00584             int idx;
00585             idx = atoi(buffer);
00586             strncpy(atom->name, get_pte_label(idx), sizeof(atom->name));
00587             atom->atomicnumber = idx;
00588             atom->mass = get_pte_mass(idx);
00589             atom->radius = get_pte_vdw_radius(idx);
00590           } else {
00591             int idx;
00592             strncpy(atom->name, buffer, sizeof(atom->name));
00593             idx = get_pte_idx(buffer);
00594             atom->atomicnumber = idx;
00595             atom->mass = get_pte_mass(idx);
00596             atom->radius = get_pte_vdw_radius(idx);
00597           }
00598           strncpy(atom->type, atom->name, sizeof(atom->type));
00599           atom->resname[0] = '\0';
00600           atom->resid = 1;
00601           atom->chain[0] = '\0';
00602           atom->segid[0] = '\0';
00603 #ifdef TEST_PLUGIN
00604           fprintf(stderr,"xsfplugin) atom %4d: %s  / mass= %f\n", i, atom->name, atom->mass);
00605 #endif
00606         }
00607         
00608         // ok. done. rewind once more and get the hell out of here.
00609         rewind(xsf->fd);
00610         return MOLFILE_SUCCESS;
00611         break;
00612         
00613         // read primitive cell info.
00614       case xsf_PRIMVEC: 
00615       {
00616         float a[3], b[3], c[3];
00617 
00618         if (xsf_read_cell(xsf->fd, a, b, c)) { // ignore unit cell info,if we cannot parse it.
00619           xsf_readbox(&(xsf->box), a, b, c);
00620           xsf_buildrotmat(xsf, a, b);
00621           // print warning, if the rotation will be significant:
00622           if ((fabs((double) a[1]) + fabs((double) a[2]) + fabs((double) b[2]))
00623               > 0.001) {
00624             fprintf(stderr, "xsfplugin) WARNING: Coordinates will be rotated to comply \n"
00625                     "xsfplugin) with VMD's conventions for periodic display...\n");
00626           }
00627           xsf_buildinvmat(xsf, a, b, c);
00628 
00629 #if defined(TEST_PLUGIN)
00630         printf("cell vectors:\n");
00631         printf("<a>: %12.8f %12.8f %12.8f\n", a[0], a[1], a[2]);
00632         printf("<b>: %12.8f %12.8f %12.8f\n", b[0], b[1], b[2]);
00633         printf("<c>: %12.8f %12.8f %12.8f\n", c[0], c[1], c[2]);
00634         printf("cell dimensions:\n");
00635         printf("a= %12.8f   b= %12.8f   c= %12.8f\n", xsf->box.A, xsf->box.B, xsf->box.C);
00636         printf("alpha= %6.2f  beta= %6.2f  gamma= %6.2f\n", xsf->box.alpha, xsf->box.beta, xsf->box.gamma);
00637         printf("reciprocal cell vectors:\n");
00638         printf("i: %12.8f %12.8f %12.8f\n", xsf->invmat[0][0], xsf->invmat[0][1], xsf->invmat[0][2]);
00639         printf("k: %12.8f %12.8f %12.8f\n", xsf->invmat[1][0], xsf->invmat[1][1], xsf->invmat[1][2]);
00640         printf("l: %12.8f %12.8f %12.8f\n", xsf->invmat[2][0], xsf->invmat[2][1], xsf->invmat[2][2]);
00641         printf("cell rotation matrix:\n");
00642         printf("x: %12.8f %12.8f %12.8f\n", xsf->rotmat[0][0], xsf->rotmat[0][1], xsf->rotmat[0][2]);
00643         printf("y: %12.8f %12.8f %12.8f\n", xsf->rotmat[1][0], xsf->rotmat[1][1], xsf->rotmat[1][2]);
00644         printf("z: %12.8f %12.8f %12.8f\n", xsf->rotmat[2][0], xsf->rotmat[2][1], xsf->rotmat[2][2]);
00645 #endif
00646         }
00647       }
00648       break;
00649 
00650       default:                  // ignore everything unknown
00651         break;
00652     }
00653   } while (! (feof(xsf->fd) || ferror(xsf->fd)));
00654 
00655   // if we reach this point, some error must have happened.
00656   return MOLFILE_ERROR;
00657 
00658 }
00659 
00660 
00661 static int read_xsf_timestep(void *v, int natoms, molfile_timestep_t *ts) {
00662   int i;
00663   
00664   xsf_t *xsf = (xsf_t *)v;
00665 
00666   // we loop until we have found, read and parsed the next.
00667   // set of atomic coordinates. we only accept ATOMS and PRIMCOORD
00668   // sections. if we happen to find a PRIMVEC section, too, we use
00669   // that to reset the cell parameters.
00670 
00671   char readbuf[1024]; // line buffer
00672   do {
00673     if (NULL == fgets(readbuf, 256, xsf->fd)) break;
00674 
00675     switch (lookup_keyword(readbuf)) {
00676 
00677       case xsf_PRIMCOORD: // number of atoms is in the next line. skip.
00678         eatline(xsf->fd);
00679         // fallthrough
00680       case xsf_ATOMS: // no specification for the number of atoms, 
00681                       // we can use the same parser for both sections.
00682 
00683         for(i=0; i<natoms; ++i) {
00684           int j, n;
00685           char *k, buffer[1024];
00686           float x, y, z;
00687           
00688           k = fgets(readbuf, 1024, xsf->fd);
00689           j=sscanf(readbuf, "%s %f %f %f", buffer, &x, &y, &z);
00690 
00691           if (k == NULL) {
00692             return MOLFILE_ERROR;
00693           } else if (j < 4) {
00694             fprintf(stderr, "xsfplugin) missing type or coordinate(s) in file '%s' for atom '%d'\n",
00695                     xsf->file_name, i+1);
00696             return MOLFILE_ERROR;
00697           } else if (j>=3) {
00698             if (ts != NULL) { 
00699               // Only save coords if we're given a timestep pointer, 
00700               // otherwise assume that VMD wants us to skip past it.
00701               float xf, yf, zf;
00702               
00703               // apply periodic boundary conditions.
00704 #ifdef TEST_PLUGIN
00705                   printf("wrap PBC: before: %12.6f %12.6f %12.6f\n", x, y, z);
00706 #endif                  
00707               switch(xsf->pbctype) {
00708                 case xsf_CRYSTAL:
00709                   xf = xsf->invmat[0][0] * x + xsf->invmat[0][1] * y + xsf->invmat[0][2] * z;
00710                   xf = xf - floor((double)xf);
00711                   yf = xsf->invmat[1][0] * x + xsf->invmat[1][1] * y + xsf->invmat[1][2] * z;
00712                   yf = yf - floor((double)yf);
00713                   zf = xsf->invmat[2][0] * x + xsf->invmat[2][1] * y + xsf->invmat[2][2] * z;
00714                   zf = zf - floor((double)zf);
00715                   x = xsf->box.cell[0][0] * xf + xsf->box.cell[0][1] * yf + xsf->box.cell[0][2] * zf;
00716                   y = xsf->box.cell[1][0] * xf + xsf->box.cell[1][1] * yf + xsf->box.cell[1][2] * zf;
00717                   z = xsf->box.cell[2][0] * xf + xsf->box.cell[2][1] * yf + xsf->box.cell[2][2] * zf;
00718                   break;
00719 
00720                 case xsf_SLAB:
00721                   xf = xsf->invmat[0][0] * x + xsf->invmat[0][1] * y + xsf->invmat[0][2] * z;
00722                   xf = xf - floor((double)xf);
00723                   yf = xsf->invmat[1][0] * x + xsf->invmat[1][1] * y + xsf->invmat[1][2] * z;
00724                   yf = yf - floor((double)yf);
00725                   zf = xsf->invmat[2][0] * x + xsf->invmat[2][1] * y + xsf->invmat[2][2] * z;
00726                   x = xsf->box.cell[0][0] * xf + xsf->box.cell[0][1] * yf + xsf->box.cell[0][2] * zf;
00727                   y = xsf->box.cell[1][0] * xf + xsf->box.cell[1][1] * yf + xsf->box.cell[1][2] * zf;
00728                   z = xsf->box.cell[2][0] * xf + xsf->box.cell[2][1] * yf + xsf->box.cell[2][2] * zf;
00729                   break;
00730                   
00731                 case xsf_POLYMER:
00732                   xf = xsf->invmat[0][0] * x + xsf->invmat[0][1] * y + xsf->invmat[0][2] * z;
00733                   xf = xf - floor((double)xf);
00734                   yf = xsf->invmat[1][0] * x + xsf->invmat[1][1] * y + xsf->invmat[1][2] * z;
00735                   zf = xsf->invmat[2][0] * x + xsf->invmat[2][1] * y + xsf->invmat[2][2] * z;
00736                   x = xsf->box.cell[0][0] * xf + xsf->box.cell[0][1] * yf + xsf->box.cell[0][2] * zf;
00737                   y = xsf->box.cell[1][0] * xf + xsf->box.cell[1][1] * yf + xsf->box.cell[1][2] * zf;
00738                   z = xsf->box.cell[2][0] * xf + xsf->box.cell[2][1] * yf + xsf->box.cell[2][2] * zf;
00739                   break;
00740                   
00741                 case xsf_MOLECULE:
00742                   xf = x;
00743                   yf = y;
00744                   zf = z;
00745                   break;
00746                   
00747                 default:
00748                   break;
00749               }
00750               
00751 #ifdef TEST_PLUGIN
00752                   printf("wrap PBC: fract:  %12.6f %12.6f %12.6f\n", xf, yf, zf);
00753                   printf("wrap PBC: after: %12.6f %12.6f %12.6f\n", x, y, z);
00754 #endif                  
00755 
00756               // In order to make the periodic display work, we need to
00757               // rotate the coordinates around the origin by the stored
00758               // rotation matrix. for xsf files the origin is not explicitely
00759               // so far, but since it is already initialized to (0,0,0) it
00760               // does no harm, to leave it in here.
00761               x -= xsf->origin[0];
00762               y -= xsf->origin[1];
00763               z -= xsf->origin[2];
00764         
00765               for (n=0; n<3; ++n) {
00766                 ts->coords[3*i + n] = (xsf->origin[n] + xsf->rotmat[n][0] * x
00767                                        + xsf->rotmat[n][1] * y + xsf->rotmat[n][2] * z);
00768               }
00769             }
00770           } else {
00771             break;
00772           }
00773         }
00774         // (re-)set unitcell dimensions
00775         if (ts != NULL) { 
00776           ts->A = xsf->box.A;
00777           ts->B = xsf->box.B;
00778           ts->C = xsf->box.C;
00779           ts->alpha = xsf->box.alpha;
00780           ts->beta  = xsf->box.beta;
00781           ts->gamma = xsf->box.gamma;
00782         }
00783   
00784         return MOLFILE_SUCCESS;
00785         break;
00786         
00787         // read and update primitive cell info.
00788       case xsf_PRIMVEC: 
00789       {
00790         float a[3], b[3], c[3];
00791         
00792         if (xsf_read_cell(xsf->fd, a, b, c)) {
00793           xsf_readbox(&(xsf->box), a,b,c);
00794           xsf_buildrotmat(xsf, a, b);
00795           // print warning, if the rotation will be significant:
00796           if ((fabs((double) a[1]) + fabs((double) a[2]) + fabs((double) b[2]))
00797               > 0.001) {
00798             fprintf(stderr, "xsfplugin) WARNING: Coordinates will be rotated to comply \n"
00799                     "xsfplugin) with VMD's conventions for periodic display...\n");
00800           }
00801           xsf_buildinvmat(xsf, a, b, c);
00802 
00803 #if defined(TEST_PLUGIN)
00804         printf("new cell vectors:\n");
00805         printf("<a>: %12.8f %12.8f %12.8f\n", a[0], a[1], a[2]);
00806         printf("<b>: %12.8f %12.8f %12.8f\n", b[0], b[1], b[2]);
00807         printf("<c>: %12.8f %12.8f %12.8f\n", c[0], c[1], c[2]);
00808         printf("new cell dimensions:\n");
00809         printf("a= %12.8f   b= %12.8f   c= %12.8f\n", xsf->box.A, xsf->box.B, xsf->box.C);
00810         printf("alpha= %6.2f  beta= %6.2f  gamma= %6.2f\n", xsf->box.alpha, xsf->box.beta, xsf->box.gamma);
00811         printf("new reciprocal cell vectors:\n");
00812         printf("i: %12.8f %12.8f %12.8f\n", xsf->invmat[0][0], xsf->invmat[0][1], xsf->invmat[0][2]);
00813         printf("k: %12.8f %12.8f %12.8f\n", xsf->invmat[1][0], xsf->invmat[1][1], xsf->invmat[1][2]);
00814         printf("l: %12.8f %12.8f %12.8f\n", xsf->invmat[2][0], xsf->invmat[2][1], xsf->invmat[2][2]);
00815         printf("new cell rotation matrix:\n");
00816         printf("x: %12.8f %12.8f %12.8f\n", xsf->rotmat[0][0], xsf->rotmat[0][1], xsf->rotmat[0][2]);
00817         printf("y: %12.8f %12.8f %12.8f\n", xsf->rotmat[1][0], xsf->rotmat[1][1], xsf->rotmat[1][2]);
00818         printf("z: %12.8f %12.8f %12.8f\n", xsf->rotmat[2][0], xsf->rotmat[2][1], xsf->rotmat[2][2]);
00819 #endif
00820         }
00821       }
00822       break;
00823 
00824       default:                  // ignore everything unknown
00825         break;
00826     }
00827   } while (! (feof(xsf->fd) || ferror(xsf->fd)));
00828 
00829   // if we reach this point, some error must have happened.
00830   return MOLFILE_ERROR;
00831 }
00832 
00833 static int read_xsf_metadata(void *v, int *nvolsets, 
00834   molfile_volumetric_t **metadata) {
00835   xsf_t *xsf = (xsf_t *)v;
00836   *nvolsets = xsf->nvolsets; 
00837   *metadata = xsf->vol;  
00838 
00839   return MOLFILE_SUCCESS;
00840 }
00841 
00842 static int read_xsf_data(void *v, int set, float *datablock, float *colorblock) {
00843   xsf_t *xsf = (xsf_t *)v;
00844   const char *block = xsf->vol[set].dataname;
00845   
00846   fprintf(stderr, "xsfplugin) trying to read xsf data set %d: %s\n", set, block);
00847   
00848   int xsize = xsf->vol[set].xsize; 
00849   int ysize = xsf->vol[set].ysize;
00850   int zsize = xsf->vol[set].zsize;
00851   int x, y, z;
00852   int n;
00853   char readbuf[1024];
00854   float dummy;
00855   
00856   // find data set ...
00857   rewind(xsf->fd);
00858   do {
00859     if (NULL == fgets(readbuf, 1024, xsf->fd)) return MOLFILE_ERROR;
00860   } while (strncmp(readbuf, block, 1024));
00861   // ... and skip five more lines to get to the beginning of the data
00862   eatline(xsf->fd);
00863   eatline(xsf->fd);
00864   eatline(xsf->fd);
00865   eatline(xsf->fd);
00866   eatline(xsf->fd);
00867   
00868   // read in the data values
00869   // XSF data grids include the points at the border on both sides,
00870   // so we have to read and skip over those on one of them.
00871   n = 0;
00872   for (z=0; z<zsize+1; z++) {
00873     for (y=0; y<ysize+1; y++) {
00874       for (x=0; x<xsize+1; x++) {
00875         if ((x>=xsize) || (y>=ysize) || (z>=zsize)) { 
00876           if (fscanf(xsf->fd, "%f", &dummy) != 1) return MOLFILE_ERROR;
00877         } else {
00878           if (fscanf(xsf->fd, "%f", &datablock[n]) != 1) return MOLFILE_ERROR;
00879           ++n;
00880         }
00881       }
00882     }
00883   }
00884   rewind(xsf->fd);
00885   return MOLFILE_SUCCESS;
00886 }
00887 
00888 static void close_xsf_read(void *v) {
00889   xsf_t *xsf = (xsf_t *)v;
00890 
00891   fclose(xsf->fd);
00892   if (xsf->vol) {
00893     delete[] xsf->vol; 
00894   }
00895   free(xsf->file_name);
00896   delete xsf;
00897 }
00898 
00899 /*
00900  * Initialization stuff here
00901  */
00902 static molfile_plugin_t plugin;
00903 
00904 VMDPLUGIN_API int VMDPLUGIN_init(void) {
00905   memset(&plugin, 0, sizeof(molfile_plugin_t));
00906   plugin.abiversion = vmdplugin_ABIVERSION;
00907   plugin.type = MOLFILE_PLUGIN_TYPE;
00908   plugin.name = "xsf";
00909   plugin.prettyname = "(Animated) XCrySDen Structure File";
00910   plugin.author = "Axel Kohlmeyer, John Stone";
00911   plugin.majorv = 0;
00912   plugin.minorv = 10;
00913   plugin.is_reentrant = VMDPLUGIN_THREADSAFE;
00914   plugin.filename_extension = "axsf,xsf";
00915   plugin.open_file_read = open_xsf_read;
00916   plugin.read_structure = read_xsf_structure;
00917   plugin.read_next_timestep =read_xsf_timestep;
00918   plugin.close_file_read = close_xsf_read;
00919   plugin.read_volumetric_metadata = read_xsf_metadata;
00920   plugin.read_volumetric_data = read_xsf_data;
00921   return VMDPLUGIN_SUCCESS;
00922 }
00923 
00924 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
00925   (*cb)(v, (vmdplugin_t *)&plugin);
00926   return VMDPLUGIN_SUCCESS;
00927 }
00928 
00929 VMDPLUGIN_API int VMDPLUGIN_fini(void) { return VMDPLUGIN_SUCCESS; }
00930 
00931 
00932 #ifdef TEST_PLUGIN
00933 
00934 int main(int argc, char *argv[]) {
00935   int natoms, optflags;
00936   void *v;
00937   int i, nvolsets, set;
00938   molfile_volumetric_t * meta;
00939 
00940 
00941   printf("got index: %d\n", lookup_keyword("     ATOMS  "));
00942   
00943   while (--argc) {
00944     ++argv;
00945     v = open_xsf_read(*argv, "xsf", &natoms);
00946     if (!v) {
00947       fprintf(stderr, "open_xsf_read failed for file %s\n", *argv);
00948       return 1;
00949     }
00950     fprintf(stderr, "open_xsf_read succeeded for file %s\n", *argv);
00951     fprintf(stderr, "input contains %d atoms\n", natoms);
00952 
00953     molfile_atom_t atoms[natoms];
00954     
00955     // try reading the structure info
00956     i = read_xsf_structure(v, &optflags, atoms);
00957     if (i) {
00958       fprintf(stderr, "read_xsf_structure failed for file %s\n", *argv);
00959       return 1;
00960     }
00961     fprintf(stderr, "read_xsf_structure succeeded for file %s\n", *argv);
00962 
00963     // try loading the EDM metadata now
00964     if (read_xsf_metadata(v, &nvolsets, &meta)) {
00965       return 1; // failed to load xsf file
00966     }
00967     fprintf(stderr, "read_xsf_metadata succeeded for file %s\n", *argv);
00968 
00969     for (set=0; set<nvolsets; set++) {
00970       printf("Loading volume set: %d\n", set);   
00971       
00972       int elements = meta[set].xsize * meta[set].ysize * meta[set].zsize;
00973       printf("   Grid Elements: %d\n", elements);
00974       printf(" Grid dimensions: X: %d Y: %d Z: %d\n", 
00975              meta[set].xsize, meta[set].ysize, meta[set].zsize);
00976 
00977       float * voldata = (float *) malloc(sizeof(float) * elements);
00978       float * coldata = NULL;
00979 
00980       if (meta[set].has_color) {
00981         coldata = (float *) malloc(sizeof(float) * elements * 3);
00982       }
00983 
00984       // try loading the data sets now
00985       if (read_xsf_data(v, set, voldata, coldata)) {
00986         return 1; // failed to load xsf file
00987       }
00988 
00989       printf("First 6 elements:\n   ");
00990       for (i=0; i<6; i++) {
00991         printf("%g, ", voldata[i]);
00992       }
00993       printf("\n"); 
00994 
00995       printf("Last 6 elements:\n   ");
00996       for (i=elements - 6; i<elements; i++) {
00997         printf("%g, ", voldata[i]);
00998       }
00999       printf("\n"); 
01000     }
01001 
01002     molfile_timestep_t ts;
01003     ts.coords = new float[3*natoms];
01004     
01005     if (read_xsf_timestep(v, natoms, &ts)) {
01006       printf("read_xsf_timestep 0: failed\n");
01007     } else {
01008       printf("read_xsf_timestep 0: success\n");
01009     }
01010     
01011     if (read_xsf_timestep(v, natoms, &ts)) {
01012       printf("read_xsf_timestep 1: failed\n");
01013     } else {
01014       printf("read_xsf_timestep 1: success\n");
01015     }
01016     
01017     close_xsf_read(v);
01018   }
01019   return 0;
01020 }
01021 
01022 #endif
01023 

Generated on Mon Jan 20 02:56:37 2020 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002