00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063 #include "largefiles.h"
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
00097 #define SMALL 1.0e-12f
00098
00099
00100 #define LINE_LEN 1024
00101
00102
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
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
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;
00175 float time;
00176 float dip2atoms;
00177 float dumx,dumy,dumz;
00178 int numfields;
00179 l_attr_t field[LAMMPS_MAX_NUM_FIELDS];
00180 inthash_t *idmap;
00181 int fieldinit;
00182 #if vmdplugin_ABIVERSION > 10
00183 molfile_timestep_metadata_t ts_meta;
00184 #endif
00185 } lammpsdata;
00186
00187 static int velinfo = 1;
00188
00189
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
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
00238 if ((!tag) || (!list))
00239 return tag;
00240
00241 ptr=list;
00242 i=pos=flag=0;
00243 len=strlen(list);
00244
00245
00246 while (pos < len) {
00247
00248 if (flag) {
00249 if (ptr[pos] == ',') {
00250 from[i] = '\0';
00251 if (0 == strcmp(tag,from)) {
00252
00253 if (strlen(to))
00254 return to;
00255
00256 flag=0;
00257 i=0;
00258 } else {
00259
00260 flag=0;
00261 i=0;
00262 }
00263 } else {
00264 from[i]=ptr[pos];
00265 if (i<30)
00266 ++i;
00267 }
00268 } else {
00269
00270 if (ptr[pos] == '=') {
00271 to[i] = '\0';
00272 flag=1;
00273 i=0;
00274 } else if (ptr[pos] == ',') {
00275 i=0;
00276 flag=0;
00277 } else {
00278 to[i]=ptr[pos];
00279 if (i<30)
00280 ++i;
00281 }
00282 }
00283 ++pos;
00284 }
00285
00286
00287 if (flag) {
00288 from[i] = '\0';
00289 if (0 == strcmp(tag,from)) {
00290
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
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;
00356 }
00357 return NULL;
00358 }
00359
00360
00361
00362
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
00401 do {
00402 if (!myFgets(buffer, LINE_LEN, fd)) return NULL;
00403 } while (strncmp(buffer, "---", 3) != 0);
00404
00405
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
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
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
00429
00430
00431
00432
00433
00434
00435
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
00449
00450
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;
00463 data->numatoms = tmp;
00464 data->coord_data = LAMMPS_COORD_NONE;
00465 myRewind(data->file);
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
00484 *optflags = MOLFILE_NOOPTIONS;
00485 data->coord_data = LAMMPS_COORD_NONE;
00486 memset(atoms, 0, data->numatoms * sizeof(molfile_atom_t));
00487
00488
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
00510 myRewind(data->file);
00511
00512
00513 do {
00514 if (!myFgets(buffer, LINE_LEN, data->file)) return MOLFILE_ERROR;
00515 } while (strncmp(buffer, "---", 3) != 0);
00516
00517
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
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
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
00546 k = strtok(data_start(fieldlist), ",] \t\n\r");
00547 i = 0;
00548 do {
00549
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
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
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
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
00712 k = strtok(fieldlist, ",] \t\r\n");
00713 atomid = i;
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;
00720 break;
00721
00722 default:
00723 break;
00724 }
00725 ++j;
00726 k = strtok(NULL, ",] \t\n\r");
00727 }
00728 if (data->dip2atoms < 0.0f) {
00729 idlist[i] = atomid;
00730 } else {
00731
00732 idlist[i] = 2*atomid;
00733 ++i;
00734 idlist[i] = 2*atomid+1;
00735 }
00736 }
00737 data->coord_data &= ~LAMMPS_COORD_UNKNOWN;
00738
00739
00740
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
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
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
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
00798 myRewind(data->file);
00799
00800
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
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
00825 memset(&thisatom, 0, sizeof(molfile_atom_t));
00826 thisatom.resid = 0;
00827 strncpy(thisatom.resname, "UNK", 4);
00828 strncpy(thisatom.chain, "",1);
00829 strncpy(thisatom.segid, "",1);
00830 atomid = i;
00831 has_element = has_mass = has_radius = 0;
00832
00833
00834 fieldlist = data_start(buffer);
00835
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;
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
00861 strncpy(thisatom.name, k, 16);
00862 thisatom.type[15] = '\0';
00863 }
00864
00865
00866
00867
00868
00869
00870
00871
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
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:
00901 default:
00902 break;
00903 }
00904 ++j;
00905 k = strtok(NULL, ",] \t\n\r");
00906 }
00907
00908
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
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
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
00962
00963
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
00982
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
00991 do {
00992 if (!myFgets(buffer, LINE_LEN, data->file)) return MOLFILE_ERROR;
00993 } while (strncmp(buffer, "---", 3) != 0);
00994
00995
00996 if (!ts) return MOLFILE_SUCCESS;
00997
00998
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
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
01021 if (strncmp(ptr, KEY_TIME, sizeof(KEY_TIME)-1) == 0)
01022 ts->physical_time = atof(buffer + i + 1);
01023 }
01024
01025
01026 numatoms = atol(buffer + i + 1);
01027 data->numatoms = numatoms;
01028
01029
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
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
01073
01074
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
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
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
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
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
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
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
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:
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
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
01293
01294
01295
01296
01297
01298
01299 int addr = 3 * j;
01300 if (data->coord_data & LAMMPS_COORD_TRICLINIC) {
01301 if (data->coord_data & LAMMPS_COORD_SCALED) {
01302
01303
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
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
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
01322
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
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
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
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
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
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 ) {
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
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 {
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
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
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
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
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, ×tep)) {
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