The ABINIT Users Mailing List ( CLOSED )

[Rob <spamrefuse@yahoo.com>]

Hi,

In case you use VMD for your visualization:

This plugin helps VMD to read directly formatted GEO and unformatted
DEN files from ABINIT.

This is an alpha-version. I might put some more work in it in the near
future to also read other ABINIT output files. I wrote this plugin
primarily for myself, but hopefully it can be useful to others too.

I have tested it on my own Linux system, on which I also run the ABINIT
code. I'm not sure how machine-independent my code is for reading
the binary DEN files.

Let me know what you think of it, if you use VMD and want to try it out.
The source code is attached.

All you have to do is:
$ gcc -O2 -Wall -I. -I$VMDBASEDIR/plugins/include -c abinitplugin.c
$ ld -shared -o abinitplugin.so abinitplugin.o

Then
  copy abinitplugin.so $VMDBASEDIR/plugins/$ARCH/molfile

Rob.



/***************************************************************************
 * RCS INFORMATION:
 *
 *      $RCSfile:  $
 *      $Author:  $       $Locker:  $             $State:  $
 *      $Revision:  $       $Date:  $
 *
 ***************************************************************************/

/*
 *  ABINIT plugin for VMD
 *  Rob Lahaye, Sungkyunkwan University, Korea 
 *  
 *  ABINIT manual   
 *  http://www.abinit.org/
 * 
 *  LINUX
 *  gcc -O2 -Wall -I. -I$VMDBASEDIR/plugins/include -c abinitplugin.c
 *  ld -shared -o abinitplugin.so abinitplugin.o
 *
 *  MACOSX
 *  gcc -O2 -Wall -I. -I$VMDBASEDIR/plugins/include -c abinitplugin.c
 *  gcc -bundle -o abinitplugin.so abinitplugin.o
 *
 *  Install
 *  copy abinitplugin.so $VMDBASEDIR/plugins/$ARCH/molfile
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include "molfile_plugin.h"
#include "periodic_table.h"
#include "unit_conversion.h"

#define LINESIZE 2048  // maximum length of a line
#define NATOM_MAX 100  // maximum number of atoms

#define DEBUG 0
#define DBGPRINT if (0) fprintf

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

typedef struct {
 char codvsn[7];
 int headform,fform;
 int bantot, date, intxc, ixc, natom, ngfft[3], nkpt, npsp,
     nspden, nspinor, nsppol, nsym, ntypat, occopt, pertcase, usepaw;
 double ecut, ecutdg, ecutsm, ecut_eff, qptn[3], rprimd[3][3], stmbias, tphysel, tsmear;
 int usewvl, *istwfk, *nband, *npwarr, *so_psp, *symafm, *symrel[3][3], *typat;
 double *kpt[3], *occ, *tnons[3], *znucltypat, *wtk;

 char title[133];
 double znuclpsp, zionpsp;
 int pspso, pspdat, pspcod, pspxc, lmn_size;

 double residm, *xred[3], etotal, fermie;

 int cplex;
} abinit_binary_header_t;


typedef struct {
  FILE *file;                   /* file pointer for reading or writing */
  char *filename;               /* file name for reading or writing */
  char filetype[3];             /* type of data file: GEOmetric or DENsity file */

  /* variables for atomic structure */
  float rprimd[3][3];           /* Real space PRIMitive translations of the unit cell */
  int natom;                    /* total number of atoms */
  int typat[NATOM_MAX];         /* number of atoms per atom type */
  molfile_atom_t *atomlist;

  /* volumetric variables for data of charge density, wavefunction, etc. */
  int nvolsets;                 /* number of volumetric datasets */
  molfile_volumetric_t *vol;    /* volume set metadata */
  
  abinit_binary_header_t *hdr;  /* header info of binary file */
} abinit_plugindata_t;


static abinit_binary_header_t *abinit_header(FILE *);


/* allocate memory for plugin data, NULLify the pointers, and initialize numbers */
static abinit_plugindata_t *abinit_plugindata_malloc()
{
  int i;

  abinit_plugindata_t *data = (abinit_plugindata_t *)malloc(sizeof(abinit_plugindata_t));

  if (!data) {
    fprintf(stderr, "\n\nABINIT plugin) ERROR: cannot allocate memory for plugin data.\n");
    return NULL;
  }

  data->file = NULL;
  data->filename = NULL;
  data->filetype[0]='\0';
  data->atomlist = NULL;
  data->vol = NULL;
  data->hdr = NULL;

  /* Initialize atom numbers to zero */
  data->natom = data->nvolsets = 0;
  for (i = 0; i < NATOM_MAX; ++i) data->typat[i] = 0;

  /* Initialize to unit vectors */
  for (i = 0; i < 3; ++i) {
    data->rprimd[i][0] = (i == 0 ? 1 : 0);
    data->rprimd[i][1] = (i == 1 ? 1 : 0);
    data->rprimd[i][2] = (i == 2 ? 1 : 0);
  }

  return data;
}


/* free up the header data */
static void abinit_headerdata_free(abinit_binary_header_t *hdr)
{
  int i;

  if (!hdr) return;

  if (hdr->istwfk) free(hdr->istwfk);
  if (hdr->nband) free(hdr->nband);
  if (hdr->npwarr) free(hdr->npwarr);
  if (hdr->so_psp) free(hdr->so_psp);
  if (hdr->symafm) free(hdr->symafm);
  for (i = 0; i < 3; ++i) {
    int j;
    for (j = 0; j < 3; ++j) if (hdr->symrel[i][j]) free(hdr->symrel[i][j]);
    if (hdr->kpt[i]) free(hdr->kpt[i]);
    if (hdr->tnons[i]) free(hdr->tnons[i]);
    if (hdr->xred[i]) free(hdr->xred[i]);
  }
  if (hdr->typat) free(hdr->typat);
  if (hdr->occ) free(hdr->occ);
  if (hdr->znucltypat) free(hdr->znucltypat);
  if (hdr->wtk) free(hdr->wtk);

  free(hdr);
  hdr = NULL;
}


/* free up the plugin data */
static void abinit_plugindata_free(abinit_plugindata_t *data)
{
  if (!data) return;

  if (data->file) fclose(data->file);
  if (data->filename) free(data->filename);
  if (data->atomlist) free(data->atomlist);
  if (data->vol) free(data->vol);

  abinit_headerdata_free(data->hdr);

  free(data);
  data = NULL;
}


/* read a non-empty line from stream; also remove comments (#...)
 * return the stripped line, or NULL when EOF is reached
 */
static char *abinit_readline(char *line, FILE *stream)
{
  char *status;

  if (!line || !stream) return NULL;

  do {
    char *cptr;

    // read one line from the stream
    status = fgets(line, LINESIZE, stream);

    // first remove comment from the line
    cptr = strchr(line, '#');
    if (cptr) *cptr = '\0';

    // next remove redundant white spaces at the end of the line
    for (cptr = &line[strlen(line) - 1]; isspace(*cptr); --cptr) *cptr = '\0';

    // continue for as long as EOF is not reached and the line is empty
  } while (status != NULL && strlen(line) == 0);

  return status;
}


