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

vaspparchgplugin.c

Go to the documentation of this file.
00001 /***************************************************************************
00002  * RCS INFORMATION:
00003  *
00004  *      $RCSfile: vaspparchgplugin.c,v $
00005  *      $Author: johns $       $Locker:  $             $State: Exp $
00006  *      $Revision: 1.2 $       $Date: 2014/10/10 14:41:01 $
00007  *
00008  ***************************************************************************/
00009 
00010 
00011 /*
00012  *  VASP plugins for VMD
00013  *  Sung Sakong, Dept. of Phys., Univsity Duisburg-Essen
00014  *  
00015  *  VASP manual   
00016  *  http://cms.mpi.univie.ac.at/vasp/
00017  * 
00018  *  LINUX
00019  *  g++ -O2 -Wall -I. -I$VMDBASEDIR/plugins/include -c vaspparchgplugin.c
00020  *  ld -shared -o vaspparchgplugin.so vaspparchgplugin.o
00021  *
00022  *  MACOSX
00023  *  c++ -O2 -Wall -I. -I$VMDBASEDIR/plugins/include -c vaspparchgplugin.c
00024  *  c++ -bundle -o vaspparchgplugin.so vaspparchgplugin.o
00025  *
00026  *  Install
00027  *  copy vaspparchgplugin.so $VMDBASEDIR/plugins/$ARCH/molfile
00028  */
00029 
00030 #include <stdio.h>
00031 #include <stdlib.h>
00032 #include <math.h>
00033 #include <string.h>
00034 
00035 #include "molfile_plugin.h"
00036 #include "vaspplugin.h"
00037 
00038 
00039 static void *open_vaspparchg_read(const char *filename, const char *filetype, int *natoms)
00040 {
00041   vasp_plugindata_t *data;
00042   char lineptr[LINESIZE];
00043   float lc;
00044   int i;
00045 
00046   /* Verify that input is OK */
00047   if (!filename || !natoms) return NULL;
00048 
00049   /* Start with undefined value; set it after successful read */
00050   *natoms = MOLFILE_NUMATOMS_UNKNOWN;
00051 
00052   data = vasp_plugindata_malloc();
00053   if (!data) return NULL;
00054 
00055   /* VASP4 is assumed in default */
00056   data->version = 4;
00057   data->file = fopen(filename, "rb");
00058   if (!data->file) {
00059     vasp_plugindata_free(data);
00060     return NULL;
00061   }
00062 
00063   data->filename = strdup(filename);
00064 
00065   /* Read system title */
00066   fgets(lineptr, LINESIZE, data->file);
00067   data->titleline = strdup(lineptr);
00068 
00069   /* Read lattice constant */
00070   fgets(lineptr, LINESIZE, data->file);
00071   lc = atof(strtok(lineptr, " "));
00072 
00073   /* Read unit cell lattice vectors and multiply by lattice constant */
00074   for(i = 0; i < 3; ++i) {
00075     float x, y, z;
00076     fgets(lineptr, LINESIZE, data->file);
00077     sscanf(lineptr, "%f %f %f", &x, &y, &z);
00078     data->cell[i][0] = x*lc;
00079     data->cell[i][1] = y*lc;
00080     data->cell[i][2] = z*lc;
00081   }
00082 
00083   /* Build rotation matrix */
00084   vasp_buildrotmat(data);
00085 
00086   /* Count number of atoms */
00087   fgets(lineptr, LINESIZE, data->file);
00088   data->numatoms = 0;
00089   for (i = 0; i < MAXATOMTYPES; ++i) {
00090     char const *tmplineptr = strdup(lineptr);
00091     char const *token = (i == 0 ? strtok(lineptr, " ") : strtok(NULL, " "));
00092     int const n = (token ? atoi(token) : -1);
00093 
00094     /* if fails to read number of atoms, then assume VASP5 */
00095     if (i == 0 && n <= 0) {
00096       data->version = 5;
00097       data->titleline =  strdup(tmplineptr);
00098       fgets(lineptr, LINESIZE, data->file);
00099       break;
00100     }else if (n <= 0) break;
00101 
00102     data->eachatom[i] = n;
00103     data->numatoms += n;
00104   }
00105 
00106   if (data->version == 5) {
00107     data->numatoms = 0;
00108     for (i = 0; i < MAXATOMTYPES; ++i) {
00109       char const *token = (i == 0 ? strtok(lineptr, " ") : strtok(NULL, " "));
00110       int const n = (token ? atoi(token) : -1);
00111       
00112       if (n <= 0) break;
00113       
00114       data->eachatom[i] = n;
00115       data->numatoms += n;
00116     }
00117   }
00118 
00119   if (data->numatoms == 0) {
00120     vasp_plugindata_free(data);
00121     fprintf(stderr, "\n\nVASP PARCHG read) ERROR: file '%s' does not contain list of atom numbers.\n", filename);
00122     return NULL;
00123   }
00124 
00125   /* Skip lines up to the grid numbers */
00126   for (i = 0; i < data->numatoms + 2; ++i) fgets(lineptr, LINESIZE, data->file);
00127 
00128   *natoms = data->numatoms;
00129 
00130   return data;
00131 }
00132 
00133 
00134 static int read_vaspparchg_metadata(void *mydata, int *nvolsets, molfile_volumetric_t **metadata)
00135 {
00136   vasp_plugindata_t *data = (vasp_plugindata_t *)mydata;
00137   char lineptr[LINESIZE];
00138   int gridx, gridy, gridz, i;
00139   char const spintext[4][20] = { "spin up+down", "spin up-down", "spin up", "spin down" };
00140 
00141   /* Verify that input is OK */
00142   if (!data || !nvolsets || !metadata) return MOLFILE_ERROR;
00143 
00144   /* Read the grid size */
00145   fgets(lineptr, LINESIZE, data->file);
00146   if (3 != sscanf(lineptr, "%d %d %d", &gridx, &gridy, &gridz)) {
00147      fprintf(stderr, "\n\nVASP PARCHG read) ERROR: file '%s' does not contain grid dimensions.\n", data->filename);
00148      return MOLFILE_ERROR;
00149   }
00150 
00151   fprintf(stderr, "\n\nVASP PARCHG read) found grid data block...\n");
00152 
00153   /* Initialize the volume set list with 4 entries:
00154    * spin up+down : always present
00155    * spin up-down / spin up /spin down : only there for spin-polarized calculations
00156    *                (the latter remain empty for non-spin-polarized calculations)
00157    */
00158   data->nvolsets = 4;
00159   data->vol = (molfile_volumetric_t *)malloc(data->nvolsets * sizeof(molfile_volumetric_t));
00160   if (!data->vol) {
00161      fprintf(stderr, "\n\nVASP PARCHG read) ERROR: Cannot allocate space for volume data.\n");
00162      return MOLFILE_ERROR;
00163   }
00164 
00165   for (i = 0; i < data->nvolsets; ++i) {
00166     molfile_volumetric_t *const set = &(data->vol[i]); /* get a handle to the current volume set meta data */
00167     int k;
00168 
00169     set->has_color = 0;
00170 
00171     /* put volume data name */
00172     sprintf(set->dataname, "Charge density (%s)", spintext[i]);
00173 
00174     set->origin[0] = set->origin[1] = set->origin[2] = 0;
00175     set->xsize = gridx + 1;
00176     set->ysize = gridy + 1;
00177     set->zsize = gridz + 1;
00178 
00179     /* Rotate unit cell vectors */
00180     for (k = 0; k < 3; ++k) {
00181       set->xaxis[k] = data->rotmat[k][0] * data->cell[0][0]
00182                 + data->rotmat[k][1] * data->cell[0][1]
00183                 + data->rotmat[k][2] * data->cell[0][2];
00184       
00185       set->yaxis[k] = data->rotmat[k][0] * data->cell[1][0] 
00186                 + data->rotmat[k][1] * data->cell[1][1]
00187                 + data->rotmat[k][2] * data->cell[1][2];
00188       
00189       set->zaxis[k] = data->rotmat[k][0] * data->cell[2][0] 
00190                 + data->rotmat[k][1] * data->cell[2][1]
00191                 + data->rotmat[k][2] * data->cell[2][2];
00192     }
00193   }
00194 
00195   *nvolsets = data->nvolsets;
00196   *metadata = data->vol;  
00197 
00198   return MOLFILE_SUCCESS;
00199 }
00200 
00201 
00202 static int read_vaspparchg_data(void *mydata, int set, float *datablock, float *colorblock)
00203 {
00204   vasp_plugindata_t *data = (vasp_plugindata_t *)mydata;
00205   char lineptr[LINESIZE];
00206   int chargedensity, error, iset, n;
00207   float volume;
00208 
00209   /* Verify that input is OK */
00210   if (!data || !datablock) return MOLFILE_ERROR;
00211   if (set >= data->nvolsets) return MOLFILE_ERROR;
00212 
00213   if (strstr(data->filename, "LOCPOT") == NULL && strstr(data->filename, "ELFCAR") == NULL) {
00214     chargedensity = 1;
00215     fprintf(stderr, "\nVASP PARCHG read) Charge density is assumed. Each value will be divided by unit cell volume.\n");
00216   } else {
00217     if (set == 1) {
00218       fprintf(stderr, "\n\nVASP PARCHG read) ERROR: ELF or local potential do not include spin difference information.\n");
00219       return MOLFILE_ERROR;
00220     }
00221     chargedensity = 0;
00222     fprintf(stderr, "\nVASP PARCHG read) ELF or local potential is assumed.\n");
00223   }
00224 
00225   volume = fabs(
00226             data->cell[0][0]*(data->cell[1][1]*data->cell[2][2]-data->cell[1][2]*data->cell[2][1])
00227           + data->cell[0][1]*(data->cell[1][2]*data->cell[2][0]-data->cell[1][0]*data->cell[2][2])
00228           + data->cell[0][2]*(data->cell[1][0]*data->cell[2][1]-data->cell[2][0]*data->cell[1][1])
00229                );
00230 
00231   /* Set file pointer to beginning of file and then skip header up to density data */
00232   rewind(data->file);
00233   for (n = 0; n < data->numatoms + data->version + 5; ++n) fgets(lineptr, LINESIZE, data->file);
00234 
00235   for(error = iset = 0; iset <= set && iset < 2 && !error; ++iset) {
00236     char const *dataname = data->vol[iset].dataname;
00237     int const xsize = data->vol[iset].xsize; 
00238     int const ysize = data->vol[iset].ysize;
00239     int const zsize = data->vol[iset].zsize;
00240     int const numberOfDatapoints = (xsize - 1) * (ysize - 1) * (zsize - 1);
00241     int ix, iy, iz;
00242 
00243     fprintf(stderr, "\nVASP PARCHG read) Patience! Reading volume set %d (%d points): %s\n", iset + 1, numberOfDatapoints, dataname);
00244 
00245     for (n = iz = 0; iz < zsize; ++iz) {
00246       for (iy = 0; iy < ysize; ++iy) {
00247         for (ix = 0; ix < xsize; ++ix, ++n) {
00248           float value;
00249           if (ix == xsize - 1) value = datablock[n - ix];
00250           else if (iy == ysize - 1) value = datablock[n - iy*xsize];
00251           else if (iz == zsize - 1) value = datablock[n - iz*ysize*xsize];
00252           else {
00253             if (1 != fscanf(data->file, "%f", &value)) return MOLFILE_ERROR;
00254             if (chargedensity) value /= volume;
00255           }
00256 
00257           /* for set == 2: spin-up   = 0.5 * set0 + 0.5 * set1
00258            * for set == 3: spin-down = 0.5 * set0 - 0.5 * set1 */
00259           if (iset == set) datablock[n] = value;
00260           else if (set >= 2 && iset == 0) datablock[n] = 0.5 * value;
00261           else if (set == 2 && iset == 1) datablock[n] += 0.5 * value;
00262           else if (set == 3 && iset == 1) datablock[n] -= 0.5 * value;
00263         }
00264       }
00265     }
00266 
00267     if(iset == 0){
00268       for (iy = 0; iy < 3; ++iy) {
00269         int ival;
00270         if (1 != fscanf(data->file, "%d", &ival)) error = 2;
00271       }
00272     }
00273 
00274     fprintf(stderr, "\nVASP PARCHG read) %s finished.\n", dataname);
00275   }
00276 
00277   if (error) fprintf(stderr, "\nVASP PARCHG read) PAW-augmentation part is incomplete, but it is ignored anyway.\n");
00278 
00279   return MOLFILE_SUCCESS;
00280 }
00281 
00282 
00283 static void close_vaspparchg_read(void *mydata)
00284 {
00285   vasp_plugindata_t *data = (vasp_plugindata_t *)mydata;
00286   vasp_plugindata_free(data);
00287 }
00288 
00289 
00290 /* registration stuff */
00291 static molfile_plugin_t plugin;
00292 
00293 int VMDPLUGIN_init(void) {
00294   memset(&plugin, 0, sizeof(molfile_plugin_t));
00295   plugin.abiversion = vmdplugin_ABIVERSION;
00296   plugin.type = MOLFILE_PLUGIN_TYPE;
00297   plugin.name = "PARCHG";
00298   plugin.prettyname = "VASP_PARCHG";
00299   plugin.author = "Sung Sakong";
00300   plugin.majorv = 0;
00301   plugin.minorv = 7;
00302   plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
00303   plugin.filename_extension = "PARCHG";
00304   plugin.open_file_read = open_vaspparchg_read;
00305   plugin.close_file_read = close_vaspparchg_read;
00306   plugin.read_volumetric_metadata = read_vaspparchg_metadata;
00307   plugin.read_volumetric_data = read_vaspparchg_data;
00308   return VMDPLUGIN_SUCCESS;
00309 }
00310 
00311 int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
00312   (*cb)(v, (vmdplugin_t *)&plugin);
00313   return VMDPLUGIN_SUCCESS;
00314 }
00315 
00316 int VMDPLUGIN_fini(void) {
00317   return VMDPLUGIN_SUCCESS;
00318 }
00319 

Generated on Fri Mar 29 03:10:19 2024 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002