ClusterLib.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
  17     Based on the algorithm in S. Stoddard JCP(1978)27:291
  18
  19 */
  20
  21 #include <Python.h>
  22 #include <stdio.h>
  23 #include <stdlib.h>
  24 #include <math.h>
  25 #include <time.h>
  26 #include <unistd.h>
  27 #include <stdarg.h>
  28 #include <string.h>
  29 #include "Numeric/arrayobject.h"
  30
  31
  32 #define EOK 0
  33 #define EMEM 1
  34 #define ESYN 2
  35 #define EFILE 3
  36
  37 #define X 0
  38 #define Y 1
  39 #define Z 2
  40
  41 #define TRUE 1
  42 #define FALSE 0
  43
  44 #define MAGIC_MIN 1e10
  45 #define MAGIC_MAX -1e10
  46
  47 #define MAX(a,b) (a>b)?a:b
  48 #define MIN(a,b) (a>b)?b:a
  49
  50 #define FLOAT float
  51 #define INT int
  52
  53 #define USAGE  "-h     help\n" \
  54                "-r     rcut\n" \
  55                "-n     use neighborlist (radius)\n" \
  56                "-l     use linked lists\n" \
  57                "-f     infile (name)\n" \
  58                "-XYZ   peridicity\n"
  59
  60 #define CLK 100000.
  61
  62 #define SQ(a) ((a)*(a))
  63
  64 /* --------------------------------------------------------------------------------------------- */
  65
  66 typedef struct s_linkedlist {
  67    unsigned int particle;
  68    struct s_linkedlist *next;
  69 } t_linkedlist;
  70
  71 /* --------------------------------------------------------------------------------------------- */
  72
  73 int Periodic[3];
  74 float Length[3];
  75 unsigned int NCellSide[3], NCell;
  76 float CellLength[3], CellLengthI[3];
  77 t_linkedlist **CellList;
  78 float min[3], max[3];
  79
  80 /* --------------------------------------------------------------------------------------------- */
  81
  82 double SqDist (double pos1, double pos2, int dim) {
  83    if (Periodic[dim] == TRUE)
  84       return SQ( fmod(pos1-pos2,Length[dim]) );
  85    else
  86       return SQ(pos1-pos2);
  87 }
  88
  89 /* --------------------------------------------------------------------------------------------- */
  90
  91 void DetermineBoxLength (unsigned int n, float *x, float *y, float *z, float *lx, float *ly, float *lz) {
  92    printf ("Determine boxlength\n");
  93    min[X] = min[Y] = min[Z] = MAGIC_MIN;
  94    max[X] = max[Y] = max[Z] = MAGIC_MAX;
  95    unsigned int i;
  96    for (i = 0; i < n; i++) {
  97       min[X] = MIN(min[X],x[i]);
  98       min[Y] = MIN(min[Y],y[i]);
  99       min[Z] = MIN(min[Z],z[i]);
 100       max[X] = MAX(max[X],x[i]);
 101       max[Y] = MAX(max[Y],y[i]);
 102       max[Z] = MAX(max[Z],z[i]);
 103    }
 104
 105    double off[3];
 106    off[X] = 0.;
 107    off[Y] = 0.;
 108    off[Z] = 0.;
 109    unsigned int ShiftFlag = FALSE;
 110    if (min[X] < 0.) {
 111       ShiftFlag = TRUE;
 112       off[X] = -min[X];
 113    }
 114
 115    if (min[Y] < 0.) {
 116       ShiftFlag = TRUE;
 117       off[Y] = -min[Y];
 118     }
 119
 120    if (min[Z] < 0.) {
 121       ShiftFlag = TRUE;
 122       off[Z] = -min[Z];
 123    }
 124
 125    if (ShiftFlag) {
 126       for (i = 0; i < n; i++) {
 127          x[i] += off[X];
 128          y[i] += off[Y];
 129          z[i] += off[Z];
 130       }
 131       unsigned int j;
 132       for (j = 0; j < 3; j++) {
 133          min[j] += off[j];
 134          max[j] += off[j];
 135       }
 136    }
 137
 138    *(lx) = max[X] - min[X];
 139    *(ly) = max[Y] - min[Y];
 140    *(lz) = max[Z] - min[Z];
 141 }
 142
 143 /* --------------------------------------------------------------------------------------------- */
 144
 145 void LinkedListInsert (unsigned int particle, t_linkedlist **list) {
 146    t_linkedlist *new = (t_linkedlist *)calloc (1, sizeof (t_linkedlist));
 147    if (!new)
 148       exit (EMEM);
 149    new->next = *list;
 150    new->particle = particle;
 151    *list = new;
 152 }
 153
 154 /* --------------------------------------------------------------------------------------------- */
 155
 156 void LinkedListPrint (t_linkedlist *list) {
 157    t_linkedlist *ptr = list;
 158    while (ptr != NULL) {
 159       ptr = ptr->next;
 160    }
 161 }
 162
 163 /* --------------------------------------------------------------------------------------------- */
 164
 165 #define CoordinateToCell(x,y,z) (int) (    floor (x*CellLengthI[X]) +    floor (y*CellLengthI[Y])*NCellSide[X] +    floor (z*CellLengthI[Z])*NCellSide[X]*NCellSide[Y]   )
 166
 167
 168 int *BuildLinkedList (int n, float *x, float *y, float *z, float RcNeighbor) {
 169    printf ("Build linked list\n");
 170    unsigned int i;
 171
 172    float CellSize = RcNeighbor;
 173    NCellSide[X] = (unsigned int) ceil (Length[X] / CellSize);
 174    NCellSide[Y] = (unsigned int) ceil (Length[Y] / CellSize);
 175    NCellSide[Z] = (unsigned int) ceil (Length[Z] / CellSize);
 176    NCell = NCellSide[X] * NCellSide[Y] * NCellSide[Z];
 177
 178    CellLength[X] = Length[X] / (float)NCellSide[X];
 179    CellLength[Y] = Length[Y] / (float)NCellSide[Y];
 180    CellLength[Z] = Length[Z] / (float)NCellSide[Z];
 181
 182    CellLengthI[X] = (float)NCellSide[X] / Length[X];
 183    CellLengthI[Y] = (float)NCellSide[Y] / Length[Y];
 184    CellLengthI[Z] = (float)NCellSide[Z] / Length[Z];
 185
 186    CellList = (t_linkedlist **) calloc (NCell, sizeof (t_linkedlist));
 187
 188    for (i = 0; i < n; i++)
 189       LinkedListInsert (i, CellList + CoordinateToCell (x[i], y[i], z[i]) );
 190
 191    return NULL;
 192 }
 193
 194 /* --------------------------------------------------------------------------------------------- */
 195
 196 int *BuildNeighborList (int n, float *x, float *y, float *z, float RcNeighbor) {
 197    printf ("Build neighbos list\n");
 198    return BuildLinkedList (n, x, y, z, RcNeighbor);
 199 }
 200
 201 /* --------------------------------------------------------------------------------------------- */
 202
 203 void MoleculeListNeighbor (int Atoms, float *x, float *y, float *z, float RcClusterSq) {
 204    printf ("Clustering: neighbor list\n");
 205    double start = clock();
 206
 207    float z0=0.0;
 208    float RcCluster = sqrt (RcClusterSq);
 209    int i, j, k, Lk, Size, *List,*Spectrum,s,n;
 210
 211    float xCell[3];
 212    int m[3];
 213    t_linkedlist *ptr;
 214
 215    List = (INT *) malloc (sizeof (INT) * Atoms);
 216    Spectrum = (INT *) malloc (sizeof (INT) * Atoms+1);
 217    if (!List || !Spectrum)
 218       exit (EMEM);
 219
 220    for (i = 0; i < Atoms; i++)
 221       List[i] = i;
 222
 223    #define NOT_YET_CLUSTERED(i) (i ==  List[i])
 224    for (i = 0; i < Atoms; i++) {
 225       if (NOT_YET_CLUSTERED(i)) {
 226          // set j
 227          j = i;
 228          do {
 229             // sum over all neighbor cells
 230             for (m[X] = -1; m[X] <= 1; m[X]++) {
 231                xCell[X] = x[j]+(float)m[X]*RcCluster;
 232                if (xCell[X] < min[X] || xCell[X] > max[X]) continue;
 233
 234                for (m[Y] = -1; m[Y] <= 1; m[Y]++) {
 235                   xCell[Y] = y[j]+(float)m[Y]*RcCluster;
 236                   if (xCell[Y] < min[Y] || xCell[Y] > max[Y]) continue;
 237
 238                   for (m[Z] = -1; m[Z] <= 1; m[Z]++) {
 239                      xCell[Z] = z[j]+(float)m[Z]*RcCluster;
 240                      if (xCell[Z] < min[Z] || xCell[Z] > max[Z]) continue;
 241
 242                      // sum over all particles in cell
 243                      for (ptr = *(CellList + CoordinateToCell (xCell[X], xCell[Y], xCell[Z])); ptr != NULL; ptr = ptr->next) {
 244                         k = ptr->particle;
 245                         //if (k <= i) continue;
 246                         if (NOT_YET_CLUSTERED(k)) {
 247                            #define RSQ(j,k) SqDist(x[j],x[k],X) + SqDist(y[j],y[k],Y) + SqDist(z[j],z[k],Z)
 248                            // FIXME: LinkedList
 249                            if (RSQ(j,k) <= RcClusterSq) { Lk = List[k]; List[k] = List[j];  List[j] = Lk; }
 250                         }
 251                      }
 252
 253                   }
 254                }
 255             }
 256             // set j
 257             j = List[j];
 258
 259          }  while (i != j);
 260       }  /* if i */
 261    }  /* for i */
 262
 263    double stop = clock();
 264
 265    for (i = 0; i < Atoms+1; i++)
 266       Spectrum[i] = 0;
 267    for (i = 0; i < Atoms; i++) {
 268       if (List[i] >= 0) {
 269          Size = 1;  j = List[i];  List[i] = -1;
 270          while (j != i) { Size++;   k = List[j];  List[j] = -1;  j = k; }
 271          //Bestimmung der Groesse des Clusters im Startteilchen i
 272          Spectrum[Size] += 1;
 273       }  /* if List */
 274    }  /* for i */
 275    printf("\nCoherent Clusters below %4.4f A\n", RcClusterSq);
 276    printf("Size Times Atoms\n");
 277    for (i = 0; i < Atoms+1; i++) // Atoms + 1 because 0 is also a possible Count
 278    {
 279       if ((s = Spectrum[i]) != 0)
 280          printf ("%6d %6d %6d\n", i, s, i*s);
 281    }
 282
 283    free (List);
 284    free (Spectrum);
 285 }
 286
 287 /* --------------------------------------------------------------------------------------------- */
 288
 289 void MoleculeList (int Atoms, float *x, float *y, float *z, float RcClusterSq) {
 290    printf ("Clustering\n");
 291    double start = clock();
 292
 293    FLOAT     Xj, Yj, Zj, RSq, z0=0.0;
 294
 295    INT        i, j, k, Lk, Size, *List,*Spectrum,s,n;
 296
 297    List = (INT *) malloc (sizeof (INT) * Atoms);
 298    Spectrum = (INT *) malloc (sizeof (INT) * Atoms+1);
 299    if (!List || !Spectrum)
 300       exit (EMEM);
 301
 302    for (i = 0; i < Atoms; i++)
 303       List[i] = i;
 304
 305    for (i = 0; i < Atoms-1; i++) {
 306       if (i == List[i]) {
 307          j = i;  Xj = x[j];  Yj = y[j];  Zj = z[j];
 308          do {
 309             for (k = i+1; k < Atoms; k++) {
 310                Lk = List[k];
 311                if (k == Lk) {
 312                   RSq = SqDist(Xj,x[k],X) + SqDist(Yj,y[k],Y) + SqDist(Zj,z[k],Z);
 313                   if (RSq <= RcClusterSq) { List[k] = List[j];  List[j] = Lk; }
 314                }  /* if k */
 315             }  /* for k */
 316             j = List[j];  Xj = x[j];  Yj = y[j];  Zj = z[j];
 317          }  while (i != j);
 318       }  /* if i */
 319    }  /* for i */
 320
 321    double stop = clock();
 322
 323    for (i = 0; i < Atoms+1; i++)
 324       Spectrum[i] = 0;
 325    for (i = 0; i < Atoms; i++) {
 326       if (List[i] >= 0) {
 327          Size = 1;  j = List[i];  List[i] = -1;
 328          while (j != i) { Size++;   k = List[j];  List[j] = -1;  j = k; }
 329          //Bestimmung der Groesse des Clusters im Startteilchen i
 330          Spectrum[Size] += 1;
 331       }  /* if List */
 332    }  /* for i */
 333    printf("\nCoherent Clusters below %4.4f A\n", RcClusterSq);
 334    printf("Size Times Atoms\n");
 335    for (i = 0; i < Atoms+1; i++) // Atoms + 1 because 0 is also a possible Count
 336    {
 337       if ((s = Spectrum[i]) != 0)
 338          printf ("%6d %6d %6d\n", i, s, i*s);
 339    }
 340
 341    free (List);
 342    free (Spectrum);
 343 }
 344
 345 /* --------------------------------------------------------------------------------------------- */
 346
 347 static PyObject * clusterdetector(PyObject *self, PyObject *args) {
 348    PyObject *inx, *iny, *inz;
 349    PyArrayObject *xx, *yy, *zz;
 350    double rcut;
 351    int list = 0;
 352
 353    if (!PyArg_ParseTuple (args, "OOOd|i", &inx, &iny, &inz, &rcut, &list))
 354       return NULL;
 355
 356    xx = (PyArrayObject *)PyArray_ContiguousFromObject (inx, PyArray_DOUBLE, 1, 0);
 357    yy = (PyArrayObject *)PyArray_ContiguousFromObject (iny, PyArray_DOUBLE, 1, 0);
 358    zz = (PyArrayObject *)PyArray_ContiguousFromObject (inz, PyArray_DOUBLE, 1, 0);
 359
 360    int n = xx->dimensions[0];
 361
 362    float *x = (float *) calloc (n, sizeof (float));
 363    float *y = (float *) calloc (n, sizeof (float));
 364    float *z = (float *) calloc (n, sizeof (float));
 365
 366    int i;
 367    for (i=0;i <n;i++) {
 368       x[i] = (float)*(double*) (xx->data+i*xx->strides[0]);
 369       y[i] = (float)*(double*) (yy->data+i*yy->strides[0]);
 370       z[i] = (float)*(double*) (zz->data+i*zz->strides[0]);
 371    }
 372
 373    DetermineBoxLength (n, x, y, z, &(Length[X]), &(Length[Y]), &(Length[Z]));
 374
 375    float RcClusterSq, RcCluster = 0.;
 376    float RcNeighbor = 0;
 377    RcNeighbor = RcCluster;
 378
 379    if (list != 0) {
 380       BuildNeighborList (n, x, y, z, RcNeighbor);
 381       MoleculeListNeighbor (n, x, y, z, RcClusterSq);
 382    } else {
 383       MoleculeList (n, x, y, z, RcClusterSq);
 384    }
 385
 386
 387    Py_DECREF (xx);
 388    Py_DECREF (yy);
 389    Py_DECREF (zz);
 390
 391    return  PyFloat_FromDouble (Length[Z]);
 392 }
 393
 394
 395 static struct PyMethodDef cluster_methods[] = {
 396    {"clusterdetector", clusterdetector, 1},
 397    {NULL, NULL}
 398 };
 399
 400 PyMODINIT_FUNC initcluster (void) {
 401    Py_InitModule ("cluster", cluster_methods);
 402    import_array ();
 403 }