/* Abinit uses and produces several types of files:
 *   regular output file (formatted/ascii text): foobar.out
 *   geometry file (formatted/ascii text): foobar_GEO
 *   charge density file (unformatted/binary): foobar_DEN
 *   wavefunction file (unformatted/binary): foobar_WFK
 *   potential file (unformatted/binary): foobar_POT (?)
 *
 * If the filetype is already set, then we only compare the
 * filetype with the given string, otherwise we first find
 * out about the filetype and then compare.
 *
 * Return value is the comparison result.
 */
static int abinit_filetype(abinit_plugindata_t *data, char *cmp)
{
  char const binary_header_tag[] = {'\x0e','\x00', '\x00', '\x00'};
  char lineptr[LINESIZE];

  if (!data || !cmp) return 0;

  // if filetype is already set, then only compare
  if (strlen(data->filetype) != 0) return (strncmp(data->filetype, cmp, 3) == 0);

  // read first 4 bytes to determine (un)formatted binary/ascii text
  rewind(data->file);
  fread(lineptr, 1, 4, data->file);
  if (memcmp(lineptr, binary_header_tag, 4) == 0) {
    // it is an unformatted (binary) file

    char codvsn[7];
    int headform, fform;

    // read code version (string of 6 bytes), headform (integer), and fform (integer)
    fread(codvsn, 1, 6, data->file);
    fread(&headform, 4, 1, data->file);
    fread(&fform, 4, 1, data->file);
   
    switch(fform) {
    case   2: // WFK Wave Function file
              strcpy(data->filetype, "WFK");
              break;
    case  52: // DEN Charge Density file
              strcpy(data->filetype, "DEN");
              break;
    case 102: // POT Potential file
              strcpy(data->filetype, "POT");
              break;
    default:  // Error
              break;
    }

  } else {
    // it is a formatted (ascii text) file

    // read first non-empty line of the file
    rewind(data->file);
    abinit_readline(lineptr, data->file);

    // GEO geometry file
    if (strstr(lineptr, " GEO file")) strcpy(data->filetype, "GEO");

    // regular output file
    if (strstr(lineptr, ".Version ")) strcpy(data->filetype, "out");
  }

  rewind(data->file);

  return (strncmp(data->filetype, cmp, 3) == 0);
}


static void *GEO_open_file_read(abinit_plugindata_t *data, int *natom)
{
  char lineptr[LINESIZE], atomname[NATOM_MAX][10];
  int i, idx;
DBGPRINT(stderr, "Enter GEO_open_file_read\n");

  // go to the line with the text 'XMOL data'
  while (abinit_readline(lineptr, data->file) != NULL) {
    if (strstr(lineptr, "XMOL data")) break;
  }
  if (!strstr(lineptr, "XMOL data")) {
    fprintf(stderr, "\n\nABINIT read) ERROR: '%s' has no 'XMOL data...' lines.\n", data->filename);
    abinit_plugindata_free(data);
    return NULL;
  }

  // next non-empty line has the number of atoms
  if (abinit_readline(lineptr, data->file) == NULL) {
    fprintf(stderr, "\n\nABINIT read) ERROR: cannot find the number of atoms in file '%s' .\n", data->filename);
    abinit_plugindata_free(data);
    return NULL;
  }
  data->natom = atoi(lineptr);
  if (data->natom <= 0 || data->natom > NATOM_MAX) {
    fprintf(stderr, "\n\nABINIT read) ERROR: file '%s' has %d number of atoms.\n", data->filename, data->natom);
    abinit_plugindata_free(data);
    return NULL;
  }

  // read through the list of atom names and skip their positions
  for (i = 0; i < NATOM_MAX; ++i) data->typat[i] = atomname[i][0] = '\0';
  for (idx = i = 0; i < data->natom; ++i) {
    int n;
    char name[10];
    if (1 != fscanf(data->file, "%s %*f %*f %*f", name)) {
      fprintf(stderr, "\n\nABINIT read) ERROR: file '%s' does not have the atom list.\n", data->filename);
      abinit_plugindata_free(data);
      return NULL;
    }

    // compare current atom name with previous read list and set typat accordingly
    for (n = 0; n < idx; ++n) if (strcmp(atomname[n], name) == 0) break;
    if (n == idx) strcpy(atomname[idx++], name);
    data->typat[i] = n + 1;

DBGPRINT(stderr, "   \"%s\": name = %s : data->typat[%d] = %d\n", data->filetype, atomname[n], i, data->typat[i]);
  }

  rewind(data->file);

  *natom = data->natom;

DBGPRINT(stderr, "Exit GEO_open_file_read\n");
  return data;
}


static int GEO_read_structure(abinit_plugindata_t *data, int *optflags, molfile_atom_t *atomlist)
{
  char lineptr[LINESIZE];
  int i, status;
DBGPRINT(stderr, "Enter GEO_read_structure\n");

  // go to the line with the 'XMOL data'
  do {
    char *line = abinit_readline(lineptr, data->file);
    status = line != NULL && !strstr(lineptr, "XMOL data");
  } while (status);

  // skip line with atom numbers
  abinit_readline(lineptr, data->file);

  // find atom types in XMOL list
  for (i = 0; i < data->natom; ++i) {
    molfile_atom_t *const atom = &(atomlist[i]);

    // required fields
    if (1 != fscanf(data->file, "%s %*f %*f %*f", atom->name)) {
      fprintf(stderr, "\n\nABINIT read) ERROR: file '%s' does not have the atom list.\n", data->filename);
      abinit_plugindata_free(data);
      return MOLFILE_ERROR;
    }
    strncpy(atom->type, atom->name, sizeof(atom->type));
    atom->resname[0] = '\0';
    atom->resid = 1;
    atom->segid[0]='\0';
    atom->chain[0]='\0';

    // Optional fields (defined in *optflags)
    atom->atomicnumber = get_pte_idx(atom->name);
    atom->mass = get_pte_mass(atom->atomicnumber);
    atom->radius = get_pte_vdw_radius(atom->atomicnumber);
DBGPRINT(stderr, "   atom %d : %d (%s)\n", i, atom->atomicnumber, atom->name);
  }

  // tell which of the optional fields in the molfile_atom_t structure are provided
  *optflags = MOLFILE_ATOMICNUMBER | MOLFILE_MASS | MOLFILE_RADIUS; 

  rewind(data->file);

DBGPRINT(stderr, "Exit GEO_read_structure\n");
  return MOLFILE_SUCCESS;
}


