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

abinitplugin.c

Go to the documentation of this file.
00001 /*
00002  *  ABINIT plugin for VMD
00003  *  Rob Lahaye, Sungkyunkwan University, Korea
00004  *  August 2010
00005  *
00006  *  LICENSE:
00007  *    You can include my code for as long as it is part of an 
00008  *  open source project. Hence: I do not allow my code to be 
00009  *  part of any closed source project.
00010  *
00011  *  ABINIT manual   
00012  *  http://www.abinit.org/
00013  * 
00014  *  LINUX
00015  *  gcc -O2 -Wall -fPIC -I. -I$VMDBASEDIR/plugins/include -c abinitplugin.c
00016  *  gcc -shared -o abinitplugin.so abinitplugin.o
00017  *
00018  *  MACOSX
00019  *  c++ -O2 -Wall -I. -I$VMDBASEDIR/plugins/include -c abinitplugin.c
00020  *  c++ -bundle -o abinitplugin.so abinitplugin.o
00021  *
00022  *  Install
00023  *  copy abinitplugin.so $VMDBASEDIR/plugins/$ARCH/molfile
00024  */
00025 
00026  /*
00027   * This plugin does NOT read the general input (*.in) and output (*.out) files of abinit;
00028   * their syntax is far too flexible and complex to be read in all its varieties.
00029   * Most (soon hopefully: all) other output files can be read by this plugin. 
00030   *
00031   * The VMD generic plugin routines jump to the appropriate routine, depending on the
00032   * selected file type:
00033   *
00034   * open_file_read / read_structure / read_next_timestep / 
00035   *            |__ GEO_open_file_read/read_structure/read_next_timestep
00036   *            |__ DEN_POT_WFK_open_file_read/read_structure/read_next_timestep
00037   *
00038   * read_volumetric_metadata / read_volumetric_data
00039   *            |__ DEN_read_volumetric_metadata/read_volumetric_data
00040   *            |__ POT_read_volumetric_metadata/read_volumetric_data
00041   *            |__ WFK_read_volumetric_metadata/read_volumetric_data (NOT YET IMPLEMENTED)
00042   *
00043   * open_file_write / write_structure / write_timestep / close_file_write
00044   *            |
00045   *            |__ the output format is in a sloppy abinit input file format; however,
00046   *                this file cannot serve as an input file to abinit, although its
00047   *                syntax is supposed to be easy to understand for a human :).
00048   */
00049 
00050 #include <stdio.h>
00051 #include <stdlib.h>
00052 #if defined(_MSC_VER)
00053 #include <io.h>
00054 #define F_OK 0
00055 // Visual C++ 2005 incorrectly displays a warning about the use of POSIX APIs
00056 // on Windows, which is supposed to be POSIX compliant...
00057 #define access _access
00058 #else
00059 #include <unistd.h>
00060 #endif
00061 #include <math.h>
00062 #include <string.h>
00063 
00064 #include "molfile_plugin.h"
00065 #include "periodic_table.h"
00066 #include "unit_conversion.h"
00067 
00068 #define LINESIZE 2048  /* maximum length of a line */
00069 #define NATOM_MAX 300  /* maximum number of atoms */
00070 
00071 #define DBGPRINT if(1) fprintf
00072 
00073 #ifndef M_PI
00074 #define M_PI 3.14159265358979323846
00075 #endif
00076 
00077 /* Structure to store a file's binary type */
00078 enum Endianness { little_endian, big_endian };
00079 typedef struct {
00080   enum Endianness endian; /* little or big endian*/
00081   int recordmarker;       /* 4 or 8 bytes */
00082 } binary_t;
00083 
00084 
00085 /* Structure to store the header data of an abinit binary file
00086  * See for example
00087  *   http://www.abinit.org/documentation/helpfiles/for-v6.0/users/abinit_help.html#6
00088  */
00089 typedef struct {
00090   binary_t bintype;
00091   char codvsn[7];
00092   int headform,fform;
00093   int bantot, date, intxc, ixc, natom, ngfft[3], nkpt, npsp,
00094       nspden, nspinor, nsppol, nsym, ntypat, occopt, pertcase, usepaw;
00095   double ecut, ecutdg, ecutsm, ecut_eff, qptn[3], rprimd[3][3], stmbias, tphysel, tsmear;
00096   int usewvl, *istwfk, *nband, *npwarr, *so_psp, *symafm, *symrel[3][3], *typat;
00097   double *kpt[3], *occ, *tnons[3], *znucltypat, *wtk;
00098 
00099   char title[133];
00100   double znuclpsp, zionpsp;
00101   int pspso, pspdat, pspcod, pspxc, lmn_size;
00102 
00103   double residm, *xred[3], etotal, fermie;
00104 
00105   int cplex;
00106 } abinit_binary_header_t;
00107 
00108 
00109 /* Structure to store important abinit data so that it can be easily shared among the routines */
00110 typedef struct {
00111   FILE *file;                  /* file pointer for reading or writing */
00112   char *filename;              /* file name for reading or writing */
00113   char filetype[4];            /* type of data file: e.g. GEOmetric or DENsity file */
00114 
00115   float rotmat[3][3];          /* rotation matrix, stored for periodic display hack */
00116 
00117   /* variables for atomic structure */
00118   float rprimd[3][3];          /* Real space PRIMitive translations of the unit cell */
00119   int natom;                   /* total number of atoms */
00120   int typat[NATOM_MAX];        /* number of atoms per atom type */
00121   molfile_atom_t *atomlist;
00122 
00123   /* volumetric variables for data of charge density, wavefunction, etc. */
00124   int nvolsets;                /* number of volumetric datasets */
00125   molfile_volumetric_t *vol;   /* volume set metadata */
00126 
00127   abinit_binary_header_t *hdr; /* header info of an abinit binary file */
00128 } abinit_plugindata_t;
00129 
00130 
00131 static int binread(void *, size_t, FILE *, binary_t);
00132 static abinit_binary_header_t *abinit_header(FILE *);
00133 
00134 
00135 /* Allocate memory for header */
00136 static abinit_binary_header_t *abinit_header_malloc()
00137 {
00138   abinit_binary_header_t *hdr = (abinit_binary_header_t *)malloc(sizeof(abinit_binary_header_t));
00139 
00140   /* zero all bytes in the structure */
00141   if (hdr) memset(hdr, 0, sizeof(abinit_binary_header_t));
00142   else fprintf(stderr, "\n\nABINIT plugin) ERROR: cannot allocate memory for header.\n");
00143 
00144   return hdr;
00145 }
00146 
00147 
00148 /* Allocate memory for plugin data */
00149 static abinit_plugindata_t *abinit_plugindata_malloc()
00150 {
00151   abinit_plugindata_t *data = (abinit_plugindata_t *)malloc(sizeof(abinit_plugindata_t));
00152 
00153   /* zero all bytes in the structure */
00154   if (data) memset(data, 0, sizeof(abinit_plugindata_t));
00155   else fprintf(stderr, "\n\nABINIT plugin) ERROR: cannot allocate memory for plugin data.\n");
00156 
00157   return data;
00158 }
00159 
00160 
00161 /* Free up the header data */
00162 static void abinit_header_free(abinit_binary_header_t *hdr)
00163 {
00164   int i;
00165 
00166   if (!hdr) return;
00167 
00168   if (hdr->istwfk) free(hdr->istwfk);
00169   if (hdr->nband) free(hdr->nband);
00170   if (hdr->npwarr) free(hdr->npwarr);
00171   if (hdr->so_psp) free(hdr->so_psp);
00172   if (hdr->symafm) free(hdr->symafm);
00173   for (i = 0; i < 3; ++i) {
00174     int j;
00175     for (j = 0; j < 3; ++j) if (hdr->symrel[i][j]) free(hdr->symrel[i][j]);
00176     if (hdr->kpt[i]) free(hdr->kpt[i]);
00177     if (hdr->tnons[i]) free(hdr->tnons[i]);
00178     if (hdr->xred[i]) free(hdr->xred[i]);
00179   }
00180   if (hdr->typat) free(hdr->typat);
00181   if (hdr->occ) free(hdr->occ);
00182   if (hdr->znucltypat) free(hdr->znucltypat);
00183   if (hdr->wtk) free(hdr->wtk);
00184 
00185   free(hdr);
00186   hdr = NULL;
00187 }
00188 
00189 
00190 /* Free up the plugin data */
00191 static void abinit_plugindata_free(abinit_plugindata_t *data)
00192 {
00193   if (!data) return;
00194 
00195   if (data->file) fclose(data->file);
00196   if (data->filename) free(data->filename);
00197   if (data->atomlist) free(data->atomlist);
00198   if (data->vol) free(data->vol);
00199 
00200   abinit_header_free(data->hdr);
00201 
00202   free(data);
00203   data = NULL;
00204 }
00205 
00206 
00207 /* This rotation matrix is determined by two vectors.
00208  * The matrix aligns the first vector along the x-axis
00209  * and places the second vector in the xy-plane.
00210  *
00211  * Input: data-structure, which contains the matrix
00212  */
00213 static void abinit_buildrotmat(abinit_plugindata_t *data)
00214 {
00215   float const *const a = data->rprimd[0];
00216   float const *const b = data->rprimd[1];
00217 
00218   /* Rotate first about y-axis and z-axis to align vector a along the x-axis.
00219    * phi  : angle between vector a and its projection on the xy-plane
00220    * theta: angle between projection of vector a on the xy-plane and the x-axis
00221    */
00222   const double len   = sqrt(a[0]*a[0] + a[1]*a[1]);
00223   const double phi   = atan2((double) a[2], (double) len);
00224   const double theta = atan2((double) a[1], (double) a[0]);
00225 
00226   const double cph = cos(phi);
00227   const double cth = cos(theta);
00228   const double sph = sin(phi);
00229   const double sth = sin(theta);
00230 
00231   /* Rotate about x-axis to place b in the xy-plane. */
00232   const double psi = atan2(-sph*cth*b[0] - sph*sth*b[1] + cph*b[2],-sth*b[0] + cth*b[1]);
00233   const double cps = cos(psi);
00234   const double sps = sin(psi);
00235 
00236   data->rotmat[0][0] =  cph * cth;
00237   data->rotmat[0][1] =  cph * sth;
00238   data->rotmat[0][2] =  sph;
00239   data->rotmat[1][0] = -sth * cps - sph * cth * sps;
00240   data->rotmat[1][1] =  cth * cps - sph * sth * sps;
00241   data->rotmat[1][2] =  cph * sps; 
00242   data->rotmat[2][0] =  sth * sps - sph * cth * cps;
00243   data->rotmat[2][1] = -cth * sps - sph * sth * cps; 
00244   data->rotmat[2][2] =  cph * cps;
00245 
00246   DBGPRINT(stderr, "   ROTATION MATRIX: %f   %f   %f\n", data->rotmat[0][0], data->rotmat[0][1], data->rotmat[0][2]);
00247   DBGPRINT(stderr, "                    %f   %f   %f\n", data->rotmat[1][0], data->rotmat[1][1], data->rotmat[1][2]);
00248   DBGPRINT(stderr, "                    %f   %f   %f\n", data->rotmat[2][0], data->rotmat[2][1], data->rotmat[2][2]);
00249 }
00250 
00251 
00252 /* Read a non-empty line from stream, remove comments (#... or !...)
00253  * and strip redundant whitespaces.
00254  *
00255  * Input: string, which can contain the line
00256  *        input stream
00257  *
00258  * Return: the stripped line, or NULL when EOF is reached
00259  */
00260 static char *abinit_readline(char *line, FILE *stream)
00261 {
00262   char *lineptr;
00263 
00264   if (!line || !stream) return NULL;
00265 
00266   do {
00267     int i;
00268     char *cptr;
00269 
00270     /* read one line from the stream */
00271     lineptr = fgets(line, LINESIZE, stream);
00272 
00273     /* first remove comment from the line */
00274     for (i = 0; i < strlen(line); ++i) {
00275         if (line[i] == '#' || line[i] == '!') {line[i] = '\0'; break;}
00276     }
00277 
00278     /* next remove redundant white spaces at the end of the line */
00279     for (cptr = &line[strlen(line) - 1]; isspace(*cptr); --cptr) *cptr = '\0';
00280 
00281     /* continue for as long as EOF is not reached and the line is empty */
00282   } while (lineptr != NULL && strlen(line) == 0);
00283 
00284   return lineptr;
00285 }
00286 
00287 
00288 /* Abinit uses and generates several types of files:
00289  *   foobar_GEO: geometry file (formatted/ascii text)
00290  *   foobar_DEN: charge density file (unformatted/binary)
00291  *   foobar_WFK: wavefunction file (unformatted/binary)
00292  *   foobar_POT: potential file (unformatted/binary)
00293  *
00294  * If the filetype is already set, then we only compare the
00295  * filetype with the given string, otherwise we first find
00296  * out about the filetype and then do the compare.
00297  *
00298  * Input: abinit-data-structure (may or may not have the filetype set)
00299           string to compare with
00300  *
00301  * Return: comparison result: 1 (equal), 0 (not equal or error)
00302  */
00303 static int abinit_filetype(abinit_plugindata_t *data, char const *cmp)
00304 {
00305   char lineptr[LINESIZE];
00306 
00307   if (!data || !cmp) return 0;
00308 
00309   /* if filetype is already set, then only compare */
00310   if (strlen(data->filetype) != 0) return (strncmp(data->filetype, cmp, 3) == 0);
00311 
00312   /* first try to read the abinit binary header */
00313   data->hdr = abinit_header(data->file);
00314   if (data->hdr) {
00315     /* header is read successfully,
00316      * therefore it must be an abinit unformatted (binary) file
00317      */
00318 
00319     switch(data->hdr->fform) {
00320     case   2: /* WFK Wave Function file */
00321               strcpy(data->filetype, "WFK");
00322               break;
00323     case  52: /* DEN Charge Density file */
00324               strcpy(data->filetype, "DEN");
00325               break;
00326     case 102: /* POT Potential file */
00327               strcpy(data->filetype, "POT");
00328               break;
00329     default:  /* Error */
00330               strcpy(data->filetype, "ERR");
00331               break;
00332     }
00333 
00334   } else {
00335     /* reading the abinit binary header failed;
00336      * we now resort to formatted (ascii text) file.
00337      */
00338 
00339     /* read the first non-empty line of the file */
00340     rewind(data->file);
00341     abinit_readline(lineptr, data->file);
00342 
00343     if (strstr(lineptr, " GEO file"))
00344       strcpy(data->filetype, "GEO"); /* GEO geometry file */
00345     else
00346       strcpy(data->filetype, "ERR"); /* Error */
00347 
00348     /* set the file pointer back to the beginning of the file */
00349     rewind(data->file);
00350   }
00351 
00352   return (strncmp(data->filetype, cmp, 3) == 0);
00353 }
00354 
00355 
00356 /* Construct a new file name by finding the first integer number
00357  * in the existing file name, and create the same file name but
00358  * with this integer number incremented.
00359  *
00360  * Note: the string 'filename' is changed into the new filename
00361  * ONLY if the file with the new filename exists! 
00362  *
00363  * Input: filename
00364  *
00365  * Return: 0 (success), 1 (no more files), 2 (error)
00366  */
00367 static int increment_filename(char *filename)
00368 {
00369   int i;
00370   char *newfilename = NULL, *endpart = NULL;
00371 
00372   DBGPRINT(stderr, "Enter increment_filename\n");
00373 
00374   DBGPRINT(stderr, "increment_filename: filename = %s \n", filename);
00375   /* search for integer in the filename starting from the end of the string */
00376   for (i = strlen(filename) - 1; i >= 0 && !newfilename; --i) {
00377 
00378     /* endpart points to the string AFTER the integer */
00379     if (!endpart && isdigit(filename[i])) endpart = strdup(filename + i + 1);
00380 
00381     /* allocate newfilename when integer is found */
00382     if (endpart && !newfilename && !isdigit(filename[i])) {
00383        newfilename = (char *)malloc(sizeof(char) * (2 + strlen(filename)));
00384        if (!newfilename) {
00385          free(endpart);
00386          return 2;
00387        }
00388 
00389        /* first copy part of the string BEFORE the integer number */
00390        strncpy(newfilename, filename, i + 1);
00391 
00392        /* second append the incremented integer number and add 'endpart' of the filename */
00393        sprintf(newfilename + i + 1, "%d%s", 1 + atoi(filename + i + 1), endpart);
00394     }
00395   }
00396 
00397   /* if the endpart is not found, then the file name does not have a number in it,
00398    * and therefore this filename is not part of a series of timesteps.
00399    * note that if this occurs, then 'newfilename' must also still be NULL !
00400    */
00401   if (!endpart) {
00402     DBGPRINT(stderr, "Exit increment_filename\n");
00403     return 1;
00404   }
00405   
00406   /* clean up */
00407   free(endpart);
00408 
00409   /* check if we can access the new file */
00410   if (access(newfilename, F_OK) != 0) {
00411     /* the new file cannot be accessed, so it does not exist */
00412     free(newfilename);
00413     DBGPRINT(stderr, "Exit increment_filename\n");
00414     return 1;
00415   } else {
00416     /* the new filename exists! Replace the old file name. */
00417     strcpy(filename, newfilename);
00418     free(newfilename);
00419     DBGPRINT(stderr, "increment_filename: filename = %s \n", filename);
00420     DBGPRINT(stderr, "Exit increment_filename\n");
00421     return 0;
00422   }
00423 }
00424 
00425 
00426 /* Geometry files are generated by ABINIT when "prtgeo 1" is given in
00427  * the input file. Note that this is NOT the abinit default!
00428  *
00429  * Input: data-structure
00430  *        reference to the number of atoms
00431  *
00432  * Return: data-structure, or NULL on error
00433  */
00434 static void *GEO_open_file_read(abinit_plugindata_t *data, int *natoms)
00435 {
00436   char lineptr[LINESIZE], atomname[NATOM_MAX][10];
00437   int i, idx;
00438 
00439   DBGPRINT(stderr, "Enter GEO_open_file_read\n");
00440 
00441   /* go to the line with the text 'XMOL data' */
00442   while (abinit_readline(lineptr, data->file) != NULL) {
00443     if (strstr(lineptr, "XMOL data")) break;
00444   }
00445   if (!strstr(lineptr, "XMOL data")) {
00446     fprintf(stderr, "\n\nABINIT read) ERROR: '%s' has no 'XMOL data...' lines.\n", data->filename);
00447     return NULL;
00448   }
00449 
00450   /* next non-empty line has the number of atoms */
00451   if (abinit_readline(lineptr, data->file) == NULL) {
00452     fprintf(stderr, "\n\nABINIT read) ERROR: cannot find the number of atoms in file '%s'.\n", data->filename);
00453     return NULL;
00454   }
00455   data->natom = atoi(lineptr);
00456   if (data->natom <= 0 || data->natom > NATOM_MAX) {
00457     fprintf(stderr, "\n\nABINIT read) ERROR: file '%s' has %d number of atoms.\n", data->filename, data->natom);
00458     return NULL;
00459   }
00460 
00461   /* read through the list of atom names and ignore their positions */
00462   for (i = 0; i < NATOM_MAX; ++i) data->typat[i] = atomname[i][0] = '\0';
00463   for (idx = i = 0; i < data->natom; ++i) {
00464     int n;
00465     char name[10];
00466     if (1 != fscanf(data->file, "%s %*f %*f %*f", name)) {
00467       fprintf(stderr, "\n\nABINIT read) ERROR: file '%s' does not have the atom list.\n", data->filename);
00468       return NULL;
00469     }
00470 
00471     /* compare current atom name with previous read list and set typat accordingly */
00472     for (n = 0; n < idx; ++n) if (strcmp(atomname[n], name) == 0) break;
00473     if (n == idx) strcpy(atomname[idx++], name);
00474     data->typat[i] = n + 1;
00475 
00476     DBGPRINT(stderr, "   \"%s\": name = %s : data->typat[%d] = %d\n", data->filetype, atomname[n], i, data->typat[i]);
00477   }
00478 
00479   rewind(data->file);
00480 
00481   *natoms = data->natom;
00482 
00483   DBGPRINT(stderr, "Exit GEO_open_file_read\n");
00484   return data;
00485 }
00486 
00487 
00488 static int GEO_read_structure(abinit_plugindata_t *data, int *optflags, molfile_atom_t *atomlist)
00489 {
00490   char lineptr[LINESIZE];
00491   int i, status;
00492 
00493   DBGPRINT(stderr, "Enter GEO_read_structure\n");
00494 
00495   /* go to the line with the 'XMOL data' */
00496   do {
00497     char *line = abinit_readline(lineptr, data->file);
00498     status = line != NULL && !strstr(lineptr, "XMOL data");
00499   } while (status);
00500 
00501   /* skip line with atom numbers */
00502   abinit_readline(lineptr, data->file);
00503 
00504   /* find atom types in XMOL list */
00505   for (i = 0; i < data->natom; ++i) {
00506     molfile_atom_t *const atom = &(atomlist[i]);
00507 
00508     /* required fields */
00509     if (1 != fscanf(data->file, "%s %*f %*f %*f", atom->name)) {
00510       fprintf(stderr, "\n\nABINIT read) ERROR: file '%s' does not have the atom list.\n", data->filename);
00511       return MOLFILE_ERROR;
00512     }
00513     strncpy(atom->type, atom->name, sizeof(atom->type));
00514     atom->resname[0] = '\0';
00515     atom->resid = 1;
00516     atom->segid[0]='\0';
00517     atom->chain[0]='\0';
00518 
00519     /* Optional fields (defined in *optflags) */
00520     atom->atomicnumber = get_pte_idx(atom->name);
00521     atom->mass = get_pte_mass(atom->atomicnumber);
00522     atom->radius = get_pte_vdw_radius(atom->atomicnumber);
00523 
00524     DBGPRINT(stderr, "   atom %d : %d (%s)\n", i, atom->atomicnumber, atom->name);
00525   }
00526 
00527   /* tell which of the optional fields in the molfile_atom_t structure are provided */
00528   *optflags = MOLFILE_ATOMICNUMBER | MOLFILE_MASS | MOLFILE_RADIUS; 
00529 
00530   rewind(data->file);
00531 
00532   DBGPRINT(stderr, "Exit GEO_read_structure\n");
00533   return MOLFILE_SUCCESS;
00534 }
00535 
00536 
00537 static int GEO_read_next_timestep(abinit_plugindata_t *data, int natoms, molfile_timestep_t *ts)
00538 {
00539   char lineptr[LINESIZE];
00540   float *a, *b, *c;
00541   int i , status;
00542 
00543   DBGPRINT(stderr, "Enter GEO_read_next_timestep\n");
00544 
00545   /* At the very first call the file pointer will be non-NULL.
00546    * All consecutive calls will have a NULL file pointer (see at the end
00547    * of this routine) and in that case "increment" the filename for the
00548    * data of the next timestep; if this file does not exist, then there
00549    * are no more timesteps.
00550    */
00551   if (!data->file) {
00552     if (increment_filename(data->filename) != 0) return MOLFILE_EOF;
00553 
00554     data->file = fopen(data->filename, "r");
00555     if (!data->file) return MOLFILE_EOF;
00556   }
00557 
00558   DBGPRINT(stderr, "GEO_read_next_timestep: filename = %s \n", data->filename);
00559 
00560   /* go to the line with 'Primitive vectors...' */
00561   do {
00562     char *line = abinit_readline(lineptr, data->file);
00563     status = ( line != NULL && !strstr(lineptr, "Primitive vectors") );
00564   } while (status);
00565 
00566   /* read unit cell vectors from file */
00567   for (i = 0; i < 3; ++i) {
00568     float length, *r = data->rprimd[i];
00569     if (3 != fscanf(data->file, "%*s %f %f %f", &r[0], &r[1], &r[2])) return MOLFILE_EOF;
00570 
00571     /* convert length units from Bohr to Angstrom */
00572     r[0] *= BOHR_TO_ANGS;
00573     r[1] *= BOHR_TO_ANGS;
00574     r[2] *= BOHR_TO_ANGS;
00575 
00576     /* lengths of the respective unit cell vector */
00577     length = sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]);
00578     switch (i) {
00579       case 0: ts->A = length; break;
00580       case 1: ts->B = length; break;
00581       case 2: ts->C = length; break;
00582     }
00583   }
00584   abinit_buildrotmat(data);
00585 
00586   /* determine angles between the vectors of the unit cell */
00587   a = data->rprimd[0];
00588   b = data->rprimd[1];
00589   c = data->rprimd[2];
00590 
00591   /* alpha: angle (in degrees) between b and c */
00592   ts->alpha = (180.0/M_PI) * acos( (b[0]*c[0] + b[1]*c[1] + b[2]*c[2]) / (ts->B*ts->C) );
00593 
00594   /* beta: angle (in degrees) between a and c */
00595   ts->beta  = (180.0/M_PI) * acos( (a[0]*c[0] + a[1]*c[1] + a[2]*c[2]) / (ts->A*ts->C) );
00596 
00597   /* gamma: angle (in degrees) between a and b */
00598   ts->gamma = (180.0/M_PI) * acos( (a[0]*b[0] + a[1]*b[1] + a[2]*b[2]) / (ts->A*ts->B) );
00599 
00600   for (i = 0; i < 9; ++i) DBGPRINT(stderr, "   data->rprimd[%d][%d] = %f %s", i%3, i/3, data->rprimd[i%3][i/3], ((i+1)%3 == 0 ? "\n" : ""));
00601 
00602   /* go to the line with 'XMOL data' */
00603   do {
00604     char *line = abinit_readline(lineptr, data->file);
00605     status = line != NULL && !strstr(lineptr, "XMOL data");
00606   } while (status);
00607 
00608   /* skip line with atom numbers */
00609   abinit_readline(lineptr, data->file);
00610 
00611   /* read coordinates (in Angstrom) of each atom */
00612   for (i = 0; i < data->natom; ++i) {
00613     float x, y, z, *coords = &ts->coords[3*i];
00614     fscanf(data->file, "%*s %f %f %f", &x, &y, &z);
00615     coords[0] = data->rotmat[0][0]*x + data->rotmat[0][1]*y + data->rotmat[0][2]*z;
00616     coords[1] = data->rotmat[1][0]*x + data->rotmat[1][1]*y + data->rotmat[1][2]*z;
00617     coords[2] = data->rotmat[2][0]*x + data->rotmat[2][1]*y + data->rotmat[2][2]*z;
00618   }
00619 
00620   /* close the file and NULLify as preparation for the next timestep call */
00621   fclose(data->file);
00622   data->file = NULL;
00623 
00624   DBGPRINT(stderr, "Exit GEO_read_next_timestep\n");
00625   return MOLFILE_SUCCESS;
00626 }
00627 
00628 
00629 /* Electron density (DEN), potential (POT), and wavefunction (WFK)
00630  * files have the structure information in the generic header,
00631  * followed by the volumetric data. For all three we can determine
00632  * the structure with the same reading routine.
00633  * The volumetric data is the similar for density and potential (DEN/POT)
00634  * files, but is different for wavefunction (WFK) files.
00635  */
00636 static void *DEN_POT_WFK_open_file_read(abinit_plugindata_t *data, int *natoms)
00637 {
00638   int i;
00639 
00640   DBGPRINT(stderr, "Enter DEN_POT_WFK_open_file_read\n");
00641 
00642   data->natom = data->hdr->natom;
00643 
00644   if (data->natom <= 0 || data->natom > NATOM_MAX) return NULL;
00645 
00646   for (i = 0; i < data->natom; ++i) data->typat[i] = data->hdr->typat[i];
00647 
00648   for (i = 0; i < data->natom; ++i) DBGPRINT(stderr, "   \"%s\": data->typat[%d] = %d\n", data->filetype, i, data->typat[i]);
00649 
00650   *natoms = data->natom;
00651 
00652   DBGPRINT(stderr, "Exit DEN_POT_WFK_open_file_read\n");
00653   return data;
00654 }
00655 
00656 
00657 static int DEN_POT_WFK_read_structure(abinit_plugindata_t *data, int *optflags, molfile_atom_t *atomlist)
00658 {
00659   int i;
00660 
00661   DBGPRINT(stderr, "Enter DEN_POT_WFK_read_structure\n");
00662 
00663   /* set the atom types, names, etc. */
00664   for (i = 0; i < data->natom; ++i) {
00665     molfile_atom_t *const atom = &(atomlist[i]);
00666 
00667     /* optional fields (defined in *optflags) */
00668     atom->atomicnumber = (int)floor(0.5 + data->hdr->znucltypat[data->hdr->typat[i] - 1]);
00669     atom->mass = get_pte_mass(atom->atomicnumber);
00670     atom->radius = get_pte_vdw_radius(atom->atomicnumber);
00671 
00672     /* required fields */
00673     strncpy(atom->name, get_pte_label(atom->atomicnumber), sizeof(atom->name));    
00674     strncpy(atom->type, atom->name, sizeof(atom->type));
00675     atom->resname[0] = '\0';
00676     atom->resid = 1;
00677     atom->segid[0]='\0';
00678     atom->chain[0]='\0';
00679 
00680     DBGPRINT(stderr, "   atom %d : %d (%s)\n", i, atom->atomicnumber, atom->name);
00681   }
00682 
00683   /* tell which of the optional fields in the molfile_atom_t structure are provided */
00684   *optflags = MOLFILE_ATOMICNUMBER | MOLFILE_MASS | MOLFILE_RADIUS; 
00685 
00686   DBGPRINT(stderr, "Exit DEN_POT_WFK_read_structure\n");
00687   return MOLFILE_SUCCESS;
00688 }
00689 
00690 
00691 static int DEN_POT_WFK_read_next_timestep(abinit_plugindata_t *data, int natoms, molfile_timestep_t *ts)
00692 {
00693   float *a, *b, *c;
00694   int i;
00695 
00696   DBGPRINT(stderr, "Enter DEN_POT_WFK_read_next_timestep\n");
00697 
00698   /* At the very first call the file pointer will be non-NULL.
00699    * All consecutive calls will have a NULL file pointer (see at the end
00700    * of this routine) and in that case "increment" the filename for the
00701    * data of the next timestep; if this file does not exist, then there
00702    * are no more timesteps.
00703    */
00704   if (!data->file) {
00705     /* alas, volumetric data cannot (yet ?) be read as timesteps */
00706     return MOLFILE_EOF;
00707 
00708     if (increment_filename(data->filename) != 0) return MOLFILE_EOF;
00709 
00710     data->file = fopen(data->filename, "r");
00711     if (!data->file) return MOLFILE_EOF;
00712 
00713     /* read the new header info from the new file */
00714     abinit_header_free(data->hdr);
00715     data->hdr = abinit_header(data->file);
00716     if (!data->hdr) return MOLFILE_EOF;
00717   }
00718 
00719   /* get unit cell vectors and convert them to Angstrom */
00720   for (i = 0; i < 3; ++i) {
00721     float length;
00722     int k;
00723     for (k = 0; k < 3; ++k) data->rprimd[i][k] = data->hdr->rprimd[i][k] * BOHR_TO_ANGS;
00724     length = sqrt(pow(data->rprimd[i][0], 2) + pow(data->rprimd[i][1], 2) + pow(data->rprimd[i][2], 2));
00725     switch (i) {
00726       case 0: ts->A = length; break;
00727       case 1: ts->B = length; break;
00728       case 2: ts->C = length; break;
00729     }    
00730   }
00731   abinit_buildrotmat(data);
00732 
00733   for (i = 0; i < 9; ++i) DBGPRINT(stderr, "   data->rprimd[%d][%d] = %f %s", i%3, i/3, data->rprimd[i%3][i/3], ((i+1)%3 == 0 ? "\n" : ""));
00734 
00735   /* determine angles between the vectors of the unit cell */
00736   a = data->rprimd[0];
00737   b = data->rprimd[1];
00738   c = data->rprimd[2];
00739 
00740   /* alpha: angle (in degrees) between b and c */
00741   ts->alpha = (180.0/M_PI) * acos( (b[0]*c[0] + b[1]*c[1] + b[2]*c[2]) / (ts->B*ts->C) );
00742 
00743   /* beta: angle (in degrees) between a and c */
00744   ts->beta  = (180.0/M_PI) * acos( (a[0]*c[0] + a[1]*c[1] + a[2]*c[2]) / (ts->A*ts->C) );
00745 
00746   /* gamma: angle (in degrees) between a and b */
00747   ts->gamma = (180.0/M_PI) * acos( (a[0]*b[0] + a[1]*b[1] + a[2]*b[2]) / (ts->A*ts->B) );
00748 
00749   /* get coordinates (in Angstrom) of each atom */
00750   for (i = 0; i < data->natom; ++i) {
00751     double **xred = data->hdr->xred;
00752     float *coords = &ts->coords[3*i];
00753     float const x = xred[0][i] * a[0] + xred[1][i] * b[0] + xred[2][i] * c[0],
00754                 y = xred[0][i] * a[1] + xred[1][i] * b[1] + xred[2][i] * c[1],
00755                 z = xred[0][i] * a[2] + xred[1][i] * b[2] + xred[2][i] * c[2];
00756 
00757     coords[0] = data->rotmat[0][0]*x + data->rotmat[0][1]*y + data->rotmat[0][2]*z;
00758     coords[1] = data->rotmat[1][0]*x + data->rotmat[1][1]*y + data->rotmat[1][2]*z;
00759     coords[2] = data->rotmat[2][0]*x + data->rotmat[2][1]*y + data->rotmat[2][2]*z;
00760   }
00761 
00762   /* close the file and NULLify as preparation for the next timestep call */
00763   fclose(data->file);
00764   data->file = NULL;
00765 
00766   DBGPRINT(stderr, "Exit DEN_POT_WFK_read_next_timestep\n");
00767   return MOLFILE_SUCCESS;
00768 }
00769 
00770 
00771 static int DEN_read_volumetric_metadata(abinit_plugindata_t *data, int *nvolsets, molfile_volumetric_t **metadata)
00772 {
00773   char const spintext[3][30] = { "Charge density spin up",
00774                                  "Charge density spin down",
00775                                  "Charge density spin up - down" };
00776   char const magntext[3][40] = { "X-projection of local magnetization",
00777                                  "Y-projection of local magnetization",
00778                                  "Z-projection of local magnetization" };
00779   int i;
00780 
00781   DBGPRINT(stderr, "Enter DEN_read_volumetric_metadata\n");
00782 
00783   /* Abinit provides the electron density, which is the negative of the
00784    * charge density. For norm-conserving pseudopotentials (usepaw = 0)
00785    * this is the pseudo-valence electron density; for a PAW calculation
00786    * (usepaw = 1) this is the pseudo-valence electron density plus the
00787    * compensation charge density.
00788    * To visualize the electron density from a PAW calculation, create
00789    * the density file with the pawprtden key word, not prtden (however,
00790    * if you want to chain this density into later calculations, you have
00791    * to use prtden.
00792    */
00793   if (data->hdr->usepaw) {
00794     fprintf(stderr, "\n\nABINIT read) WARNING: be sure that you have used \"pawprtden 1\"\n");
00795     fprintf(stderr,     "                      in order to visualize the electron density!\n\n");
00796   }
00797 
00798   /* Initialize the volume set list
00799    * Total charge density (spin up+down) : always present
00800    *    spin up / spin down / spin up-down \ only there for
00801    *    X / Y / Z local magnetization      /  spin-polarized calculations
00802    *                    (the latter two are absent for non-spin-polarized calculations)
00803    */
00804   data->nvolsets = (data->hdr->nspden == 1 ? 1 : 4);
00805   data->vol = (molfile_volumetric_t *)malloc(data->nvolsets * sizeof(molfile_volumetric_t));
00806   if (!data->vol) {
00807     fprintf(stderr, "\n\nABINIT read) ERROR: cannot allocate space for volumetric data.\n");
00808     return MOLFILE_ERROR;
00809   }
00810 
00811   /* get unit cell vectors and convert to Angstrom */
00812   for (i = 0; i < 3; ++i) {
00813     int k;
00814     for (k = 0; k < 3; ++k) data->rprimd[i][k] = data->hdr->rprimd[i][k] * BOHR_TO_ANGS;
00815   }
00816   abinit_buildrotmat(data);
00817 
00818   for (i = 0; i < 9; ++i) DBGPRINT(stderr, "   data->rprimd[%d][%d] = %f %s", i%3, i/3, data->rprimd[i%3][i/3], ((i+1)%3 == 0 ? "\n" : ""));
00819 
00820   for (i = 0; i < data->nvolsets; ++i) {
00821 
00822     int k;
00823 
00824     /* handle to the current volume set meta data */
00825     molfile_volumetric_t *const set = &(data->vol[i]);
00826 
00827     /* set the volume data name */
00828     if (i == 0) {
00829         /* the first data set always is the total charge density */
00830         sprintf(set->dataname, "Total charge density");
00831     } else if (data->hdr->nspden <= 2) {
00832         /* datasets contain spin polarized (up/down) densities */
00833         sprintf(set->dataname, "%s", spintext[i-1]);
00834     } else if (data->hdr->nspden == 4) {
00835         /* subsequent datasets contain magnetization data */
00836         sprintf(set->dataname, "%s", magntext[i-1]);
00837     } else {
00838         /* we should never get here.... */
00839         sprintf(set->dataname, "%s", "ERROR: no datasets available");
00840     }
00841 
00842     /* rotate unit cell vectors */
00843     for (k = 0 ; k < 3; ++k) {
00844       set->xaxis[k] = data->rotmat[k][0] * data->rprimd[0][0]
00845                     + data->rotmat[k][1] * data->rprimd[0][1]
00846                     + data->rotmat[k][2] * data->rprimd[0][2];
00847 
00848       set->yaxis[k] = data->rotmat[k][0] * data->rprimd[1][0] 
00849                     + data->rotmat[k][1] * data->rprimd[1][1]
00850                     + data->rotmat[k][2] * data->rprimd[1][2];
00851 
00852       set->zaxis[k] = data->rotmat[k][0] * data->rprimd[2][0] 
00853                     + data->rotmat[k][1] * data->rprimd[2][1]
00854                     + data->rotmat[k][2] * data->rprimd[2][2];
00855     }
00856     DBGPRINT(stderr, "   set->xaxis[%d] set->yaxis[%d] set->zaxis[%d]\n", k, k, k);
00857     for (k = 0 ; k < 3; ++k) DBGPRINT(stderr, "   %f         %f        %f\n", set->xaxis[k], set->yaxis[k], set->zaxis[k]);
00858 
00859     /* Add one more to the grid size and later fill the extra voxel
00860      * with the same value from the beginning of the row to make
00861      * the volumetric data smooth across the cell boundaries.
00862      */
00863     set->xsize = data->hdr->ngfft[0] + 1;
00864     set->ysize = data->hdr->ngfft[1] + 1;
00865     set->zsize = data->hdr->ngfft[2] + 1;
00866 
00867     set->has_color = 0;
00868     set->origin[0] = set->origin[1] = set->origin[2] = 0;
00869   }
00870 
00871   *nvolsets = data->nvolsets;
00872   *metadata = data->vol;  
00873 
00874   DBGPRINT(stderr, "Exit DEN_read_volumetric_metadata.\n");
00875   return MOLFILE_SUCCESS;
00876 }
00877 
00878 
00879 static int DEN_read_volumetric_data(abinit_plugindata_t *data, int set, float *datablock, float *colorblock)
00880 {
00881   double const density_conversion = 1.0/pow(BOHR_TO_ANGS, 3);
00882   int iset;
00883 
00884   DBGPRINT(stderr, "Enter DEN_read_volumetric_data\n");
00885 
00886   /* this should never happen, nevertheless check it to prevent disasters */
00887   if (set >= data->nvolsets) return MOLFILE_ERROR;
00888 
00889   /* The density output file contains the header plus:
00890    * do ispden=1,nspden
00891    * write(unit) (rhor(ir),ir=1,cplex*ngfft(1)*ngfft(2)*ngfft(3))
00892    * enddo
00893    */
00894   for (iset = 0; iset <= set && iset < data->hdr->nspden; ++iset) {
00895     int const xsize = data->vol[iset].xsize; 
00896     int const ysize = data->vol[iset].ysize;
00897     int const zsize = data->vol[iset].zsize;
00898 
00899     char recordmarker[10];
00900     int n, ix, iy, iz;
00901 
00902     /* Fill the datablock with the density or magnetization data.
00903      * Note that for each 'iset'-loop the datablock is overwritten,
00904      * so that only the last datablock readings remain.
00905      */
00906     for (n = iz = 0; iz < zsize; ++iz) {
00907       for (iy = 0; iy < ysize; ++iy) {
00908         for (ix = 0; ix < xsize; ++ix, ++n) {
00909           double value;
00910 
00911           /* The datablock grid is one voxel larger in each direction in order
00912            * to smoothen the density at the cell boundaries; hence fill the extra
00913            * voxel with the same value as the one at the beginning of that row.
00914            */
00915           if (ix == xsize - 1) value = datablock[n - ix];
00916           else if (iy == ysize - 1) value = datablock[n - iy*xsize];
00917           else if (iz == zsize - 1) value = datablock[n - iz*ysize*xsize];
00918           else if (data->hdr->cplex == 1) {
00919             /* get the volumetric data and convert */
00920             binread(&value, 8, data->file, data->hdr->bintype);
00921             value *= density_conversion;
00922           } else if (data->hdr->cplex == 2) {
00923             /* get volumetric data as a complex number and convert */
00924             double a, b;
00925             binread(&a, 8, data->file, data->hdr->bintype);
00926             binread(&b, 8, data->file, data->hdr->bintype);
00927             value = sqrt(a*a + b*b) * density_conversion;
00928           } else {
00929             /* we should never get here */
00930             return MOLFILE_ERROR;
00931           }
00932 
00933           /* fill the datablock according to the data provided */
00934           if (data->hdr->nspden <= 2) {
00935             /* data is charge density as one or two sets: total and spin up */
00936 
00937             if (set == 0 || set == 1) {
00938               /* first two sets are always directly provided as total and spin up */
00939               datablock[n] = value;
00940             } else if (set == 2) {
00941               /* third set we calculate from first two sets: down = total - up */
00942               datablock[n] = (iset == 0 ? value : datablock[n] - value);
00943             } else if (set == 3) {
00944               /* fourth set we calculate from first two sets: up - down = 2*up - total */
00945               datablock[n] = (iset == 0 ? -value : datablock[n] + 2*value);
00946             } else {
00947               /* we should never get here */
00948               return MOLFILE_ERROR;
00949             }
00950 
00951           } else if (data->hdr->nspden == 4) {
00952             /* data is magnetization as four sets: total charge density and X-Y-Z-magnetization */
00953             datablock[n] = value;
00954 
00955           } else {
00956             /* we should never get here...*/
00957             return MOLFILE_ERROR;
00958           }
00959 
00960         }
00961       }
00962     }
00963 
00964     /* skip the recordmarker bytes twice */
00965     fread(recordmarker, 1, data->hdr->bintype.recordmarker, data->file);
00966     fread(recordmarker, 1, data->hdr->bintype.recordmarker, data->file);
00967   }
00968 
00969   DBGPRINT(stderr, "Exit DEN_read_volumetric_data\n");
00970   return MOLFILE_SUCCESS;
00971 }
00972 
00973 
00974 static int POT_read_volumetric_metadata(abinit_plugindata_t *data, int *nvolsets, molfile_volumetric_t **metadata)
00975 {
00976   int i;
00977 
00978   DBGPRINT(stderr, "Enter POT_read_volumetric_metadata\n");
00979 
00980   data->nvolsets = data->hdr->nspden;
00981 
00982   data->vol = (molfile_volumetric_t *)malloc(data->nvolsets * sizeof(molfile_volumetric_t));
00983   if (!data->vol) {
00984     fprintf(stderr, "\n\nABINIT read) ERROR: cannot allocate space for volumetric data.\n");
00985     return MOLFILE_ERROR;
00986   }
00987 
00988   /* get unit cell vectors and convert to Angstrom */
00989   for (i = 0; i < 3; ++i) {
00990     int k;
00991     for (k = 0; k < 3; ++k) data->rprimd[i][k] = data->hdr->rprimd[i][k] * BOHR_TO_ANGS;
00992   }
00993   abinit_buildrotmat(data);
00994 
00995   for (i = 0; i < 9; ++i) DBGPRINT(stderr, "   data->rprimd[%d][%d] = %f %s", i%3, i/3, data->rprimd[i%3][i/3], ((i+1)%3 == 0 ? "\n" : ""));
00996 
00997   for (i = 0; i < data->nvolsets; ++i) {
00998 
00999     int k;
01000 
01001     /* handle to the current volume set meta data */
01002     molfile_volumetric_t *const set = &(data->vol[i]);
01003 
01004     if (data->nvolsets == 1) strcpy(set->dataname, "Total potential");
01005     else if (data->nvolsets == 2) {
01006         if (i == 0) strcpy(set->dataname, "Spin up potential");
01007         if (i == 1) strcpy(set->dataname, "Spin down potential");
01008     } else if (data->nvolsets == 4) {
01009         if (i == 0) strcpy(set->dataname, "Spin up-up potential");
01010         if (i == 1) strcpy(set->dataname, "Spin down-down potential");
01011         if (i == 2) strcpy(set->dataname, "Real part of spin up-down potential");
01012         if (i == 3) strcpy(set->dataname, "Imaginary part of spin up-down potential");
01013     }
01014 
01015     /* rotate unit cell vectors */
01016     for (k = 0 ; k < 3; ++k) {
01017       set->xaxis[k] = data->rotmat[k][0] * data->rprimd[0][0]
01018                     + data->rotmat[k][1] * data->rprimd[0][1]
01019                     + data->rotmat[k][2] * data->rprimd[0][2];
01020 
01021       set->yaxis[k] = data->rotmat[k][0] * data->rprimd[1][0] 
01022                     + data->rotmat[k][1] * data->rprimd[1][1]
01023                     + data->rotmat[k][2] * data->rprimd[1][2];
01024 
01025       set->zaxis[k] = data->rotmat[k][0] * data->rprimd[2][0] 
01026                     + data->rotmat[k][1] * data->rprimd[2][1]
01027                     + data->rotmat[k][2] * data->rprimd[2][2];
01028     }
01029     DBGPRINT(stderr, "   set->xaxis[%d] set->yaxis[%d] set->zaxis[%d]\n", k, k, k);
01030     for (k = 0 ; k < 3; ++k) DBGPRINT(stderr, "   %f         %f        %f\n", set->xaxis[k], set->yaxis[k], set->zaxis[k]);
01031 
01032     /* Add one more to the grid size and fill the extra voxel
01033      * with the same value at the beginning of the row to make
01034      * the volumetric data smooth across the cell boundaries.
01035      */
01036     set->xsize = data->hdr->ngfft[0] + 1;
01037     set->ysize = data->hdr->ngfft[1] + 1;
01038     set->zsize = data->hdr->ngfft[2] + 1;
01039 
01040     set->has_color = 0;
01041     set->origin[0] = set->origin[1] = set->origin[2] = 0;
01042   }
01043 
01044   *nvolsets = data->nvolsets;
01045   *metadata = data->vol;  
01046 
01047   DBGPRINT(stderr, "Exit POT_read_volumetric_metadata.\n");
01048   return MOLFILE_SUCCESS;
01049 }
01050 
01051 
01052 static int POT_read_volumetric_data(abinit_plugindata_t *data, int set, float *datablock, float *colorblock)
01053 {
01054   int n, iset;
01055 
01056   DBGPRINT(stderr, "Enter POT_read_volumetric_data\n");
01057 
01058   /* this should never happen, nevertheless check it to prevent disasters */
01059   if (set >= data->nvolsets) return MOLFILE_ERROR;
01060 
01061   /* The potential files contain the header plus:
01062    * do ispden=1,nspden
01063    * write(unit) (potential(ir),ir=1,cplex*ngfft(1)*ngfft(2)*ngfft(3))
01064    * enddo
01065    */
01066   for (n = iset = 0; iset <= set; ++iset) {
01067     int const xsize = data->vol[iset].xsize; 
01068     int const ysize = data->vol[iset].ysize;
01069     int const zsize = data->vol[iset].zsize;
01070 
01071     char recordmarker[10];
01072     int ix, iy, iz;
01073 
01074     /* Fill the datablock with the density data.
01075      * Note that for each 'iset'-loop the datablock is overwritten,
01076      * so that only the last datablock readings remain.
01077      */
01078     for (n = iz = 0; iz < zsize; ++iz) {
01079       for (iy = 0; iy < ysize; ++iy) {
01080         for (ix = 0; ix < xsize; ++ix, ++n) {
01081           double value;
01082 
01083           /* The datablock grid is one voxel larger in each direction in order
01084            * to smoothen the density at the cell boundaries; hence fill the extra
01085            * voxel with the same value as the one at the beginning of that row.
01086            */
01087           if (ix == xsize - 1) value = datablock[n - ix];
01088           else if (iy == ysize - 1) value = datablock[n - iy*xsize];
01089           else if (iz == zsize - 1) value = datablock[n - iz*ysize*xsize];
01090           else if (data->hdr->cplex == 1) {
01091             /* get the volumetric data and convert */
01092             binread(&value, 8, data->file, data->hdr->bintype);
01093             value *= HARTREE_TO_EV;
01094           } else if (data->hdr->cplex == 2) {
01095             /* get volumetric data as a complex number and convert */
01096             double a, b;
01097             binread(&a, 8, data->file, data->hdr->bintype);
01098             binread(&b, 8, data->file, data->hdr->bintype);
01099             value = sqrt(a*a + b*b) * HARTREE_TO_EV;
01100           } else return MOLFILE_ERROR;
01101 
01102           datablock[n] = value;
01103         }
01104       }
01105     }
01106 
01107     /* THIS HAS NOT BEEN VERIFIED YET: skip the recordmarker bytes twice */
01108     fread(recordmarker, 1, data->hdr->bintype.recordmarker, data->file);
01109     fread(recordmarker, 1, data->hdr->bintype.recordmarker, data->file);
01110   }
01111 
01112   DBGPRINT(stderr, "Exit POT_read_volumetric_data\n");
01113   return MOLFILE_SUCCESS;
01114 }
01115 
01116 
01117 static int WFK_read_volumetric_metadata(abinit_plugindata_t *data, int *nvolsets, molfile_volumetric_t **metadata)
01118 {
01119   /* net yet implemented */
01120   DBGPRINT(stderr, "Enter/Exit WFK_read_volumetric_metadata\n");
01121   fprintf(stderr, "\n\nABINIT read) WARNING: loading WFK is NOT YET IMPLEMENTED!\n");
01122   return MOLFILE_ERROR;
01123 }
01124 
01125 
01126 static int WFK_read_volumetric_data(abinit_plugindata_t *data, int set, float *datablock, float *colorblock)
01127 {
01128   /* net yet implemented */
01129   DBGPRINT(stderr, "Enter/Exit WFK_read_volumetric_data: NOT YET IMPLEMENTED!\n");
01130   fprintf(stderr, "\n\nABINIT read) WARNING: loading WFK is NOT YET IMPLEMENTED!\n");
01131   return MOLFILE_ERROR;
01132 }
01133 
01134 /* ===================================
01135  * Generic vmd routines
01136  * ===================================
01137  */
01138 
01139 
01140 /* VMD calls this one just once to verify access to the file. */
01141 static void *open_file_read(const char *filename, const char *filetype, int *natoms)
01142 {
01143   void *result = NULL;
01144   abinit_plugindata_t *data;
01145 
01146   DBGPRINT(stderr, "Enter open_file_read\n");
01147 
01148   /* verify that the input variables are OK */
01149   if (!filename || !natoms) return NULL;
01150 
01151   /* start with undefined value and set it after successful read */
01152   *natoms = MOLFILE_NUMATOMS_UNKNOWN;
01153 
01154   /* allocate memory for the abinit data structure */
01155   data = abinit_plugindata_malloc();
01156   if (!data) return NULL;
01157 
01158   /* allocate memory for filename (add an extra 10 bytes for some more flexibility) */
01159   data->filename = (char *)malloc( sizeof(char) * (strlen(filename) + 10) );
01160 
01161   /* open the file for reading */
01162   data->file = fopen(filename, "rb");
01163 
01164   if (!data->file || !data->filename) {
01165     abinit_plugindata_free(data);
01166     return NULL;
01167   }
01168   strcpy(data->filename, filename);
01169 
01170   if (abinit_filetype(data, "GEO"))
01171     result = GEO_open_file_read(data, natoms);
01172   else if (abinit_filetype(data, "DEN") || abinit_filetype(data, "POT") || abinit_filetype(data, "WFK"))
01173     result = DEN_POT_WFK_open_file_read(data, natoms);
01174 
01175   if (result == NULL) abinit_plugindata_free(data);
01176 
01177   DBGPRINT(stderr, "Exit open_file_read\n");
01178   return result;
01179 }
01180 
01181 
01182 /* VMD calls this once to find out about the atom types in the structure */
01183 static int read_structure(void *mydata, int *optflags, molfile_atom_t *atomlist)
01184 {
01185   int result = MOLFILE_ERROR;
01186   abinit_plugindata_t *data = (abinit_plugindata_t *)mydata;
01187 
01188   DBGPRINT(stderr, "Enter read_structure\n");
01189 
01190   if (!data || !optflags || !atomlist) return MOLFILE_ERROR;
01191 
01192   if (abinit_filetype(data, "GEO"))
01193     result = GEO_read_structure(data, optflags, atomlist);
01194   else if (abinit_filetype(data, "DEN") || abinit_filetype(data, "POT") || abinit_filetype(data, "WFK"))
01195     result = DEN_POT_WFK_read_structure(data, optflags, atomlist);
01196 
01197   DBGPRINT(stderr, "Exit read_structure\n");
01198   return result;
01199 }
01200 
01201 
01202 /* VMD keeps calling for a next timestep, until it gets End-Of-File here */
01203 static int read_next_timestep(void *mydata, int natoms, molfile_timestep_t *ts)
01204 {
01205   int result = MOLFILE_EOF;
01206   abinit_plugindata_t *data = (abinit_plugindata_t *)mydata;
01207 
01208   DBGPRINT(stderr, "Enter read_next_timestep\n");
01209 
01210   /* Save coordinatess only if we are given a timestep pointer.
01211    * Otherwise assume that VMD wants us to skip past it.
01212    */
01213   if (!ts || !data) return MOLFILE_EOF;
01214 
01215   /* Double check that the number of atoms are correct */
01216   if (natoms != data->natom) return MOLFILE_EOF;
01217 
01218   if (abinit_filetype(data, "GEO"))
01219     result = GEO_read_next_timestep(data, natoms, ts);
01220   else if (abinit_filetype(data, "DEN") || abinit_filetype(data, "POT") || abinit_filetype(data, "WFK"))
01221     result = DEN_POT_WFK_read_next_timestep(data, natoms, ts);
01222 
01223   DBGPRINT(stderr, "Exit read_next_timestep\n");
01224   return result;
01225 }
01226 
01227 
01228 static void close_file_read(void *mydata)
01229 {
01230   abinit_plugindata_t *data = (abinit_plugindata_t *)mydata;
01231 
01232   DBGPRINT(stderr, "Enter close_read\n");
01233 
01234   abinit_plugindata_free(data);
01235 
01236   DBGPRINT(stderr, "Exit close_read\n");
01237 }
01238 
01239 
01240 
01241 /* This writes the basic structure data in the syntax of the abinit input format. */
01242 static void *open_file_write(const char *filename, const char *filetype, int natoms)
01243 {
01244   abinit_plugindata_t *data = abinit_plugindata_malloc();
01245 
01246   DBGPRINT(stderr, "Enter open_file_write\n");
01247 
01248   if (!data) return NULL;
01249 
01250   /* allocate memory for filename (add an extra 10 bytes for some more flexibility) */
01251   data->filename = (char *)malloc( sizeof(char) * (strlen(filename) + 10) );
01252 
01253   /* open the file for writing */
01254   data->file = fopen(filename, "w");
01255 
01256   if (!data->filename || !data->file) {
01257     abinit_plugindata_free(data);
01258     fprintf(stderr, "ABINIT write) ERROR: unable to open file '%s' for writing\n", filename);
01259     return NULL;
01260   }
01261   strcpy(data->filename, filename);
01262 
01263   data->natom = natoms;
01264 
01265   DBGPRINT(stderr, "Exit open_file_write\n");
01266   return data;
01267 }
01268 
01269 
01270 static int write_structure(void *mydata, int optflags, const molfile_atom_t *atoms)
01271 {
01272   abinit_plugindata_t *data = (abinit_plugindata_t *)mydata;
01273   int i, znucl[NATOM_MAX], ntypat;
01274 
01275   DBGPRINT(stderr, "Enter write_structure\n");
01276 
01277   if (!data || !atoms) return MOLFILE_ERROR;
01278 
01279   for (i = 0; i < NATOM_MAX; ++i) znucl[i] = 0;
01280 
01281 
01282   for (ntypat = i = 0; i < data->natom; ++i) {
01283     int const idx = get_pte_idx(atoms[i].type);
01284     int k;
01285 
01286     /* check if we already have this atom's idx in the list */
01287     for (k = 0; k < ntypat; ++k) if (idx == znucl[k]) break;
01288 
01289     /* if it is not yet in the list, increment ntypat */
01290     if (k == ntypat) ntypat++;
01291 
01292     znucl[k] = idx;
01293     data->typat[i] = k + 1;
01294   }
01295 
01296   /* write header with info */
01297   fprintf(data->file, "# Format below is in a sloppy ABINIT style.\n");
01298   fprintf(data->file, "# See http://www.abinit.org/ for the meaning of the keywords used here.\n\n");
01299 
01300   /* write the atom data to file */
01301   fprintf(data->file, "# Definition of the atom types\nntypat %d\nznucl ", ntypat);
01302   for (i = 0; i < ntypat; ++i) fprintf(data->file, " %d", znucl[i]);
01303   fprintf(data->file, "\n\n");
01304 
01305   fprintf(data->file, "# Definition of the atoms\nnatom %d\ntypat ", data->natom);
01306   for (i = 0; i < data->natom; ++i) fprintf(data->file, " %d", data->typat[i]);
01307   fprintf(data->file, "\n\n");
01308 
01309   DBGPRINT(stderr, "Exit write_structure\n");
01310   return MOLFILE_SUCCESS;
01311 }
01312 
01313 
01314 static int write_timestep(void *mydata, const molfile_timestep_t *ts)
01315 {
01316   abinit_plugindata_t *data = (abinit_plugindata_t *)mydata; 
01317   int i;
01318 
01319   DBGPRINT(stderr, "Enter write_timestep\n");
01320 
01321   if (!data || !ts) return MOLFILE_ERROR;
01322 
01323   fprintf(data->file, "# Definition of the unit cell in Bohr\n");
01324   fprintf(data->file, "acell %f %f %f\n", ts->A * ANGS_TO_BOHR, ts->B * ANGS_TO_BOHR, ts->C * ANGS_TO_BOHR);
01325   fprintf(data->file, "angdeg %f %f %f\n\n", ts->alpha, ts->beta, ts->gamma);
01326 
01327   fprintf(data->file, "# location of the atoms in Bohr\nxcart ");
01328   for (i = 0; i < data->natom; ++i) {
01329     float const rx = ts->coords[3*i    ]  * ANGS_TO_BOHR,
01330                 ry = ts->coords[3*i + 1]  * ANGS_TO_BOHR,
01331                 rz = ts->coords[3*i + 2]  * ANGS_TO_BOHR;
01332     fprintf(data->file, "%s%17.12f %17.12f %17.12f\n", (i != 0 ? "      " : ""), rx, ry, rz);
01333   }
01334   fprintf(data->file, "\n\n");
01335 
01336   DBGPRINT(stderr, "Exit write_timestep\n");
01337   return MOLFILE_SUCCESS;
01338 }
01339 
01340 
01341 static void close_file_write(void *mydata)
01342 {
01343   abinit_plugindata_t *data = (abinit_plugindata_t *)mydata;
01344   DBGPRINT(stderr, "Enter close_file_write\n");
01345 
01346   abinit_plugindata_free(data);
01347 
01348   DBGPRINT(stderr, "Exit close_file_write\n");
01349 }
01350 
01351 
01352 static int read_volumetric_metadata(void *mydata, int *nvolsets, molfile_volumetric_t **metadata)
01353 {
01354   int result = MOLFILE_ERROR;
01355   abinit_plugindata_t *data = (abinit_plugindata_t *)mydata;
01356 
01357   DBGPRINT(stderr, "Enter read_volumetric_metadata\n");
01358 
01359   if (!data || !nvolsets || !metadata) return MOLFILE_ERROR;
01360 
01361   if (abinit_filetype(data, "DEN"))
01362     result = DEN_read_volumetric_metadata(data, nvolsets, metadata);
01363   else if (abinit_filetype(data, "POT"))
01364     result = POT_read_volumetric_metadata(data, nvolsets, metadata);
01365   else if (abinit_filetype(data, "WFK"))
01366     result = WFK_read_volumetric_metadata(data, nvolsets, metadata);
01367 
01368   DBGPRINT(stderr, "Exit read_volumetric_metadata\n");
01369   return result;
01370 }
01371 
01372 static int read_volumetric_data(void *mydata, int set, float *datablock, float *colorblock)
01373 {
01374   int result = MOLFILE_ERROR;
01375   abinit_plugindata_t *data = (abinit_plugindata_t *)mydata;
01376 
01377   DBGPRINT(stderr, "Enter read_volumetric_data\n");
01378 
01379   if (!data || !datablock) return MOLFILE_ERROR;
01380 
01381   if (abinit_filetype(data, "DEN"))
01382     result = DEN_read_volumetric_data(data, set, datablock, colorblock);
01383   else if (abinit_filetype(data, "POT"))
01384     result = POT_read_volumetric_data(data, set, datablock, colorblock);
01385   else if (abinit_filetype(data, "WFK"))
01386     result = WFK_read_volumetric_data(data, set, datablock, colorblock);
01387 
01388   DBGPRINT(stderr, "Exit read_volumetric_data\n");
01389   return result;
01390 }
01391 
01392 
01393 /* ===================================
01394  * Registration stuff
01395  * ===================================
01396  */
01397 
01398 static molfile_plugin_t abinitplugin;
01399 
01400 VMDPLUGIN_API int VMDPLUGIN_init() {
01401   /* zero all bytes in the plugin prior to setting the fields */
01402   memset(&abinitplugin, 0, sizeof(molfile_plugin_t));
01403 
01404   /* see molfile_plugin.h for details on the fields */
01405   /* header  */
01406   abinitplugin.abiversion   = vmdplugin_ABIVERSION; /* the ABI for the base plugin type */
01407   abinitplugin.type         = MOLFILE_PLUGIN_TYPE;  /* string descriptor of the plugin type */
01408   abinitplugin.name         = "ABINIT";             /* name for the plugin */
01409   abinitplugin.prettyname   = "ABINIT";             /* name in filetype list */
01410   abinitplugin.author       = "Rob Lahaye";         /* string identifier */
01411   abinitplugin.majorv       = 0;                    /* major version */
01412   abinitplugin.minorv       = 4;                    /* minor version */
01413   abinitplugin.is_reentrant = VMDPLUGIN_THREADSAFE; /* can this library be run concurrently with itself? */
01414 
01415   /* the rest of the plugin */
01416   abinitplugin.filename_extension       = "*|*_GEO|*_DEN|*_WFK|*_POT|*_VHA|*_VHXC|*_VXC"; /* filename extension for this file type */
01417   abinitplugin.open_file_read           = open_file_read;           /* try to open the file for reading */
01418   abinitplugin.read_structure           = read_structure;           /* read molecular structure from the given file handle */
01419   abinitplugin.read_bonds               = 0;                        /* read bond information for the molecule */
01420   abinitplugin.read_next_timestep       = read_next_timestep;       /* read the next timestep from the file */
01421   abinitplugin.close_file_read          = close_file_read;          /* close the file and release all data */
01422   abinitplugin.open_file_write          = open_file_write;          /* open a coordinate file for writing */
01423   abinitplugin.write_structure          = write_structure;          /* write a timestep to the coordinate file */
01424   abinitplugin.write_timestep           = write_timestep;           /* write a timestep to the coordinate file */
01425   abinitplugin.close_file_write         = close_file_write;         /* close the file and release all data */
01426   abinitplugin.read_volumetric_metadata = read_volumetric_metadata; /* retrieve metadata pertaining to volumetric datasets in this file */
01427   abinitplugin.read_volumetric_data     = read_volumetric_data;     /* read the specified volumetric data set into the space pointed to by datablock */
01428   abinitplugin.read_rawgraphics         = 0;                        /* read raw graphics data stored in this file */
01429   abinitplugin.read_molecule_metadata   = 0;                        /* read molecule metadata */
01430   abinitplugin.write_bonds              = 0;                        /* write bond information for the molecule */
01431 
01432   return VMDPLUGIN_SUCCESS;
01433 }
01434 
01435 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
01436 
01437   (*cb)(v, (vmdplugin_t *)&abinitplugin);
01438   return VMDPLUGIN_SUCCESS;
01439 }
01440 
01441 VMDPLUGIN_API int VMDPLUGIN_fini() {
01442   return VMDPLUGIN_SUCCESS;
01443 }
01444 
01445 
01446 /* ===================================
01447  * Procedures to read the binary data
01448  * ===================================
01449  */
01450 
01451 
01452 /* Do a binary read of 'size' bytes from stream.
01453  * Store the bytes into the pointer according to the match
01454  * between the endianness of the file and the system reading it.
01455  * Return the number of bytes succesfully read.
01456  */
01457 static int binread(void *ptr, size_t size, FILE *stream, binary_t bintype)
01458 {
01459   char *data = (char *)ptr;
01460   char *storage = malloc(size * sizeof(char));
01461   int const result = fread(storage, 1, size, stream);
01462   unsigned int i = 1;
01463   enum Endianness myEndian = ( *(char *)&i == 1 ? little_endian : big_endian );
01464 
01465   for (i = 0; i < size; ++i) {
01466     int const index = (myEndian == bintype.endian ? i : size - i - 1);
01467     data[i] = storage[index];
01468   }
01469 
01470   free(storage);
01471   return result;
01472 }
01473 
01474 
01475 /* Read the full header of a binary ABINIT file until the point
01476  * where the data starts.
01477  * Return a pointer to the header struct, or NULL if the reading fails.
01478  */ 
01479 static abinit_binary_header_t *abinit_header(FILE *fp)
01480 {
01481   int const debug = 0;
01482 
01483   /* the pattern of the first 4 or 8 bytes of a binary abinit file determines the
01484    * endianness and recordmarker length of the file.
01485    */
01486   char const Big_Endian_4_pattern[]    = {'\x00','\x00','\x00','\x0e'},
01487              Big_Endian_8_pattern[]    = {'\x00','\x00','\x00','\x00','\x00','\x00','\x00','\x0e'}, 
01488              Little_Endian_4_pattern[] = {'\x0e','\x00','\x00','\x00'},
01489              Little_Endian_8_pattern[] = {'\x0e','\x00','\x00','\x00','\x00','\x00','\x00','\x00'};
01490 
01491   char skip[1024];
01492   int i, bc = 0;
01493 
01494   abinit_binary_header_t *hdr = abinit_header_malloc();
01495   if (!hdr) return NULL;
01496 
01497   /* be sure to start from the beginning of the file */
01498   rewind(fp);
01499 
01500   if (debug) fprintf(stderr, "START OF BINARY FILE DEBUG INFO\n");
01501 
01502   /* Determine endianness and record length of the binary file.
01503    * Note that this pattern is followed by a version string, which has no null termination.
01504    * Therefore the pattern is always followed by a non-null byte.
01505    */
01506   fread(skip, 1, 8, fp);
01507   if      (memcmp(skip, Big_Endian_4_pattern, 4) == 0)    {hdr->bintype.endian = big_endian;    hdr->bintype.recordmarker = 4;}
01508   else if (memcmp(skip, Big_Endian_8_pattern, 8) == 0)    {hdr->bintype.endian = big_endian;    hdr->bintype.recordmarker = 8;}
01509   else if (memcmp(skip, Little_Endian_8_pattern, 8) == 0) {hdr->bintype.endian = little_endian; hdr->bintype.recordmarker = 8;}
01510   else if (memcmp(skip, Little_Endian_4_pattern, 4) == 0) {hdr->bintype.endian = little_endian; hdr->bintype.recordmarker = 4;}
01511   else {abinit_header_free(hdr); return NULL;}
01512 
01513   if (debug) fprintf(stderr, "Binary file is in %s with a record-marker of %d bytes\n",
01514                                (hdr->bintype.endian == big_endian ? "Big Endian" : "Little Endian"), hdr->bintype.recordmarker);
01515 
01516   /* start reading again from the beginning of the file */
01517   rewind(fp);
01518 
01519   /*The header
01520    * The wavefunction files, density files, and potential files all begin with the same records,
01521    * called the "header".
01522    *
01523    *character*6 :: codvsn
01524    *integer :: headform,fform
01525    *integer :: bantot,date,intxc,ixc,natom,ngfft(3),nkpt,npsp,
01526    *nspden,nspinor,nsppol,nsym,ntypat,occopt,pertcase,usepaw
01527    *integer :: usewvl, cplex, nspden
01528    *double precision :: acell(3),ecut,ecutdg,ecutsm,ecut_eff,qptn(3),rprimd(3,3),stmbias,tphysel,tsmear
01529    *integer :: istwfk(nkpt),nband(nkpt*nsppol),npwarr(nkpt),so_psp(npsp),&
01530    *& symafm(nsym),symrel(3,3,nsym),typat(natom),nrhoijsel(nspden),rhoijselect(*,nspden)
01531    *double precision :: kpt(3,nkpt),occ(bantot),tnons(3,nsym),znucltypat(ntypat),wtk(nkpt)
01532    *character*132 :: title
01533    *double precision :: znuclpsp,zionpsp
01534    *integer :: pspso,pspdat,pspcod,pspxc,lmax,lloc,mmax=integers
01535    *double precision :: residm,xred(3,natom),etotal,fermie,rhoij(*,nspden)
01536    *
01537    * write(unit=header) codvsn,headform,fform
01538    * write(unit=header) bantot,date,intxc,ixc,natom,ngfft(1:3),&
01539    *& nkpt,nspden,nspinor,nsppol,nsym,npsp,ntypat,occopt,pertcase,usepaw,&
01540    *& ecut,ecutdg,ecutsm,ecut_eff,qptn(1:3),rprimd(1:3,1:3),stmbias,tphysel,tsmear,usewvl
01541 
01542    * write(unit=header) istwfk(1:nkpt),nband(1:nkpt*nsppol),&
01543    *& npwarr(1:nkpt),so_psp(1:npsp),symafm(1:nsym),symrel(1:3,1:3,1:nsym),typat(1:natom),&
01544    *& kpt(1:3,1:nkpt),occ(1:bantot),tnons(1:3,1:nsym),znucltypat(1:ntypat),wtk(1:nkpt)
01545    * do ipsp=1,npsp
01546    *! (npsp lines, 1 for each pseudopotential ; npsp=ntypat, except if alchemical pseudo-atoms)
01547    *  write(unit=unit) title,znuclpsp,zionpsp,pspso,pspdat,pspcod,pspxc,lmn_size
01548    * enddo
01549    *!(in case of usepaw==0, final record: residm, coordinates, total energy, Fermi energy)
01550    * write(unit=unit) residm,xred(1:3,1:natom),etotal,fermie
01551    *!(in case of usepaw==1, there are some additional records)
01552    * if (usepaw==1)then
01553    *  write(unit=unit)( pawrhoij(iatom)%nrhoijsel(1:nspden),iatom=1,natom), cplex, nspden
01554    *  write(unit=unit)((pawrhoij(iatom)%rhoijselect(1:      nrhoijsel(ispden),ispden),ispden=1,nspden),iatom=1,natom),&
01555    *&                 ((pawrhoij(iatom)%rhoijp     (1:cplex*nrhoijsel(ispden),ispden),ispden=1,nspden),iatom=1,natom)
01556    * endif
01557    */
01558 
01559   /* skip first the recordmarker bytes */
01560   bc += fread(skip, 1, hdr->bintype.recordmarker, fp);
01561 
01562   /* code version */
01563   bc += fread(hdr->codvsn, sizeof(char), 6, fp); hdr->codvsn[6] = '\0';
01564   if (debug) fprintf(stderr, "codvsn = '%s' (code version)\n", hdr->codvsn); 
01565 
01566   /* format of the header */
01567   bc += binread(&hdr->headform, 4, fp, hdr->bintype);
01568   if (debug) fprintf(stderr, "headform = '%d' (format of the header)\n", hdr->headform);
01569 
01570   /* specification for data type: 2 for wf; 52 for density; 102 for potential */
01571   bc += binread(&hdr->fform, 4, fp, hdr->bintype);
01572   if (debug) fprintf(stderr, "fform = '%d' (2 for wf; 52 for density; 102 for potential)\n", hdr->fform);
01573 
01574   /* skip the recordmarker bytes two times */
01575   bc += fread(skip, 1, 2 * hdr->bintype.recordmarker, fp);
01576 
01577   /* total number of bands (sum of nband on all kpts and spins) */
01578   bc += binread(&hdr->bantot, 4, fp, hdr->bintype);
01579   if (debug) fprintf(stderr, "bantot = '%d' (sum of nband on all kpts and spins)\n", hdr->bantot);
01580 
01581   /* starting date */
01582   bc += binread(&hdr->date, 4, fp, hdr->bintype);
01583   if (debug) fprintf(stderr, "date = '%d' (starting date)\n", hdr->date);
01584 
01585   bc += binread(&hdr->intxc, 4, fp, hdr->bintype);
01586   if (debug) fprintf(stderr, "intxc = '%d'\n", hdr->intxc);
01587 
01588   bc += binread(&hdr->ixc, 4, fp, hdr->bintype);
01589   if (debug) fprintf(stderr, "ixc = '%d'\n", hdr->ixc);
01590 
01591   bc += binread(&hdr->natom, 4, fp, hdr->bintype);
01592   if (debug) fprintf(stderr, "natom = '%d'\n", hdr->natom);
01593 
01594   /* failsafe */
01595   if (hdr->natom <= 0) {
01596     fprintf(stderr, "ABINIT read) ERROR Binary Header: natom = %d is wrong!",  hdr->natom);
01597     abinit_header_free(hdr);
01598     return NULL;
01599   }
01600 
01601   bc += binread(&hdr->ngfft[0], 4, fp, hdr->bintype);
01602   if (debug) fprintf(stderr, "ngfft[0] = '%d'\n", hdr->ngfft[0]);
01603   bc += binread(&hdr->ngfft[1], 4, fp, hdr->bintype);
01604   if (debug) fprintf(stderr, "ngfft[1] = '%d'\n", hdr->ngfft[1]);
01605   bc += binread(&hdr->ngfft[2], 4, fp, hdr->bintype);
01606   if (debug) fprintf(stderr, "ngfft[2] = '%d'\n", hdr->ngfft[2]);
01607 
01608   bc += binread(&hdr->nkpt, 4, fp, hdr->bintype);
01609   if (debug) fprintf(stderr, "nkpt = '%d'\n", hdr->nkpt);
01610 
01611   bc += binread(&hdr->nspden, 4, fp, hdr->bintype);
01612   if (debug) fprintf(stderr, "nspden = '%d'\n", hdr->nspden);
01613 
01614   /* failsafe */
01615   if (hdr->nspden != 1 && hdr->nspden != 2 && hdr->nspden != 4) {
01616     fprintf(stderr, "ABINIT read) ERROR Binary Header: nspden = %d is wrong!",  hdr->nspden);
01617     abinit_header_free(hdr);
01618     return NULL;
01619   }
01620 
01621   bc += binread(&hdr->nspinor, 4, fp, hdr->bintype);
01622   if (debug) fprintf(stderr, "nspinor = '%d'\n", hdr->nspinor);
01623 
01624   bc += binread(&hdr->nsppol, 4, fp, hdr->bintype);
01625   if (debug) fprintf(stderr, "nsppol = '%d'\n", hdr->nsppol);
01626 
01627   bc += binread(&hdr->nsym, 4, fp, hdr->bintype);
01628   if (debug) fprintf(stderr, "nsym = '%d'\n", hdr->nsym);
01629 
01630   bc += binread(&hdr->npsp, 4, fp, hdr->bintype);
01631   if (debug) fprintf(stderr, "npsp = '%d'\n", hdr->npsp);
01632 
01633   bc += binread(&hdr->ntypat, 4, fp, hdr->bintype);
01634   if (debug) fprintf(stderr, "ntypat = '%d'\n", hdr->ntypat);
01635 
01636   bc += binread(&hdr->occopt, 4, fp, hdr->bintype);
01637   if (debug) fprintf(stderr, "occopt = '%d'\n", hdr->occopt);
01638 
01639   /* the index of the perturbation, 0 if GS calculation */
01640   bc += binread(&hdr->pertcase, 4, fp, hdr->bintype);
01641   if (debug) fprintf(stderr, "pertcase = '%d' (the index of the perturbation, 0 if GS calculation)\n", hdr->pertcase);
01642 
01643   bc += binread(&hdr->usepaw, 4, fp, hdr->bintype);
01644   if (debug) fprintf(stderr, "usepaw = '%d' (0=norm-conserving psps, 1=paw)\n", hdr->usepaw);
01645 
01646   /* failsafe */
01647   if (hdr->usepaw != 0 && hdr->usepaw != 1) {
01648     fprintf(stderr, "ABINIT read) ERROR Binary Header: usepaw = %d is wrong!", hdr->usepaw);
01649     abinit_header_free(hdr);
01650     return NULL;
01651   }
01652 
01653   bc += binread(&hdr->ecut, 8, fp, hdr->bintype);
01654   if (debug) fprintf(stderr, "ecut = '%g'\n", hdr->ecut);
01655 
01656   /* input variable (ecut for NC psps, pawecutdg for paw) */
01657   bc += binread(&hdr->ecutdg, 8, fp, hdr->bintype);
01658   if (debug) fprintf(stderr, "ecutdg = '%g' (ecut for NC psps, pawecutdg for paw)\n", hdr->ecutdg);
01659 
01660   bc += binread(&hdr->ecutsm, 8, fp, hdr->bintype);
01661   if (debug) fprintf(stderr, "ecutsm = '%g'\n", hdr->ecutsm);
01662 
01663   /* ecut*dilatmx**2 (dilatmx is an input variable) */
01664   bc += binread(&hdr->ecut_eff, 8, fp, hdr->bintype);
01665   if (debug) fprintf(stderr, "ecut_eff = '%g' (ecut*dilatmx**2 [dilatmx is an input variable])\n", hdr->ecut_eff);
01666 
01667   bc += binread(&hdr->qptn[0], 8, fp, hdr->bintype);
01668   if (debug) fprintf(stderr, "qptn[0] = '%g'\n", hdr->qptn[0]);
01669   bc += binread(&hdr->qptn[1], 8, fp, hdr->bintype);
01670   if (debug) fprintf(stderr, "qptn[1] = '%g'\n", hdr->qptn[1]);
01671   bc += binread(&hdr->qptn[2], 8, fp, hdr->bintype);
01672   if (debug) fprintf(stderr, "qptn[2] = '%g'\n", hdr->qptn[2]);
01673 
01674   for (i = 0; i < 3; ++i) {
01675     bc += binread(&hdr->rprimd[i][0], 8, fp, hdr->bintype);
01676     if (debug) fprintf(stderr, "rprimd[%d][0] = '%g'\n", i, hdr->rprimd[i][0]);
01677     bc += binread(&hdr->rprimd[i][1], 8, fp, hdr->bintype);
01678     if (debug) fprintf(stderr, "rprimd[%d][1] = '%g'\n", i, hdr->rprimd[i][1]);
01679     bc += binread(&hdr->rprimd[i][2], 8, fp, hdr->bintype);
01680     if (debug) fprintf(stderr, "rprimd[%d][2] = '%g'\n", i, hdr->rprimd[i][2]);
01681   }
01682 
01683   bc += binread(&hdr->stmbias, 8, fp, hdr->bintype);
01684   if (debug) fprintf(stderr, "stmbias = '%g'\n", hdr->stmbias);
01685 
01686   bc += binread(&hdr->tphysel, 8, fp, hdr->bintype);
01687   if (debug) fprintf(stderr, "tphysel = '%g'\n", hdr->tphysel);
01688 
01689   bc += binread(&hdr->tsmear, 8, fp, hdr->bintype);
01690   if (debug) fprintf(stderr, "tsmear = '%g'\n", hdr->tsmear);
01691 
01692   bc += binread(&hdr->usewvl, 4, fp, hdr->bintype);
01693   if (debug) fprintf(stderr, "usewvl = '%d'\n", hdr->usewvl);
01694 
01695   /* failsafe */
01696   if (hdr->usewvl != 0 && hdr->usewvl != 1) {
01697     fprintf(stderr, "ABINIT read) ERROR Binary Header: usewvl = %d is wrong!",  hdr->usewvl);
01698     abinit_header_free(hdr);
01699     return NULL;
01700   }
01701 
01702   /* skip the recordmarker bytes two times */
01703   bc += fread(skip, 1, 2 * hdr->bintype.recordmarker, fp);
01704 
01705   hdr->istwfk = (int *)malloc(sizeof(int) * hdr->nkpt);
01706   hdr->nband  = (int *)malloc(sizeof(int) * hdr->nkpt * hdr->nsppol);
01707   hdr->npwarr = (int *)malloc(sizeof(int) * hdr->nkpt);
01708   hdr->so_psp = (int *)malloc(sizeof(int) * hdr->npsp);
01709   hdr->symafm = (int *)malloc(sizeof(int) * hdr->nsym);
01710   hdr->typat  = (int *)malloc(sizeof(int) * hdr->natom);
01711 
01712   if (!hdr->istwfk || !hdr->nband || !hdr->npwarr || !hdr->so_psp || !hdr->symafm || !hdr->typat) {
01713     abinit_header_free(hdr);
01714     return NULL;
01715   }
01716   for (i = 0; i < 3; ++i) {
01717     int j;
01718     for (j = 0; j < 3; ++j) {
01719       hdr->symrel[i][j] = (int *)malloc(sizeof(int) * hdr->nsym);
01720       if (!hdr->symrel[i][j]) {abinit_header_free(hdr); return NULL;}
01721     }
01722   }
01723 
01724   for (i = 0; i < hdr->nkpt; ++i) {
01725     bc += binread(&hdr->istwfk[i], 4, fp, hdr->bintype);
01726     if (debug) fprintf(stderr, "istwfk[%d] = '%d'\n", i, hdr->istwfk[i]);
01727   }
01728 
01729   for (i = 0; i < hdr->nkpt * hdr->nsppol; ++i) {
01730     bc += binread(&hdr->nband[i], 4, fp, hdr->bintype);
01731     if (debug) fprintf(stderr, "nband[%d] = '%d'\n", i, hdr->nband[i]);
01732   }
01733 
01734   for (i = 0; i < hdr->nkpt; ++i) {
01735     bc += binread(&hdr->npwarr[i], 4, fp, hdr->bintype);
01736     if (debug) fprintf(stderr, "npwarr[%d] = '%d'\n", i, hdr->npwarr[i]);
01737   }
01738 
01739   for (i = 0; i < hdr->npsp; ++i) {
01740     bc += binread(&hdr->so_psp[i], 4, fp, hdr->bintype);
01741     if (debug) fprintf(stderr, "so_psp[%d] = '%d'\n", i, hdr->so_psp[i]);
01742   }
01743 
01744   for (i = 0; i < hdr->nsym; ++i) {
01745     bc += binread(&hdr->symafm[i], 4, fp, hdr->bintype);
01746     if (debug) fprintf(stderr, "symafm[%d] = '%d'\n", i, hdr->symafm[i]);
01747   }
01748 
01749   for (i = 0; i < hdr->nsym; ++i) {
01750     int j;
01751     for (j = 0; j < 3; ++j) {
01752       int k;     
01753       for (k = 0; k < 3; ++k) {
01754         bc += binread(&hdr->symrel[k][j][i], 4, fp, hdr->bintype);
01755         if (debug) fprintf(stderr, "symrel[%d][%d][%2d]= '%2d'", k, j, i, hdr->symrel[k][j][i]);
01756       }
01757       if (debug) fprintf(stderr, "\n");
01758     }
01759     if (debug) fprintf(stderr, "\n");
01760   }
01761 
01762   for (i = 0; i < hdr->natom; ++i) {
01763     bc += binread(&hdr->typat[i], 4, fp, hdr->bintype);
01764     if (debug) fprintf(stderr, "typat[%d] = '%d'\n", i, hdr->typat[i]);
01765   }
01766 
01767   for (i = 0; i < 3; ++i) {
01768     hdr->kpt[i] = (double *)malloc(sizeof(double) * hdr->nkpt);
01769     hdr->tnons[i] = (double *)malloc(sizeof(double) * hdr->nsym);
01770     if (!hdr->kpt[i] || !hdr->tnons[i]) {abinit_header_free(hdr); return NULL;}
01771   }
01772 
01773   hdr->occ = (double *)malloc(sizeof(double) * hdr->bantot);
01774   hdr->znucltypat = (double *)malloc(sizeof(double) * hdr->ntypat);
01775   hdr->wtk = (double *)malloc(sizeof(double) * hdr->nkpt);
01776   if (!hdr->occ || !hdr->znucltypat || !hdr->wtk) {abinit_header_free(hdr); return NULL;}
01777 
01778   for (i = 0; i < hdr->nkpt; ++i) {
01779     int j;
01780     for (j = 0; j < 3; ++j) {
01781       bc += binread(&hdr->kpt[j][i], 8, fp, hdr->bintype);
01782       if (debug) fprintf(stderr, "kpt[%d][%2d] = '%g'\n", j, i, hdr->kpt[j][i]);
01783     }
01784   }
01785 
01786   for (i = 0; i < hdr->bantot; ++i) {
01787     bc += binread(&hdr->occ[i], 8, fp, hdr->bintype);
01788     if (debug) fprintf(stderr, "occ[%d] = '%g'\n", i, hdr->occ[i]);
01789   }
01790 
01791   for (i = 0; i < hdr->nsym; ++i) {
01792     int j;
01793     for (j = 0; j < 3; ++j) {
01794       bc += binread(&hdr->tnons[j][i], 8, fp, hdr->bintype);
01795       if (debug) fprintf(stderr, "tnons[%d][%2d] = '%g' ", j, i, hdr->tnons[j][i]);
01796     }
01797     if (debug) fprintf(stderr, "\n");
01798   }
01799 
01800   for (i = 0; i < hdr->ntypat; ++i) {
01801     bc += binread(&hdr->znucltypat[i], 8, fp, hdr->bintype);
01802     if (debug) fprintf(stderr, "znucltypat[%d] = '%g'\n", i, hdr->znucltypat[i]);
01803   }
01804 
01805   for (i = 0; i < hdr->nkpt; ++i) {
01806     bc += binread(&hdr->wtk[i], 8, fp, hdr->bintype);
01807     if (debug) fprintf(stderr, "wtk[%d] = '%g'\n", i, hdr->wtk[i]);
01808   }
01809 
01810   for (i = 0; i < hdr->npsp; ++i) {
01811 
01812     /* skip the recordmarker bytes two times */
01813     bc += fread(skip, 1, 2 * hdr->bintype.recordmarker, fp);
01814 
01815     bc += fread(hdr->title, sizeof(char), 132, fp); hdr->title[132] = '\0';
01816     if (debug) fprintf(stderr, "title = '%s'\n", hdr->title);
01817 
01818     bc += binread(&hdr->znuclpsp, 8, fp, hdr->bintype);
01819     if (debug) fprintf(stderr, "znuclpsp = '%g'\n", hdr->znuclpsp);
01820 
01821     bc += binread(&hdr->zionpsp, 8, fp, hdr->bintype);
01822     if (debug) fprintf(stderr, "zionpsp = '%g'\n", hdr->zionpsp);
01823 
01824     bc += binread(&hdr->pspso, 4, fp, hdr->bintype);
01825     if (debug) fprintf(stderr, "pspso = '%d'\n", hdr->pspso);
01826 
01827     bc += binread(&hdr->pspdat, 4, fp, hdr->bintype);
01828     if (debug) fprintf(stderr, "pspdat = '%d'\n", hdr->pspdat);
01829 
01830     bc += binread(&hdr->pspcod, 4, fp, hdr->bintype);
01831     if (debug) fprintf(stderr, "pspcod = '%d'\n", hdr->pspcod);
01832 
01833     bc += binread(&hdr->pspxc, 4, fp, hdr->bintype);
01834     if (debug) fprintf(stderr, "pspxc = '%d'\n", hdr->pspxc);
01835 
01836     bc += binread(&hdr->lmn_size, 4, fp, hdr->bintype);
01837     if (debug) fprintf(stderr, "lmn_size = '%d'\n", hdr->lmn_size);
01838   }
01839 
01840   /* skip the recordmarker bytes two times */
01841   bc += fread(skip, 1, 2 * hdr->bintype.recordmarker, fp);
01842 
01843   for (i = 0; i < 3; ++i) {
01844     hdr->xred[i] = (double *)malloc(sizeof(double) * hdr->natom);
01845     if (!hdr->xred[i]) {abinit_header_free(hdr); return NULL;}
01846   }
01847 
01848   bc += binread(&hdr->residm, 8, fp, hdr->bintype);
01849   if (debug) fprintf(stderr, "residm = '%g'\n", hdr->residm);
01850 
01851   for (i = 0; i < hdr->natom; ++i) {
01852     int j;
01853     for (j = 0; j < 3; ++j) {
01854       bc += binread(&hdr->xred[j][i], 8, fp, hdr->bintype);
01855       if (debug) fprintf(stderr, "xred[%d][%d] = '%g'\n", j, i, hdr->xred[j][i]);
01856 
01857       /* failsafe */
01858       if (hdr->xred[j][i] < -1 || hdr->xred[j][i] > 1) {
01859         fprintf(stderr, "Binary Header Error: hdr->xred[%d][%d] = %g; something must be wrong!", j, i, hdr->xred[j][i]);
01860         {abinit_header_free(hdr); return NULL;}
01861       }
01862     }
01863   }
01864 
01865   bc += binread(&hdr->etotal, 8, fp, hdr->bintype);
01866   if (debug) fprintf(stderr, "etotal = '%g'\n", hdr->etotal);
01867 
01868   bc += binread(&hdr->fermie, 8, fp, hdr->bintype);
01869   if (debug) fprintf(stderr, "fermie = '%g'\n", hdr->fermie);
01870 
01871   if (hdr->usepaw == 1) {
01872     struct pawrhoij_t {
01873       int *nrhoijsel;
01874       int **rhoijselect;
01875       double **rhoij;
01876     };
01877     struct pawrhoij_t *pawrhoij = malloc(hdr->natom * sizeof(struct pawrhoij_t));
01878 
01879     for (i = 0; i < hdr->natom; ++i) {
01880       int j;
01881       pawrhoij[i].nrhoijsel = (int *)malloc(sizeof(int) * hdr->nspden);
01882       if (!pawrhoij[i].nrhoijsel) {abinit_header_free(hdr); return NULL;}
01883       for (j = 0; j < hdr->nspden; ++j){
01884         bc += binread(&pawrhoij[i].nrhoijsel[j], 4, fp, hdr->bintype);
01885         if (debug) fprintf(stderr, "pawrhoij[%d].nrhoijsel[%d] = '%d'\n", i, j, pawrhoij[i].nrhoijsel[j]);
01886       }
01887       pawrhoij[i].rhoijselect = (int **)malloc(sizeof(int) * hdr->nspden);
01888       pawrhoij[i].rhoij = (double **)malloc(sizeof(double) * hdr->nspden);
01889       if (!pawrhoij[i].rhoijselect || !pawrhoij[i].rhoij) {abinit_header_free(hdr); return NULL;}
01890       for (j = 0; j < hdr->nspden; ++j) {
01891         int k;
01892         pawrhoij[i].rhoijselect[j] = (int *)malloc(sizeof(int) * pawrhoij[i].nrhoijsel[j]);
01893         pawrhoij[i].rhoij[j] = (double *)malloc(sizeof(double) * pawrhoij[i].nrhoijsel[j]);
01894         if (!pawrhoij[i].rhoijselect[j]) {abinit_header_free(hdr); return NULL;}
01895         for (k = 0; k < pawrhoij[i].nrhoijsel[j]; ++k) {
01896           bc += binread(&pawrhoij[i].rhoijselect[j][k], 4, fp, hdr->bintype);
01897          if (debug) fprintf(stderr, "pawrhoij[%d].rhoijselect[%d][%d] = '%d'\n", i, j, k, pawrhoij[i].rhoijselect[j][k]);
01898         }
01899         for (k = 0; k < pawrhoij[i].nrhoijsel[j]; ++k) {
01900           bc += binread(&pawrhoij[i].rhoij[j][k], 8, fp, hdr->bintype);
01901          if (debug) fprintf(stderr, "pawrhoij[%d].rhoij[%d][%d] = '%g'\n", i, j, k, pawrhoij[i].rhoij[j][k]);
01902         }
01903       }
01904       for (j = 0; j < hdr->nspden; ++j) {
01905         free(pawrhoij[i].rhoijselect[j]);
01906         free(pawrhoij[i].rhoij[j]);
01907       }
01908       free(pawrhoij[i].rhoijselect);
01909       free(pawrhoij[i].rhoij);
01910       free(pawrhoij[i].nrhoijsel);
01911     }
01912     free(pawrhoij);
01913   } /* end of "if (usepaw == 1)" */
01914 
01915   /* cplex depends on other variables:
01916    * In GS calculations, cplex = 1.
01917    * In response function calculations (non-zero pertcase), cplex = 1 at the gamma point (qpt=(0,0,0))
01918    * and cplex = 2 if qpt/=(0,0,0).
01919    */
01920   hdr->cplex = 1;
01921   if (hdr->pertcase != 0 && fabs(hdr->qptn[0]) > 10e-6 && fabs(hdr->qptn[1]) > 10e-6 && fabs(hdr->qptn[2]) > 10e-6) hdr->cplex = 2;
01922 
01923 
01924   /* skip the recordmarker bytes two times */
01925   bc += fread(skip, 1, 2 * hdr->bintype.recordmarker, fp);
01926 
01927   if (debug) fprintf(stderr, "END OF BINARY FILE DEBUG INFO\n");
01928 
01929   return hdr;
01930 }

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