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

lammpsplugin.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: lammpsplugin.c,v $
00013  *      $Author: johns $       $Locker:  $             $State: Exp $
00014  *      $Revision: 1.49 $       $Date: 2016/11/28 05:01:54 $
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 (boxxlo, boxxhi)
00026  *      %f %f (boxylo, boxyhi)
00027  *      %f %f (boxzlo, boxzhi)
00028  *    ITEM: ATOMS
00029  *      %d %d %f %f %f  (atomid, atomtype, x, y, z)
00030  *      ...
00031  * newer LAMMPS versions have instead
00032  *    ITEM: ATOMS id x y z
00033  *      %d %d %f %f %f  (atomid, atomtype, x, y, z)
00034  *      ...
00035  * also triclinic boxes are possible:
00036  *    ITEM: BOX BOUNDS xy xz yz
00037  *      %f %f %f (boxxlo, boxxhi, xy)
00038  *      %f %f %f (boxylo, boxyhi, xz)
00039  *      %f %f %f (boxzlo, boxzhi, yz)
00040  * 
00041  * as of 11 April 2011 box bounds always include periodicity settings.
00042  *    ITEM: BOX BOUNDS pp pp pp xy xz yz
00043  * instead of p (periodic) also f (fixed), s (shrinkwrap) and m (shrinkwrap 
00044  *         with minimum) are possible boundaries.
00045  *
00046  * SPECIAL NOTICE: these are box _boundaries_ not lengths.
00047  *                 the dimensions of the box still need to 
00048  *                 be calculated from them and xy,xz,yz.
00049  *
00050  * the newer format allows to handle custom dumps with velocities
00051  * and other features that are not yet in VMD and the molfile API.
00052  */
00053 
00054 #include "largefiles.h"   /* platform dependent 64-bit file I/O defines */
00055 
00056 #include <stdio.h>
00057 #include <stdlib.h>
00058 #include <string.h>
00059 #include <ctype.h>
00060 #include <math.h>
00061 #include "molfile_plugin.h"
00062 
00063 #include "periodic_table.h"
00064 
00065 #define THISPLUGIN plugin
00066 #include "vmdconio.h"
00067 
00068 #define VMDPLUGIN_STATIC
00069 #include "hash.h"
00070 #include "inthash.h"
00071 
00072 #ifndef LAMMPS_DEBUG
00073 #define LAMMPS_DEBUG 0
00074 #endif
00075 
00076 #ifndef M_PI_2
00077 #define M_PI_2 1.57079632679489661922
00078 #endif
00079 
00080 #ifndef MIN
00081 #define MIN(A,B) ((A) < (B) ? (A) : (B))
00082 #endif
00083 #ifndef MAX
00084 #define MAX(A,B) ((A) > (B) ? (A) : (B))
00085 #endif
00086 
00087 /* small magnitude floating point number */
00088 #define SMALL 1.0e-12f
00089 
00090 /* maximum supported length of line for line buffers */
00091 #define LINE_LEN 1024
00092 
00093 /* lammps item keywords */
00094 #define KEY_ATOMS "NUMBER OF ATOMS"
00095 #define KEY_BOX   "BOX BOUNDS"
00096 #define KEY_DATA  "ATOMS"
00097 #define KEY_TSTEP "TIMESTEP"
00098 
00099 /* lammps coordinate data styles */
00100 #define LAMMPS_COORD_NONE       0x000U
00101 #define LAMMPS_COORD_WRAPPED    0x001U
00102 #define LAMMPS_COORD_SCALED     0x002U
00103 #define LAMMPS_COORD_IMAGES     0x004U
00104 #define LAMMPS_COORD_UNWRAPPED  0x008U
00105 #define LAMMPS_COORD_UNKNOWN    0x010U
00106 #define LAMMPS_COORD_VELOCITIES 0x020U
00107 #define LAMMPS_COORD_FORCES     0x040U
00108 #define LAMMPS_COORD_DIPOLE     0x080U
00109 #define LAMMPS_COORD_TRICLINIC  0x100U
00110 
00112 #define LAMMPS_MAX_NUM_FIELDS 64
00113 enum lammps_attribute {
00114   LAMMPS_FIELD_UNKNOWN=0,
00115   LAMMPS_FIELD_ATOMID, LAMMPS_FIELD_MOLID,  LAMMPS_FIELD_TYPE,
00116   LAMMPS_FIELD_POSX,   LAMMPS_FIELD_POSY,   LAMMPS_FIELD_POSZ, 
00117   LAMMPS_FIELD_POSXS,  LAMMPS_FIELD_POSYS,  LAMMPS_FIELD_POSZS,
00118   LAMMPS_FIELD_POSXU,  LAMMPS_FIELD_POSYU,  LAMMPS_FIELD_POSZU,
00119   LAMMPS_FIELD_POSXSU, LAMMPS_FIELD_POSYSU, LAMMPS_FIELD_POSZSU,
00120   LAMMPS_FIELD_IMGX,   LAMMPS_FIELD_IMGY,   LAMMPS_FIELD_IMGZ,
00121   LAMMPS_FIELD_VELX,   LAMMPS_FIELD_VELY,   LAMMPS_FIELD_VELZ,
00122   LAMMPS_FIELD_FORX,   LAMMPS_FIELD_FORY,   LAMMPS_FIELD_FORZ,
00123   LAMMPS_FIELD_CHARGE, LAMMPS_FIELD_RADIUS, LAMMPS_FIELD_DIAMETER,
00124   LAMMPS_FIELD_ELEMENT,LAMMPS_FIELD_MASS,   LAMMPS_FIELD_QUATW,
00125   LAMMPS_FIELD_QUATI,  LAMMPS_FIELD_QUATJ,  LAMMPS_FIELD_QUATK,
00126   LAMMPS_FIELD_MUX,    LAMMPS_FIELD_MUY,    LAMMPS_FIELD_MUZ,
00127   LAMMPS_FIELD_USER0,  LAMMPS_FIELD_USER1,  LAMMPS_FIELD_USER2,
00128   LAMMPS_FIELD_USER3,  LAMMPS_FIELD_USER4,  LAMMPS_FIELD_USER5,
00129   LAMMPS_FIELD_USER6,  LAMMPS_FIELD_USER7,  LAMMPS_FIELD_USER8,
00130   LAMMPS_FILED_USER9
00131 };
00132 
00133 typedef enum lammps_attribute l_attr_t;
00134 
00135 /* for transparent reading of .gz files */
00136 #ifdef _USE_ZLIB
00137 #include <zlib.h>
00138 #define FileDesc gzFile
00139 #define myFgets(buf,size,fd) gzgets(fd,buf,size)
00140 #define myFprintf gzprintf
00141 #define myFopen gzopen
00142 #define myFclose gzclose
00143 #define myRewind gzrewind
00144 #else
00145 #define FileDesc FILE*
00146 #define myFprintf fprintf
00147 #define myFopen fopen
00148 #define myFclose fclose
00149 #define myFgets(buf,size,fd) fgets(buf,size,fd)
00150 #define myRewind rewind
00151 #endif
00152 
00153 typedef struct {
00154   FileDesc file;
00155   FILE *fp;
00156   char *file_name;
00157   int *atomtypes;
00158   int numatoms;
00159   int maxatoms;
00160   int nstep;
00161   unsigned int coord_data; /* indicate type of coordinate data   */
00162   float dip2atoms;         /* scaling factor for dipole to atom data */
00163   float dumx,dumy,dumz;    /* location of dummy/disabled atoms */
00164   int numfields;           /* number of data fields present */
00165   l_attr_t field[LAMMPS_MAX_NUM_FIELDS]; /* type of data fields in dumps */
00166   inthash_t *idmap;        /* for keeping track of atomids */
00167   int fieldinit;           /* whether the field mapping was initialized */
00168 #if vmdplugin_ABIVERSION > 10
00169   molfile_timestep_metadata_t ts_meta;
00170 #endif
00171 } lammpsdata;
00172 
00173 /* merge sort for integer atom id map: merge function */
00174 static void id_merge(int *output, int *left, int nl, int *right, int nr)
00175 {
00176     int i,l,r;
00177     i = l = r = 0;
00178     while ((l < nl) && (r < nr)) {
00179         if (left[l] < right[r])
00180             output[i++] = left[l++];
00181         else
00182             output[i++] = right[r++];
00183     }
00184     while (l < nl)
00185         output[i++] = left[l++];
00186     while (r < nr)
00187         output[i++] = right[r++];
00188 }
00189 
00190 /* bottom up merge sort for integer atom id map: main function */
00191 static void id_sort(int *idmap, int num)
00192 {
00193     int *hold;
00194     int i,j,k;
00195     
00196     hold = (int *)malloc(num*sizeof(int));
00197     if (hold == NULL) return;
00198 
00199     for (i=1; i < num; i *=2) {
00200         memcpy(hold,idmap,num*sizeof(int));
00201         for (j=0; j < (num - i); j += 2*i) {
00202             k =(j+2*i > num) ? num-j-i : i;
00203             id_merge(idmap+j, hold+j, i, hold+j+i, k);
00204         }
00205     }
00206     free((void *)hold);
00207 }
00208 
00215 static const char *remap_field(const char *tag, const char *list) 
00216 {
00217   int i, pos, len, flag;
00218   const char *ptr;
00219   static char to[32], from[32];
00220 
00221   /* no point in trying to match against NULL pointers */
00222   if ((!tag) || (!list))
00223     return tag;
00224 
00225   ptr=list;
00226   i=pos=flag=0;
00227   len=strlen(list);
00228 
00229   /* loop over whole string */
00230   while (pos < len) {
00231 
00232     if (flag) { /* case 1: determine the "from" string */
00233       if (ptr[pos] == ',') { /* end of value */
00234         from[i] = '\0';
00235         if (0 == strcmp(tag,from)) { 
00236           /* only return a token if it is non-NULL */
00237           if (strlen(to))
00238             return to;
00239 
00240           flag=0;
00241           i=0;
00242         } else {
00243           /* try next key */
00244           flag=0;
00245           i=0;
00246         }
00247       } else { /* copy into "from" */
00248         from[i]=ptr[pos];
00249         if (i<30)
00250           ++i;
00251       }
00252     } else { /* case 2: determine the "to" string */
00253 
00254       if (ptr[pos] == '=') { /* end of the key */
00255         to[i] = '\0';
00256         flag=1;
00257         i=0;
00258       } else if (ptr[pos] == ',') { /* incomplete entry. reset "to".*/
00259         i=0;
00260         flag=0;
00261       } else { /* copy into "to" */
00262         to[i]=ptr[pos];
00263         if (i<30)
00264           ++i;
00265       }
00266     }
00267     ++pos;
00268   }
00269 
00270   /* we reached end of the list */
00271   if (flag) {
00272     from[i] = '\0';
00273     if (0 == strcmp(tag,from)) { 
00274       /* only return a token if it is non-NULL */
00275       if (strlen(to))
00276         return to;
00277     }
00278   }
00279 
00280   return tag;
00281 }
00282 
00288 static char* find_next_item(FileDesc fd, char* linebuf, int buflen) {
00289   char* ptr;
00290 
00291   while(myFgets(linebuf, buflen, fd)) {
00292 
00293     /* strip of leading whitespace */
00294     ptr = linebuf;
00295     while (ptr && (*ptr == ' ' || *ptr == '\t'))
00296       ++ptr;
00297 
00298     /* check if this is an "item" */
00299     if(0 == strncmp(ptr, "ITEM:", 5)) {
00300       ptr += 5;
00301       return ptr;
00302     }
00303   }
00304 
00305   return NULL;
00306 }
00307 
00316 static char *find_item_keyword(FileDesc fd, const char* keyword,
00317                                char *linebuf, int buflen) {
00318   char *ptr;
00319   int len;
00320   
00321   while(1) {
00322     ptr = find_next_item(fd, linebuf, buflen);
00323 
00324     if (ptr == NULL) 
00325       break;
00326     
00327     while (ptr && (*ptr == ' ' || *ptr == '\t'))
00328       ++ptr;
00329 
00330 #if LAMMPS_DEBUG
00331     fprintf(stderr, "text=%s/%s", keyword, ptr);
00332 #endif
00333     len = strlen(keyword);
00334     if (0 == strncmp(ptr, keyword, len) ) {
00335       ptr += len;
00336       if (*ptr == '\0' || *ptr == ' ' || *ptr == '\n' || *ptr == '\r') {
00337 #if LAMMPS_DEBUG
00338         fprintf(stderr, "return=%s", ptr);
00339 #endif
00340         return ptr;
00341       } else continue; /* keyword was not an exact match, try again. */
00342     }
00343   }
00344 #if LAMMPS_DEBUG
00345   fprintf(stderr, "return='NULL'\n");
00346 #endif
00347   return NULL;
00348 }
00349 
00350  
00351 static void *open_lammps_read(const char *filename, const char *filetype, 
00352                            int *natoms) {
00353   FileDesc fd;
00354   lammpsdata *data;
00355   char buffer[LINE_LEN];
00356   char *ptr;
00357   const char *envvar;
00358   long tmp, maxatoms;
00359 
00360   fd = myFopen(filename, "rb");
00361   if (!fd) return NULL;
00362  
00363   data = (lammpsdata *)calloc(1, sizeof(lammpsdata));
00364   data->file = fd;
00365   data->file_name = strdup(filename);
00366   data->dip2atoms = -1.0;
00367   data->fieldinit = 0;
00368   *natoms = 0;
00369   maxatoms = 0;
00370   
00371   ptr = find_item_keyword(data->file, KEY_ATOMS,  buffer, LINE_LEN);
00372   if (ptr == NULL) {
00373     vmdcon_printf(VMDCON_ERROR, "lammpsplugin) Unable to find '%s' item.\n",
00374                   KEY_ATOMS);
00375     return NULL;
00376   }
00377 
00378   if (!myFgets(buffer, LINE_LEN, data->file)) {
00379     vmdcon_printf(VMDCON_ERROR, "lammpsplugin) dump file '%s' should "
00380                   "have the number of atoms after line ITEM: %s\n", 
00381                   filename, KEY_ATOMS);
00382     return NULL;
00383   }
00384 
00385   tmp = atol(buffer);
00386   /* we currently only support 32-bit integer atom numbers */
00387   if (tmp > 0x7FFFFFFF) {
00388     vmdcon_printf(VMDCON_ERROR, "lammpsplugin) dump file '%s' contains "
00389                   "%ld atoms which is more than what this plugin supports.\n", 
00390                   filename, tmp);
00391     return NULL;
00392   }
00393 
00394   /* hack to allow trajectories with varying number of atoms.
00395    * we initialize the structure to the value of LAMMPSMAXATOMS
00396    * if that environment variable exists and has a larger value
00397    * than the number of actual atoms in the first frame of the
00398    * dump file. This way we can provision additional space and
00399    * automatically support systems where atoms are lost.
00400    * coordinates of atoms that are no longer present are set to
00401    * either 0.0 or the content of LAMMPSDUMMYPOS.
00402    */
00403   envvar=getenv("LAMMPSMAXATOMS");
00404   if (envvar) maxatoms = atol(envvar);
00405   data->dumx = data->dumy = data->dumz = 0.0f;
00406   envvar=getenv("LAMMPSDUMMYPOS");
00407   if (envvar) sscanf(envvar,"%f%f%f", &(data->dumx),
00408                      &(data->dumy), &(data->dumz));
00409 
00410   if (maxatoms > tmp) {
00411     vmdcon_printf(VMDCON_INFO, "lammpsplugin) provisioning space for up to "
00412                   "%ld atoms.\n", maxatoms);
00413   } else maxatoms = tmp;
00414   *natoms = maxatoms;
00415 
00416   /* hack to allow displaying dipoles as two atoms. 
00417    * the presence of the environment variable toggles
00418    * feature. its value is a scaling factor.
00419    */
00420   envvar=getenv("LAMMPSDIPOLE2ATOMS");
00421   if (envvar) {
00422     data->dip2atoms = (float) strtod(envvar,NULL);
00423     maxatoms *= 2;
00424     tmp *=2;
00425   }
00426   *natoms = maxatoms;
00427  
00428   data->maxatoms = maxatoms;  /* size of per-atom storage */
00429   data->numatoms = tmp;       /* number of atoms initialized atoms */
00430   data->coord_data = LAMMPS_COORD_NONE;  
00431   myRewind(data->file); /* prepare for first read_timestep call */
00432  
00433   return data;
00434 }
00435 
00436 
00437 static int read_lammps_structure(void *mydata, int *optflags, 
00438                                  molfile_atom_t *atoms) {
00439   int i, j;
00440   char buffer[LINE_LEN];
00441   lammpsdata *data = (lammpsdata *)mydata;
00442   int atomid, atomtype, needhash;
00443   float x, y, z;
00444   char *fieldlist;
00445   molfile_atom_t thisatom;
00446   const char *k;
00447   int *idlist=NULL;
00448   
00449   /* clear atom info. */
00450   *optflags = MOLFILE_NOOPTIONS; 
00451   data->coord_data = LAMMPS_COORD_NONE;
00452   memset(atoms, 0, data->numatoms * sizeof(molfile_atom_t));
00453 
00454   /*  fake info for dummy atoms */
00455   strcpy(thisatom.name,"@");
00456   strcpy(thisatom.type,"X");
00457   strcpy(thisatom.resname,"@@@");
00458   strcpy(thisatom.segid,"@@@");
00459   strcpy(thisatom.chain,"@");
00460   thisatom.resid = -1;
00461   thisatom.occupancy = -1.0;
00462   thisatom.bfactor = -1.0;
00463   thisatom.mass = 0.0;
00464   thisatom.charge = 0.0;
00465   thisatom.radius = 0.0;
00466   thisatom.atomicnumber = 0;
00467   for (i=data->numatoms; i < data->maxatoms; ++i)
00468     memcpy(atoms+i, &thisatom, sizeof(molfile_atom_t)); 
00469 
00470 #if vmdplugin_ABIVERSION > 10
00471   data->ts_meta.count = -1;
00472   data->ts_meta.has_velocities = 0;
00473 #endif
00474  
00475   /* go to the beginning of the file */
00476   myRewind(data->file); /* prepare for first read_timestep call */
00477 
00478   /* find the boundary box info to determine if triclinic or not. */
00479   fieldlist = find_item_keyword(data->file, KEY_BOX, buffer, LINE_LEN);
00480   if (fieldlist == NULL) {
00481     vmdcon_printf(VMDCON_ERROR, "lammpsplugin) Could not find box boundaries in timestep.\n");
00482     return MOLFILE_ERROR;
00483   }
00484   k = myFgets(buffer, LINE_LEN, data->file);
00485   if (k == NULL) {
00486     vmdcon_printf(VMDCON_ERROR, "lammpsplugin) Could not find box boundaries in timestep.\n");
00487     return MOLFILE_ERROR;
00488   }
00489 
00490   j = sscanf(buffer, "%f%f%f", &x, &y, &z);
00491   if (j == 3) {
00492     data->coord_data |= LAMMPS_COORD_TRICLINIC;
00493   } else if (j < 2) {
00494     vmdcon_printf(VMDCON_ERROR, "lammpsplugin) Could not find box boundaries in timestep.\n");
00495     return MOLFILE_ERROR;
00496   }
00497 
00498   /* find the sections with atoms */
00499   fieldlist = find_item_keyword(data->file, KEY_DATA, buffer, LINE_LEN);
00500   if (fieldlist == NULL) {
00501     vmdcon_printf(VMDCON_ERROR, "lammpsplugin) Couldn't find data to "
00502                   "read structure from file '%s'.\n", data->file_name);
00503     return MOLFILE_ERROR;
00504   }
00505 
00506 #if LAMMPS_DEBUG  
00507   fprintf(stderr,"fieldlist for atoms: %s", fieldlist);
00508 #if 0  /* simulate old style trajectory */
00509   fieldlist = strdup("\n");
00510 #endif
00511 #endif
00512 
00513   /* parse list of fields */
00514   i = 0;
00515   k = strtok(fieldlist, " \t\n\r");
00516   if (k == NULL) {
00517     /* assume old style lammps trajectory  */
00518     vmdcon_printf(VMDCON_WARN, "lammpsplugin) Found old style trajectory. "
00519                   "assuming data is ordered "
00520                   "'id type x|xs|xu y|ys|yu z|zs|zu [...]'.\n");
00521     data->coord_data |= LAMMPS_COORD_UNKNOWN;
00522   } else {
00523     /* try to identify supported output types */
00524     do {
00525       /* hack to allow re-mapping of arbitrary fields. */
00526       const char *envvar;
00527       envvar=getenv("LAMMPSREMAPFIELDS");
00528       if (envvar)
00529         k=remap_field(k,envvar);
00530       
00531       if (0 == strcmp(k, "id")) {
00532         data->field[i] = LAMMPS_FIELD_ATOMID;
00533       } else if (0 == strcmp(k, "mol")) {
00534         data->field[i] = LAMMPS_FIELD_MOLID;
00535       } else if (0 == strcmp(k, "type")) {
00536         data->field[i] = LAMMPS_FIELD_TYPE;
00537       } else if (0 == strcmp(k, "x")) {
00538         data->field[i] = LAMMPS_FIELD_POSX;
00539         data->coord_data |= LAMMPS_COORD_WRAPPED;
00540       } else if (0 == strcmp(k, "y")) {
00541         data->field[i] = LAMMPS_FIELD_POSY;
00542         data->coord_data |= LAMMPS_COORD_WRAPPED;
00543       } else if (0 == strcmp(k, "z")) {
00544         data->field[i] = LAMMPS_FIELD_POSZ;
00545         data->coord_data |= LAMMPS_COORD_WRAPPED;
00546       } else if (0 == strcmp(k, "xs")) {
00547         data->field[i] = LAMMPS_FIELD_POSXS;
00548         data->coord_data |= LAMMPS_COORD_SCALED;
00549       } else if (0 == strcmp(k, "ys")) {
00550         data->field[i] = LAMMPS_FIELD_POSYS;
00551         data->coord_data |= LAMMPS_COORD_SCALED;
00552       } else if (0 == strcmp(k, "zs")) {
00553         data->field[i] = LAMMPS_FIELD_POSZS;
00554         data->coord_data |= LAMMPS_COORD_SCALED;
00555       } else if (0 == strcmp(k, "xu")) {
00556         data->field[i] = LAMMPS_FIELD_POSXU;
00557         data->coord_data |= LAMMPS_COORD_UNWRAPPED;
00558       } else if (0 == strcmp(k, "yu")) {
00559         data->field[i] = LAMMPS_FIELD_POSYU;
00560         data->coord_data |= LAMMPS_COORD_UNWRAPPED;
00561       } else if (0 == strcmp(k, "zu")) {
00562         data->field[i] = LAMMPS_FIELD_POSZU;
00563         data->coord_data |= LAMMPS_COORD_UNWRAPPED;
00564       } else if (0 == strcmp(k, "xus")) {
00565         data->field[i] = LAMMPS_FIELD_POSXU;
00566         data->coord_data |= LAMMPS_COORD_UNWRAPPED;
00567         data->coord_data |= LAMMPS_COORD_SCALED;
00568       } else if (0 == strcmp(k, "yus")) {
00569         data->field[i] = LAMMPS_FIELD_POSYU;
00570         data->coord_data |= LAMMPS_COORD_UNWRAPPED;
00571         data->coord_data |= LAMMPS_COORD_SCALED;
00572       } else if (0 == strcmp(k, "zus")) {
00573         data->field[i] = LAMMPS_FIELD_POSZU;
00574         data->coord_data |= LAMMPS_COORD_UNWRAPPED;
00575         data->coord_data |= LAMMPS_COORD_SCALED;
00576       } else if (0 == strcmp(k, "ix")) {
00577         data->field[i] = LAMMPS_FIELD_IMGX;
00578         data->coord_data |= LAMMPS_COORD_IMAGES;
00579       } else if (0 == strcmp(k, "iy")) {
00580         data->field[i] = LAMMPS_FIELD_IMGY;
00581         data->coord_data |= LAMMPS_COORD_IMAGES;
00582       } else if (0 == strcmp(k, "iz")) {
00583         data->field[i] = LAMMPS_FIELD_IMGZ;
00584         data->coord_data |= LAMMPS_COORD_IMAGES;
00585       } else if (0 == strcmp(k, "vx")) {
00586         data->field[i] = LAMMPS_FIELD_VELX;
00587 #if vmdplugin_ABIVERSION > 10
00588         data->coord_data |= LAMMPS_COORD_VELOCITIES;
00589         data->ts_meta.has_velocities = 1;
00590 #endif
00591       } else if (0 == strcmp(k, "vy")) {
00592         data->field[i] = LAMMPS_FIELD_VELY;
00593 #if vmdplugin_ABIVERSION > 10
00594         data->coord_data |= LAMMPS_COORD_VELOCITIES;
00595         data->ts_meta.has_velocities = 1;
00596 #endif
00597       } else if (0 == strcmp(k, "vz")) {
00598         data->field[i] = LAMMPS_FIELD_VELZ;
00599 #if vmdplugin_ABIVERSION > 10
00600         data->coord_data |= LAMMPS_COORD_VELOCITIES;
00601         data->ts_meta.has_velocities = 1;
00602 #endif
00603       } else if (0 == strcmp(k, "fx")) {
00604         data->field[i] = LAMMPS_FIELD_FORX;
00605         data->coord_data |= LAMMPS_COORD_FORCES;
00606       } else if (0 == strcmp(k, "fy")) {
00607         data->field[i] = LAMMPS_FIELD_FORY;
00608         data->coord_data |= LAMMPS_COORD_FORCES;
00609       } else if (0 == strcmp(k, "fz")) {
00610         data->field[i] = LAMMPS_FIELD_FORZ;
00611         data->coord_data |= LAMMPS_COORD_FORCES;
00612       } else if (0 == strcmp(k, "q")) {
00613         data->field[i] = LAMMPS_FIELD_CHARGE;
00614         *optflags |= MOLFILE_CHARGE; 
00615       } else if (0 == strcmp(k, "radius")) {
00616         data->field[i] = LAMMPS_FIELD_RADIUS;
00617         *optflags |= MOLFILE_RADIUS; 
00618       } else if (0 == strcmp(k, "diameter")) {
00619         data->field[i] = LAMMPS_FIELD_RADIUS;
00620         *optflags |= MOLFILE_RADIUS; 
00621       } else if (0 == strcmp(k, "element")) {
00622         data->field[i] = LAMMPS_FIELD_ELEMENT;
00623         *optflags |= MOLFILE_ATOMICNUMBER; 
00624         *optflags |= MOLFILE_MASS; 
00625         *optflags |= MOLFILE_RADIUS; 
00626       } else if (0 == strcmp(k, "mass")) {
00627         data->field[i] = LAMMPS_FIELD_MASS;
00628         *optflags |= MOLFILE_MASS; 
00629       } else if (0 == strcmp(k, "mux")) {
00630         data->field[i] = LAMMPS_FIELD_MUX;
00631         data->coord_data |= LAMMPS_COORD_DIPOLE;
00632       } else if (0 == strcmp(k, "muy")) {
00633         data->field[i] = LAMMPS_FIELD_MUY;
00634         data->coord_data |= LAMMPS_COORD_DIPOLE;
00635       } else if (0 == strcmp(k, "muz")) {
00636         data->field[i] = LAMMPS_FIELD_MUZ;
00637         data->coord_data |= LAMMPS_COORD_DIPOLE;
00638       } else {
00639         data->field[i] = LAMMPS_FIELD_UNKNOWN;
00640       }
00641       ++i;
00642       data->numfields = i;
00643       k = strtok(NULL," \t\n\r");
00644     } while ((k != NULL) && (i < LAMMPS_MAX_NUM_FIELDS));
00645   
00646     vmdcon_printf(VMDCON_INFO, "lammpsplugin) New style dump with %d data "
00647                   "fields. Coordinate data flags: 0x%02x\n",
00648                   data->numfields, data->coord_data);
00649     
00650     if ( !(data->coord_data & LAMMPS_COORD_DIPOLE) && (data->dip2atoms >= 0.0f)) {
00651       vmdcon_printf(VMDCON_ERROR, "lammpsplugin) conversion of dipoles to "
00652                     "two atoms requested, but no dipole data found\n");
00653       free(idlist);
00654       return MOLFILE_ERROR;
00655     }
00656   }
00657 
00658   idlist = (int *)malloc(data->numatoms * sizeof(int));
00659 
00660   /* read and parse ATOMS data section to build idlist */
00661   for(i=0; i<data->numatoms; i++) {
00662     k = myFgets(buffer, LINE_LEN, data->file);
00663 
00664     if (k == NULL) { 
00665       vmdcon_printf(VMDCON_ERROR, "lammpsplugin) Error while reading "
00666                     "structure from lammps dump file '%s': atom missing in "
00667                     "the first timestep\n", data->file_name);
00668       vmdcon_printf(VMDCON_ERROR, "lammpsplugin) expecting '%d' atoms, "
00669                     "found only '%d'\n", data->numatoms, i+1);
00670       free(idlist);
00671       return MOLFILE_ERROR;
00672     }
00673 
00674     /* if we have an old-style trajectory we have to guess what is there.
00675      * this chunk of code should only be executed once. LAMMPS_COORD_UNKNOWN
00676      * will be kept set until the very end or when we find that one position
00677      * is outside the box. */
00678     if (data->coord_data == LAMMPS_COORD_UNKNOWN) {
00679       int ix, iy, iz;
00680       j = sscanf(buffer, "%d%d%f%f%f%d%d%d", &atomid, &atomtype, 
00681                  &x, &y, &z, &ix, &iy, &iz);
00682       if (j > 4) {  /* assume id type xs ys zs .... format */
00683         data->coord_data |= LAMMPS_COORD_SCALED;
00684         data->numfields = 5;
00685         data->field[0] = LAMMPS_FIELD_ATOMID;
00686         data->field[1] = LAMMPS_FIELD_TYPE;
00687         data->field[2] = LAMMPS_FIELD_POSXS;
00688         data->field[3] = LAMMPS_FIELD_POSYS;
00689         data->field[4] = LAMMPS_FIELD_POSZS;
00690       } else {
00691         vmdcon_printf(VMDCON_ERROR, "lammpsplugin) Error while reading "
00692                       "structure from lammps dump file '%s'. Unsupported "
00693                       "dump file format.\n", data->file_name);
00694         free(idlist);
00695         return MOLFILE_ERROR;
00696       }
00697     }
00698 
00699     atomid = i; /* default, if no atomid given */
00700 
00701     /* parse the line of data and find the atom id */
00702     j = 0;
00703     k = strtok(buffer, " \t\n\r");
00704     while ((k != NULL) && (j < data->numfields)) {
00705       switch (data->field[j]) {
00706 
00707         case LAMMPS_FIELD_ATOMID:
00708           atomid = atoi(k) - 1; /* convert to 0 based list */
00709           break;
00710           
00711         default:
00712           break; /* ignore the rest */
00713       }
00714       ++j;
00715       k = strtok(NULL, " \t\n\r");
00716     }
00717     if (data->dip2atoms < 0.0f) {
00718       idlist[i] = atomid;
00719     } else {
00720       /* increment for fake second atom */
00721       idlist[i] =  2*atomid;
00722       ++i;
00723       idlist[i] =  2*atomid+1;
00724     }
00725 
00726     /* for old-style files, we have to use some heuristics to determine
00727      * if we have scaled or unscaled (absolute coordinates).
00728      * we assume scaled unless proven differently, and we assume unwrapped
00729      * unless we have images present. */
00730     if ( (data->coord_data & LAMMPS_COORD_UNKNOWN) != 0) {
00731       x=y=z=0.0f;
00732       j = sscanf(buffer, "%*d%*d%f%f%f", &x, &y, &z);
00733       if ((x<-0.1) || (x>1.1) || (y<-0.1) || (y>1.1) 
00734           || (z<-0.1) || (x>1.1)) {
00735         data->coord_data &= ~LAMMPS_COORD_UNKNOWN;
00736         if ((data->coord_data & LAMMPS_COORD_IMAGES) != 0) {
00737           data->coord_data |= LAMMPS_COORD_WRAPPED;
00738           data->field[2] = LAMMPS_FIELD_POSX;
00739           data->field[3] = LAMMPS_FIELD_POSY;
00740           data->field[4] = LAMMPS_FIELD_POSZ;
00741         } else {
00742           data->coord_data |= LAMMPS_COORD_UNWRAPPED;
00743           data->field[2] = LAMMPS_FIELD_POSXU;
00744           data->field[3] = LAMMPS_FIELD_POSYU;
00745           data->field[4] = LAMMPS_FIELD_POSZU;
00746         }
00747       }
00748     }
00749   }
00750   data->coord_data &= ~LAMMPS_COORD_UNKNOWN;
00751 
00752   /* pick coordinate type that we want to read and disable the rest. 
00753      we want unwrapped > wrapped > scaled. */
00754   if (data->coord_data & LAMMPS_COORD_UNWRAPPED) {
00755     data->coord_data &= ~(LAMMPS_COORD_WRAPPED|LAMMPS_COORD_SCALED
00756                           |LAMMPS_COORD_IMAGES);
00757   } else if (data->coord_data & LAMMPS_COORD_WRAPPED) {
00758     data->coord_data &= ~LAMMPS_COORD_SCALED;
00759   } else if (!(data->coord_data & LAMMPS_COORD_SCALED)) {
00760     /* we don't have any proper coordinates, not even scaled: bail out. */
00761     vmdcon_printf(VMDCON_ERROR, "lammpsplugin) No usable coordinate data "
00762                   "found in lammps dump file '%s'.\n", data->file_name);
00763     return MOLFILE_ERROR;
00764   }
00765   
00766   if (data->coord_data & LAMMPS_COORD_SCALED) {
00767     vmdcon_printf(VMDCON_INFO, "lammpsplugin) Reconstructing atomic "
00768                   "coordinates from fractional coordinates and box vectors.\n");
00769   } else {
00770     vmdcon_printf(VMDCON_INFO, "lammpsplugin) Using absolute atomic "
00771                   "coordinates directly.\n");
00772   }
00773   if (data->coord_data & LAMMPS_COORD_DIPOLE) {
00774     vmdcon_printf(VMDCON_INFO, "lammpsplugin) Detected dipole vector data.\n");
00775   }
00776   
00777   if (data->coord_data & LAMMPS_COORD_TRICLINIC) {
00778     vmdcon_printf(VMDCON_INFO, "lammpsplugin) Detected triclinic box.\n");
00779   }
00780   
00781 #if vmdplugin_ABIVERSION > 10
00782   if (data->coord_data & LAMMPS_COORD_VELOCITIES) {
00783     vmdcon_printf(VMDCON_INFO, "lammpsplugin) Importing atomic velocities.\n");
00784   }
00785 #endif
00786 
00787   /* sort list of atomids and figure out if we need the hash table */
00788   id_sort(idlist, data->numatoms);
00789   needhash=0;
00790   for (i=0; i < data->numatoms; ++i)
00791     if (idlist[i] != i) needhash=1;
00792 
00793   /* set up an integer hash to keep a sorted atom id map */
00794   if (needhash) {
00795     vmdcon_printf(VMDCON_INFO, "lammpsplugin) Using hash table to track "
00796                   "atom identities.\n");
00797     data->idmap = (inthash_t *)calloc(1, sizeof(inthash_t));
00798     inthash_init(data->idmap, data->numatoms);
00799     for (i=0; i < data->numatoms; ++i) {
00800       if (inthash_insert(data->idmap, idlist[i], i) != HASH_FAIL) {
00801         vmdcon_printf(VMDCON_ERROR, "lammpsplugin) Duplicate atomid %d or "
00802                       "unsupported dump file format.\n", idlist[i]);
00803         free(idlist);
00804         return MOLFILE_ERROR;
00805       }
00806     }
00807   } else {
00808     data->idmap = NULL;
00809   }
00810   free(idlist);
00811   
00812 
00813   /* now go back to the beginning of the file to parse it properly. */
00814   myRewind(data->file); 
00815 
00816   /* find the sections with atoms */
00817   fieldlist = find_item_keyword(data->file, KEY_DATA, buffer, LINE_LEN);
00818   if (fieldlist == NULL) {
00819     vmdcon_printf(VMDCON_ERROR, "lammpsplugin) Couldn't find data to "
00820                   "read structure from file '%s'.\n", data->file_name);
00821     return MOLFILE_ERROR;
00822   }
00823 
00824   /* read and parse ATOMS data section */
00825   for(i=0; i<data->numatoms; i++) {
00826     int has_element, has_mass, has_radius;
00827     
00828     k = myFgets(buffer, LINE_LEN, data->file);
00829 
00830     if (k == NULL) { 
00831       vmdcon_printf(VMDCON_ERROR, "lammpsplugin) Error while reading "
00832                     "structure from lammps dump file '%s': atom missing in "
00833                     "the first timestep\n", data->file_name);
00834       vmdcon_printf(VMDCON_ERROR, "lammpsplugin) expecting '%d' atoms, "
00835                     "found only '%d'\n", data->numatoms, i+1);
00836       free(idlist);
00837       return MOLFILE_ERROR;
00838     }
00839 
00840     /* some defaults */
00841     memset(&thisatom, 0, sizeof(molfile_atom_t)); 
00842     thisatom.resid = 0; /* mapped to MolID, if present */
00843     strncpy(thisatom.resname, "UNK", 4);
00844     strncpy(thisatom.chain, "",1);
00845     strncpy(thisatom.segid, "",1);
00846     atomid = i; /* needed if there is no atomid in a custom dump. */
00847     has_element = has_mass = has_radius = 0;
00848 
00849     /* parse the line of data */
00850     j = 0;
00851     k = strtok(buffer, " \t\n\r");
00852     while ((k != NULL) && (j < data->numfields)) {
00853       int idx;
00854 
00855       switch (data->field[j]) {
00856 
00857         case LAMMPS_FIELD_ATOMID:
00858           atomid = atoi(k) - 1; /* convert to 0 based list */
00859           break;
00860 
00861         case LAMMPS_FIELD_ELEMENT:
00862           strncpy(thisatom.name, k, 16);
00863           thisatom.name[15] = '\0';
00864           idx = get_pte_idx(k);
00865           thisatom.atomicnumber = idx;
00866           has_element = 1;
00867           break;
00868 
00869         case LAMMPS_FIELD_TYPE:
00870           strncpy(thisatom.type, k, 16); 
00871           if (has_element == 0) {
00872             /* element label has preference for name */
00873             strncpy(thisatom.name, k, 16);
00874             thisatom.type[15] = '\0';
00875           }
00876           /* WARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNING *\
00877            * Don't try using the atomid as name. This will waste a _lot_ of  *
00878            * memory due to names being stored in a string hashtable.         *
00879            * VMD currently cannot handle changing atomids anyways. We thus   *
00880            * use a hash table to track atom ids. Within VMD the atomids are  *
00881            * then lost, but atoms can be identified uniquely via 'serial'    *
00882            * or 'index' atom properties.
00883            * WARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNING*/
00884           break;
00885           
00886         case LAMMPS_FIELD_MOLID:
00887           thisatom.resid = atoi(k);
00888           break;
00889 
00890         case LAMMPS_FIELD_CHARGE:
00891           thisatom.charge = atof(k);
00892           break;
00893 
00894         case LAMMPS_FIELD_MASS:
00895           thisatom.mass = atof(k);
00896           has_mass = 1;
00897           break;
00898 
00899         case LAMMPS_FIELD_RADIUS:
00900           thisatom.radius = atof(k);
00901           has_radius = 1;
00902           break;
00903 
00904         case LAMMPS_FIELD_DIAMETER:
00905           /* radius has preference over diameter */
00906           if (has_radius == 0) {
00907             thisatom.radius = 0.5*atof(k);
00908             has_radius = 1;
00909           }
00910           break;
00911 
00912         case LAMMPS_FIELD_UNKNOWN: /* fallthrough */
00913         default:
00914           break;                /* do nothing */
00915       }
00916       ++j;
00917       k = strtok(NULL, " \t\n\r");
00918     }
00919 
00920     /* guess missing data from element name */
00921     if (has_element) {
00922       int idx = get_pte_idx(thisatom.name);
00923       if (!has_mass)   thisatom.mass   = get_pte_mass(idx);
00924       if (!has_radius) thisatom.radius = get_pte_vdw_radius(idx);
00925     }
00926 
00927     /* find position of this atom in the global list and copy its data */
00928     if (data->dip2atoms < 0.0f) {
00929       if (data->idmap != NULL) {
00930         j = inthash_lookup(data->idmap, atomid);
00931       } else {
00932         j = atomid;
00933       }
00934       memcpy(atoms+j, &thisatom, sizeof(molfile_atom_t)); 
00935     } else {
00936       if (data->idmap != NULL) {
00937         j = inthash_lookup(data->idmap, 2*atomid);
00938       } else {
00939         j = 2*atomid;
00940       }
00941       memcpy(atoms+j, &thisatom, sizeof(molfile_atom_t)); 
00942       /* the fake second atom has the same data */
00943       ++i;
00944       if (data->idmap != NULL) {
00945         j = inthash_lookup(data->idmap, 2*atomid+1);
00946       } else {
00947         j = 2*atomid+1;
00948       }
00949       memcpy(atoms+j, &thisatom, sizeof(molfile_atom_t)); 
00950     }
00951   }
00952 
00953   myRewind(data->file);
00954   data->fieldinit = 1;
00955   return MOLFILE_SUCCESS;
00956 }
00957 
00958 #if vmdplugin_ABIVERSION > 10
00959 /***********************************************************/
00960 static int read_timestep_metadata(void *mydata,
00961                                   molfile_timestep_metadata_t *meta) {
00962   lammpsdata *data = (lammpsdata *)mydata;
00963   
00964   meta->count = -1;
00965   meta->has_velocities = data->ts_meta.has_velocities;
00966   if (meta->has_velocities) {
00967     vmdcon_printf(VMDCON_INFO, "lammpsplugin) Importing velocities from "
00968                       "custom LAMMPS dump file.\n");
00969   }
00970   return MOLFILE_SUCCESS;
00971 }
00972 #endif
00973 
00974 /* convert cosine of angle to degrees. 
00975    bracket within -1.0 <= x <= 1.0 to avoid NaNs
00976    due to rounding errors. */
00977 static float cosangle2deg(double val)
00978 {
00979   if (val < -1.0) val=-1.0;
00980   if (val >  1.0) val= 1.0;
00981   return (float) (90.0 - asin(val)*90.0/M_PI_2);
00982 }
00983 
00984 static int read_lammps_timestep(void *mydata, int natoms, molfile_timestep_t *ts) {
00985   int i, j;
00986   char buffer[LINE_LEN];
00987   float x, y, z, vx, vy, vz;
00988   int atomid, numres, numatoms;
00989   float xlo, xhi, ylo, yhi, zlo, zhi, xy, xz, yz, ylohi;
00990   float xlo_bound, xhi_bound, ylo_bound, yhi_bound, zlo_bound, zhi_bound;
00991 
00992   lammpsdata *data = (lammpsdata *)mydata;
00993   /* we need to read/parse the structure information,
00994    * even if we only want to read coordinates later */
00995   if (data->fieldinit == 0) {
00996     molfile_atom_t *atoms;
00997     atoms = (molfile_atom_t *)malloc(natoms*sizeof(molfile_atom_t));
00998     read_lammps_structure(mydata,&natoms,atoms);
00999     free(atoms);
01000   }
01001 
01002   /* check if there is another time step in the file. */
01003   if (NULL == find_item_keyword(data->file, KEY_TSTEP, buffer, LINE_LEN)) 
01004     return MOLFILE_ERROR;
01005  
01006   /* check if we should read or skip this step. */
01007   if (!ts) return MOLFILE_SUCCESS;
01008 
01009   /* search for the number of atoms in the timestep */
01010   if (NULL==find_item_keyword(data->file, KEY_ATOMS, buffer, LINE_LEN)) {
01011     vmdcon_printf(VMDCON_ERROR, "lammpsplugin) Unable to find item: %s for "
01012                   "current timestep in file %s.\n", KEY_ATOMS, data->file_name);
01013     return MOLFILE_ERROR;
01014   }
01015 
01016   if (!myFgets(buffer, LINE_LEN, data->file)) {
01017     vmdcon_printf(VMDCON_ERROR, "lammpsplugin) Premature EOF for %s.\n", data->file_name);
01018     return MOLFILE_ERROR;
01019   }
01020 
01021   /* check if we have sufficient storage for atoms in this frame */
01022   numatoms = atoi(buffer);
01023   data->numatoms = numatoms;
01024   if (data->dip2atoms < 0.0f) {
01025     if (natoms < numatoms) {
01026       vmdcon_printf(VMDCON_ERROR, "lammpsplugin) Too many atoms in timestep."
01027                     " %d vs. %d\n", numatoms, natoms);
01028       return MOLFILE_ERROR;
01029     }
01030   } else {
01031     if (natoms/2 < numatoms) {
01032       vmdcon_printf(VMDCON_ERROR, "lammpsplugin) Too many atoms in timestep."
01033                     " %d vs. %d\n", numatoms, natoms/2);
01034       return MOLFILE_ERROR;
01035     }
01036   }
01037 
01038   /* now read the boundary box of the timestep */
01039   if (NULL == find_item_keyword(data->file, KEY_BOX, buffer, LINE_LEN)) {
01040     vmdcon_printf(VMDCON_ERROR, "lammpsplugin) Could not find box boundaries in timestep.\n");
01041     return MOLFILE_ERROR;
01042   }
01043 
01044   if (NULL == myFgets(buffer, LINE_LEN, data->file)) return MOLFILE_ERROR;
01045   numres = sscanf(buffer,"%f%f%f", &xlo_bound, &xhi_bound, &xy);
01046 
01047   if (NULL == myFgets(buffer, LINE_LEN, data->file)) return MOLFILE_ERROR;
01048   numres += sscanf(buffer,"%f%f%f", &ylo_bound, &yhi_bound, &xz);
01049 
01050   if (NULL == myFgets(buffer, LINE_LEN, data->file)) return MOLFILE_ERROR;
01051   numres += sscanf(buffer,"%f%f%f", &zlo_bound, &zhi_bound, &yz);
01052 
01053   xlo = xlo_bound;
01054   xhi = xhi_bound;
01055   ylo = ylo_bound;
01056   yhi = yhi_bound;
01057   zlo = zlo_bound;
01058   zhi = zhi_bound;
01059 
01060   if (data->coord_data & LAMMPS_COORD_TRICLINIC) {
01061     float xdelta;
01062 
01063     /* triclinic box. conveniently, LAMMPS makes the same assumptions
01064        about the orientation of the simulation cell than VMD. so
01065        hopefully no coordinate rotations or translations will be required */
01066 
01067     if (numres != 9) {
01068       vmdcon_printf(VMDCON_ERROR, "lammpsplugin) Inconsistent triclinic box specifications.\n");
01069       return MOLFILE_ERROR;
01070     }
01071 
01072     /* adjust box bounds */
01073     xdelta = MIN(0.0f,xy);
01074     xdelta = MIN(xdelta,xz);
01075     xdelta = MIN(xdelta,xy+xz);
01076     xlo -= xdelta;
01077 
01078     xdelta = MAX(0.0f,xy);
01079     xdelta = MAX(xdelta,xz);
01080     xdelta = MAX(xdelta,xy+xz);
01081     xhi -= xdelta;
01082 
01083     ylo -= MIN(0.0f,yz);
01084     yhi -= MAX(0.0f,yz);
01085     
01086   } else {
01087     if (numres != 6) {
01088       vmdcon_printf(VMDCON_ERROR, "lammpsplugin) Inconsistent orthorhombic box specifications.\n");
01089       return MOLFILE_ERROR;
01090     }
01091     xy = 0.0f;
01092     xz = 0.0f;
01093     yz = 0.0f;
01094   }
01095 
01096   /* convert bounding box info to real box lengths */
01097   ts->A = xhi-xlo;
01098   ylohi = yhi-ylo;
01099   ts->B = sqrt(ylohi*ylohi + xy*xy);
01100   ts->C = sqrt((zhi-zlo)*(zhi-zlo) + xz*xz + yz*yz);
01101 
01102   /* compute angles from box lengths and tilt factors */
01103   ts->alpha = cosangle2deg((xy*xz + ylohi*yz)/(ts->B * ts->C));
01104   ts->beta  = cosangle2deg(xz/ts->C);
01105   ts->gamma = cosangle2deg(xy/ts->B);
01106 
01107   /* read the coordinates */
01108   if (NULL == find_item_keyword(data->file, KEY_DATA, buffer, LINE_LEN)) {
01109     vmdcon_printf(VMDCON_ERROR, "lammpsplugin) could not find atom data for timestep.\n");
01110     return MOLFILE_ERROR;
01111   }
01112 
01113   /* initialize all coordinates to dummy location */
01114   for (i=0; i<natoms; i++) {
01115     ts->coords[3*i+0] = data->dumx;
01116     ts->coords[3*i+1] = data->dumy;
01117     ts->coords[3*i+2] = data->dumz;
01118   }
01119 
01120 #if vmdplugin_ABIVERSION > 10
01121   /* initialize velocities to zero if present */
01122   if (ts->velocities != NULL) memset(ts->velocities,0,3*natoms*sizeof(float));
01123 #endif
01124 
01125   for (i=0; i<numatoms; i++) {
01126     float ix, iy, iz, mux, muy, muz;
01127     char *k;
01128     
01129     k = myFgets(buffer, LINE_LEN, data->file);
01130 
01131     if (k == NULL) { 
01132       vmdcon_printf(VMDCON_ERROR, "lammpsplugin) Error while reading "
01133                     "data from lammps dump file '%s'.\n", data->file_name);
01134       vmdcon_printf(VMDCON_ERROR, "lammpsplugin) expecting '%d' atoms, "
01135                     "found only '%d'\n", numatoms, i+1);
01136       return MOLFILE_ERROR;
01137     }
01138 
01139     x=y=z=ix=iy=iz=vx=vy=vz=mux=muy=muz=0.0f;
01140     atomid=i;
01141     
01142     /* parse the line of data */
01143     j = 0;
01144     k = strtok(buffer, " \t\n\r");
01145     while ((k != NULL) && (j < data->numfields)) {
01146       switch (data->field[j]) {
01147 
01148         case LAMMPS_FIELD_ATOMID:
01149           atomid = atoi(k) - 1;
01150           break;
01151 
01152         case LAMMPS_FIELD_POSX:
01153           if (data->coord_data & LAMMPS_COORD_WRAPPED)
01154             x = atof(k);
01155           break;
01156 
01157         case LAMMPS_FIELD_POSY:
01158           if (data->coord_data & LAMMPS_COORD_WRAPPED)
01159             y = atof(k);
01160           break;
01161 
01162         case LAMMPS_FIELD_POSZ:
01163           if (data->coord_data & LAMMPS_COORD_WRAPPED)
01164             z = atof(k);
01165           break;
01166 
01167         case LAMMPS_FIELD_POSXU:
01168           if (data->coord_data & LAMMPS_COORD_UNWRAPPED)
01169             x = atof(k);
01170           break;
01171 
01172         case LAMMPS_FIELD_POSYU:
01173           if (data->coord_data & LAMMPS_COORD_UNWRAPPED)
01174             y = atof(k);
01175           break;
01176 
01177         case LAMMPS_FIELD_POSZU:
01178           if (data->coord_data & LAMMPS_COORD_UNWRAPPED)
01179             z = atof(k);
01180           break;
01181 
01182         case LAMMPS_FIELD_POSXS:
01183           if (data->coord_data & LAMMPS_COORD_SCALED)
01184             x = atof(k);
01185           break;
01186 
01187         case LAMMPS_FIELD_POSYS:
01188           if (data->coord_data & LAMMPS_COORD_SCALED)
01189             y = atof(k);
01190           break;
01191 
01192         case LAMMPS_FIELD_POSZS:
01193           if (data->coord_data & LAMMPS_COORD_SCALED)
01194             z = atof(k);
01195           break;
01196 
01197         case LAMMPS_FIELD_IMGX:
01198           if (data->coord_data & LAMMPS_COORD_IMAGES)
01199             ix = atof(k);
01200           break;
01201 
01202         case LAMMPS_FIELD_IMGY:
01203           if (data->coord_data & LAMMPS_COORD_IMAGES)
01204             iy = atof(k);
01205           break;
01206 
01207         case LAMMPS_FIELD_IMGZ:
01208           if (data->coord_data & LAMMPS_COORD_IMAGES)
01209             iz = atof(k);
01210           break;
01211 
01212         case LAMMPS_FIELD_MUX:
01213           if (data->coord_data & LAMMPS_COORD_DIPOLE)
01214             mux = atof(k);
01215           break;
01216 
01217         case LAMMPS_FIELD_MUY:
01218           if (data->coord_data & LAMMPS_COORD_DIPOLE)
01219             muy = atof(k);
01220           break;
01221 
01222         case LAMMPS_FIELD_MUZ:
01223           if (data->coord_data & LAMMPS_COORD_DIPOLE)
01224             muz = atof(k);
01225           break;
01226 
01227 #if vmdplugin_ABIVERSION > 10
01228         case LAMMPS_FIELD_VELX:
01229           vx = atof(k);
01230           break;
01231 
01232         case LAMMPS_FIELD_VELY:
01233           vy = atof(k);
01234           break;
01235 
01236         case LAMMPS_FIELD_VELZ:
01237           vz = atof(k);
01238           break;
01239 #endif
01240 
01241         default: /* do nothing */
01242           break;
01243       }
01244 
01245       ++j;
01246       k = strtok(NULL, " \t\n\r");
01247     } 
01248     
01249     if (data->dip2atoms < 0.0f) {
01250       if (data->idmap != NULL) {
01251         j = inthash_lookup(data->idmap, atomid);
01252       } else {
01253         j = atomid;
01254       }
01255     } else {
01256       if (data->idmap != NULL) {
01257         j = inthash_lookup(data->idmap, 2*atomid);
01258       } else {
01259         j = 2*atomid;
01260       }
01261     }
01262     
01263     if (data->idmap && (j == HASH_FAIL)) {
01264       /* we have space in the hash table. add this atom */
01265       if (inthash_entries(data->idmap) < data->maxatoms) {
01266         if (data->dip2atoms < 0.0f) {
01267           inthash_insert(data->idmap,atomid,i);
01268           j = inthash_lookup(data->idmap,atomid);
01269         } else {
01270           inthash_insert(data->idmap,2*atomid,2*i);
01271           inthash_insert(data->idmap,2*atomid+1,2*i+1);
01272           j = inthash_lookup(data->idmap, 2*atomid);
01273         }
01274       } else j = -1;
01275     }
01276 
01277     if ((j < 0) || (j >= data->maxatoms)) {
01278       vmdcon_printf(VMDCON_WARN, "lammpsplugin) ignoring out of range "
01279                     "atom #%d with id %d\n", i, atomid);
01280     } else {
01281       /* copy coordinates. we may have coordinates in different
01282        * formats available in custom dumps. those have been checked
01283        * before and we prefer to use unwrapped > wrapped > scaled 
01284        * in this order. in the second two cases, we also apply image
01285        * shifts, if that data is available. unnecessary or unsupported
01286        * combinations of flags have been cleared based on data in 
01287        * the first frame in read_lammps_structure(). */
01288       int addr = 3 * j;
01289       if (data->coord_data & LAMMPS_COORD_TRICLINIC) {
01290         if (data->coord_data & LAMMPS_COORD_SCALED) {
01291           /* we have fractional coordinates, so they need 
01292            * to be scaled accordingly. */
01293           ts->coords[addr    ] = xlo + x * ts->A + y * xy + z * xz;
01294           ts->coords[addr + 1] = ylo + y * ylohi + z * yz;
01295           ts->coords[addr + 2] = zlo + z * (zhi-zlo);
01296         } else {
01297           /* ... but they can also be absolute values */
01298           ts->coords[addr    ] = x;
01299           ts->coords[addr + 1] = y;
01300           ts->coords[addr + 2] = z;
01301         }
01302         if (data->coord_data & LAMMPS_COORD_IMAGES) {
01303           /* we have image counter data to unwrap coordinates. */
01304           ts->coords[addr    ] += ix * ts->A + iy * xy + iz * xz;
01305           ts->coords[addr + 1] += iy * ylohi + iz * yz;
01306           ts->coords[addr + 2] += iz * (zhi-zlo);
01307         }
01308       } else {
01309         if (data->coord_data & LAMMPS_COORD_SCALED) {
01310           /* we have fractional coordinates, so they need 
01311            * to be scaled by a/b/c etc. */
01312           ts->coords[addr    ] = xlo + x * ts->A;
01313           ts->coords[addr + 1] = ylo + y * ts->B;
01314           ts->coords[addr + 2] = zlo + z * ts->C;
01315         } else {
01316           /* ... but they can also be absolute values */
01317           ts->coords[addr    ] = x;
01318           ts->coords[addr + 1] = y;
01319           ts->coords[addr + 2] = z;
01320         }
01321         if (data->coord_data & LAMMPS_COORD_IMAGES) {
01322           /* we have image counter data to unwrap coordinates. */
01323           ts->coords[addr    ] += ix * ts->A;
01324           ts->coords[addr + 1] += iy * ts->B;
01325           ts->coords[addr + 2] += iz * ts->C;
01326         }
01327       }
01328         
01329 #if vmdplugin_ABIVERSION > 10
01330       if (ts->velocities != NULL) {
01331         ts->velocities[addr    ] = vx;
01332         ts->velocities[addr + 1] = vy;
01333         ts->velocities[addr + 2] = vz;
01334       }
01335 #endif
01336 
01337       /* translate and copy atom positions for dipole data */
01338       if (data->dip2atoms >= 0.0f) {
01339         const float sf = data->dip2atoms;
01340         
01341         x = ts->coords[addr    ];
01342         y = ts->coords[addr + 1];
01343         z = ts->coords[addr + 2];
01344         ts->coords[addr    ] = x + sf*mux;
01345         ts->coords[addr + 1] = y + sf*muy;
01346         ts->coords[addr + 2] = z + sf*muz;
01347 
01348         if (data->idmap != NULL) {
01349           j = inthash_lookup(data->idmap, 2*atomid + 1);
01350         } else {
01351           j = 2*atomid + 1;
01352         }
01353         addr = 3*j;
01354         ts->coords[addr    ] = x - sf*mux;
01355         ts->coords[addr + 1] = y - sf*muy;
01356         ts->coords[addr + 2] = z - sf*muz;
01357 #if vmdplugin_ABIVERSION > 10
01358         if (ts->velocities != NULL) {
01359           ts->velocities[addr    ] = vx;
01360           ts->velocities[addr + 1] = vy;
01361           ts->velocities[addr + 2] = vz;
01362         }
01363 #endif
01364       }
01365     }
01366   }
01367 
01368   return MOLFILE_SUCCESS;
01369 }
01370     
01371 static void close_lammps_read(void *mydata) {
01372   lammpsdata *data = (lammpsdata *)mydata;
01373   myFclose(data->file);
01374   free(data->file_name);
01375 #if LAMMPS_DEBUG
01376   if (data->idmap != NULL) 
01377     fprintf(stderr, "inthash stats: %s\n", inthash_stats(data->idmap));
01378 #endif
01379   if (data->idmap != NULL) {
01380     inthash_destroy(data->idmap);
01381     free(data->idmap);
01382   }
01383   free(data);
01384 }
01385 
01386 static void *open_lammps_write(const char *filename, const char *filetype, 
01387                            int natoms) {
01388   FILE *fp;
01389   lammpsdata *data;
01390 
01391   fp = fopen(filename, "w");
01392   if (!fp) { 
01393     vmdcon_printf(VMDCON_ERROR, "Error) Unable to open lammpstrj file %s for writing\n",
01394             filename);
01395     return NULL;
01396   }
01397   
01398   data = (lammpsdata *)malloc(sizeof(lammpsdata));
01399   data->numatoms = natoms;
01400   data->fp = fp;
01401   data->file_name = strdup(filename);
01402   data->nstep = 0;
01403   return data;
01404 }
01405 
01406 static int write_lammps_structure(void *mydata, int optflags, 
01407                                const molfile_atom_t *atoms) {
01408   lammpsdata *data = (lammpsdata *)mydata;
01409   int i, j;
01410   hash_t atomtypehash;
01411 
01412   hash_init(&atomtypehash,128);
01413 
01414   /* generate 1 based lookup table for atom types */
01415   for (i=0, j=1; i < data->numatoms; i++)
01416     if (hash_insert(&atomtypehash, atoms[i].type, j) == HASH_FAIL)
01417       j++;
01418   
01419   data->atomtypes = (int *) malloc(data->numatoms * sizeof(int));
01420 
01421   for (i=0; i < data->numatoms ; i++)
01422     data->atomtypes[i] = hash_lookup(&atomtypehash, atoms[i].type);
01423 
01424   hash_destroy(&atomtypehash);
01425   
01426   return MOLFILE_SUCCESS;
01427 }
01428 
01429 static int write_lammps_timestep(void *mydata, const molfile_timestep_t *ts) {
01430   lammpsdata *data = (lammpsdata *)mydata; 
01431   const float *pos;
01432   float xmin[3], xmax[3], xcen[3];
01433   int i, tric, pbcx, pbcy, pbcz;
01434 
01435   fprintf(data->fp, "ITEM: TIMESTEP\n");
01436   fprintf(data->fp, "%d\n", data->nstep);
01437   fprintf(data->fp, "ITEM: NUMBER OF ATOMS\n");
01438   fprintf(data->fp, "%d\n", data->numatoms);
01439 
01440   pos = ts->coords;
01441 
01442   xmax[0] = xmax[1] = xmax[2] = -1.0e30f;
01443   xmin[0] = xmin[1] = xmin[2] =  1.0e30f;
01444   tric = pbcx = pbcy = pbcz = 0;
01445 
01446 #if defined(_MSC_VER)
01447   if ((fabs(ts->alpha - 90.0f) > SMALL) ||
01448       (fabs(ts->beta  - 90.0f) > SMALL) ||
01449       (fabs(ts->gamma - 90.0f) > SMALL)) 
01450     tric = 1;
01451     if (fabs(ts->A > SMALL)) pbcx = 1;
01452     if (fabs(ts->B > SMALL)) pbcy = 1;
01453     if (fabs(ts->C > SMALL)) pbcz = 1;
01454 #else
01455   if ((fabsf(ts->alpha - 90.0f) > SMALL) ||
01456       (fabsf(ts->beta  - 90.0f) > SMALL) ||
01457       (fabsf(ts->gamma - 90.0f) > SMALL))
01458     tric = 1;
01459     if (fabsf(ts->A > SMALL)) pbcx = 1;
01460     if (fabsf(ts->B > SMALL)) pbcy = 1;
01461     if (fabsf(ts->C > SMALL)) pbcz = 1;
01462 #endif  
01463 
01464   /* find min/max coordinates to approximate lo/hi values */
01465   for (i = 0; i < data->numatoms; ++i) {
01466     xmax[0] = (pos[0] > xmax[0]) ? pos[0] : xmax[0];
01467     xmax[1] = (pos[1] > xmax[1]) ? pos[1] : xmax[1];
01468     xmax[2] = (pos[2] > xmax[2]) ? pos[2] : xmax[2];
01469     xmin[0] = (pos[0] < xmin[0]) ? pos[0] : xmin[0];
01470     xmin[1] = (pos[1] < xmin[1]) ? pos[1] : xmin[1];
01471     xmin[2] = (pos[2] < xmin[2]) ? pos[2] : xmin[2];
01472     pos += 3;
01473   }
01474   xcen[0] = 0.5f * (xmax[0] + xmin[0]);
01475   xcen[1] = 0.5f * (xmax[1] + xmin[1]);
01476   xcen[2] = 0.5f * (xmax[2] + xmin[2]);
01477 
01478   pos = ts->coords;
01479 
01480   if (tric == 0 ) {  /* orthogonal box */
01481 
01482     if (pbcx) xmax[0] = xcen[0] + 0.5f*ts->A;
01483     if (pbcx) xmin[0] = xcen[0] - 0.5f*ts->A;
01484     if (pbcy) xmax[1] = xcen[1] + 0.5f*ts->B;
01485     if (pbcy) xmin[1] = xcen[1] - 0.5f*ts->B;
01486     if (pbcz) xmax[2] = xcen[2] + 0.5f*ts->C;
01487     if (pbcz) xmin[2] = xcen[2] - 0.5f*ts->C;
01488 
01489     /* flag using PBC when box length exists, else shrinkwrap BC */
01490     fprintf(data->fp, "ITEM: BOX BOUNDS %s %s %s\n", pbcx ? "pp" : "ss",
01491             pbcy ? "pp" : "ss", pbcz ? "pp" : "ss");
01492     fprintf(data->fp, "%g %g\n", xmin[0], xmax[0]);
01493     fprintf(data->fp, "%g %g\n", xmin[1], xmax[1] );
01494     fprintf(data->fp, "%g %g\n", xmin[2], xmax[2] );
01495 
01496   } else { /* triclinic box */
01497 
01498     double lx, ly, lz, xy, xz, yz, xbnd;
01499 
01500     lx = ts->A;
01501     xy = ts->B * cos(ts->gamma/90.0*M_PI_2);
01502     xz = ts->C * cos(ts->beta/90.0*M_PI_2);
01503     ly = sqrt(ts->B*ts->B - xy*xy);
01504     if (fabs(ly) > SMALL) 
01505       yz = (ts->B*ts->C*cos(ts->alpha/90.0*M_PI_2) - xy*xz) / ly;
01506     else
01507       yz = 0.0;
01508     lz = sqrt(ts->C*ts->C - xz*xz - yz*yz);
01509 
01510     if (pbcx) xmax[0] = xcen[0] + 0.5f*lx;
01511     if (pbcx) xmin[0] = xcen[0] - 0.5f*lx;
01512     if (pbcy) xmax[1] = xcen[1] + 0.5f*ly;
01513     if (pbcy) xmin[1] = xcen[1] - 0.5f*ly;
01514     if (pbcz) xmax[2] = xcen[2] + 0.5f*lz;
01515     if (pbcz) xmin[2] = xcen[2] - 0.5f*lz;
01516 
01517     /* go from length to boundary */
01518 
01519     xbnd = 0.0;
01520     xbnd = (xy > xbnd) ? xy : xbnd;
01521     xbnd = (xz > xbnd) ? xz : xbnd;
01522     xbnd = (xy+xz > xbnd) ? (xy + xz) : xbnd;
01523     xmax[0] += xbnd;
01524 
01525     xbnd = 0.0;
01526     xbnd = (xy < xbnd) ? xy : xbnd;
01527     xbnd = (xz < xbnd) ? xz : xbnd;
01528     xbnd = (xy+xz < xbnd) ? (xy + xz) : xbnd;
01529     xmin[0] += xbnd;
01530     
01531     xbnd = 0.0;
01532     xbnd = (yz > xbnd) ? yz : xbnd;
01533     xmax[1] += xbnd;
01534 
01535     xbnd = 0.0;
01536     xbnd = (yz < xbnd) ? yz : xbnd;
01537     xmin[1] += xbnd;
01538 
01539     /* flag using PBC when box length exists, else shrinkwrap BC */
01540     fprintf(data->fp, "ITEM: BOX BOUNDS %s %s %s xy xz yz\n", pbcx ? "pp" : "ss",
01541             pbcy ? "pp" : "ss", pbcz ? "pp" : "ss");
01542     fprintf(data->fp, "%g %g %g\n", xmin[0], xmax[0], xy);
01543     fprintf(data->fp, "%g %g %g\n", xmin[1], xmax[1], xz);
01544     fprintf(data->fp, "%g %g %g\n", xmin[2], xmax[2], yz);
01545   }
01546   
01547   /* coordinates are written as unwrapped coordinates */
01548   fprintf(data->fp, "ITEM: ATOMS id type xu yu zu\n");
01549   for (i = 0; i < data->numatoms; ++i) {
01550     fprintf(data->fp, " %d %d %g %g %g\n", 
01551             i+1, data->atomtypes[i], pos[0], pos[1], pos[2]);
01552     pos += 3;
01553   }
01554 
01555   data->nstep ++;
01556   return MOLFILE_SUCCESS;
01557 }
01558 
01559 
01560 static void close_lammps_write(void *mydata) {
01561   lammpsdata *data = (lammpsdata *)mydata;
01562 
01563   fclose(data->fp);
01564   free(data->atomtypes);
01565   free(data->file_name);
01566   free(data);
01567 }
01568 
01569 
01570 /* registration stuff */
01571 static molfile_plugin_t plugin;
01572 
01573 VMDPLUGIN_API int VMDPLUGIN_init() {
01574   memset(&plugin, 0, sizeof(molfile_plugin_t));
01575   plugin.abiversion = vmdplugin_ABIVERSION;
01576   plugin.type = MOLFILE_PLUGIN_TYPE;
01577   plugin.name = "lammpstrj";
01578   plugin.prettyname = "LAMMPS Trajectory";
01579   plugin.author = "Marco Kalweit, Axel Kohlmeyer, Lutz Maibaum, John Stone";
01580   plugin.majorv = 0;
01581   plugin.minorv = 22;
01582   plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
01583 #ifdef _USE_ZLIB
01584   plugin.filename_extension = "lammpstrj,lammpstrj.gz";
01585 #else
01586   plugin.filename_extension = "lammpstrj";
01587 #endif
01588   plugin.open_file_read = open_lammps_read;
01589   plugin.read_structure = read_lammps_structure;
01590   plugin.read_next_timestep = read_lammps_timestep;
01591 #if vmdplugin_ABIVERSION > 10
01592   plugin.read_timestep_metadata    = read_timestep_metadata;
01593 #endif
01594   plugin.close_file_read = close_lammps_read;
01595   plugin.open_file_write = open_lammps_write;
01596   plugin.write_structure = write_lammps_structure;
01597   plugin.write_timestep = write_lammps_timestep;
01598   plugin.close_file_write = close_lammps_write;
01599 
01600   return VMDPLUGIN_SUCCESS;
01601 }
01602 
01603 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
01604   (*cb)(v, (vmdplugin_t *)&plugin);
01605   return VMDPLUGIN_SUCCESS;
01606 }
01607 
01608 VMDPLUGIN_API int VMDPLUGIN_fini() {
01609   return VMDPLUGIN_SUCCESS;
01610 }
01611 
01612 
01613 #ifdef TEST_PLUGIN
01614 
01615 int main(int argc, char *argv[]) {
01616   molfile_timestep_t timestep;
01617   molfile_atom_t *atoms = NULL;
01618   void *v;
01619   int natoms;
01620   int i, j, opts;
01621 #if vmdplugin_ABIVERSION > 10
01622   molfile_timestep_metadata_t ts_meta;
01623 #endif
01624 
01625   while (--argc >=0) {
01626     ++argv;
01627     v = open_lammps_read(*argv, "lammps", &natoms);
01628     if (!v) {
01629       fprintf(stderr, "open_lammps_read failed for file %s\n", *argv);
01630       return 1;
01631     }
01632     fprintf(stderr, "open_lammps_read succeeded for file %s\n", *argv);
01633     fprintf(stderr, "number of atoms: %d\n", natoms);
01634 
01635     timestep.coords = (float *)malloc(3*sizeof(float)*natoms);
01636     atoms = (molfile_atom_t *)malloc(sizeof(molfile_atom_t)*natoms);
01637     if (read_lammps_structure(v, &opts, atoms) == MOLFILE_ERROR) {
01638       close_lammps_read(v);
01639       continue;
01640     }
01641       
01642     fprintf(stderr, "read_lammps_structure: options=0x%08x\n", opts);
01643 #if 0
01644     for (i=0; i<natoms; ++i) {
01645       fprintf(stderr, "atom %09d: name=%s, type=%s, resname=%s, resid=%d, segid=%s, chain=%s\n",
01646                       i, atoms[i].name, atoms[i].type, atoms[i].resname, atoms[i].resid,
01647                       atoms[i].segid, atoms[i].chain);
01648     }
01649 #endif
01650 #if vmdplugin_ABIVERSION > 10
01651     read_timestep_metadata(v,&ts_meta);
01652     if (ts_meta.has_velocities) {
01653       fprintf(stderr, "found timestep velocities metadata.\n");
01654     }
01655     timestep.velocities = (float *) malloc(3*natoms*sizeof(float));
01656 #endif
01657     j = 0;
01658     while (!read_lammps_timestep(v, natoms, &timestep)) {
01659       for (i=0; i<10; ++i) {
01660         fprintf(stderr, "atom %09d: type=%s, resid=%d, "
01661                       "x/y/z = %.3f %.3f %.3f "
01662 #if vmdplugin_ABIVERSION > 10
01663                       "vx/vy/vz = %.3f %.3f %.3f "
01664 #endif
01665                       "\n",
01666                       i, atoms[i].type, atoms[i].resid, 
01667                       timestep.coords[3*i], timestep.coords[3*i+1], 
01668                       timestep.coords[3*i+2]
01669 #if vmdplugin_ABIVERSION > 10
01670                       ,timestep.velocities[3*i], timestep.velocities[3*i+1], 
01671                       timestep.velocities[3*i+2]
01672 #endif
01673                       );    
01674       }
01675       j++;
01676     }
01677     fprintf(stderr, "ended read_next_timestep on frame %d\n", j);
01678 
01679     close_lammps_read(v);
01680   }
01681 #if vmdplugin_ABIVERSION > 10
01682   free(timestep.velocities);
01683 #endif
01684   free(timestep.coords);
01685   free(atoms);
01686   return 0;
01687 }
01688 
01689 #endif

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