static int GEO_read_next_timestep(abinit_plugindata_t *data, int natom, molfile_timestep_t *ts)
{
  char lineptr[LINESIZE];
  float *a, *b, *c;
  int i , status;
DBGPRINT(stderr, "Enter GEO_read_next_timestep\n");

  /* At the very first call the file pointer will be non-NULL, but all
   * consecutive calls will have a NULL file pointer (see at the end of
   * this routine); in that case find the new next timestep file; if this
   * file does not exist, then there are no more timesteps.
   */
  if (!data->file) {
    /* Construct a new file name by finding the first integer number
     * in the current file name, and create the same file name but
     * with this integer number incremented.
     */
    char *newfilename = NULL, *endpart = NULL;
    for (i = strlen(data->filename) - 1; i >= 0 && !newfilename; --i) {
      if (!endpart && isdigit(data->filename[i])) endpart = strdup(data->filename + i + 1);
      if (endpart && !newfilename && !isdigit(data->filename[i])) {
         newfilename = (char *)malloc(sizeof(char) * (1 + strlen(data->filename)));
         if (!newfilename) {
           free(endpart);
           return MOLFILE_EOF;
         }

         strncpy(newfilename, data->filename, i + 1);
         sprintf(newfilename + i + 1, "%d%s", 1 + atoi(data->filename + i + 1), endpart);
      }
    }

    // if the file name does not have a number in it, then there's no further timestep
    if (!endpart) {
       free(newfilename);
       return MOLFILE_EOF;
    }

    // clean up
    free(endpart);

    // try to open the new file; if it fails, we stop the timestep calls
    data->file = fopen(newfilename, "r");
    if (!data->file) {
       free(newfilename);
       return MOLFILE_EOF;
    }

    // the new file name exists! Substitute it for the old file name.
    free(data->filename);
    data->filename = newfilename;
  }

  // go to the line with the 'Primitive vectors...'
  do {
    char *line = abinit_readline(lineptr, data->file);
    status = ( line != NULL && !strstr(lineptr, "Primitive vectors") );
  } while (status);

  // read unit cell vectors from file
  for (i = 0; i < 3; ++i) {
    float length, *r = data->rprimd[i];
    if (3 != fscanf(data->file, "%*s %f %f %f", &r[0], &r[1], &r[2])) return MOLFILE_EOF;

    // convert length units from Bohr to Angstrom
    r[0] *= BOHR_TO_ANGS;
    r[1] *= BOHR_TO_ANGS;
    r[2] *= BOHR_TO_ANGS;

    // lengths of the respective unit cell vector
    length = sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]);
    switch (i) {
      case 0: ts->A = length; break;
      case 1: ts->B = length; break;
      case 2: ts->C = length; break;
    }
  }

  // vectors of the unit cell
  a = data->rprimd[0];
  b = data->rprimd[1];
  c = data->rprimd[2];

  // alpha: angle (in degrees) between b and c
  ts->alpha = (180.0/M_PI) * acos( (b[0]*c[0] + b[1]*c[1] + b[2]*c[2]) / (ts->B*ts->C) );

  // beta: angle (in degrees) between a and c
  ts->beta  = (180.0/M_PI) * acos( (a[0]*c[0] + a[1]*c[1] + a[2]*c[2]) / (ts->A*ts->C) );

  // gamma: angle (in degrees) between between a and b
  ts->gamma = (180.0/M_PI) * acos( (a[0]*b[0] + a[1]*b[1] + a[2]*b[2]) / (ts->A*ts->B) );

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" : ""));

  // go to the line with the 'XMOL data'
  do {
    char *line = abinit_readline(lineptr, data->file);
    status = line != NULL && !strstr(lineptr, "XMOL data");
  } while (status);

  // skip line with atom numbers
  abinit_readline(lineptr, data->file);

  // read coordinates (in Angstrom) of each atom
  for (i = 0; i < data->natom; ++i) {
    float *coords = &ts->coords[3*i];
    fscanf(data->file, "%*s %f %f %f", &coords[0], &coords[1], &coords[2]);
  }

  // close the file and NULLify as preparation for the next timestep call
  fclose(data->file);
  data->file = NULL;

DBGPRINT(stderr, "Exit GEO_read_next_timestep\n");
  return MOLFILE_SUCCESS;
}


static void *DEN_open_file_read(abinit_plugindata_t *data, int *natom)
{
  int i;
DBGPRINT(stderr, "Enter DEN_open_file_read\n");

  // read all information from the header of the file
  data->hdr = abinit_header(data->file);

  if (!data->hdr) {
    abinit_plugindata_free(data);
    return NULL;
  }

  data->natom = data->hdr->natom;

  if (data->natom <= 0 || data->natom > NATOM_MAX) {
    abinit_plugindata_free(data);
    return NULL;
  }

  for (i = 0; i < data->natom; ++i) data->typat[i] = data->hdr->typat[i];

for (i = 0; i < data->natom; ++i) DBGPRINT(stderr, "   \"%s\": data->typat[%d] = %d\n", data->filetype, i, data->typat[i]);

  *natom = data->natom;

DBGPRINT(stderr, "Exit DEN_open_file_read\n");
  return data;
}


static int DEN_read_structure(abinit_plugindata_t *data, int *optflags, molfile_atom_t *atomlist)
{
  int i;
DBGPRINT(stderr, "Enter DEN_read_structure\n");

  // set the atom types, names, etc.
  for (i = 0; i < data->natom; ++i) {
    molfile_atom_t *const atom = &(atomlist[i]);
    
    // optional fields (defined in *optflags)
    atom->atomicnumber = (int)floor(0.5 + data->hdr->znucltypat[data->hdr->typat[i] - 1]);
    atom->mass = get_pte_mass(atom->atomicnumber);
    atom->radius = get_pte_vdw_radius(atom->atomicnumber);

    // required fields
    strncpy(atom->name, get_pte_label(atom->atomicnumber), sizeof(atom->name));    
    strncpy(atom->type, atom->name, sizeof(atom->type));
    atom->resname[0] = '\0';
    atom->resid = 1;
    atom->segid[0]='\0';
    atom->chain[0]='\0';
DBGPRINT(stderr, "   atom %d : %d (%s)\n", i, atom->atomicnumber, atom->name);
  }

  // tell which of the optional fields in the molfile_atom_t structure are provided
  *optflags = MOLFILE_ATOMICNUMBER | MOLFILE_MASS | MOLFILE_RADIUS; 

DBGPRINT(stderr, "Exit DEN_read_structure\n");
  return MOLFILE_SUCCESS;
}


