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

pqrplugin.c

Go to the documentation of this 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: pqrplugin.c,v $
00013  *      $Author: johns $       $Locker:  $             $State: Exp $
00014  *      $Revision: 1.21 $       $Date: 2016/11/28 05:01:54 $
00015  *
00016  ***************************************************************************/
00017 
00018 /*
00019  * PQR is be a space-delimited variant of the PDB format, with atom radius
00020  * and charge information replacing the B-factor and occupancy columns, and
00021  * containing only ATOM records.
00022  *
00023  * Most of this code is lifted from pdbplugin and readpdb.h
00024  */
00025 
00026 #include "molfile_plugin.h"
00027 #include <stdio.h>
00028 #include <stdlib.h>
00029 #include <string.h>
00030 
00031 #define PQR_RECORD_LENGTH 80
00032 
00033 /*  record type defines */
00034 enum {PQR_REMARK, PQR_ATOM, PQR_UNKNOWN, PQR_END, PQR_EOF, PQR_ERROR, PQR_CRYST1};
00035 
00036 /* 
00037  * read the next record from the specified pqr file, and put the string
00038  * found in the given string pointer (the caller must provide adequate (81
00039  * chars) buffer space); return the type of record found
00040  */
00041 static int read_pqr_record(FILE *f, char *retStr) {
00042 
00043   char inbuf[PQR_RECORD_LENGTH+2];
00044   int recType = PQR_UNKNOWN;
00045  
00046   /* XXX This PQR record reading code breaks with files that use
00047    * Mac or DOS style line breaks with ctrl-M characters.  We need
00048    * to replace the use of fgets() and comparisons against \n with
00049    * code that properly handles the other cases.
00050    */
00051  
00052   /*  read the next line  */
00053   if(inbuf != fgets(inbuf, PQR_RECORD_LENGTH+1, f)) {
00054     retStr[0] = '\0';
00055     if (feof(f)) {
00056       recType = PQR_EOF;
00057     }
00058     else {
00059       recType = PQR_ERROR;
00060     }
00061   } 
00062   else {
00063     /*  remove the newline character, if there is one */
00064     if(inbuf[strlen(inbuf)-1] == '\n')
00065       inbuf[strlen(inbuf)-1] = '\0';
00066 
00067     if (!strncmp(inbuf, "ATOM ",  5) || !strncmp(inbuf, "HETATM", 6)) {
00068       /* Note that by only comparing 5 chars for "ATOM " rather than 6,     */
00069       /* we allow PQR files containing > 99,999 atoms generated by AMBER    */
00070       /* to load which would otherwise fail.  Not needed for HETATM since   */
00071       /* those aren't going to show up in files produced for/by MD engines. */
00072       strcpy(retStr,inbuf);
00073       recType = PQR_ATOM;
00074     } else if (!strncmp(inbuf, "CRYST1", 6)) {
00075       recType = PQR_CRYST1;
00076       strcpy(retStr,inbuf);
00077     } else if (!strncmp(inbuf, "END", 3)) {
00078       strcpy(retStr,inbuf);
00079       recType = PQR_END;
00080     } else {
00081       retStr[0] = '\0';
00082       recType = PQR_UNKNOWN;
00083     }
00084   }
00085 
00086   /* read the '\r', if there was one */
00087   {
00088     int ch = fgetc(f);
00089     if (ch != '\r') {
00090       ungetc(ch, f);
00091     }
00092   }
00093   
00094   return recType;
00095 }
00096 
00097 /* Extract the alpha/beta/gamma a/b/c unit cell info from a CRYST1 record */
00098 static void get_pqr_cryst1(const char *record, 
00099                            float *alpha, float *beta, float *gamma, 
00100                            float *a, float *b, float *c) {
00101   char tmp[PQR_RECORD_LENGTH+3]; /* space for line + cr + lf + NUL */
00102   char ch, *s;
00103   memset(tmp, 0, sizeof(tmp));
00104   strncpy(tmp, record, PQR_RECORD_LENGTH);
00105 
00106   s = tmp+6 ;          ch = tmp[15]; tmp[15] = 0;
00107   *a = (float) atof(s);
00108   s = tmp+15; *s = ch; ch = tmp[24]; tmp[24] = 0;
00109   *b = (float) atof(s);
00110   s = tmp+24; *s = ch; ch = tmp[33]; tmp[33] = 0;
00111   *c = (float) atof(s);
00112   s = tmp+33; *s = ch; ch = tmp[40]; tmp[40] = 0;
00113   *alpha = (float) atof(s);
00114   s = tmp+40; *s = ch; ch = tmp[47]; tmp[47] = 0;
00115   *beta = (float) atof(s);
00116   s = tmp+47; *s = ch; ch = tmp[54]; tmp[54] = 0;
00117   *gamma = (float) atof(s);
00118 }
00119 
00120 /*
00121  * Read the fields in an ATOM line. Return 0 on success, nonzero otherwise.
00122  */
00123 static int get_pqr_fields(char *record, char *name, char *resname, char *chain,
00124     char *segname, char *resid,
00125     float *x, float *y, float *z, float *charge, float *radius) {
00126 
00127   /*
00128    * PQR files are like pdb files, but not quite.  They're fixed format
00129    * for the the identifier data up to resid, but after that it's free-format
00130    * for the floating point data.
00131    */
00132   strncpy(name, record+12, 4);    name[4] = 0;
00133   strncpy(resname, record+17, 4); resname[4] = 0;
00134   strncpy(chain, record+21, 1);   chain[1] = 0;
00135   strncpy(resid, record+22, 4);   resid[4] = 0;
00136   /* XXX what to do about the segid? */
00137   segname[0] = 0;
00138 
00139   if (sscanf(record+30, "%f%f%f%f%f", x, y, z, charge, radius) != 5) {
00140     return 1;
00141   }
00142   /* success */
00143   return 0;
00144 }
00145 
00146 /*
00147  * Read the coordinates in an ATOM line. Return 0 on success, nonzero otherwise.
00148  */
00149 static int get_pqr_coordinates(char *record, float *x, float *y, float *z) {
00150   /* skip over the fixed format part and only parse the coordinates */
00151   if (sscanf(record+30, "%f%f%f", x, y, z) != 3) {
00152     return 1;
00153   } 
00154   return 0;
00155 }
00156 
00157 /*
00158  * Return a default radius for the given atom.
00159  */
00160 static float get_atom_radius(molfile_atom_t *atom) {
00161   char *name, *type;
00162   name = atom->name;
00163   type = atom->type;
00164 
00165   /* XXX - yeah, this needs to be fixed */
00166   return 1.0;
00167 }
00168 
00169 /*
00170  * Append a PQR ATOM record to the given file.
00171  */
00172 static int write_pqr_atom(FILE *fd, int index, const molfile_atom_t *atom, 
00173                    float x, float y, float z) {
00174   int rc;
00175 
00176   rc = fprintf(fd, "ATOM  %5d %-4s %s %5d    %8.3f %8.3f %8.3f %.3f %.3f\n",
00177                index, atom->name, atom->resname, atom->resid, 
00178                x, y, z, atom->charge, atom->radius);
00179 
00180   return (rc > 0);
00181 }
00182 
00183 
00184 /*
00185  * API functions start here
00186  */
00187 
00188 typedef struct {
00189   FILE *fd;
00190   int natoms;
00191   molfile_atom_t *atomlist;
00192 } pqrdata;
00193 
00194 static void *open_pqr_read(const char *filepath, const char *filetype, 
00195     int *natoms) {
00196   FILE *fd;
00197   pqrdata *pqr;
00198   char pqr_record[PQR_RECORD_LENGTH+2];
00199   int record_type;
00200 
00201   fd = fopen(filepath, "r");
00202   if (!fd) 
00203     return NULL;
00204   pqr = (pqrdata *)malloc(sizeof(pqrdata));
00205   pqr->fd = fd;
00206 
00207   *natoms=0;
00208   do {
00209     record_type = read_pqr_record(pqr->fd, pqr_record);
00210     if (record_type == PQR_ATOM) {
00211       *natoms += 1;
00212     } else if (record_type == PQR_ERROR) {
00213       printf("pqrplugin) Error reading file after opening.\n");
00214       free(pqr);
00215       return NULL;
00216     }
00217   } while ((record_type != PQR_EOF) && (record_type != PQR_END));
00218 
00219   /* If no atoms were found, this is probably not a PQR file! */
00220   if (!*natoms) {
00221     printf("pqrplugin) file '%s' contains no atoms.\n", filepath);
00222     free(pqr);
00223     return NULL;
00224   }
00225 
00226   rewind(pqr->fd);
00227   pqr->natoms = *natoms;
00228   return pqr; 
00229 }
00230 
00231 /* 
00232  * Read atom structure, but not coordinate information.
00233  */
00234 static int read_pqr_structure(void *mydata, int *optflags, 
00235     molfile_atom_t *atoms) { 
00236   pqrdata *pqr = (pqrdata *)mydata;
00237   molfile_atom_t *atom;
00238   char pqr_record[PQR_RECORD_LENGTH + 1];
00239   int i, record_type;
00240   char ridstr[8];
00241   float newpos;
00242   long fpos = ftell(pqr->fd);
00243  
00244   *optflags = MOLFILE_CHARGE | MOLFILE_RADIUS;
00245 
00246   i = 0;
00247   do {
00248     record_type = read_pqr_record(pqr->fd, pqr_record);
00249     if (((record_type == PQR_EOF) || (record_type == PQR_END)) 
00250         && (i < pqr->natoms)) {
00251       printf("pqrplugin) unexpected end-of-file while reading structure.\n");
00252 printf("XXX: %d of %d \n", i, pqr->natoms);
00253       return MOLFILE_ERROR;
00254     } else if (record_type == PQR_ERROR) {
00255       printf("pqrplugin) error reading atom coordinates.\n");
00256       return MOLFILE_ERROR;
00257     } else if (record_type == PQR_ATOM) {
00258       if (i >= pqr->natoms) {
00259         printf("pqrplugin) too many atoms.\n");
00260         return MOLFILE_ERROR;
00261       }
00262       atom = atoms+i;
00263       get_pqr_fields(pqr_record, atom->name, atom->resname, atom->chain, 
00264           atom->segid, ridstr,
00265           &newpos, &newpos, &newpos, &atom->charge, &atom->radius);
00266       atom->resid = atoi(ridstr);
00267       strcpy(atom->type, atom->name);
00268       i++;
00269     }
00270   } while((record_type != PQR_EOF) && (record_type != PQR_END));
00271 
00272   fseek(pqr->fd, fpos, SEEK_SET);
00273 
00274   return MOLFILE_SUCCESS;
00275 }
00276 
00277 /* 
00278  * Read atom coordinates. PQR files contain only a single "timestep".
00279  */
00280 static int read_pqr_timestep(void *v, int natoms, molfile_timestep_t *ts) {
00281   pqrdata *pqr = (pqrdata *)v;
00282   char pqr_record[PQR_RECORD_LENGTH+2];
00283   int record_type, i;
00284   float *x, *y, *z;
00285 
00286   if (pqr->natoms == 0) 
00287     return MOLFILE_ERROR; /* EOF */
00288 
00289   if (ts != NULL) {
00290     x = ts->coords;
00291     y = x+1;
00292     z = x+2;
00293   } else {
00294     x = y = z = NULL;
00295   }
00296 
00297   i = 0;
00298   do {
00299     record_type = read_pqr_record(pqr->fd, pqr_record);
00300     if (((record_type == PQR_EOF) || (record_type == PQR_END)) 
00301         && (i < pqr->natoms)) {
00302       if (i == 0) {
00303         /* don't emit an error if there're no more timesteps */
00304         return MOLFILE_EOF;
00305       } else {
00306         printf("pqrplugin) unexpected end-of-file while reading timestep.\n");
00307         return MOLFILE_ERROR;
00308       }
00309     } else if (record_type == PQR_ERROR) {
00310       printf("pqrplugin) error reading atom coordinates.\n");
00311       return MOLFILE_ERROR;
00312     } else if (record_type == PQR_CRYST1) {
00313       if (ts) {
00314         get_pqr_cryst1(pqr_record, &ts->alpha, &ts->beta, &ts->gamma,
00315                                &ts->A, &ts->B, &ts->C);
00316       }
00317     } else if (record_type == PQR_ATOM) {
00318       /* just get the coordinates, and store them */
00319       if (ts != NULL) {
00320         get_pqr_coordinates(pqr_record, x, y, z);
00321         x += 3;
00322         y += 3;
00323         z += 3;
00324       } 
00325       i++;
00326     }
00327   } while(i < pqr->natoms);
00328 
00329   return MOLFILE_SUCCESS;
00330 }
00331 
00332 static void close_pqr_read(void *v) { 
00333   pqrdata *pqr = (pqrdata *)v;
00334   if (pqr->fd) {
00335     fclose(pqr->fd);
00336     pqr->fd = 0;
00337   }
00338   free(pqr);
00339 }
00340 
00341 static void *open_pqr_write(const char *path, const char *filetype, 
00342     int natoms) {
00343   FILE *fd;
00344   pqrdata *pqr;
00345   fd = fopen(path, "w");
00346   if (!fd) {
00347     printf("pqrplugin) unable to open file %s for writing\n", path);
00348     return NULL;
00349   }
00350 
00351   pqr = (pqrdata *)malloc(sizeof(pqrdata));
00352   pqr->fd = fd;
00353   pqr->natoms = natoms; 
00354   pqr->atomlist = NULL;
00355   return pqr;
00356 }
00357  
00358 static int write_pqr_structure(void *v, int optflags, 
00359     const molfile_atom_t *atoms) {
00360   int i;
00361   pqrdata *pqr = (pqrdata *)v;
00362   int natoms = pqr->natoms;
00363   pqr->atomlist = (molfile_atom_t *)malloc(natoms*sizeof(molfile_atom_t));
00364 
00365   memcpy(pqr->atomlist, atoms, natoms*sizeof(molfile_atom_t));
00366 
00367   /* If charge and radius aren't given, we assign defaultvalues. */
00368   if (!(optflags & MOLFILE_CHARGE)) {
00369     printf("pqrplugin) Warning no atom charges available, assigning zero\n");
00370     for (i=0; i<natoms; i++) pqr->atomlist[i].charge = 0.0f;
00371   }
00372   if (!(optflags & MOLFILE_RADIUS)) {
00373     printf("pqrplugin) Warning no atom radii available, assigning radii of 1.0\n");
00374     for (i=0; i<natoms; i++) pqr->atomlist[i].radius = 
00375       get_atom_radius(&(pqr->atomlist[i]));
00376   }
00377 
00378   return MOLFILE_SUCCESS;
00379 }
00380 
00381 static void write_pqr_cryst1(FILE *fd, const molfile_timestep_t *ts) {
00382   fprintf(fd, "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1           1\n", 
00383     ts->A, ts->B, ts->C, ts->alpha, ts->beta, ts->gamma);
00384 }
00385 
00386 static int write_pqr_timestep(void *v, const molfile_timestep_t *ts) {
00387   pqrdata *pqr = (pqrdata *)v; 
00388   const molfile_atom_t *atom;
00389   const float *pos;
00390   int i;
00391 
00392   if (pqr->natoms == 0)
00393     return MOLFILE_SUCCESS;
00394 
00395   write_pqr_cryst1(pqr->fd, ts);
00396   atom = pqr->atomlist;
00397   pos = ts->coords;
00398   for (i=0; i<pqr->natoms; i++) {
00399     if (!write_pqr_atom(pqr->fd, i+1, atom, pos[0], pos[1], pos[2])) {
00400       printf("pqrplugin) error writing atom %d; file may be incomplete.\n", i+1);
00401       return MOLFILE_ERROR;
00402     }
00403     atom++;
00404     pos += 3;
00405   }
00406 
00407   fprintf(pqr->fd, "END\n");
00408 
00409   return MOLFILE_SUCCESS;
00410 }
00411  
00412 static void close_pqr_write(void *v) {
00413   pqrdata *pqr = (pqrdata *)v; 
00414   fclose(pqr->fd);
00415   free(pqr->atomlist);
00416   free(pqr);
00417 }
00418 
00419 /*
00420  * Initialization stuff down here
00421  */
00422 
00423 static molfile_plugin_t plugin;
00424 
00425 VMDPLUGIN_API int VMDPLUGIN_init() {
00426   memset(&plugin, 0, sizeof(molfile_plugin_t));
00427   plugin.abiversion = vmdplugin_ABIVERSION;
00428   plugin.type = MOLFILE_PLUGIN_TYPE;
00429   plugin.name = "pqr";
00430   plugin.prettyname = "PQR";
00431   plugin.author = "Eamon Caddigan";
00432   plugin.majorv = 0;
00433   plugin.minorv = 6;
00434   plugin.is_reentrant = VMDPLUGIN_THREADSAFE;
00435   plugin.filename_extension = "pqr";
00436   plugin.open_file_read = open_pqr_read;
00437   plugin.read_structure = read_pqr_structure;
00438   plugin.read_next_timestep = read_pqr_timestep;
00439   plugin.close_file_read = close_pqr_read;
00440   plugin.open_file_write = open_pqr_write;
00441   plugin.write_structure = write_pqr_structure;
00442   plugin.write_timestep = write_pqr_timestep;
00443   plugin.close_file_write = close_pqr_write;
00444   return VMDPLUGIN_SUCCESS;
00445 }
00446 
00447 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
00448   (*cb)(v, (vmdplugin_t *)&plugin);
00449   return VMDPLUGIN_SUCCESS;
00450 }
00451 
00452 VMDPLUGIN_API int VMDPLUGIN_fini() {
00453   return VMDPLUGIN_SUCCESS;
00454 }
00455 

Generated on Sat Aug 15 03:07:01 2020 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002