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

fs4plugin.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: fs4plugin.C,v $
00013  *      $Author: johns $       $Locker:  $             $State: Exp $
00014  *      $Revision: 1.23 $       $Date: 2016/11/28 05:01:54 $
00015  *
00016  ***************************************************************************/
00017 
00018 /* 
00019  * fsfour format density maps
00020  *
00021  * More info for this format can be found at 
00022  * <http://www.csb.yale.edu/userguides/datamanip/phases/FSFOUR.html>
00023  *
00024  * Old versions of the cns2fsfour and ccp2fsfour utilities produced a
00025  * slightly broken variant of this format, omitting the data that isn't used
00026  * by XTalView. This plugin currently reads these files, but the cell
00027  * dimensions will be incorrect.
00028  *
00029  */
00030 
00031 #include <stdlib.h>
00032 #include <stdio.h>
00033 #include <ctype.h>
00034 #include <math.h>
00035 #include <string.h>
00036 
00037 #if defined(_AIX)
00038 #include <strings.h>
00039 #endif
00040 
00041 #include "molfile_plugin.h"
00042 #include "endianswap.h"
00043 #include "fortread.h"
00044 
00045 #ifndef M_PI
00046 #define M_PI 3.14159265358979323846
00047 #endif
00048 
00049 typedef struct {
00050   FILE *fd;
00051   int nsets;
00052   int swap;
00053   int f, m, s, x, y, z; // indecies mapping fast-medium-slow to x-y-z
00054   float scale;
00055   molfile_volumetric_t *vol;
00056 } fs4_t;
00057 
00058 
00059 static void *open_fs4_read(const char *filepath, const char *filetype,
00060     int *natoms) {
00061   fs4_t *fs4;
00062   FILE *fd;
00063   float header[32], scale, fmsCellSize[3], alpha, beta, gamma, z1, z2, z3;
00064   int dataBegin, blocksize, geom[16], fmsGridSize[3], norn, swap=0;
00065 
00066   fd = fopen(filepath, "rb");
00067   if (!fd) {
00068     fprintf(stderr, "fs4plugin) Error opening file.\n");
00069     return NULL;
00070   }
00071 
00072   // Use the first four-byte integer in the file to determine the file's
00073   // byte-order
00074   fread(&dataBegin, sizeof(int), 1, fd);
00075   if (dataBegin > 255) {
00076     // check if the bytes need to be swapped
00077     swap4_aligned(&dataBegin, 1);
00078     if (dataBegin <= 255) {
00079       swap = 1;
00080     } else {
00081       fprintf(stderr, "fs4plugin) Cannot read file: header block is too large.\n");
00082       return NULL;
00083     }
00084   }
00085 
00086   // Read the header
00087   rewind(fd);
00088   blocksize = fortread_4(header, 32, swap, fd);
00089 
00090   // Handle files produced by old versions of (cns|ccp)2fsfour
00091   if (blocksize == 28) {
00092     printf("fs4plugin) Recognized %s cns2fsfour map.\n", 
00093         swap ? "opposite-endian" : "same-endian");
00094 
00095     // Read the geometry block
00096     blocksize = fortread_4(geom, 16, swap, fd);
00097     if (blocksize != 7) {
00098       fprintf(stderr, "fs4plugin) Incorrect size for geometry block.\n");
00099       return NULL;
00100     }
00101 
00102     fmsGridSize[0] = geom[0];
00103     fmsGridSize[1] = geom[1];
00104     fmsGridSize[2] = geom[2];
00105     norn = geom[4];
00106 
00107     // Warn about assumptions
00108     printf("fs4plugin) Warning: file does not contain unit cell lengths or angles.\n");
00109 
00110     scale = 50.0;
00111     fmsCellSize[0] = 1.0;
00112     fmsCellSize[1] = 1.0;
00113     fmsCellSize[2] = 1.0;
00114     alpha = 90.0;
00115     beta = 90.0;
00116     gamma = 90.0;
00117   }
00118 
00119   // Handle standard fsfour files
00120   else if (blocksize == 31) {
00121     printf("fs4plugin) Recognize standard fsfour map.\n");
00122 
00123     fmsCellSize[0] = header[21];
00124     fmsCellSize[1] = header[22];
00125     fmsCellSize[2] = header[23];
00126     alpha = header[24];
00127     beta = header[25];
00128     gamma = header[26];
00129     
00130     // Skip the symmetry block if one present
00131     blocksize = fortread_4(geom, 16, swap, fd);
00132     if (blocksize == 9) {
00133       printf("fs4plugin) Skipping symmetry block.\n");
00134       blocksize = fortread_4(geom, 16, swap, fd);
00135     }
00136 
00137     // Read the geometry block
00138     if (blocksize != 13) {
00139       fprintf(stderr, "fs4plugin) Incorrect size for geometry block.\n");
00140       return NULL;
00141     }
00142 
00143     fmsGridSize[0] = geom[0];
00144     fmsGridSize[1] = geom[1];
00145     fmsGridSize[2] = geom[2];
00146 
00147     scale = *((float *) geom + 3);
00148     if (scale == 0) {
00149       scale = 50;
00150     }
00151 
00152     norn = geom[4];
00153     if ((norn < 0) || (norn > 2)) {
00154       fprintf(stderr, "fs4plugin) norn out of range.\n");
00155       return NULL;
00156     }
00157   }
00158 
00159   // Unrecognized format
00160   else {
00161     fprintf(stderr, "fs4plugin) Unrecognized map format.\n");
00162     return NULL;
00163   }
00164 
00165   // Convert degrees to radians
00166   alpha *= M_PI / 180.0;
00167   beta  *= M_PI / 180.0;
00168   gamma *= M_PI / 180.0;
00169 
00170   // Warn about assumptions
00171   printf("fs4plugin) Warning: file does not contain molecule center.\nCentering at <0, 0, 0>\n");
00172 
00173   // Allocate and initialize the fs4 structure
00174   fs4 = new fs4_t;
00175   fs4->fd = fd;
00176   fs4->vol = NULL;
00177   *natoms = MOLFILE_NUMATOMS_NONE;
00178   fs4->nsets = 1; // this file contains only one data set
00179   fs4->swap = swap;
00180   fs4->scale = scale;
00181   if (norn == 0) {
00182     // x fast, z medium, y slow
00183     fs4->x = 0;
00184     fs4->y = 2;
00185     fs4->z = 1;
00186     fs4->f = 0;
00187     fs4->m = 2;
00188     fs4->s = 1;
00189   }
00190   else if (norn == 1) {
00191     // y fast, z medium, x slow
00192     fs4->x = 2;
00193     fs4->y = 0;
00194     fs4->z = 1;
00195     fs4->f = 1;
00196     fs4->m = 2;
00197     fs4->s = 0;
00198   }
00199   else { // norn ==2
00200     // x fast, y medium, z slow
00201     fs4->x = 0;
00202     fs4->y = 1;
00203     fs4->z = 2;
00204     fs4->f = 0;
00205     fs4->m = 1;
00206     fs4->s = 2;
00207   }
00208 
00209   fs4->vol = new molfile_volumetric_t[1];
00210   strcpy(fs4->vol[0].dataname, "Fsfour Electron Density Map");
00211 
00212   // Place the origin at {0 0 0}
00213   fs4->vol[0].origin[0] = 0.0;
00214   fs4->vol[0].origin[1] = 0.0;
00215   fs4->vol[0].origin[2] = 0.0;
00216 
00217   fs4->vol[0].xaxis[0] = fmsCellSize[0];
00218   fs4->vol[0].xaxis[1] = 0.0;
00219   fs4->vol[0].xaxis[2] = 0.0;
00220 
00221   fs4->vol[0].yaxis[0] = cos(gamma) * fmsCellSize[1];
00222   fs4->vol[0].yaxis[1] = sin(gamma) * fmsCellSize[1];
00223   fs4->vol[0].yaxis[2] = 0.0;
00224  
00225   z1 = cos(beta);
00226   z2 = (cos(alpha) - cos(beta)*cos(gamma)) / sin(gamma);
00227   z3 = sqrt(1.0 - z1*z1 - z2*z2);
00228   fs4->vol[0].zaxis[0] = z1 * fmsCellSize[2];
00229   fs4->vol[0].zaxis[1] = z2 * fmsCellSize[2];
00230   fs4->vol[0].zaxis[2] = z3 * fmsCellSize[2];
00231 
00232   fs4->vol[0].xsize = fmsGridSize[fs4->x];
00233   fs4->vol[0].ysize = fmsGridSize[fs4->y];
00234   fs4->vol[0].zsize = fmsGridSize[fs4->z];
00235 
00236   fs4->vol[0].has_color = 0;
00237 
00238   return fs4;
00239 }
00240 
00241 
00242 static int read_fs4_metadata(void *v, int *nsets, 
00243   molfile_volumetric_t **metadata) {
00244   fs4_t *fs4 = (fs4_t *)v;
00245   *nsets = fs4->nsets; 
00246   *metadata = fs4->vol;  
00247 
00248   return MOLFILE_SUCCESS;
00249 }
00250 
00251 
00252 static int read_fs4_data(void *v, int set, float *dstBlock, 
00253                          float *colorblock) {
00254   fs4_t *fs4 = (fs4_t *) v;
00255   int *srcBlock, index;
00256   int col, row, plane, colSize, rowSize, planeSize;
00257   int xyzGridSize[3], xyzIndexIncrement[3];
00258 
00259   xyzGridSize[0] = fs4->vol[0].xsize;
00260   xyzGridSize[1] = fs4->vol[0].ysize;
00261   xyzGridSize[2] = fs4->vol[0].zsize;
00262   xyzIndexIncrement[0] = 1;
00263   xyzIndexIncrement[1] = xyzGridSize[0];
00264   xyzIndexIncrement[2] = xyzGridSize[0] * xyzGridSize[1];
00265 
00266   colSize = xyzGridSize[fs4->f];
00267   rowSize = xyzGridSize[fs4->m];
00268   planeSize = xyzGridSize[fs4->s];
00269 
00270   srcBlock = new int[colSize];
00271 
00272   index = 0;
00273   for (plane = 0; plane < planeSize; plane++) {
00274     for (row = 0; row < rowSize; row++) {
00275 
00276       // Read one row of data
00277       if (fortread_4(srcBlock, colSize, fs4->swap, fs4->fd) != colSize) {
00278         fprintf(stderr, "fs4plugin) Error reading data: incorrect record size.\n");
00279         delete [] srcBlock;
00280         return MOLFILE_ERROR;
00281       }
00282 
00283       for (col = 0; col < colSize; col++) {
00284         dstBlock[index] = (float) srcBlock[col] / fs4->scale;
00285         index += xyzIndexIncrement[fs4->f];
00286       }
00287 
00288       index += xyzIndexIncrement[fs4->m] - colSize*xyzIndexIncrement[fs4->f];
00289     } // end for (row)
00290 
00291     index += xyzIndexIncrement[fs4->s] - rowSize*xyzIndexIncrement[fs4->m];
00292   } // end for (plane)
00293 
00294   delete [] srcBlock;
00295   return MOLFILE_SUCCESS;
00296 }
00297 
00298 static void close_fs4_read(void *v) {
00299   fs4_t *fs4 = (fs4_t *)v;
00300 
00301   fclose(fs4->fd);
00302   if (fs4->vol != NULL)
00303     delete [] fs4->vol; 
00304   delete fs4;
00305 }
00306 
00307 /*
00308  * Initialization stuff here
00309  */
00310 static molfile_plugin_t plugin;
00311 
00312 VMDPLUGIN_API int VMDPLUGIN_init(void) { 
00313   memset(&plugin, 0, sizeof(molfile_plugin_t));
00314   plugin.abiversion = vmdplugin_ABIVERSION;
00315   plugin.type = MOLFILE_PLUGIN_TYPE;
00316   plugin.name = "fs";
00317   plugin.prettyname = "FS4 Density Map";
00318   plugin.author = "Eamon Caddigan";
00319   plugin.majorv = 0;
00320   plugin.minorv = 6;
00321   plugin.is_reentrant = VMDPLUGIN_THREADSAFE;
00322   plugin.filename_extension = "fs,fs4";
00323   plugin.open_file_read = open_fs4_read;
00324   plugin.read_volumetric_metadata = read_fs4_metadata;
00325   plugin.read_volumetric_data = read_fs4_data;
00326   plugin.close_file_read = close_fs4_read;
00327   return VMDPLUGIN_SUCCESS; 
00328 }
00329 
00330 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
00331   (*cb)(v, (vmdplugin_t *)&plugin);
00332   return VMDPLUGIN_SUCCESS;
00333 }
00334 
00335 VMDPLUGIN_API int VMDPLUGIN_fini(void) { return VMDPLUGIN_SUCCESS; }
00336 
00337 #ifdef TEST_FS4_PLUGIN
00338 
00339 int main(int argc, char **argv) {
00340   fs4_t *fs4;
00341   int natoms;
00342   char *filetype;
00343   float *datablock;
00344   
00345   while (--argc) {
00346     ++argv;
00347 
00348     fs4 = (fs4_t*) open_fs4_read(*argv, "fs", &natoms);
00349     if (!fs4) {
00350       fprintf(stderr, "open_fs4_read failed for file %s\n", *argv);
00351       return 1;
00352     }
00353 
00354     printf("a:\t%f\nb:\t%f\nc:\t%f\n",
00355            fs4->vol[0].xaxis[0], fs4->vol[0].yaxis[1], fs4->vol[0].zaxis[2]);
00356     printf("ncol:\t%d\nnrow:\t%d\nnplane:\t%d\n", 
00357            fs4->vol[0].xsize, fs4->vol[0].ysize, fs4->vol[0].zsize);
00358 
00359     datablock = new float[fs4->vol[0].xsize * fs4->vol[0].ysize * 
00360                           fs4->vol[0].zsize];
00361 
00362     if (read_fs4_data(fs4, 0, datablock, NULL) != 0) {
00363       fprintf(stderr, "read_fs4_data failed for file %s\n", *argv);
00364       return 1;
00365     }
00366 
00367     delete [] datablock;
00368   }
00369 
00370   return 0;
00371 }
00372 
00373 # endif
00374 

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