Cluster.c

   1 /*
   2     This program is free software: you can redistribute it and/or modify
   3     it under the terms of the GNU General Public License as published by
   4     the Free Software Foundation, either version 3 of the License, or
   5     (at your option) any later version.
   6
   7     This program is distributed in the hope that it will be useful,
   8     but WITHOUT ANY WARRANTY; without even the implied warranty of
   9     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  10     GNU General Public License for more details.
  11
  12     You should have received a copy of the GNU General Public License
  13     along with this program.  If not, see <http://www.gnu.org/licenses/>.
  14
  15     (C)opyright   Gerolf Ziegenhain 2008 <gerolf@ziegenhain.com>
  16                   Christian Anders 2008 <anders@physik.uni-kl.de>
  17
  18     Based on the algorithm in S. Stoddard JCP(1978)27:291
  19
  20 */
  21
  22
  23 #include <stdio.h>
  24 #include <stdlib.h>
  25 #include <math.h>
  26 #include <time.h>
  27 #include <unistd.h>
  28 #include <stdarg.h>
  29 #include <string.h>
  30 #include <string.h>
  31
  32 #include "H5LT.h"
  33 #include "H5TA.h"
  34 #include <hdf5.h>
  35 #undef MAX
  36 #undef MIN
  37
  38
  39 int getopt(int argc, char * const argv[], const char *optstring);
  40 extern char *optarg;
  41 extern int optind, opterr, optopt;
  42 FILE *popen(const char *command, const char *type);
  43
  44
  45 #define EOK 0
  46 #define EMEM 1
  47 #define ESYN 2
  48 #define EFILE 3
  49
  50 #define X 0
  51 #define Y 1
  52 #define Z 2
  53
  54 #define TRUE 1
  55 #define FALSE 0
  56
  57 #define MAGIC_MIN 1e10
  58 #define MAGIC_MAX -1e10
  59
  60 #define MAX(a,b) (a>b)?a:b
  61 #define MIN(a,b) (a>b)?b:a
  62
  63 #define FLOAT float
  64 #define INT int
  65
  66 #define USAGE  "-h     help\n" \
  67                "-r     rcut\n" \
  68                "-n     use neighborlist (radius)\n" \
  69                "-l     use linked lists\n" \
  70                "-f     infile (name)\n" \
  71                "-XYZ   peridicity\n"
  72
  73 #define CLK 100000.
  74
  75 #define SQ(a) ((a)*(a))
  76 #define PROGRESS(n) if (n % 1000 == 0) llog (" %d", n)
  77
  78 /* --------------------------------------------------------------------------------------------- */
  79
  80 typedef struct s_linkedlist {
  81    unsigned int particle;
  82    struct s_linkedlist *next;
  83 } t_linkedlist;
  84
  85 /* --------------------------------------------------------------------------------------------- */
  86
  87 int Periodic[3];
  88 float Length[3];
  89 unsigned int NCellSide[3], NCell;
  90 float CellLength[3], CellLengthI[3];
  91 t_linkedlist **CellList;
  92 float min[3], max[3];
  93
  94 /* --------------------------------------------------------------------------------------------- */
  95
  96 void llog (const char *format, ...) {
  97    va_list arg;
  98    va_start (arg, format);
  99    vfprintf (stderr, format, arg);
 100    fflush (stderr);
 101    va_end (arg);
 102 }
 103
 104 /* --------------------------------------------------------------------------------------------- */
 105
 106 double SqDist (double pos1, double pos2, int dim) {
 107    if (Periodic[dim] == TRUE)
 108       return SQ( fmod(pos1-pos2,Length[dim]) );
 109    else
 110       return SQ(pos1-pos2);
 111 }
 112
 113 /* --------------------------------------------------------------------------------------------- */
 114
 115 void DetermineBoxLength (unsigned int n, float *x, float *y, float *z, float *lx, float *ly, float *lz) {
 116    llog ("\n\tDetermining box length: ");
 117    min[X] = min[Y] = min[Z] = MAGIC_MIN;
 118    max[X] = max[Y] = max[Z] = MAGIC_MAX;
 119    unsigned int i;
 120    for (i = 0; i < n; i++) {
 121       min[X] = MIN(min[X],x[i]);
 122       min[Y] = MIN(min[Y],y[i]);
 123       min[Z] = MIN(min[Z],z[i]);
 124       max[X] = MAX(max[X],x[i]);
 125       max[Y] = MAX(max[Y],y[i]);
 126       max[Z] = MAX(max[Z],z[i]);
 127    }
 128
 129    double off[3];
 130    off[X] = 0.;
 131    off[Y] = 0.;
 132    off[Z] = 0.;
 133    unsigned int ShiftFlag = FALSE;
 134    if (min[X] < 0.) {
 135       ShiftFlag = TRUE;
 136       off[X] = -min[X];
 137    }
 138
 139    if (min[Y] < 0.) {
 140       ShiftFlag = TRUE;
 141       off[Y] = -min[Y];
 142     }
 143
 144    if (min[Z] < 0.) {
 145       ShiftFlag = TRUE;
 146       off[Z] = -min[Z];
 147    }
 148
 149    if (ShiftFlag) {
 150       llog ("\n\tShifting atoms by %f %f %f", off[X], off[Y], off[Z]);
 151       for (i = 0; i < n; i++) {
 152          x[i] += off[X];
 153          y[i] += off[Y];
 154          z[i] += off[Z];
 155       }
 156       unsigned int j;
 157       for (j = 0; j < 3; j++) {
 158          min[j] += off[j];
 159          max[j] += off[j];
 160       }
 161    }
 162
 163    *(lx) = max[X] - min[X];
 164    *(ly) = max[Y] - min[Y];
 165    *(lz) = max[Z] - min[Z];
 166    llog ("\n\tBox dimensions: %f %f %f %f %f %f", min[X], max[X], min[Y], max[Y], min[Z], max[Z]);
 167 }
 168
 169 /* --------------------------------------------------------------------------------------------- */
 170
 171 void LinkedListInsert (unsigned int particle, t_linkedlist **list) {
 172    t_linkedlist *new = (t_linkedlist *)calloc (1, sizeof (t_linkedlist));
 173    if (!new)
 174       exit (EMEM);
 175    new->next = *list;
 176    new->particle = particle;
 177    *list = new;
 178 }
 179
 180 /* --------------------------------------------------------------------------------------------- */
 181
 182 void LinkedListPrint (t_linkedlist *list) {
 183    t_linkedlist *ptr = list;
 184    while (ptr != NULL) {
 185       llog ("%ld ", ptr->particle);
 186       ptr = ptr->next;
 187    }
 188 }
 189
 190 /* --------------------------------------------------------------------------------------------- */
 191
 192 #define CoordinateToCell(x,y,z) (int) (    floor (x*CellLengthI[X]) +    floor (y*CellLengthI[Y])*NCellSide[X] +    floor (z*CellLengthI[Z])*NCellSide[X]*NCellSide[Y]   )
 193
 194
 195 int *BuildLinkedList (int n, float *x, float *y, float *z, float RcNeighbor) {
 196    llog ("\nBuilding cell list...");
 197    unsigned int i;
 198
 199    float CellSize = RcNeighbor;
 200    NCellSide[X] = (unsigned int) ceil (Length[X] / CellSize);
 201    NCellSide[Y] = (unsigned int) ceil (Length[Y] / CellSize);
 202    NCellSide[Z] = (unsigned int) ceil (Length[Z] / CellSize);
 203    NCell = NCellSide[X] * NCellSide[Y] * NCellSide[Z];
 204
 205    CellLength[X] = Length[X] / (float)NCellSide[X];
 206    CellLength[Y] = Length[Y] / (float)NCellSide[Y];
 207    CellLength[Z] = Length[Z] / (float)NCellSide[Z];
 208
 209    CellLengthI[X] = (float)NCellSide[X] / Length[X];
 210    CellLengthI[Y] = (float)NCellSide[Y] / Length[Y];
 211    CellLengthI[Z] = (float)NCellSide[Z] / Length[Z];
 212
 213    llog ("\n\tNcells: %d %d %d -> %d    CellSize: %f   CellLength: %f %f %f",
 214       NCellSide[X], NCellSide[Y], NCellSide[Z], NCell, CellSize, CellLength[X], CellLength[Y], CellLength[Z]);
 215
 216    CellList = (t_linkedlist **) calloc (NCell, sizeof (t_linkedlist));
 217    llog ("\n\tMemory: %fMB", (sizeof (t_linkedlist)*NCell+n*sizeof (unsigned int))/1024./1024.);
 218
 219    for (i = 0; i < n; i++)
 220       LinkedListInsert (i, CellList + CoordinateToCell (x[i], y[i], z[i]) );
 221
 222    return NULL;
 223 }
 224
 225 /* --------------------------------------------------------------------------------------------- */
 226
 227 int *BuildNeighborList (int n, float *x, float *y, float *z, float RcNeighbor) {
 228    llog ("\nBuilding neighbor list...");
 229    return BuildLinkedList (n, x, y, z, RcNeighbor);
 230 }
 231
 232 /* --------------------------------------------------------------------------------------------- */
 233
 234 void MoleculeListNeighborLinkedList (int Atoms, float *x, float *y, float *z, float RcClusterSq) {
 235    llog ("\nBegin clustering...");
 236    double start = clock();
 237
 238    float z0=0.0;
 239    float RcCluster = sqrt (RcClusterSq);
 240    int i, j, k, Lk, Size, s, n;
 241
 242    float xCell[3];
 243    int m[3];
 244    t_linkedlist *ptr;
 245
 246    t_linkedlist **List = (t_linkedlist **) calloc (Atoms, sizeof (t_linkedlist));
 247    int *Spectrum = (int *) calloc (Atoms+1, sizeof (int));
 248    if (!List || !Spectrum)
 249       exit (EMEM);
 250    llog ("\n\t-> %fMB memory\n", (double) (sizeof (t_linkedlist)*Atoms + sizeof (int)*(Atoms+1))/1024./1024. );
 251
 252    for (i = 0; i < Atoms; i++) {
 253       List[i]->particle = i;
 254       List[i]->next = List[i];
 255    }
 256    llog ("\n\tstart");
 257    #define NOT_YET_CLUSTERED(i) (List[i] == List[i]->next)
 258    for (i = 0; i < Atoms; i++) {
 259       PROGRESS(i);
 260       if (NOT_YET_CLUSTERED(i)) {
 261          // set j
 262          j = i;
 263          do {
 264             // sum over all neighbor cells
 265             for (m[X] = -1; m[X] <= 1; m[X]++) {
 266                xCell[X] = x[j]+(float)m[X]*RcCluster;
 267                if (xCell[X] < min[X] || xCell[X] > max[X]) continue;
 268
 269                for (m[Y] = -1; m[Y] <= 1; m[Y]++) {
 270                   xCell[Y] = y[j]+(float)m[Y]*RcCluster;
 271                   if (xCell[Y] < min[Y] || xCell[Y] > max[Y]) continue;
 272
 273                   for (m[Z] = -1; m[Z] <= 1; m[Z]++) {
 274                      xCell[Z] = z[j]+(float)m[Z]*RcCluster;
 275                      if (xCell[Z] < min[Z] || xCell[Z] > max[Z]) continue;
 276
 277                      // sum over all particles in cell
 278                      for (ptr = *(CellList + CoordinateToCell (xCell[X], xCell[Y], xCell[Z])); ptr != NULL; ptr = ptr->next) {
 279                         k = ptr->particle;
 280                         //if (k <= i) continue;
 281                         if (NOT_YET_CLUSTERED(k)) {
 282                            #define RSQ(j,k) SqDist(x[j],x[k],X) + SqDist(y[j],y[k],Y) + SqDist(z[j],z[k],Z)
 283                            // FIXME: LinkedList
 284                            if (RSQ(j,k) <= RcClusterSq) {
 285                               List[k]->next = List[j]->next;
 286                               List[j]->next = List[k]->next;
 287                            }
 288                         }
 289                      }
 290
 291                   }
 292                }
 293             }
 294             // set j
 295             j = List[j]->next->particle;
 296
 297             PROGRESS(j);
 298          }  while (i != j);
 299       }  /* if i */
 300    }  /* for i */
 301
 302    double stop = clock();
 303    llog ("\n\tfinished in %.1f seconds", (stop-start)/CLK);
 304 /*
 305    for (i = 0; i < Atoms+1; i++)
 306       Spectrum[i] = 0;
 307    for (i = 0; i < Atoms; i++) {
 308       if (List[i] >= 0) {
 309          Size = 1;  j = List[i];  List[i] = -1;
 310          while (j != i) { Size++;   k = List[j];  List[j] = -1;  j = k; }
 311          //Bestimmung der Groesse des Clusters im Startteilchen i
 312          Spectrum[Size] += 1;
 313       }
 314    }
 315    printf("\nCoherent Clusters below %4.4f A\n", RcClusterSq);
 316    printf("Size Times Atoms\n");
 317    for (i = 0; i < Atoms+1; i++) // Atoms + 1 because 0 is also a possible Count
 318    {
 319       if ((s = Spectrum[i]) != 0)
 320          printf ("%6d %6d %6d\n", i, s, i*s);
 321    }
 322 */
 323    free (List);
 324    free (Spectrum);
 325    #undef NOT_YET_CLUSTERED
 326 }
 327
 328 /* --------------------------------------------------------------------------------------------- */
 329
 330 void MoleculeListNeighbor (int Atoms, float *x, float *y, float *z, float RcClusterSq) {
 331    llog ("\nBegin clustering...");
 332    double start = clock();
 333
 334    float z0=0.0;
 335    float RcCluster = sqrt (RcClusterSq);
 336    int i, j, k, Lk, Size, *List,*Spectrum,s,n;
 337
 338    float xCell[3];
 339    int m[3];
 340    t_linkedlist *ptr;
 341
 342    List = (INT *) malloc (sizeof (INT) * Atoms);
 343    Spectrum = (INT *) malloc (sizeof (INT) * Atoms+1);
 344    if (!List || !Spectrum)
 345       exit (EMEM);
 346    llog ("\n\t-> %fMB memory\n", (double) (sizeof (INT)*(Atoms*2+1))/1024./1024. );
 347
 348    for (i = 0; i < Atoms; i++)
 349       List[i] = i;
 350
 351    #define NOT_YET_CLUSTERED(i) (i ==  List[i])
 352    for (i = 0; i < Atoms; i++) {
 353       PROGRESS(i);
 354       if (NOT_YET_CLUSTERED(i)) {
 355          // set j
 356          j = i;
 357          do {
 358             // sum over all neighbor cells
 359             for (m[X] = -1; m[X] <= 1; m[X]++) {
 360                xCell[X] = x[j]+(float)m[X]*RcCluster;
 361                if (xCell[X] < min[X] || xCell[X] > max[X]) continue;
 362
 363                for (m[Y] = -1; m[Y] <= 1; m[Y]++) {
 364                   xCell[Y] = y[j]+(float)m[Y]*RcCluster;
 365                   if (xCell[Y] < min[Y] || xCell[Y] > max[Y]) continue;
 366
 367                   for (m[Z] = -1; m[Z] <= 1; m[Z]++) {
 368                      xCell[Z] = z[j]+(float)m[Z]*RcCluster;
 369                      if (xCell[Z] < min[Z] || xCell[Z] > max[Z]) continue;
 370
 371                      // sum over all particles in cell
 372                      for (ptr = *(CellList + CoordinateToCell (xCell[X], xCell[Y], xCell[Z])); ptr != NULL; ptr = ptr->next) {
 373                         k = ptr->particle;
 374                         //if (k <= i) continue;
 375                         if (NOT_YET_CLUSTERED(k)) {
 376                            #define RSQ(j,k) SqDist(x[j],x[k],X) + SqDist(y[j],y[k],Y) + SqDist(z[j],z[k],Z)
 377                            // FIXME: LinkedList
 378                            if (RSQ(j,k) <= RcClusterSq) { Lk = List[k]; List[k] = List[j];  List[j] = Lk; }
 379                         }
 380                      }
 381
 382                   }
 383                }
 384             }
 385             // set j
 386             j = List[j];
 387
 388             PROGRESS(j);
 389          }  while (i != j);
 390       }  /* if i */
 391    }  /* for i */
 392
 393    double stop = clock();
 394    llog ("\n\tfinished in %.1f seconds", (stop-start)/CLK);
 395
 396    for (i = 0; i < Atoms+1; i++)
 397       Spectrum[i] = 0;
 398    for (i = 0; i < Atoms; i++) {
 399       if (List[i] >= 0) {
 400          Size = 1;  j = List[i];  List[i] = -1;
 401          while (j != i) { Size++;   k = List[j];  List[j] = -1;  j = k; }
 402          //Bestimmung der Groesse des Clusters im Startteilchen i
 403          Spectrum[Size] += 1;
 404       }  /* if List */
 405    }  /* for i */
 406    printf("\nCoherent Clusters below %4.4f A\n", RcClusterSq);
 407    printf("Size Times Atoms\n");
 408    for (i = 0; i < Atoms+1; i++) // Atoms + 1 because 0 is also a possible Count
 409    {
 410       if ((s = Spectrum[i]) != 0)
 411          printf ("%6d %6d %6d\n", i, s, i*s);
 412    }
 413
 414    free (List);
 415    free (Spectrum);
 416 }
 417
 418 /* --------------------------------------------------------------------------------------------- */
 419
 420 void MoleculeList (int Atoms, float *x, float *y, float *z, float RcClusterSq) {
 421    llog ("\nBegin clustering...");
 422    double start = clock();
 423
 424    FLOAT     Xj, Yj, Zj, RSq, z0=0.0;
 425
 426    INT        i, j, k, Lk, Size, *List,*Spectrum,s,n;
 427
 428    List = (INT *) malloc (sizeof (INT) * Atoms);
 429    Spectrum = (INT *) malloc (sizeof (INT) * Atoms+1);
 430    if (!List || !Spectrum)
 431       exit (EMEM);
 432    llog ("\n\t-> %fMB memory\n", (double) (sizeof (INT)*(Atoms*2+1))/1024./1024. );
 433
 434    for (i = 0; i < Atoms; i++)
 435       List[i] = i;
 436
 437    for (i = 0; i < Atoms-1; i++) {
 438       PROGRESS(i);
 439       if (i == List[i]) {
 440          j = i;  Xj = x[j];  Yj = y[j];  Zj = z[j];
 441          do {
 442             for (k = i+1; k < Atoms; k++) {
 443                Lk = List[k];
 444                if (k == Lk) {
 445                   RSq = SqDist(Xj,x[k],X) + SqDist(Yj,y[k],Y) + SqDist(Zj,z[k],Z);
 446                   if (RSq <= RcClusterSq) { List[k] = List[j];  List[j] = Lk; }
 447                }  /* if k */
 448             }  /* for k */
 449             j = List[j];  Xj = x[j];  Yj = y[j];  Zj = z[j];
 450             PROGRESS(j);
 451          }  while (i != j);
 452       }  /* if i */
 453    }  /* for i */
 454
 455    double stop = clock();
 456    llog ("\n\tfinished in %.1f seconds", (stop-start)/CLK);
 457
 458    for (i = 0; i < Atoms+1; i++)
 459       Spectrum[i] = 0;
 460    for (i = 0; i < Atoms; i++) {
 461       if (List[i] >= 0) {
 462          Size = 1;  j = List[i];  List[i] = -1;
 463          while (j != i) { Size++;   k = List[j];  List[j] = -1;  j = k; }
 464          //Bestimmung der Groesse des Clusters im Startteilchen i
 465          Spectrum[Size] += 1;
 466       }  /* if List */
 467    }  /* for i */
 468    printf("\nCoherent Clusters below %4.4f A\n", RcClusterSq);
 469    printf("Size Times Atoms\n");
 470    for (i = 0; i < Atoms+1; i++) // Atoms + 1 because 0 is also a possible Count
 471    {
 472       if ((s = Spectrum[i]) != 0)
 473          printf ("%6d %6d %6d\n", i, s, i*s);
 474    }
 475
 476    free (List);
 477    free (Spectrum);
 478 }
 479
 480 /* --------------------------------------------------------------------------------------------- */
 481
 482 void ParseArgs (int argc, char *argv[], float *rcut, unsigned int *Neighborlist, float *RcNeighbor,
 483       unsigned int *LinkedList, unsigned int *FileFlag, char **filename) {
 484    int opt;
 485    while ((opt = getopt(argc, argv, "h:r:nlf:XYZ")) != -1) {
 486       switch (opt) {
 487          case 'h': printf (USAGE); break;
 488          case 'r': *rcut = (float) atof (optarg); break;
 489          case 'n': *Neighborlist = TRUE; /**RcNeighbor = (float) atof (optarg); */break;
 490          case 'l': *LinkedList = TRUE; break;
 491          case 'f': *filename = optarg; *FileFlag = TRUE; break;
 492          case 'X': Periodic[X] = TRUE; break;
 493          case 'Y': Periodic[Y] = TRUE; break;
 494          case 'Z': Periodic[Z] = TRUE; break;
 495          default:  printf (USAGE); break;
 496       }
 497    }
 498    llog ("\n\tPeriodicity: %d %d %d", Periodic[X], Periodic[Y], Periodic[Z]);
 499 }
 500
 501 /* --------------------------------------------------------------------------------------------- */
 502
 503 void ReadStdin (unsigned int *n, float **x, float **y, float **z) {
 504    unsigned int i, ret;
 505    double start = clock();
 506    llog ("\nReading from stdin...");
 507    llog ("\n\tExpecting (x, y, z)");
 508    *n = 0;
 509    *x = NULL;
 510    *y = NULL;
 511    *z = NULL;
 512    float tmp;
 513    do {
 514       *n += 1;
 515       PROGRESS(*n);
 516       i = *n-1;
 517       *x = (float *) realloc (*x, *n*sizeof (float));
 518       *y = (float *) realloc (*y, *n*sizeof (float));
 519       *z = (float *) realloc (*z, *n*sizeof (float));
 520       if (!*x || !*y || !*z)
 521          exit (EMEM);
 522       ret = scanf ("%f %f %f\n", &((*x)[i]), &((*y)[i]), &((*z)[i]));
 523       fflush(stdin);
 524    } while (ret != EOF);
 525    *n -= 1;
 526    double stop = clock();
 527    llog ("\n\tread %d atoms -> memory %gMB in %.1f seconds", *n, (3* *n *sizeof (float))/(1024.*1024.), (stop-start)/CLK);
 528 }
 529
 530 /* --------------------------------------------------------------------------------------------- */
 531
 532 void ReadFile (char *filename, unsigned int *n, float **x, float **y, float **z) {
 533    unsigned int i, ret;
 534    double start = clock ();
 535    llog ("\nReading from %s...", filename);
 536    llog ("\n\tExpecting (x, y, z)");
 537    *n = 0;
 538    *x = NULL;
 539    *y = NULL;
 540    *z = NULL;
 541    float tmp;
 542    FILE *in = fopen (filename, "r");
 543
 544    char command[2048];
 545    sprintf (command, "cat %s | awk '{n++}END{print n}'", filename);
 546    llog ("\n\tUsing %s", command);
 547    FILE *pipe = popen (command, "r");
 548    if (!pipe)
 549       exit (EFILE);
 550    char pos[1024];
 551    fgets (pos, sizeof(pos), pipe);
 552    fclose (pipe);
 553    unsigned int lines = (unsigned int) atoi (pos);
 554    llog ("\n\tWilling to read %d lines", lines);
 555
 556    if (!in)
 557       exit (EFILE);
 558    do {
 559       *n += 1;
 560       PROGRESS(*n);
 561       i = *n-1;
 562       *x = (float *) realloc (*x, *n*sizeof (float));
 563       *y = (float *) realloc (*y, *n*sizeof (float));
 564       *z = (float *) realloc (*z, *n*sizeof (float));
 565       if (!*x || !*y || !*z)
 566          exit (EMEM);
 567 //      ret = fscanf (in, "%f %f %f %f\n", &((*x)[i]), &((*y)[i]), &((*z)[i]), &tmp);
 568       ret = fscanf (in, "%f %f %f\n", &((*x)[i]), &((*y)[i]), &((*z)[i]));
 569    } while (*n-1 < lines);//ret != EOF);
 570    llog ("\nbreak with %d\n", ret);
 571    *n -= 1;
 572    fclose (in);
 573    double stop = clock();
 574    llog ("\n\tread %d atoms -> memory %gMB in %.1f seconds", *n, (3* *n *sizeof (float))/(1024.*1024.), (stop-start)/CLK);
 575 }
 576
 577 /* --------------------------------------------------------------------------------------------- */
 578
 579 void ReadHDF (char *filename, unsigned int *n, float **x, float **y, float **z) {
 580    llog ("\nReading hdf file %s", filename);
 581    char *src = "/Particles/All";
 582    hid_t hdf = H5Fopen (filename, H5F_ACC_RDONLY, H5P_DEFAULT);
 583    hid_t particles = H5Dopen (hdf, src);
 584    hid_t type = H5Dget_type (particles);
 585    hid_t size = H5Dget_space (particles);
 586    hid_t plist = H5Dget_create_plist (particles);
 587
 588    hsize_t nfields, nrecords;
 589    H5TBget_table_info ( hdf, src, &nfields, &nrecords );
 590    llog ("\thas %d fields with %d records\n", nfields, nrecords);
 591
 592    char **field_names = (char **)calloc (nfields, sizeof (char*));
 593    int i;
 594    for (i = 0; i < nfields; i++) {
 595          field_names[i] = (char *)calloc (nfields, sizeof (char[256]));
 596    }
 597    size_t *field_sizes = (size_t *)calloc (nfields, sizeof (size_t)),
 598           *field_offsets = (size_t *)calloc (nfields, sizeof (size_t)),
 599           *type_size = (size_t *)calloc (nfields, sizeof (size_t));
 600
 601    H5TBget_field_info (hdf, src, field_names, field_sizes, field_offsets, type_size);
 602    for (i = 0; i < nfields; i++)
 603       llog ("\t-> %s {%d}\n", field_names[i], field_sizes[i], field_offsets[i], type_size[i]);
 604
 605    size_t total = 0;
 606    for (i = 0; i < nfields; i++)
 607       total += field_sizes[i];
 608    llog ("\ntotal size:%d\n", total);
 609
 610    void *dst_buf = NULL;//calloc (10, total);
 611
 612    hsize_t start = 0, nn = 1;
 613    char *field = "x";
 614    int fidx = -1;
 615    for (i = 0; i < nfields; i++) {
 616       if (strcmp (field_names[i], field) == 0) {
 617          fidx = i;
 618       }
 619    }
 620 //   H5TBread_records (hdf, src, start, nn, type_size[0], field_offsets, field_sizes, dst_buf);
 621    dst_buf = calloc (nn, sizeof(float));
 622    H5TBread_fields_name (hdf, src, field, start, 0, type_size[0],  field_offsets, field_sizes, dst_buf);
 623
 624    void *ptr = dst_buf;
 625 //   for (i = 0; i < fidx; i++)
 626   //    ptr += field_sizes[i];
 627
 628    float *variable1 = (float *) ptr;
 629
 630    llog ("\ndone\n");
 631    //printf ("\nok data:\n");
 632    //for (i = 0; i < nfields; i++)
 633    //   printf ("%f ", dst_buf[i]);
 634    hid_t x1 = H5Tcreate (H5T_COMPOUND, sizeof (float));
 635 //   H5Tinsert (x1, "x", HOFFSET(x1,
 636 //   H5Dread (particles,
 637
 638 //   H5TBread_table (hdf, "/Particles/All", type_size[0], field_offsets, field_sizes, dst_buf);
 639
 640 //   HD5read (particles, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, dst_buf);
 641
 642    H5Dclose (particles);
 643    H5Fclose (hdf);
 644    exit (0);
 645 }
 646
 647 /* --------------------------------------------------------------------------------------------- */
 648
 649 int main (int argc, char *argv[]) {
 650    unsigned int n;
 651    float *x, *y, *z;
 652
 653    float RcClusterSq, RcCluster = 0.;
 654    float RcNeighbor = 0;
 655    unsigned int Neighborlist = FALSE,
 656       LinkedList = FALSE;
 657
 658    unsigned int FileFlag = FALSE;
 659    char *filename[1024];
 660    Periodic[X] = Periodic[Y] = Periodic [Z] = FALSE;
 661
 662    ParseArgs (argc, argv, &RcCluster, &Neighborlist, &RcNeighbor, &LinkedList, &FileFlag, filename);
 663
 664    if (RcCluster == 0.)
 665       exit (ESYN);
 666    else
 667       RcClusterSq = SQ(RcCluster);
 668
 669    if (!FileFlag)
 670       ReadStdin (&n, &x, &y, &z);
 671    else
 672       ReadFile (*filename, &n, &x, &y, &z);
 673
 674    DetermineBoxLength (n, x, y, z, &(Length[X]), &(Length[Y]), &(Length[Z]));
 675
 676    if (Neighborlist) {
 677       RcNeighbor = RcCluster;
 678       BuildNeighborList (n, x, y, z, RcNeighbor);
 679       if (LinkedList)
 680          MoleculeListNeighborLinkedList (n, x, y, z, RcClusterSq);
 681       else
 682          MoleculeListNeighbor (n, x, y, z, RcClusterSq);
 683    } else
 684       MoleculeList (n, x, y, z, RcClusterSq);
 685
 686    return EOK;
 687 }