pbeqplugin.C Source File

00001 /***************************************************************************
00002  *cr
00003  *cr            (C) Copyright 1995-2016 The Board of Trustees of the
00004  *cr                        University of Illinois
00005  *cr                         All Rights Reserved
00006  *cr
00007  ***************************************************************************/
00008 
00009 /***************************************************************************
00010  * RCS INFORMATION:
00011  *
00012  *      $RCSfile: pbeqplugin.C,v $
00013  *      $Author: johns $       $Locker:  $             $State: Exp $
00014  *      $Revision: 1.7 $       $Date: 2016/11/28 05:01:54 $
00015  *
00016  ***************************************************************************/
00017 
00018 /* 
00019  * "unformatted" binary potential map, created by CHARMM PBEQ module
00020  *
00021  * Format (fortran): 
00022  *      INTEGER UNIT
00023  *      INTEGER NCLX,NCLY,NCLZ
00024  *      REAL*8  DCEL
00025  *      REAL*8  XBCEN,YBCEN,ZBCEN,EPSW,EPSP,CONC,TMEMB,ZMEMB,EPSM
00026  *      REAL*4  PHI(*)
00027  *C
00028  *C     Local variable
00029  *      INTEGER I
00030  *
00031  *      WRITE(UNIT) NCLX,NCLY,NCLZ,DCEL,XBCEN,YBCEN,ZBCEN
00032  *      WRITE(UNIT) EPSW,EPSP,CONC,TMEMB,ZMEMB,EPSM
00033  *      WRITE(UNIT)(PHI(I),I=1,NCLX*NCLY*NCLZ)
00034  *
00035  */
00036 
00037 #include <stdlib.h>
00038 #include <stdio.h>
00039 #include <ctype.h>
00040 #include <math.h>
00041 #include <string.h>
00042 
00043 #if defined(_AIX)
00044 #include <strings.h>
00045 #endif
00046 
00047 #include "molfile_plugin.h"
00048 #include "endianswap.h"
00049 
00050 typedef struct {
00051   FILE *fd;
00052   int nsets;
00053   int ndata;
00054   int nclx;
00055   int ncly;
00056   int nclz;
00057   int swap;
00058   molfile_volumetric_t *vol;
00059 } pbeq_t;
00060 
00061 
00062 static void *open_pbeq_read(const char *filepath, const char *filetype,
00063     int *natoms) {
00064   FILE *fd;
00065   pbeq_t *pbeq;
00066   int nclx, ncly, nclz;
00067   int trash, length;
00068   double dcel; 
00069   double xbcen, ybcen, zbcen;
00070   double epsw, epsp, conc, tmemb, zmemb, epsm;
00071   int swap=0; 
00072 
00073   fd = fopen(filepath, "rb");
00074   if (!fd) {
00075     printf("pbeqplugin) Error opening file %s.\n", filepath);
00076     return NULL;
00077   }
00078 
00079   // skip first Fortran length record for
00080   // WRITE(UNIT) NCLX,NCLY,NCLZ,DCEL,XBCEN,YBCEN,ZBCEN
00081   if (fread(&length, 4, 1, fd) != 1)
00082     return NULL;
00083 
00084   if (fread(&nclx, 4, 1, fd) != 1)
00085     return NULL;
00086   if (fread(&ncly, 4, 1, fd) != 1)
00087     return NULL;
00088   if (fread(&nclz, 4, 1, fd) != 1)
00089     return NULL;
00090 
00091   // test endianness first
00092   if (length != 44) {
00093     swap = 1;
00094     swap4_aligned(&length, 1);
00095     if (length != 44) {
00096       printf("pbeqplugin) length record != 44, unrecognized format (length: %d)\n", length);
00097       return NULL;
00098     }
00099   }
00100 
00101   if (swap) {
00102     swap4_aligned(&nclx, 1);
00103     swap4_aligned(&ncly, 1);
00104     swap4_aligned(&nclz, 1);
00105   } 
00106 
00107   // this is a risky strategy for detecting endianness, 
00108   // but it works so far, and the charmm potential maps don't
00109   // have version numbers, magic numbers, or anything else that we
00110   // might otherwise use for this purpose.
00111   if ((nclx > 4000 && ncly > 4000 && nclz > 4000) ||
00112       (nclx * ncly * nclz < 0)) {
00113     printf("pbeqplugin) inconclusive byte ordering, bailing out\n");
00114     return NULL;
00115   }
00116 
00117   // read the rest of the header
00118   if (fread(&dcel, 8, 1, fd) != 1) 
00119     return NULL;
00120   if (fread(&xbcen, 8, 1, fd) != 1) 
00121     return NULL;
00122   if (fread(&ybcen, 8, 1, fd) != 1) 
00123     return NULL;
00124   if (fread(&zbcen, 8, 1, fd) != 1) 
00125     return NULL;
00126 
00127   // skip second Fortran length record for
00128   // WRITE(UNIT) NCLX,NCLY,NCLZ,DCEL,XBCEN,YBCEN,ZBCEN
00129   if (fread(&trash, 4, 1, fd) != 1)
00130     return NULL;
00131 
00132   // skip first Fortran length record for 
00133   // WRITE(UNIT) EPSW,EPSP,CONC,TMEMB,ZMEMB,EPSM
00134   if (fread(&trash, 4, 1, fd) != 1)
00135     return NULL;
00136 
00137   if (fread(&epsw, 8, 1, fd) != 1) 
00138     return NULL;
00139   if (fread(&epsp, 8, 1, fd) != 1) 
00140     return NULL;
00141   if (fread(&conc, 8, 1, fd) != 1) 
00142     return NULL;
00143   if (fread(&tmemb, 8, 1, fd) != 1) 
00144     return NULL;
00145   if (fread(&zmemb, 8, 1, fd) != 1) 
00146     return NULL;
00147   if (fread(&epsm, 8, 1, fd) != 1) 
00148     return NULL;
00149 
00150   // skip second Fortran length record for 
00151   // WRITE(UNIT) EPSW,EPSP,CONC,TMEMB,ZMEMB,EPSM
00152   if (fread(&trash, 4, 1, fd) != 1)
00153     return NULL;
00154 
00155   // byte swap the header data if necessary
00156   if (swap) {
00157     swap8_aligned(&dcel, 1);
00158     swap8_aligned(&xbcen, 1);
00159     swap8_aligned(&ybcen, 1);
00160     swap8_aligned(&zbcen, 1);
00161     swap8_aligned(&epsw, 1);
00162     swap8_aligned(&epsp, 1);
00163     swap8_aligned(&conc, 1);
00164     swap8_aligned(&tmemb, 1);
00165     swap8_aligned(&zmemb, 1);
00166     swap8_aligned(&epsm, 1);
00167   }
00168 
00169 #if 0
00170   // print header info for debugging of early versions
00171   printf("pbeqplugin) nclx:%d nxly:%d nclz:%d\n", nclx, ncly, nclz);
00172   printf("pbeqplugin) dcel: %f\n", dcel); 
00173   printf("pbeqplugin) x/y/zbcen: %g %g %g\n", xbcen, ybcen, zbcen); 
00174   printf("pbeqplugin) epsw/p: %g %g  conc: %g  tmemb: %g zmemb: %g  epsm: %g\n",
00175          epsw, epsp, conc, tmemb, zmemb, epsm);
00176 #endif
00177 
00178   /* Allocate and initialize the pbeq structure */
00179   pbeq = new pbeq_t;
00180   pbeq->fd = fd;
00181   pbeq->vol = NULL;
00182   *natoms = MOLFILE_NUMATOMS_NONE;
00183   pbeq->nsets = 1; /* this file contains only one data set */
00184   pbeq->ndata = nclx*ncly*nclz;
00185   pbeq->nclx = nclx;
00186   pbeq->ncly = ncly;
00187   pbeq->nclz = nclz;
00188   pbeq->swap = swap;
00189 
00190   pbeq->vol = new molfile_volumetric_t[1];
00191   strcpy(pbeq->vol[0].dataname, "CHARMM PBEQ Potential Map");
00192 
00193   pbeq->vol[0].origin[0] = -0.5*((nclx-1) * dcel) + xbcen;
00194   pbeq->vol[0].origin[1] = -0.5*((ncly-1) * dcel) + ybcen;
00195   pbeq->vol[0].origin[2] = -0.5*((nclz-1) * dcel) + zbcen;
00196 
00197   // print origin info, for debuggin of early versions
00198   printf("pbeqplugin) box LL origin: %g %g %g\n", 
00199          pbeq->vol[0].origin[0],
00200          pbeq->vol[0].origin[1],
00201          pbeq->vol[0].origin[2]);
00202 
00203   pbeq->vol[0].xaxis[0] = (nclx-1) * dcel;
00204   pbeq->vol[0].xaxis[1] = 0;
00205   pbeq->vol[0].xaxis[2] = 0;
00206 
00207   pbeq->vol[0].yaxis[0] = 0;
00208   pbeq->vol[0].yaxis[1] = (ncly-1) * dcel;
00209   pbeq->vol[0].yaxis[2] = 0;
00210 
00211   pbeq->vol[0].zaxis[0] = 0;
00212   pbeq->vol[0].zaxis[1] = 0;
00213   pbeq->vol[0].zaxis[2] = (nclz-1) * dcel;
00214 
00215   pbeq->vol[0].xsize = nclx;
00216   pbeq->vol[0].ysize = ncly;
00217   pbeq->vol[0].zsize = nclz;
00218 
00219   pbeq->vol[0].has_color = 0;
00220 
00221   return pbeq;
00222 }
00223 
00224 static int read_pbeq_metadata(void *v, int *nsets, 
00225   molfile_volumetric_t **metadata) {
00226   pbeq_t *pbeq = (pbeq_t *)v;
00227   *nsets = pbeq->nsets; 
00228   *metadata = pbeq->vol;  
00229 
00230   return MOLFILE_SUCCESS;
00231 }
00232 
00233 static int read_pbeq_data(void *v, int set, float *datablock,
00234                          float *colorblock) {
00235   pbeq_t *pbeq = (pbeq_t *)v;
00236   int ndata = pbeq->ndata;
00237   int nclx = pbeq->nclx;
00238   int ncly = pbeq->ncly;
00239   int nclz = pbeq->nclz;
00240   FILE *fd = pbeq->fd;
00241   int trash;
00242 
00243   // skip first Fortran length record for 
00244   // WRITE(UNIT)(PHI(I),I=1,NCLX*NCLY*NCLZ)
00245   if (fread(&trash, 4, 1, fd) != 1)
00246     return MOLFILE_ERROR;
00247 
00248   /* Read the densities. Order for file is z fast, y medium, x slow */
00249   int x, y, z;
00250   int count=0;
00251   for (x=0; x<nclx; x++) {
00252     for (y=0; y<ncly; y++) {
00253       for (z=0; z<nclz; z++) {
00254         int addr = z*nclx*ncly + y*nclx + x;
00255         if (fread(datablock + addr, 4, 1, fd) != 1) {
00256           printf("pbeqplugin) Error reading potential map cell: %d,%d,%d\n", x, y, z);
00257           printf("pbeqplugin) offset: %d\n", (int) ftell(fd));
00258           return MOLFILE_ERROR;
00259         }
00260         count++;
00261       }
00262     }
00263   }
00264 
00265 
00266 #if 0
00267   // skip last Fortran length record for 
00268   // WRITE(UNIT)(PHI(I),I=1,NCLX*NCLY*NCLZ)
00269   if (fread(&trash, 4, 1, fd) != 1)
00270     return NULL;
00271 
00272   // print debugging info for early versions
00273   printf("pbeqplugin) read %d phi values from block\n", count);
00274   for (x=0; x<1000000; x++) {
00275     int trash;
00276     if (fread(&trash, 4, 1, fd) != 1) {
00277       printf("pbeqplugin) read %d extra phi values past the end of the block\n", x);
00278       break;
00279     }
00280   }   
00281 #endif
00282 
00283   if (pbeq->swap) {
00284     swap4_aligned(datablock, ndata);
00285   }
00286 
00287   return MOLFILE_SUCCESS;
00288 }
00289 
00290 static void close_pbeq_read(void *v) {
00291   pbeq_t *pbeq = (pbeq_t *)v;
00292 
00293   fclose(pbeq->fd);
00294   if (pbeq->vol != NULL)
00295     delete [] pbeq->vol; 
00296   delete pbeq;
00297 }
00298 
00299 /*
00300  * Initialization stuff here
00301  */
00302 static molfile_plugin_t plugin;
00303 
00304 VMDPLUGIN_API int VMDPLUGIN_init(void) { 
00305   memset(&plugin, 0, sizeof(molfile_plugin_t));
00306   plugin.abiversion = vmdplugin_ABIVERSION;
00307   plugin.type = MOLFILE_PLUGIN_TYPE;
00308   plugin.name = "pbeq";
00309   plugin.prettyname = "CHARMM PBEQ Binary Potential Map";
00310   plugin.author = "John Stone";
00311   plugin.majorv = 0;
00312   plugin.minorv = 4;
00313   plugin.is_reentrant = VMDPLUGIN_THREADSAFE;
00314   plugin.filename_extension = "pbeq, phi80";
00315   plugin.open_file_read = open_pbeq_read;
00316   plugin.read_volumetric_metadata = read_pbeq_metadata;
00317   plugin.read_volumetric_data = read_pbeq_data;
00318   plugin.close_file_read = close_pbeq_read;
00319   return VMDPLUGIN_SUCCESS; 
00320 }
00321 
00322 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
00323   (*cb)(v, (vmdplugin_t *)&plugin);
00324   return VMDPLUGIN_SUCCESS;
00325 }
00326 
00327 VMDPLUGIN_API int VMDPLUGIN_fini(void) { return VMDPLUGIN_SUCCESS; }
00328