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

lammpsyamlplugin.c

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

Generated on Thu Nov 7 03:09:41 2024 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002