static int DEN_read_next_timestep(abinit_plugindata_t *data, int natom, molfile_timestep_t *ts)
{
  float *a, *b, *c;
  int i;
DBGPRINT(stderr, "Enter DEN_read_next_timestep\n");

  /* At the very first call the file pointer will be non-NULL, but all
   * consecutive calls will have a NULL file pointer (see at the end of
   * this routine); in that case find the new next timestep file; if this
   * file does not exist, then there are no more timesteps.
   */
  if (!data->file) {
    // alas, volumetric data cannot be read as timesteps
    return MOLFILE_EOF;
    
    /* Construct a new file name by finding the first integer number in
     * the current file name, and create the same file name but with
     * this integer number incremented.
     */
    char *newfilename = NULL, *endpart = NULL;
    for (i = strlen(data->filename) - 1; i >= 0; --i) {
      if (!endpart && isdigit(data->filename[i])) endpart = strdup(data->filename + i + 1);
      if (endpart && !newfilename && !isdigit(data->filename[i])) {
         newfilename = (char *)malloc(sizeof(char) * (1 + strlen(data->filename)));
         if (!newfilename) {
           free(endpart);
           return MOLFILE_EOF;
         }

         strncpy(newfilename, data->filename, i + 1);
         sprintf(newfilename + i + 1, "%d%s", 1 + atoi(data->filename + i + 1), endpart);
      }
    }

    // if the file name does not have a number in it, then there's no further timestep
    if (!endpart) {
       free(newfilename);
       return MOLFILE_EOF;
    }

    // clean up
    free(endpart);

    // try to open the new file; if it fails, we stop the timestep calls
    data->file = fopen(newfilename, "r");
    if (!data->file) {
       free(newfilename);
       return MOLFILE_EOF;
    }

    // the new file name exists! Substitute it for the old file name.
    free(data->filename);
    data->filename = newfilename;

    // read the new header info from the new file
    abinit_headerdata_free(data->hdr);
    data->hdr = abinit_header(data->file);
  }

  // get unit cell vectors and convert them to Angstrom
  for (i = 0; i < 3; ++i) {
    float length;
    data->rprimd[i][0] = data->hdr->rprimd[i][0] * BOHR_TO_ANGS;
    data->rprimd[i][1] = data->hdr->rprimd[i][1] * BOHR_TO_ANGS;
    data->rprimd[i][2] = data->hdr->rprimd[i][2] * BOHR_TO_ANGS;
    length = sqrt(pow(data->rprimd[i][0], 2) + pow(data->rprimd[i][1], 2) + pow(data->rprimd[i][2], 2));
    switch (i) {
      case 0: ts->A = length; break;
      case 1: ts->B = length; break;
      case 2: ts->C = length; break;
    }    
  }

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" : ""));

  // vectors of the unit cell
  a = data->rprimd[0];
  b = data->rprimd[1];
  c = data->rprimd[2];

  // alpha: angle (in degrees) between b and c
  ts->alpha = (180.0/M_PI) * acos( (b[0]*c[0] + b[1]*c[1] + b[2]*c[2]) / (ts->B*ts->C) );

  // beta: angle (in degrees) between a and c
  ts->beta  = (180.0/M_PI) * acos( (a[0]*c[0] + a[1]*c[1] + a[2]*c[2]) / (ts->A*ts->C) );

  // gamma: angle (in degrees) between between a and b
  ts->gamma = (180.0/M_PI) * acos( (a[0]*b[0] + a[1]*b[1] + a[2]*b[2]) / (ts->A*ts->B) );

  // get coordinates (in Angstrom) of each atom
  for (i = 0; i < data->natom; ++i) {
    float *coords = &ts->coords[3*i];
    double **xred = data->hdr->xred;
    coords[0] = xred[0][i] * a[0] + xred[1][i] * b[0] + xred[2][i] * c[0];
    coords[1] = xred[0][i] * a[1] + xred[1][i] * b[1] + xred[2][i] * c[1];
    coords[2] = xred[0][i] * a[2] + xred[1][i] * b[2] + xred[2][i] * c[2];
  }

  // close the file and NULLify as preparation for the next timestep call
  fclose(data->file);
  data->file = NULL;

DBGPRINT(stderr, "Exit DEN_read_next_timestep\n");
  return MOLFILE_SUCCESS;
}


static int DEN_read_volumetric_metadata(abinit_plugindata_t *data, int *nvolsets, molfile_volumetric_t **metadata)
{
  int i;
DBGPRINT(stderr, "Enter DEN_read_volumetric_metadata\n");

  data->nvolsets = data->hdr->nspden;

  data->vol = (molfile_volumetric_t *)malloc(data->nvolsets * sizeof(molfile_volumetric_t));
  if (!data->vol) {
    fprintf(stderr, "\n\nABINIT DEN read) ERROR: Cannot allocate space for volume data.\n");
    abinit_plugindata_free(data);
    return MOLFILE_ERROR;
  }

  // get unit cell vectors and convert to Angstrom
  for (i = 0; i < 3; ++i) {
    int k;
    for (k = 0; k < 3; ++k) data->rprimd[i][k] = data->hdr->rprimd[i][k] * BOHR_TO_ANGS;
  }
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" : ""));

  for (i = 0; i < data->nvolsets; ++i) {
    /* get a handle to the current volume set meta data */
    molfile_volumetric_t *const set = &(data->vol[i]);
    int k;

    if (i == 0) {
      sprintf(set->dataname, "Total charge density");
    } else if (data->nvolsets == 2) {
      sprintf(set->dataname, "Spin-up charge density");
    } else if (data->nvolsets == 4) {
      if (i == 1) sprintf(set->dataname, "X-projection of local magnetization");
      if (i == 2) sprintf(set->dataname, "Y-projection of local magnetization");
      if (i == 3) sprintf(set->dataname, "Z-projection of local magnetization");
    }

    for (k = 0 ; k < 3; ++k) {
      set->xaxis[k] = data->rprimd[0][k];
      set->yaxis[k] = data->rprimd[1][k];
      set->zaxis[k] = data->rprimd[2][k];
    }

    /* add one more to the grid size and fill the extra voxel with the same
     * value at the beginning of the row to make the volumetric data smooth
     * across the cell boundary
     */
    set->xsize = data->hdr->ngfft[0] + 1;
    set->ysize = data->hdr->ngfft[1] + 1;
    set->zsize = data->hdr->ngfft[2] + 1;

    set->has_color = 0;
    set->origin[0] = set->origin[1] = set->origin[2] = 0;
  }

  *nvolsets = data->nvolsets;
  *metadata = data->vol;  

DBGPRINT(stderr, "Exit DEN_read_volumetric_metadata.\n");
  return MOLFILE_SUCCESS;
}


static int DEN_read_volumetric_data(abinit_plugindata_t *data, int set, float *datablock, float *colorblock)
{
  int n, iset;
DBGPRINT(stderr, "Enter DEN_read_volumetric_data\n");

  if (set >= data->nvolsets) return MOLFILE_ERROR;

  /* Abinit provides the electron density, which is the negative of the
   * charge density. For norm-conserving pseudopotentials (usepaw = 0)
   * this is the pseudo-valence electron density; for a PAW calculation
   * (usepaw = 1) this is the pseudo-valence electron density plus the
   * compensation charge density.
   * To visualize the electron density from a PAW calculation, create
   * the density file with the pawprtden key word, not prtden (however,
   * if you want to chain this density into later calculations, you have
   * to use prtden.
   */
  if (data->hdr->usepaw) fprintf(stderr, "WARNING: Be sure that you have used \"pawprtden 1\" in order to visualize the electron density!\n");

  for (n = iset = 0; iset <= set; ++iset) {
    int const xsize = data->vol[iset].xsize; 
    int const ysize = data->vol[iset].ysize;
    int const zsize = data->vol[iset].zsize;

    int ix, iy, iz;

    /* Fill the datablock with the density data.
     * Note that for each 'iset' the datablock is overwritten,
     * so that only the last datablock setting remains.
     */
    for (n = iz = 0; iz < zsize; ++iz) {
      for (iy = 0; iy < ysize; ++iy) {
        for (ix = 0; ix < xsize; ++ix, ++n) {
          double value;

          /* The datablock grid is one voxel larger in each direction in order
	   * to smoothen the density at the cell borders; hence fill the extra
	   * voxel with the same value as the one at the beginning of that row.
	   */
          if (ix == xsize - 1) value = datablock[n - ix];
	  else if (iy == ysize - 1) value = datablock[n - iy*xsize];
	  else if (iz == zsize - 1) value = datablock[n - iz*ysize*xsize];
	  else {
            // get charge density (in electrons/Bohr^3) and convert to Angstrom
            fread(&value, 1, 8, data->file);
	    value /= pow(BOHR_TO_ANGS, 3);
	  }

          datablock[n] = value;
	}
      }
    }
  }
  
DBGPRINT(stderr, "Exit DEN_read_volumetric_data\n");
  return MOLFILE_SUCCESS;
}


