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

Generated on Fri Dec 6 03:09:47 2024 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002