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

vaspchgcarplugin.c

Go to the documentation of this file.
00001 /***************************************************************************
00002  * RCS INFORMATION:
00003  *
00004  *      $RCSfile: vaspchgcarplugin.c,v $
00005  *      $Author: johns $       $Locker:  $             $State: Exp $
00006  *      $Revision: 1.13 $       $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 vaspchgcarplugin.c
00020  *  ld -shared -o vaspchgcarplugin.so vaspchgcarplugin.o
00021  *
00022  *  MACOSX
00023  *  c++ -O2 -Wall -I. -I$VMDBASEDIR/plugins/include -c vaspchgcarplugin.c
00024  *  c++ -bundle -o vaspchgcarplugin.so vaspchgcarplugin.o
00025  *
00026  *  Install
00027  *  copy vaspchgcarplugin.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_vaspchgcar_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 CHGCAR 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_vaspchgcar_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 CHGCAR read) ERROR: file '%s' does not contain grid dimensions.\n", data->filename);
00148      return MOLFILE_ERROR;
00149   }
00150 
00151   fprintf(stderr, "\n\nVASP CHGCAR 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 CHGCAR 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 static int read_vaspchgcar_data(void *mydata, int set, float *datablock, float *colorblock)
00202 {
00203   vasp_plugindata_t *data = (vasp_plugindata_t *)mydata;
00204   char lineptr[LINESIZE];
00205   int chargedensity, error, iset, n;
00206   float volume;
00207 
00208   /* Verify that input is OK */
00209   if (!data || !datablock) return MOLFILE_ERROR;
00210   if (set >= data->nvolsets) return MOLFILE_ERROR;
00211 
00212   if (strstr(data->filename, "LOCPOT") == NULL && strstr(data->filename, "ELFCAR") == NULL) {
00213     chargedensity = 1;
00214     fprintf(stderr, "\nVASP CHGCAR read) Charge density is assumed. Each value will be divided by unit cell volume.\n");
00215   } else {
00216     if (set == 1) {
00217       fprintf(stderr, "\n\nVASP CHGCAR read) ERROR: ELF or local potential do not include spin difference information.\n");
00218       return MOLFILE_ERROR;
00219     }
00220     chargedensity = 0;
00221     fprintf(stderr, "\nVASP CHGCAR read) ELF or local potential is assumed.\n");
00222   }
00223 
00224   volume = fabs(
00225             data->cell[0][0]*(data->cell[1][1]*data->cell[2][2]-data->cell[1][2]*data->cell[2][1])
00226           + data->cell[0][1]*(data->cell[1][2]*data->cell[2][0]-data->cell[1][0]*data->cell[2][2])
00227           + data->cell[0][2]*(data->cell[1][0]*data->cell[2][1]-data->cell[2][0]*data->cell[1][1])
00228                );
00229 
00230   /* Set file pointer to beginning of file and then skip header up to density data */
00231   rewind(data->file);
00232   for (n = 0; n < data->numatoms + data->version + 5; ++n) fgets(lineptr, LINESIZE, data->file);
00233 
00234   for(error = iset = 0; iset <= set && iset < 2 && !error; ++iset) {
00235     char const *dataname = data->vol[iset].dataname;
00236     int const xsize = data->vol[iset].xsize; 
00237     int const ysize = data->vol[iset].ysize;
00238     int const zsize = data->vol[iset].zsize;
00239     int const numberOfDatapoints = (xsize - 1) * (ysize - 1) * (zsize - 1);
00240     int ix, iy, iz;
00241 
00242     fprintf(stderr, "\nVASP CHGCAR read) Patience! Reading volume set %d (%d points): %s\n", iset + 1, numberOfDatapoints, dataname);
00243 
00244     for (n = iz = 0; iz < zsize; ++iz) {
00245       for (iy = 0; iy < ysize; ++iy) {
00246         for (ix = 0; ix < xsize; ++ix, ++n) {
00247           float value;
00248           if (ix == xsize - 1) value = datablock[n - ix];
00249           else if (iy == ysize - 1) value = datablock[n - iy*xsize];
00250           else if (iz == zsize - 1) value = datablock[n - iz*ysize*xsize];
00251           else {
00252             if (1 != fscanf(data->file, "%f", &value)) return MOLFILE_ERROR;
00253             if (chargedensity) value /= volume;
00254           }
00255 
00256           /* for set == 2: spin-up   = 0.5 * set0 + 0.5 * set1
00257            * for set == 3: spin-down = 0.5 * set0 - 0.5 * set1 */
00258           if (iset == set) datablock[n] = value;
00259           else if (set >= 2 && iset == 0) datablock[n] = 0.5 * value;
00260           else if (set == 2 && iset == 1) datablock[n] += 0.5 * value;
00261           else if (set == 3 && iset == 1) datablock[n] -= 0.5 * value;
00262         }
00263       }
00264     }
00265 
00266     /* Skip paw-augmentation part 
00267      * augmentation parts are found only in CHGCAR not in PARCHG
00268      * I have no good idea for that */
00269     for (iy = 0; iy < data->numatoms; ++iy) {
00270       int numaug;
00271       if (1 != fscanf(data->file, "%*s %*s %*d %d", &numaug)) error = 1;
00272       
00273       for (ix = 0; ix < numaug && !error; ++ix) {
00274         float val;
00275         if (1 != fscanf(data->file, "%f", &val)) error = 2;
00276       }
00277       fgets(lineptr, LINESIZE, data->file);
00278     }
00279 
00280     /* After the charge density data there are floating numbers (times of atoms)
00281      * and a line with three grid integers, all of which should be skipped. */
00282     if(iset==0){
00283       for (iy = 0; iy < data->numatoms; ++iy) {
00284         float val;
00285         if (1 != fscanf(data->file, "%f", &val)) error = 2;
00286       }
00287       for (iy = 0; iy < 3; ++iy) {
00288         int ival;
00289         if (1 != fscanf(data->file, "%d", &ival)) error = 2;
00290       }
00291     }
00292 
00293     fprintf(stderr, "\nVASP CHGCAR read) %s finished.\n", dataname);
00294   }
00295 
00296   if (error) fprintf(stderr, "\nVASP CHGCAR read) PAW-augmentation part is incomplete, but it is ignored anyway.\n");
00297 
00298   return MOLFILE_SUCCESS;
00299 }
00300 
00301 
00302 static void close_vaspchgcar_read(void *mydata)
00303 {
00304   vasp_plugindata_t *data = (vasp_plugindata_t *)mydata;
00305   vasp_plugindata_free(data);
00306 }
00307 
00308 
00309 /* registration stuff */
00310 static molfile_plugin_t plugin;
00311 
00312 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 = "CHGCAR";
00317   plugin.prettyname = "VASP_CHGCAR";
00318   plugin.author = "Sung Sakong";
00319   plugin.majorv = 0;
00320   plugin.minorv = 7;
00321   plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
00322   plugin.filename_extension = "CHGCAR";
00323   plugin.open_file_read = open_vaspchgcar_read;
00324   plugin.close_file_read = close_vaspchgcar_read;
00325   plugin.read_volumetric_metadata = read_vaspchgcar_metadata;
00326   plugin.read_volumetric_data = read_vaspchgcar_data;
00327   return VMDPLUGIN_SUCCESS;
00328 }
00329 
00330 int VMDPLUGIN_fini(void) {
00331   return VMDPLUGIN_SUCCESS;
00332 }
00333 
00334 int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
00335   (*cb)(v, (vmdplugin_t *)&plugin);
00336   return VMDPLUGIN_SUCCESS;
00337 }

Generated on Thu Mar 28 03:08:20 2024 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002