// VMD calls this one just once to verify access to the file
static void *open_file_read(const char *filename, const char *filetype, int *natom)
{
  void *result = NULL;
  abinit_plugindata_t *data;
DBGPRINT(stderr, "Enter open_file_read\n");

  // verify that the input variables are OK
  if (!filename || !natom) return NULL;

  // start with undefined value and set it after successful read
  *natom = MOLFILE_NUMATOMS_UNKNOWN;

  // allocate memory for the abinit data structure
  data = abinit_plugindata_malloc();
  if (!data) return NULL;

  // open the file for reading
  data->file = fopen(filename, "r");
  data->filename = strdup(filename);
  if (!data->file || !data->filename) {
    abinit_plugindata_free(data);
    return NULL;
  }

  if (abinit_filetype(data, "GEO")) result = GEO_open_file_read(data, natom);
  if (abinit_filetype(data, "DEN")) result = DEN_open_file_read(data, natom);

DBGPRINT(stderr, "Exit open_file_read\n");
  return result;
}


// VMD calls this once to find out about the atom types in the structure
static int read_structure(void *mydata, int *optflags, molfile_atom_t *atomlist)
{
  int result = MOLFILE_ERROR;
  abinit_plugindata_t *data = (abinit_plugindata_t *)mydata;
DBGPRINT(stderr, "Enter read_structure\n");
 
  if (!data || !optflags || !atomlist) return MOLFILE_ERROR;

  if (abinit_filetype(data, "GEO")) result = GEO_read_structure(data, optflags, atomlist);
  if (abinit_filetype(data, "DEN")) result = DEN_read_structure(data, optflags, atomlist);

DBGPRINT(stderr, "Exit read_structure\n");
  return result;
}


// VMD keeps calling for a next timestep, until it gets End-Of-File here
static int read_next_timestep(void *mydata, int natom, molfile_timestep_t *ts)
{
  int result = MOLFILE_EOF;
  abinit_plugindata_t *data = (abinit_plugindata_t *)mydata;
DBGPRINT(stderr, "Enter read_next_timestep\n");

  /* Save coordinatess only if we are given a timestep pointer.
   * Otherwise assume that VMD wants us to skip past it.
   */
  if (!ts || !data) return MOLFILE_EOF;

  // Double check that the number of atoms are correct
  if (natom != data->natom) return MOLFILE_EOF;

  if (abinit_filetype(data, "GEO")) result = GEO_read_next_timestep(data, natom, ts);
  if (abinit_filetype(data, "DEN")) result = DEN_read_next_timestep(data, natom, ts);

DBGPRINT(stderr, "Exit vmd_read_next_timestep\n");
  return result;
}


static void close_file_read(void *mydata)
{
DBGPRINT(stderr, "Enter vmd_close_read\n");
  abinit_plugindata_t *data = (abinit_plugindata_t *)mydata;
  abinit_plugindata_free(data);
DBGPRINT(stderr, "Exit vmd_close_read\n");
}


static void *open_file_write(const char *filename, const char *filetype, int natoms)
{
  abinit_plugindata_t *data = abinit_plugindata_malloc();

  if (!data) return NULL;

  data->filename = strdup(filename);
  data->file = fopen(filename, "w");
  if (!data->filename || !data->file) {
    abinit_plugindata_free(data);
    fprintf(stderr, "ABINIT write) ERROR: Unable to open file '%s' for writing\n", filename);
    return NULL;
  }
  fclose(data->file);

  data->natom = natoms;

  return data;
}


static int write_structure(void *mydata, int optflags, const molfile_atom_t *atoms)
{
  abinit_plugindata_t *data = (abinit_plugindata_t *)mydata;

  if (!data || !atoms) return MOLFILE_ERROR;

  data->atomlist = (molfile_atom_t *)malloc(sizeof(molfile_atom_t) * data->natom);
  if (!data->atomlist) return MOLFILE_ERROR;

  memcpy(data->atomlist, atoms, sizeof(molfile_atom_t) * data->natom);

  return MOLFILE_SUCCESS;
}


static int write_timestep(void *mydata, const molfile_timestep_t *ts)
{
  return MOLFILE_ERROR;

//  abinit_plugindata_t *data = (abinit_plugindata_t *)mydata; 

  return MOLFILE_SUCCESS;
}


static void close_file_write(void *mydata)
{
  abinit_plugindata_t *data = (abinit_plugindata_t *)mydata;
  abinit_plugindata_free(data);
}

static int read_volumetric_metadata(void *mydata, int *nvolsets, molfile_volumetric_t **metadata)
{
  int result = MOLFILE_ERROR;
  abinit_plugindata_t *data = (abinit_plugindata_t *)mydata;
DBGPRINT(stderr, "Enter vmd_read_volumetric_metadata\n");

  if (!data || !nvolsets || !metadata) return MOLFILE_ERROR;

  if (abinit_filetype(data, "DEN")) result = DEN_read_volumetric_metadata(data, nvolsets, metadata);

DBGPRINT(stderr, "Exit vmd_read_volumetric_metadata\n");
  return result;
}

static int read_volumetric_data(void *mydata, int set, float *datablock, float *colorblock)
{
  int result = MOLFILE_ERROR;
  abinit_plugindata_t *data = (abinit_plugindata_t *)mydata;
DBGPRINT(stderr, "Enter vmd_read_volumetric_data\n");

  if (!data || !datablock) return MOLFILE_ERROR;

  if (strncmp(data->filetype, "DEN", 3) == 0) result = DEN_read_volumetric_data(data, set, datablock, colorblock);

DBGPRINT(stderr, "Exit vmd_read_volumetric_data\n");
  return result;
}


/* registration stuff */
static molfile_plugin_t abinitplugin = {
  vmdplugin_ABIVERSION,     /* the ABI for the base plugin type */
  MOLFILE_PLUGIN_TYPE,      /* type */
  "ABINIT",                 /* filetype */
  "ABINIT",                 /* pretty name in filetype list */
  "Rob Lahaye",             /* author */
  0,                        /* major version */
  1,                        /* minor version */
  VMDPLUGIN_THREADSAFE,     /* is reentrant */

  "*|*_GEO|*_DEN",          /* filename_extension */
  open_file_read,           /* open_file_read */
  read_structure,           /* read_structure */
  0,                        /* read_bonds */
  read_next_timestep,       /* read_next_timestep */
  close_file_read,          /* close_file_read */
  open_file_write,          /* open_file_write */
  write_structure,          /* write_structure */
  write_timestep,           /* write_timestep */
  close_file_write,         /* close_file_write */
  read_volumetric_metadata, /* read_volumetric_metadata */
  read_volumetric_data,     /* read_volumetric_data */
  0,                        /* read_rawgraphics */
  0,                        /* read_molecule_metadata */
  0                         /* write_bonds */
};

