lammpsplugin.c Source File

00001 /***************************************************************************
00002  *cr
00003  *cr            (C) Copyright 1995-2006 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: lammpsplugin.c,v $
00013  *      $Author: johns $       $Locker:  $             $State: Exp $
00014  *      $Revision: 1.21 $       $Date: 2008/02/29 20:13:22 $
00015  *
00016  ***************************************************************************/
00017 
00018 /*
00019  *  LAMMPS atom style dump file format:
00020  *    ITEM: TIMESTEP
00021  *      %d (timestep number)
00022  *    ITEM: NUMBER OF ATOMS
00023  *      %d (number of atoms)
00024  *    ITEM: BOX BOUNDS
00025  *      %f %f (alpha, a)
00026  *      %f %f (beta,  b)
00027  *      %f %f (gamma, c)
00028  *    ITEM: ATOMS
00029  *      %d %d %f %f %f  (atomid, atomtype, x, y, z)
00030  *      ...
00031  *
00032  */
00033 
00034 #include "largefiles.h"   /* platform dependent 64-bit file I/O defines */
00035 
00036 #include <stdio.h>
00037 #include <stdlib.h>
00038 #include <string.h>
00039 #include <ctype.h>
00040 #include <math.h>
00041 #include "molfile_plugin.h"
00042 
00043 #ifndef M_PI_2
00044 #define M_PI_2 1.57079632679489661922
00045 #endif
00046 
00047 /* for transparent reading of .gz files */
00048 #ifdef _USE_ZLIB
00049 #include <zlib.h>
00050 #define FileDesc gzFile
00051 #define myFgets(buf,size,fd) gzgets(fd,buf,size)
00052 #define myFopen gzopen
00053 #define myClose gzclose
00054 #define myRewind gzrewind
00055 #else
00056 #define FileDesc FILE*
00057 #define myFopen fopen
00058 #define myClose fclose
00059 #define myFgets(buf,size,fd) fgets(buf,size,fd)
00060 #define myRewind rewind
00061 #endif
00062 
00063 typedef struct {
00064   FileDesc file;
00065   int numatoms;
00066   int scaled_data;
00067   char *file_name;
00068 } lammpsdata;
00069 
00070 static char* find_next_item(FileDesc fd, char* szItem) {
00071   char szLine[1024];
00072   char* pcolon;
00073   int nReturn;
00074 
00075   while(myFgets(szLine, 1024, fd)) {
00076     pcolon = strrchr(szLine, ':');
00077     if(pcolon) {
00078       *pcolon = 0;
00079       if(0 == strcmp(szLine, "ITEM")) {
00080         *pcolon = ':';
00081         pcolon+=2;
00082         strcpy(szItem, pcolon);
00083         nReturn = (int)strcspn(szItem, "\n\r"); /* remove the carriage return */
00084         if(nReturn)
00085           szItem[nReturn] = 0;
00086         return szItem;
00087       }
00088     }
00089   }
00090 
00091   return NULL;
00092 }
00093 
00094 
00095 static int find_item(FileDesc fd, char* szItem) {
00096   char szBuffer[1024];
00097   while(find_next_item(fd, szBuffer)) {
00098     if(0 == strcmp(szBuffer, szItem))
00099       return 1;
00100   }
00101   return 0;
00102 }
00103 
00104  
00105 static void *open_lammps_read(const char *filename, const char *filetype, 
00106                            int *natoms) {
00107   FileDesc fd;
00108   lammpsdata *data;
00109   char szBuffer[1024];
00110 
00111   fd = myFopen(filename, "rb");
00112   if (!fd) return NULL;
00113  
00114   data = (lammpsdata *)malloc(sizeof(lammpsdata));
00115   data->file = fd;
00116   data->file_name = strdup(filename);
00117   *natoms = 0;
00118 
00119   if(!find_item(data->file, "NUMBER OF ATOMS")) {
00120     fprintf(stderr, "lammpsplugin) Unable to find 'NUMBER OF ATOMS' item\n");
00121     return NULL;
00122   }
00123 
00124   if(!myFgets(szBuffer, 1024, data->file)) {
00125     fprintf(stderr, "lammpsplugin) dump file '%s' should have the number of atoms after line ITEM: NUMBER OF ATOMS.\n", filename);
00126     return NULL;
00127   }
00128   *natoms = atoi(szBuffer);
00129  
00130   data->numatoms=*natoms;
00131   data->scaled_data=1;
00132   myRewind(data->file); /* prepare for first read_timestep call */
00133  
00134   return data;
00135 }
00136 
00137 
00138 static int read_lammps_structure(void *mydata, int *optflags, molfile_atom_t *atoms) {
00139   int i, j, idx;
00140   char szBuffer[1024];
00141   molfile_atom_t *atom;
00142   lammpsdata *data = (lammpsdata *)mydata;
00143   int atomid;
00144   char atom_type[1025];
00145   float x, y, z;
00146   char *k;
00147 
00148   *optflags = MOLFILE_NOOPTIONS; /* no optional data */
00149   memset(atoms, 0, data->numatoms * sizeof(molfile_atom_t)); /* clear atoms structure */
00150  
00151   /* go to the beginning of the file */
00152   myRewind(data->file); /* prepare for first read_timestep call */
00153  
00154   /* find the sections with atoms */
00155   if(!find_item(data->file, "ATOMS")) {
00156     fprintf(stderr, "lammpsplugin) couldn't find atoms to read structure from file '%s'.\n", data->file_name);
00157     return MOLFILE_ERROR;
00158   }
00159 
00160   /* now read the atoms */
00161   for(i=0;i<data->numatoms;i++) {
00162     k = myFgets(szBuffer, 1024, data->file);
00163 
00164     /* Read in atom type, X, Y, Z, skipping any remaining data fields */
00165     j = 0;
00166     if (k != NULL)
00167         j = sscanf(szBuffer, "%d %s %f %f %f", &atomid, atom_type, &x, &y, &z);
00168   
00169     if (k == NULL) {
00170       fprintf(stderr, "lammpsplugin) Error while reading structure from lammps dump file '%s': atom missing in the first timestep\n",data->file_name);
00171       fprintf(stderr, "lammpsplugin) expecting '%d' atoms, found only '%d'\n",data->numatoms,i+1);
00172       return MOLFILE_ERROR;
00173     } else if (j < 5) {
00174       fprintf(stderr, "lammpsplugin) Error while reading structure from lammps dump file '%s': missing type or coordinate(s) for atom '%d'\n", data->file_name, i+1);
00175       return MOLFILE_ERROR;
00176     }
00177 
00178     if ((atomid <= 0) || (atomid > data->numatoms)) {
00179       fprintf(stderr, "lammpsplugin) Error while reading structure from lammps dump file '%s': atomid %d is not in valid range [1, %d]\n",data->file_name, atomid, data->numatoms);
00180       return MOLFILE_ERROR;
00181     }
00182 
00183     idx = atomid - 1;
00184     strncpy(atoms[idx].type, atom_type, 16); /* type is char[16]. */
00185     atoms[idx].type[15] = '\0';              /* force null termination */
00186     /* WARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNING *
00187      * Don't use the atomid as name. This will only waste a _lot_ of memory.  *
00188      * The _identical_ information is available via the 'serial' or 'index'   *
00189      * atom properties (numbers starting from 1 or 0, respectively).          *
00190      * WARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNING */
00191     strncpy(atoms[idx].name, atom_type, 16); /* name is char[16]. */
00192     atoms[idx].name[15] = '\0';              /* force null termination */
00193     strncpy(atoms[idx].resname, "UNK", 4);
00194     atoms[idx].resid = 1;
00195     strncpy(atoms[idx].chain, "",1);
00196     strncpy(atoms[idx].segid, "",1);
00197     /* skip to the end of line */
00198 
00199     /* coordinates can be either scaled (fractions of the box)
00200        or unscaled (absolute coordinates) */
00201     if ((x<-0.1) || (x>1.1) || (y<-0.1) || (y>1.1) || (z<-0.1) || (x>1.1)) {
00202         data->scaled_data=0;
00203     }
00204   }
00205 
00206   myRewind(data->file);
00207   return MOLFILE_SUCCESS;
00208 }
00209 
00210 
00211 
00212 static int read_lammps_timestep(void *mydata, int natoms, 
00213                                 molfile_timestep_t *ts) {
00214   int i, j;
00215   char atom_name[1025], szBuffer[1025];
00216   float x, y, z;
00217   int atomid;
00218   float lo[3],hi[3], l[3];
00219 
00220   lammpsdata *data = (lammpsdata *)mydata;
00221 
00222   /* check if there is antother time steps */
00223   if(!find_item(data->file, "TIMESTEP")) return MOLFILE_ERROR;
00224  
00225   /* check if we actually have to read something */
00226   if(!ts) return MOLFILE_SUCCESS;
00227 
00228   /* check the number of atoms in the timestep */
00229   if(!find_item(data->file, "NUMBER OF ATOMS")) {
00230     fprintf(stderr, "lammpsplugin) lammps dump file '%s', unable to find item: NUMBER OF ATOMS for current timestep.\n", data->file_name);
00231     return MOLFILE_ERROR;
00232   }
00233 
00234   if(!myFgets(szBuffer, 1024, data->file)) {
00235     fprintf(stderr, "lammpsplugin) lammps dump file '%s', should have the number of atoms after line ITEM: NUMBER OF ATOMS.\n", data->file_name);
00236     return MOLFILE_ERROR;
00237   }
00238 
00239   if(natoms != atoi(szBuffer)) {
00240     fprintf(stderr, "lammpsplugin) lammps dump file '%s', wrong number of atoms in timestep.\n", data->file_name);
00241     return MOLFILE_ERROR;
00242   }
00243 
00244   /* read the boundary box of the timestep */
00245   if(!find_item(data->file, "BOX BOUNDS")) {
00246     fprintf(stderr, "lammpsplugin) lammps dump file '%s', could not find ITEM: BOX BOUNDS.\n", data->file_name);
00247     return MOLFILE_ERROR;
00248   }
00249 
00250   if (NULL == myFgets(szBuffer, 1024, data->file))  return MOLFILE_ERROR;
00251   sscanf(szBuffer,"%f %f", lo, hi);
00252   l[0] = hi[0] - lo[0];
00253 
00254   if (NULL == myFgets(szBuffer, 1024, data->file))  return MOLFILE_ERROR;
00255   sscanf(szBuffer,"%f %f", lo+1, hi+1);
00256   l[1] = hi[1] - lo[1];
00257 
00258   if (NULL == myFgets(szBuffer, 1024, data->file))  return MOLFILE_ERROR;
00259   sscanf(szBuffer,"%f %f", lo+2, hi+2);
00260   l[2] = hi[2] - lo[2];
00261  
00262   ts->A = l[0];
00263   ts->B = l[1];
00264   ts->C = l[2];
00265   /* XXX: non-orthogonal boxes are now possible, too. */
00266   ts->alpha = 90.0;
00267   ts->beta = 90.0;
00268   ts->gamma = 90.0;
00269 
00270   /* read the coordinates */
00271   if (!find_item(data->file, "ATOMS")) {
00272     fprintf(stderr, "lammpsplugin) lammps dump file '%s', could not find ITEM: ATOMS.\n", data->file_name);
00273     return MOLFILE_ERROR;
00274   }
00275 
00276   for (i=0; i<natoms; i++) {
00277     if(!myFgets(szBuffer, 1024, data->file)) {
00278       /* there is an error or end of file, both cases should not occur */
00279       fprintf(stderr, "lammpsplugin) lammps dump file '%s', unexpected end of file or error while reading coordinates for atom %d.\n", data->file_name, i);
00280       return MOLFILE_ERROR;
00281     }
00282 
00283     /* Read in atom type, X, Y, Z, skipping any remaining data fields */
00284     j = sscanf(szBuffer, "%d %s %f %f %f", &atomid, atom_name, &x, &y, &z);
00285     if (j != 5) {
00286       fprintf(stderr, "lammpsplugin) lammps dump file '%s', error while parsing coordinates for atom %d.\n", data->file_name, i);
00287       return MOLFILE_ERROR;
00288     }
00289 
00290 
00291     if (atomid > 0 && atomid <= natoms) {
00292       int addr = 3 * (atomid - 1);
00293       if(data->scaled_data) {
00294           /* LAMMPS coordinates are fractional unit cell coords 
00295            * by default, so they need to be scaled by a/b/c etc. */
00296           ts->coords[addr    ] = lo[0] + x * l[0];
00297           ts->coords[addr + 1] = lo[1] + y * l[1];
00298           ts->coords[addr + 2] = lo[2] + z * l[2];
00299       } else {
00300           /* ... but they can also be absolute values */
00301           ts->coords[addr    ] = x;
00302           ts->coords[addr + 1] = y;
00303           ts->coords[addr + 2] = z;
00304       }
00305     } else {
00306       fprintf(stderr, "lammpsplugin) ignoring illegal atom index %d\n", atomid);
00307     }
00308   }
00309 
00310   return MOLFILE_SUCCESS;
00311 }
00312     
00313 static void close_lammps_read(void *mydata) {
00314   lammpsdata *data = (lammpsdata *)mydata;
00315   myClose(data->file);
00316   free(data->file_name);
00317   free(data);
00318 }
00319 
00320 
00321 /* registration stuff */
00322 static molfile_plugin_t plugin;
00323 
00324 VMDPLUGIN_API int VMDPLUGIN_init() {
00325   memset(&plugin, 0, sizeof(molfile_plugin_t));
00326   plugin.abiversion = vmdplugin_ABIVERSION;
00327   plugin.type = MOLFILE_PLUGIN_TYPE;
00328   plugin.name = "lammpstrj";
00329   plugin.prettyname = "LAMMPS Trajectory";
00330   plugin.author = "Marco Kalweit, Axel Kohlmeyer, Lutz Maibaum, John Stone";
00331   plugin.majorv = 0;
00332   plugin.minorv = 8;
00333   plugin.is_reentrant = VMDPLUGIN_THREADSAFE;
00334 #ifdef _USE_ZLIB
00335   plugin.filename_extension = "lammpstrj,lammpstrj.gz";
00336 #else
00337   plugin.filename_extension = "lammpstrj";
00338 #endif
00339   plugin.open_file_read = open_lammps_read;
00340   plugin.read_structure = read_lammps_structure;
00341   plugin.read_next_timestep = read_lammps_timestep;
00342   plugin.close_file_read = close_lammps_read;
00343   return VMDPLUGIN_SUCCESS;
00344 }
00345 
00346 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
00347   (*cb)(v, (vmdplugin_t *)&plugin);
00348   return VMDPLUGIN_SUCCESS;
00349 }
00350 
00351 VMDPLUGIN_API int VMDPLUGIN_fini() {
00352   return VMDPLUGIN_SUCCESS;
00353 }
00354 
00355 
00356 #ifdef TEST_PLUGIN
00357 
00358 int main(int argc, char *argv[]) {
00359   molfile_timestep_t timestep;
00360   molfile_atom_t *atoms;
00361   void *v;
00362   int natoms;
00363   int i, nsets, set, opts;
00364 
00365   while (--argc) {
00366     ++argv;
00367     v = open_lammps_read(*argv, "lammps", &natoms);
00368     if (!v) {
00369       fprintf(stderr, "open_lammps_read failed for file %s\n", *argv);
00370       return 1;
00371     }
00372     fprintf(stderr, "open_lammps_read succeeded for file %s\n", *argv);
00373     fprintf(stderr, "number of atoms: %d\n", natoms);
00374 
00375     timestep.coords = (float *)malloc(3*sizeof(float)*natoms);
00376     atoms = (molfile_atom_t *)malloc(sizeof(molfile_atom_t)*natoms);
00377     read_lammps_structure(v, &opts, atoms);
00378     fprintf(stderr, "read_lammps_structure: options=0x%08x\n", opts);
00379 #if 0
00380     for (i=0; i<natoms; ++i) {
00381       fprintf(stderr, "atom %09d: name=%s, type=%s, resname=%s, resid=%d, segid=%s, chain=%s\n",
00382                       i, atoms[i].name, atoms[i].type, atoms[i].resname, atoms[i].resid,
00383                       atoms[i].segid, atoms[i].chain);
00384     }
00385 #endif
00386     i = 0;
00387     while (!read_lammps_timestep(v, natoms, &timestep)) {
00388       i++;
00389     }
00390     fprintf(stderr, "ended read_next_timestep on frame %d\n", i);
00391 
00392     close_lammps_read(v);
00393   }
00394   return 0;
00395 }
00396 
00397 #endif
00398