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

carplugin.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: carplugin.c,v $
00013  *      $Author: johns $       $Locker:  $             $State: Exp $
00014  *      $Revision: 1.18 $       $Date: 2016/11/28 05:01:53 $
00015  *
00016  ***************************************************************************/
00017 
00018 /*
00019  * Plugin for Insight II/Discover car (cartesian coordinate file) file format.
00020  * 
00021  * TODO: 
00022  * + 2D PBC info (probably just as a simplified 3D PBC)
00023  * + HELIX info (not sure how this will be handled)
00024  *
00025  * File Header:
00026  *   !BIOSYM archive 3
00027  *   HELIX 
00028  *     (if HELIX information is present)
00029  *   PBC=ON, PBC=OFF, PBC=2D 
00030  *     (one of these three choices must be present. PBC cannot be "ON" if
00031  *     HELIX information is present.)
00032  * Coordinate Header:
00033  *   Title (this line may be blank, but must be present)
00034  *     1-64 title for the system
00035  *     65-80 energy
00036  *   ! DATE day month date time year
00037  *     (day, month, date, time, year are optional)
00038  *   PBC information if PBC=ON:
00039  *     1-3 PBC
00040  *     4-13 a cell vector a in angstroms
00041  *     14-23 b cell vector b in angstroms
00042  *     24-33 c cell vector c in angstroms
00043  *     34-43 alpha cell angle alpha in degrees
00044  *     44-53 beta cell angle beta in degrees
00045  *     54-63 gamma cell angle gamma in degrees
00046  *     64-80 space group name
00047  *   PBC information if PBC=2D:
00048  *     1-3 PBC
00049  *     4-13 k plane vector k in angstroms
00050  *     14-23 l plane vector l in angstroms
00051  *     24-33 gamma plane angle gamma in degrees
00052  *     34-50 plane group name
00053  * Molecule Data:
00054  *   If helix info is present:
00055  *     1-5 HELIX
00056  *     6-15 sigma in degrees
00057  *     16-25 d in angstroms
00058  *     26-35 kappa angle between l axis and helix axis in degrees
00059  *     36-45 lambda angle between k axis and helix axis in degrees
00060  *     46-55 Tk fractional position of helix axis along k axis
00061  *     56-65 Tl fractional position of helix axis along l axis
00062  *   Atom data:
00063  *     1-5 atom name
00064  *     7-20 x Cartesian coordinate of atom in angstroms
00065  *     22-35 y Cartesian coordinate of atom in angstroms
00066  *     37-50 z Cartesian coordinate of atom in angstroms
00067  *     52-55 type of residue containing atom
00068  *     57-63 residue sequence name relative to beginning of current
00069  *     molecule, left justified
00070  *     64-70 potential type of atom left justified
00071  *     72-73 element symbol
00072  *     75-80 partial charge on atom
00073  *   Final line for a given molecule:
00074  *     1-3 end
00075  * Final line for the entire molecular system input:
00076  *     1-3 end
00077  *
00078  */
00079 
00080 #include <stdio.h>
00081 #include <stdlib.h>
00082 #include <string.h>
00083 #include <ctype.h>
00084 #include "molfile_plugin.h"
00085 
00086 #define LINESIZE 1024
00087 
00088 /* possible values for pbc */
00089 enum {PBC_ON, PBC_OFF, PBC_2D};
00090 
00091 typedef struct {
00092   FILE *file;
00093   int numatoms, pbc, helix, eof;
00094   long coord_location;
00095   molfile_atom_t *atomlist;
00096 } cardata;
00097 
00098 /* Parse a line contianing atom data from a car file and store it in the
00099  * atom structure. Returns 1 on success, 0 on failure.
00100  */
00101 static int read_car_structure_line(molfile_atom_t *atom, char *line) {
00102   char name[LINESIZE], type[LINESIZE];
00103   int resid;
00104   float charge;
00105 
00106   if (sscanf(line, "%s %*f %*f %*f %*s %d %*s %s %f", name, &resid, type, &charge) 
00107       != 4)
00108     return 0;
00109 
00110   /* check the length of the name and type strings before copying. */
00111   if ( (strlen(name) >= 8) || (strlen(type) >= 8) )
00112     return 0;
00113   strcpy(atom->name, name);
00114   strcpy(atom->type, type);
00115 
00116   /* Check that the resid won't overflow the resname string, then copy it
00117    * over
00118    */
00119   if (resid > 9999999)
00120     atom->resname[0] = '\0';
00121   else
00122     sprintf(atom->resname, "%d", resid);
00123 
00124   atom->resid = resid;
00125 
00126   atom->chain[0] = '\0';
00127   atom->segid[0] = '\0';
00128 
00129   atom->charge = charge;
00130 
00131   return 1;
00132 }
00133 
00134 /* Parse a line contianing atom coordinates from a car file and store them
00135  * in the array 'coords' in XYZ order. Returns 1 on success, 0 on failure.
00136  */
00137 static int read_car_coordinates(float *coords, char *line) {
00138   float x, y, z;
00139 
00140   if (sscanf(line, "%*s %f %f %f %*s %*d %*s %*s %*f", &x, &y, &z) 
00141       != 3)
00142     return 0;
00143 
00144   coords[0] = x;
00145   coords[1] = y;
00146   coords[2] = z;
00147 
00148   return 1;
00149 }
00150 
00151 static void *open_car_read(const char *filename, const char *filetype, 
00152                            int *natoms) {
00153   FILE *fd;
00154   cardata *data;
00155   char line[LINESIZE];
00156 
00157   fd = fopen(filename, "rb"); if (!fd) return NULL;
00158   
00159   data = (cardata *)malloc(sizeof(cardata));
00160   data->eof = 0;
00161   data->file = fd;
00162 
00163   /* First line: "!BIOSYM archive n", where n indicates the file format */
00164   fgets(line, LINESIZE, fd);
00165   if (strncmp(line, "!BIOSYM archive", 15)) {
00166     fprintf(stderr, "ERROR) badly formatted/missing header.\n");
00167     return NULL;
00168   }
00169 
00170   /* Second line: "HELIX" if helix information is present. Followed by PBC
00171    * info on the next line: "PBC=ON|OFF|2D". 
00172    * If HELIX is present, PBC cannot be "ON": return an error if this happens. 
00173    */
00174   fgets(line, LINESIZE, fd);
00175   if (!strncmp(line, "HELIX", 5)) {
00176     data->helix = 1;
00177     fgets(line, LINESIZE, fd);
00178     fprintf(stdout, "WARNING) ignoring helix information.\n");
00179   }
00180   else {
00181     data->helix = 0;
00182   }
00183 
00184   if (!strncmp(line, "PBC=ON", 6)) {
00185     data->pbc = PBC_ON;
00186   }
00187   else if (!strncmp(line, "PBC=OFF", 7)) {
00188     data->pbc = PBC_OFF;
00189   }
00190   else if (!strncmp(line, "PBC=2D", 6)) {
00191     data->pbc = PBC_2D;
00192     fprintf(stdout, "WARNING) ignoring 2D PBC information.\n");
00193   }
00194   else {
00195     fprintf(stderr, "ERROR) badly formatted/missing PBC info.\n");
00196     return NULL;
00197   }
00198 
00199   if (data->helix && (data->pbc == PBC_ON)) {
00200     fprintf(stderr, "ERROR) car file contains helix and 3D PBC information.");
00201     return NULL;
00202   }
00203 
00204   /* Next line: title/energy for the system. Skipped. */
00205   fgets(line, LINESIZE, fd);
00206 
00207   /* Next line: "!DATE [day month date time year]". */
00208   fgets(line, LINESIZE, fd);
00209   if (strncmp(line, "!DATE", 5)) {
00210     fprintf(stderr, "ERROR) badly formatted/missing date.\n");
00211     return NULL;
00212   }
00213     
00214   /* Store the location of the beginning of the PBC/coordinate data. */
00215   data->coord_location = ftell(fd);
00216 
00217   /* Skip the PBC and HELIX entries, if present */
00218   if (data->pbc != PBC_OFF) 
00219     fgets(line, LINESIZE, fd);
00220   if (data->helix)
00221     fgets(line, LINESIZE, fd);
00222 
00223   /* Count the atoms in all molecules*/
00224   data->numatoms = 0;
00225   fgets(line, LINESIZE, fd);
00226   while(strncmp(line, "end", 3)) {
00227     while(strncmp(line, "end", 3)) {
00228       data->numatoms++;
00229       fgets(line, LINESIZE, fd);
00230 
00231       if (feof(fd)) {
00232         fprintf(stderr, "ERROR) unexpected end-of-file.\n");
00233         return NULL;
00234       }
00235       if (ferror(fd)) {
00236         fprintf(stderr, "ERROR) error reading car file.\n");
00237         return NULL;
00238       }
00239     }
00240     fgets(line, LINESIZE, fd);
00241   }
00242   *natoms = data->numatoms;
00243 
00244   return data;
00245 }
00246 
00247 static int read_car_structure(void *mydata, int *optflags, 
00248                               molfile_atom_t *atoms) {
00249   int mol_num;
00250   char line[LINESIZE];
00251   molfile_atom_t *atom;
00252   cardata *data = (cardata *)mydata;
00253 
00254   *optflags = MOLFILE_CHARGE; /* car files contain partial charges */
00255 
00256   /* move to the beginning of the atom data in the file, skipping any PBC or
00257    * HELIX information that may be present
00258    */
00259   fseek(data->file, data->coord_location, SEEK_SET);
00260   if (data->pbc != PBC_OFF) 
00261     fgets(line, LINESIZE, data->file);
00262   if (data->helix)
00263     fgets(line, LINESIZE, data->file);
00264 
00265   mol_num = 0;
00266   atom = atoms;
00267   fgets(line, LINESIZE, data->file);
00268   /* Loop through all molecules */
00269   while(strncmp(line, "end", 3)) {
00270     /* Read the structure for each molecule */
00271     while(strncmp(line, "end", 3)) {
00272       if (!read_car_structure_line(atom, line)) {
00273         fprintf(stderr, "ERROR) badly formatted structure line:\n%s\n", line);
00274         return MOLFILE_ERROR;
00275       }
00276 
00277       /* XXX - this code can easily create buffer overflows
00278                and cause data corruption or segfaults:
00279 
00280          XXX - use the chain name to identify different molecules
00281          sprintf(atom->chain, "%d", mol_num);
00282 
00283          we avoid the buffer overflow and use typical capital 
00284          letters for the chain identifiers instead.
00285       */
00286       sprintf(atom->chain, "%c", (char) (((int) 'A') + (mol_num % 26))); 
00287 
00288       atom++;
00289       
00290       fgets(line, LINESIZE, data->file);
00291       if (feof(data->file)) {
00292         fprintf(stderr, "ERROR) unexpected end-of-file while reading structure.\n");
00293         return MOLFILE_ERROR;
00294       }
00295       if (ferror(data->file)) {
00296         fprintf(stderr, "ERROR) error reading car file while reading structure.\n");
00297         return MOLFILE_ERROR;
00298       }
00299     }
00300     fgets(line, LINESIZE, data->file);
00301     mol_num++;
00302   }
00303 
00304   return MOLFILE_SUCCESS;
00305 }
00306 
00307 static int read_car_timestep(void *mydata, int natoms, molfile_timestep_t *ts) {
00308   char line[LINESIZE];
00309   cardata *data = (cardata *)mydata;
00310   float *coords = NULL;
00311 
00312   /* return if coordinates have been read */
00313   if (data->eof)
00314     return MOLFILE_EOF;
00315 
00316   /* move to the beginning of the atom data in the file */
00317   fseek(data->file, data->coord_location, SEEK_SET);
00318 
00319   /* read the PBC record if present */
00320   if (data->pbc == PBC_ON) {
00321     fgets(line, LINESIZE, data->file);
00322 
00323     if (ts) {
00324       if ( sscanf(line, "PBC %f %f %f %f %f %f %*s", 
00325                   &ts->A, &ts->B, &ts->C, 
00326                   &ts->alpha, &ts->beta, &ts->gamma) != 6 ) {
00327         fprintf(stderr, "ERROR) badly formatted PBC line:\n%s\n", line);
00328         return MOLFILE_ERROR;
00329       }
00330     }
00331   }
00332   else if (data->pbc == PBC_2D) {
00333     /* XXX - Ignore 2D PBC information */
00334     fgets(line, LINESIZE, data->file);
00335   }
00336 
00337   /* skip the helix record if present */
00338   if (data->helix)
00339     fgets(line, LINESIZE, data->file);
00340 
00341   if (ts)
00342     coords = ts->coords;
00343 
00344   fgets(line, LINESIZE, data->file);
00345   /* Loop through all molecules */
00346   while(strncmp(line, "end", 3)) {
00347     /* Read the coordinates for each molecule */
00348     while(strncmp(line, "end", 3)) {
00349       /* only save coords if we're given a timestep pointer. */
00350       if (ts) {
00351         if (!read_car_coordinates(coords, line)) {
00352           fprintf(stderr, "ERROR) badly formatted coordinate line:\n%s\n", line);
00353           return MOLFILE_ERROR;
00354         }
00355         coords += 3;
00356       }
00357 
00358       fgets(line, LINESIZE, data->file);
00359       if (feof(data->file)) {
00360         fprintf(stderr, "ERROR) unexpected end-of-file while reading coordinates.\n");
00361         return MOLFILE_ERROR;
00362       }
00363       if (ferror(data->file)) {
00364         fprintf(stderr, "ERROR) file error while reading coordinates.\n");
00365         return MOLFILE_ERROR;
00366       }
00367     }
00368     fgets(line, LINESIZE, data->file);
00369   }
00370 
00371   /* set eof since this file contians only one "timestep" */
00372   data->eof = 1;
00373 
00374   return MOLFILE_SUCCESS;
00375 }
00376     
00377 static void close_car_read(void *mydata) {
00378   if (mydata) {
00379     cardata *data = (cardata *)mydata;
00380     fclose(data->file);
00381     free(data);
00382   }
00383 }
00384 
00385 
00386 /* registration stuff */
00387 static molfile_plugin_t plugin;
00388 
00389 VMDPLUGIN_API int VMDPLUGIN_init() {
00390   memset(&plugin, 0, sizeof(molfile_plugin_t));
00391   plugin.abiversion = vmdplugin_ABIVERSION;
00392   plugin.type = MOLFILE_PLUGIN_TYPE;
00393   plugin.name = "car";
00394   plugin.prettyname = "InsightII car";
00395   plugin.author = "Eamon Caddigan";
00396   plugin.majorv = 0;
00397   plugin.minorv = 5;
00398   plugin.is_reentrant = VMDPLUGIN_THREADSAFE;
00399   plugin.filename_extension = "car";
00400   plugin.open_file_read = open_car_read;
00401   plugin.read_structure = read_car_structure;
00402   plugin.read_next_timestep = read_car_timestep;
00403   plugin.close_file_read = close_car_read;
00404   return VMDPLUGIN_SUCCESS;
00405 }
00406 
00407 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
00408   (*cb)(v, (vmdplugin_t *)&plugin);
00409   return VMDPLUGIN_SUCCESS;
00410 }
00411 
00412 VMDPLUGIN_API int VMDPLUGIN_fini() {
00413   return VMDPLUGIN_SUCCESS;
00414 }
00415 

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