int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {

/* This doesn't work; will have to find out later why.
  molfile_plugin_t abinitplugin;

  // see molfile_plugin.h for details on the fields
  // header 
  abinitplugin.abiversion = vmdplugin_ABIVERSION; // the ABI for the base plugin type
  abinitplugin.type       = MOLFILE_PLUGIN_TYPE; // string descriptor of the plugin type
  abinitplugin.name       = "ABINIT"; // name for the plugin
  abinitplugin.prettyname = "ABINIT"; // name in filetype list
  abinitplugin.author     = "Rob Lahaye\nSungkyunkwan University\n2009"; // string identifier
  abinitplugin.majorv     = 0; // major version
  abinitplugin.minorv     = 1; // minor version
  abinitplugin.is_reentrant = VMDPLUGIN_THREADSAFE; // can this library be run concurrently with itself?

  abinitplugin.filename_extension = "*|*_GEO|*_DEN"; // filename extension for this file type
  abinitplugin.open_file_read = open_file_read; // try to open the file for reading
  abinitplugin.read_structure = read_structure; // read molecular structure from the given file handle
  abinitplugin.read_bonds = 0; // read bond information for the molecule
  abinitplugin.read_next_timestep = read_next_timestep; // read the next timestep from the file
  abinitplugin.close_file_read = close_file_read; // close the file and release all data
  abinitplugin.open_file_write = open_file_write; // open a coordinate file for writing
  abinitplugin.write_structure = write_structure; // write a timestep to the coordinate file
  abinitplugin.write_timestep = write_timestep; // write a timestep to the coordinate file
  abinitplugin.close_file_write = close_file_write; // close the file and release all data
  abinitplugin.read_volumetric_metadata = read_volumetric_metadata; // retrieve metadata pertaining to volumetric datasets in this file
  abinitplugin.read_volumetric_data = read_volumetric_data; // read the specified volumetric data set into the space pointed to by datablock
  abinitplugin.read_rawgraphics = 0; // read raw graphics data stored in this file
  abinitplugin.read_molecule_metadata = 0; // read molecule metadata
  abinitplugin.write_bonds = 0; // write bond information for the molecule
*/  

  (*cb)(v, (vmdplugin_t *)&abinitplugin);
  return VMDPLUGIN_SUCCESS;
}

int VMDPLUGIN_init() {
  return VMDPLUGIN_SUCCESS;
}

int VMDPLUGIN_fini() {
  return VMDPLUGIN_SUCCESS;
}

/* Read the full header of a binary ABINIT file until the point
 * where the data starts. Return a pointer to the header struct,
 * or NULL if the reading fails.
 */ 
static abinit_binary_header_t *abinit_header(FILE *fp)
{
 int const debug = 0;
 int i, bc = 0;

 char skips[1024];
 int skipi;
 abinit_binary_header_t *hdr = (abinit_binary_header_t *)malloc(sizeof(abinit_binary_header_t));

 if (!hdr) return NULL;

 // be sure to start from the very beginning of the file
 rewind(fp);

 // skip first 4 bytes
 bc += fread(skips, 1, 4, fp);

 // code version
 bc += fread(hdr->codvsn, sizeof(char), 6, fp); hdr->codvsn[6] = '\0';
 if (debug) fprintf(stderr, "codvsn = '%s' (code version)\n", hdr->codvsn); 

 // format of the header
 bc += fread(&hdr->headform, 4, 1, fp);
 if (debug) fprintf(stderr, "headform = '%d' (format of the header)\n", hdr->headform);

 // specification for data type: 2 for wf; 52 for density; 102 for potential
 bc += fread(&hdr->fform, 4, 1, fp);
 if (debug) fprintf(stderr, "fform = '%d' (2 for wf; 52 for density; 102 for potential)\n", hdr->fform);

 bc += fread(&skipi, 1, 4, fp);
 if (debug) fprintf(stderr, "skip 4 bytes = '%d'\n", skipi);

 bc += fread(&skipi, 1, 4, fp);
 if (debug) fprintf(stderr, "skip 4 bytes = '%d'\n", skipi);

 // total number of bands (sum of nband on all kpts and spins)
 bc += fread(&hdr->bantot, 4, 1, fp);
 if (debug) fprintf(stderr, "bantot = '%d' (sum of nband on all kpts and spins)\n", hdr->bantot);

 // starting date
 bc += fread(&hdr->date, 4, 1, fp);
 if (debug) fprintf(stderr, "date = '%d' (starting date)\n", hdr->date);

 bc += fread(&hdr->intxc, 4, 1, fp);
 if (debug) fprintf(stderr, "intxc = '%d'\n", hdr->intxc);

 bc += fread(&hdr->ixc, 4, 1, fp);
 if (debug) fprintf(stderr, "ixc = '%d'\n", hdr->ixc);

 bc += fread(&hdr->natom, 4, 1, fp);
 if (debug) fprintf(stderr, "natom = '%d'\n", hdr->natom);

 // failsafe
 if (hdr->natom <= 0) {
   fprintf(stderr, "Binary Header Error: natom = %d; something must be wrong!",  hdr->natom);
   return NULL;
 }

 bc += fread(&hdr->ngfft[0], 4, 1, fp);
 if (debug) fprintf(stderr, "ngfft[0] = '%d'\n", hdr->ngfft[0]);
 bc += fread(&hdr->ngfft[1], 4, 1, fp);
 if (debug) fprintf(stderr, "ngfft[1] = '%d'\n", hdr->ngfft[1]);
 bc += fread(&hdr->ngfft[2], 4, 1, fp);
 if (debug) fprintf(stderr, "ngfft[2] = '%d'\n", hdr->ngfft[2]);

 bc += fread(&hdr->nkpt, 4, 1, fp);
 if (debug) fprintf(stderr, "nkpt = '%d'\n", hdr->nkpt);

 bc += fread(&hdr->nspden, 4, 1, fp);
 if (debug) fprintf(stderr, "nspden = '%d'\n", hdr->nspden);

 // failsafe
 if (hdr->nspden != 1 && hdr->nspden != 2 && hdr->nspden != 4) {
   fprintf(stderr, "Binary Header Error: nspden = %d; something must be wrong!",  hdr->nspden);
   return NULL;
 }

 bc += fread(&hdr->nspinor, 4, 1, fp);
 if (debug) fprintf(stderr, "nspinor = '%d'\n", hdr->nspinor);

 bc += fread(&hdr->nsppol, 4, 1, fp);
 if (debug) fprintf(stderr, "nsppol = '%d'\n", hdr->nsppol);

 bc += fread(&hdr->nsym, 4, 1, fp);
 if (debug) fprintf(stderr, "nsym = '%d'\n", hdr->nsym);

 bc += fread(&hdr->npsp, 4, 1, fp);
 if (debug) fprintf(stderr, "npsp = '%d'\n", hdr->npsp);

 bc += fread(&hdr->ntypat, 4, 1, fp);
 if (debug) fprintf(stderr, "ntypat = '%d'\n", hdr->ntypat);

 bc += fread(&hdr->occopt, 4, 1, fp);
 if (debug) fprintf(stderr, "occopt = '%d'\n", hdr->occopt);

 // the index of the perturbation, 0 if GS calculation
 bc += fread(&hdr->pertcase, 4, 1, fp);
 if (debug) fprintf(stderr, "pertcase = '%d' (the index of the perturbation, 0 if GS calculation)\n", hdr->pertcase);

 bc += fread(&hdr->usepaw, 4, 1, fp);
 if (debug) fprintf(stderr, "usepaw = '%d' (0=norm-conserving psps, 1=paw)\n", hdr->usepaw);

 // failsafe
 if (hdr->usepaw != 0 && hdr->usepaw != 1) {
   fprintf(stderr, "Binary Header Error: usepaw = %d; something must be wrong!",  hdr->usepaw);
   return NULL;
 }

 bc += fread(&hdr->ecut, 8, 1, fp);
 if (debug) fprintf(stderr, "ecut = '%g'\n", hdr->ecut);

 // input variable (ecut for NC psps, pawecutdg for paw)
 bc += fread(&hdr->ecutdg, 8, 1, fp);
 if (debug) fprintf(stderr, "ecutdg = '%g' (ecut for NC psps, pawecutdg for paw)\n", hdr->ecutdg);

 bc += fread(&hdr->ecutsm, 8, 1, fp);
 if (debug) fprintf(stderr, "ecutsm = '%g'\n", hdr->ecutsm);

 // ecut*dilatmx**2 (dilatmx is an input variable)
 bc += fread(&hdr->ecut_eff, 8, 1, fp);
 if (debug) fprintf(stderr, "ecut_eff = '%g' (ecut*dilatmx**2 [dilatmx is an input variable])\n", hdr->ecut_eff);

 bc += fread(&hdr->qptn[0], 8, 1, fp);
 if (debug) fprintf(stderr, "qptn[0] = '%g'\n", hdr->qptn[0]);
 bc += fread(&hdr->qptn[1], 8, 1, fp);
 if (debug) fprintf(stderr, "qptn[1] = '%g'\n", hdr->qptn[1]);
 bc += fread(&hdr->qptn[2], 8, 1, fp);
 if (debug) fprintf(stderr, "qptn[2] = '%g'\n", hdr->qptn[2]);

 for (i = 0; i < 3; ++i) {
   bc += fread(&hdr->rprimd[i][0], 8, 1, fp);
   if (debug) fprintf(stderr, "rprimd[%d][0] = '%g'\n", i, hdr->rprimd[i][0]);
   bc += fread(&hdr->rprimd[i][1], 8, 1, fp);
   if (debug) fprintf(stderr, "rprimd[%d][1] = '%g'\n", i, hdr->rprimd[i][1]);
   bc += fread(&hdr->rprimd[i][2], 8, 1, fp);
   if (debug) fprintf(stderr, "rprimd[%d][2] = '%g'\n", i, hdr->rprimd[i][2]);
 }

 bc += fread(&hdr->stmbias, 8, 1, fp);
 if (debug) fprintf(stderr, "stmbias = '%g'\n", hdr->stmbias);

 bc += fread(&hdr->tphysel, 8, 1, fp);
 if (debug) fprintf(stderr, "tphysel = '%g'\n", hdr->tphysel);

 bc += fread(&hdr->tsmear, 8, 1, fp);
 if (debug) fprintf(stderr, "tsmear = '%g'\n", hdr->tsmear);

 bc += fread(&hdr->usewvl, 1, 4, fp);
 if (debug) fprintf(stderr, "usewvl = '%d'\n", hdr->usewvl);

 // failsafe
 if (hdr->usewvl != 0 && hdr->usewvl != 1) {
   fprintf(stderr, "Binary Header Error: usewvl = %d; something must be wrong!",  hdr->usewvl);
   return NULL;
 }

 bc += fread(&skipi, 1, 4, fp);
 if (debug) fprintf(stderr, "skip 4 bytes = '%d'\n", skipi);

 bc += fread(&skipi, 1, 4, fp);
 if (debug) fprintf(stderr, "skip 4 bytes = '%d'\n", skipi);

 hdr->istwfk = (int *)malloc(sizeof(int) * hdr->nkpt);
 hdr->nband  = (int *)malloc(sizeof(int) * hdr->nkpt * hdr->nsppol);
 hdr->npwarr = (int *)malloc(sizeof(int) * hdr->nkpt);
 hdr->so_psp = (int *)malloc(sizeof(int) * hdr->npsp);
 hdr->symafm = (int *)malloc(sizeof(int) * hdr->nsym);
 hdr->typat  = (int *)malloc(sizeof(int) * hdr->natom);

 if (!hdr->istwfk || !hdr->nband || !hdr->npwarr || !hdr->so_psp || !hdr->symafm || !hdr->typat) return NULL;
 for (i = 0; i < 3; ++i) {
   int j;
   for (j = 0; j < 3; ++j) {
     hdr->symrel[i][j] = (int *)malloc(sizeof(int) * hdr->nsym);
     if (!hdr->symrel[i][j]) return NULL;
   }
 }

 for (i = 0; i < hdr->nkpt; ++i) {
   bc += fread(&hdr->istwfk[i], 4, 1, fp);
   if (debug) fprintf(stderr, "istwfk[%d] = '%d'\n", i, hdr->istwfk[i]);
 }

 for (i = 0; i < hdr->nkpt * hdr->nsppol; ++i) {
   bc += fread(&hdr->nband[i], 4, 1, fp);
   if (debug) fprintf(stderr, "nband[%d] = '%d'\n", i, hdr->nband[i]);
 }

 for (i = 0; i < hdr->nkpt; ++i) {
   bc += fread(&hdr->npwarr[i], 4, 1, fp);
   if (debug) fprintf(stderr, "npwarr[%d] = '%d'\n", i, hdr->npwarr[i]);
 }

 for (i = 0; i < hdr->npsp; ++i) {
   bc += fread(&hdr->so_psp[i], 4, 1, fp);
   if (debug) fprintf(stderr, "so_psp[%d] = '%d'\n", i, hdr->so_psp[i]);
 }

 for (i = 0; i < hdr->nsym; ++i) {
   bc += fread(&hdr->symafm[i], 4, 1, fp);
   if (debug) fprintf(stderr, "symafm[%d] = '%d'\n", i, hdr->symafm[i]);
 }

 for (i = 0; i < hdr->nsym; ++i) {
   int j;
   for (j = 0; j < 3; ++j) {
     int k;     
     for (k = 0; k < 3; ++k) {
       bc += fread(&hdr->symrel[k][j][i], 4, 1, fp);
       if (debug) fprintf(stderr, "symrel[%d][%d][%2d]= '%2d'", k, j, i, hdr->symrel[k][j][i]);
     }
     if (debug) fprintf(stderr, "\n");
   }
   if (debug) fprintf(stderr, "\n");
 }

 for (i = 0; i < hdr->natom; ++i) {
   bc += fread(&hdr->typat[i], 4, 1, fp);
   if (debug) fprintf(stderr, "typat[%d] = '%d'\n", i, hdr->typat[i]);
 }

 for (i = 0; i < 3; ++i) {
   hdr->kpt[i] = (double *)malloc(sizeof(double) * hdr->nkpt);
   hdr->tnons[i] = (double *)malloc(sizeof(double) * hdr->nsym);
   if (!hdr->kpt[i] || !hdr->tnons[i]) return NULL;
 }

 hdr->occ = (double *)malloc(sizeof(double) * hdr->bantot);
 hdr->znucltypat = (double *)malloc(sizeof(double) * hdr->ntypat);
 hdr->wtk = (double *)malloc(sizeof(double) * hdr->nkpt);
 if (!hdr->occ || !hdr->znucltypat || !hdr->wtk) return NULL;

 for (i = 0; i < hdr->nkpt; ++i) {
   int j;
   for (j = 0; j < 3; ++j) {
     bc += fread(&hdr->kpt[j][i], 8, 1, fp);
     if (debug) fprintf(stderr, "kpt[%d][%2d] = '%g'\n", j, i, hdr->kpt[j][i]);
   }
 }

 for (i = 0; i < hdr->bantot; ++i) {
   bc += fread(&hdr->occ[i], 8, 1, fp);
   if (debug) fprintf(stderr, "occ[%d] = '%g'\n", i, hdr->occ[i]);
 }

 for (i = 0; i < hdr->nsym; ++i) {
   int j;
   for (j = 0; j < 3; ++j) {
     bc += fread(&hdr->tnons[j][i], 8, 1, fp);
     if (debug) fprintf(stderr, "tnons[%d][%2d] = '%g'", j, i, hdr->tnons[j][i]);
   }
   if (debug) fprintf(stderr, "\n");
 }

 for (i = 0; i < hdr->ntypat; ++i) {
   bc += fread(&hdr->znucltypat[i], 8, 1, fp);
   if (debug) fprintf(stderr, "znucltypat[%d] = '%g'\n", i, hdr->znucltypat[i]);
 }

 for (i = 0; i < hdr->nkpt; ++i) {
   bc += fread(&hdr->wtk[i], 8, 1, fp);
   if (debug) fprintf(stderr, "wtk[%d] = '%g'\n", i, hdr->wtk[i]);
 }

 for (i = 0; i < hdr->npsp; ++i) {

   bc += fread(&skipi, 1, 4, fp);
   if (debug) fprintf(stderr, "skip 4 bytes = '%d'\n", skipi);

   bc += fread(&skipi, 1, 4, fp);
   if (debug) fprintf(stderr, "skip 4 bytes = '%d'\n", skipi);

   bc += fread(hdr->title, 132, 1, fp); hdr->title[132] = '\0';
   if (debug) fprintf(stderr, "title = '%s'\n", hdr->title);

   bc += fread(&hdr->znuclpsp, 1, 8, fp);
   if (debug) fprintf(stderr, "znuclpsp = '%g'\n", hdr->znuclpsp);

   bc += fread(&hdr->zionpsp, 1, 8, fp);
   if (debug) fprintf(stderr, "zionpsp = '%g'\n", hdr->zionpsp);

   bc += fread(&hdr->pspso, 1, 4, fp);
   if (debug) fprintf(stderr, "pspso = '%d'\n", hdr->pspso);

   bc += fread(&hdr->pspdat, 1, 4, fp);
   if (debug) fprintf(stderr, "pspdat = '%d'\n", hdr->pspdat);

   bc += fread(&hdr->pspcod, 1, 4, fp);
   if (debug) fprintf(stderr, "pspcod = '%d'\n", hdr->pspcod);

   bc += fread(&hdr->pspxc, 1, 4, fp);
   if (debug) fprintf(stderr, "pspxc = '%d'\n", hdr->pspxc);

   bc += fread(&hdr->lmn_size, 1, 4, fp);
   if (debug) fprintf(stderr, "lmn_size = '%d'\n", hdr->lmn_size);
 }

 bc += fread(&skipi, 1, 4, fp);
 if (debug) fprintf(stderr, "skip = '%d'\n", skipi);

 bc += fread(&skipi, 1, 4, fp);
 if (debug) fprintf(stderr, "skip = '%d'\n", skipi);

 for (i = 0; i < 3; ++i) {
   hdr->xred[i] = (double *)malloc(sizeof(double) * hdr->natom);
   if (!hdr->xred[i]) return NULL;
 }

 bc += fread(&hdr->residm, 1, 8, fp);
 if (debug) fprintf(stderr, "residm = '%g'\n", hdr->residm);

 for (i = 0; i < hdr->natom; ++i) {
   int j;
   for (j = 0; j < 3; ++j) {
     bc += fread(&hdr->xred[j][i], 1, 8, fp);
     if (debug) fprintf(stderr, "xred[%d][%d] = '%g'\n", j, i, hdr->xred[j][i]);

     // failsafe
     if (hdr->xred[j][i] < -1 || hdr->xred[j][i] > 1) {
       fprintf(stderr, "Binary Header Error: hdr->xred[%d][%d] = %g; something must be wrong!", j, i, hdr->xred[j][i]);
       return NULL;
     }
   }
 }
 
 bc += fread(&hdr->etotal, 1, 8, fp);
 if (debug) fprintf(stderr, "etotal = '%g'\n", hdr->etotal);

 bc += fread(&hdr->fermie, 1, 8, fp);
 if (debug) fprintf(stderr, "fermie = '%g'\n", hdr->fermie);

 if (hdr->usepaw == 1) {
   struct PAWRHOIJ {
     int *nrhoijsel;
     int **rhoijselect;
     double **rhoij;
   } pawrhoij[hdr->natom];

   for (i = 0; i < hdr->natom; ++i) {
     int j;
     pawrhoij[i].nrhoijsel = (int *)malloc(sizeof(int) * hdr->nspden);
     if (!pawrhoij[i].nrhoijsel) return NULL;
     for (j = 0; j < hdr->nspden; ++j){
       bc += fread(&pawrhoij[i].nrhoijsel[j], 4, 1, fp);
       if (debug) fprintf(stderr, "pawrhoij[%d].nrhoijsel[%d] = '%d'\n", i, j, pawrhoij[i].nrhoijsel[j]);
     }
     pawrhoij[i].rhoijselect = (int **)malloc(sizeof(int) * hdr->nspden);
     pawrhoij[i].rhoij = (double **)malloc(sizeof(double) * hdr->nspden);
     if (!pawrhoij[i].rhoijselect || !pawrhoij[i].rhoij) return NULL;
     for (j = 0; j < hdr->nspden; ++j) {
       int k;
       pawrhoij[i].rhoijselect[j] = (int *)malloc(sizeof(int) * pawrhoij[i].nrhoijsel[j]);
       pawrhoij[i].rhoij[j] = (double *)malloc(sizeof(double) * pawrhoij[i].nrhoijsel[j]);
       if (!pawrhoij[i].rhoijselect[j]) return NULL;
       for (k = 0; k < pawrhoij[i].nrhoijsel[j]; ++k) {
         bc += fread(&pawrhoij[i].rhoijselect[j][k], 4, 1, fp);
	 if (debug) fprintf(stderr, "pawrhoij[%d].rhoijselect[%d][%d] = '%d'\n", i, j, k, pawrhoij[i].rhoijselect[j][k]);
       }
       for (k = 0; k < pawrhoij[i].nrhoijsel[j]; ++k) {
         bc += fread(&pawrhoij[i].rhoij[j][k], 8, 1, fp);
	 if (debug) fprintf(stderr, "pawrhoij[%d].rhoij[%d][%d] = '%g'\n", i, j, k, pawrhoij[i].rhoij[j][k]);
       }
     }
     for (j = 0; j < hdr->nspden; ++j) {
       free(pawrhoij[i].rhoijselect[j]);
       free(pawrhoij[i].rhoij[j]);
     }
     free(pawrhoij[i].rhoijselect);
     free(pawrhoij[i].rhoij);
     free(pawrhoij[i].nrhoijsel);
   }
 } // end of "if (usepaw == 1)"

 bc += fread(&skipi, 1, 4, fp);
 if (debug) fprintf(stderr, "skip 4 bytes = '%d'\n", skipi);

 bc += fread(&skipi, 1, 4, fp);
 if (debug) fprintf(stderr, "skip 4 bytes = '%d'\n", skipi);

 return hdr;
}